language-Modula2-0.1: examples/Modula-2_Libraries/andrea-m2/lib/generic/largenumb.mod
IMPLEMENTATION MODULE LargeNumbers;
(* based on ExNumLib in "Modula-2 Library Modules: A Programmers Reference"
by Robert D. Walker, Tab Books Inc, 1988 *)
(* V3.0, J.Andrea, Nov.14/92 - character operations working *)
(* V2.1, J.Andrea, Nov.13/92 - use lookup table for character operations *)
(* V2.0, J.Andrea, Nov.13/92 - store mantissa as characters not integers *)
(* V1.0, J.Andrea, Nov.13/92 - add tests for non-initialized numbers *)
(* V0.0, implemented and some fixes by J.Andrea, Nov.12/92 *)
FROM StringOperations IMPORT Length, Append;
FROM Conversions IMPORT IntToString;
CONST
exp_max_len = 7; (* how many digits could an exponent hold *)
max_real_digits = 7; (* precision of a REAL type floating point number *)
an_value = 12321; (* an unlikely value in garbage memory space *)
VAR
max_digits :INTEGER;
i :INTEGER;
c, d :CHAR;
digit_to_int :ARRAY [0C..11C] OF INTEGER; (* range chr(0) to chr(9) *)
int_to_digit :ARRAY [0..9] OF CHAR;
char_to_digit :ARRAY ['0'..'9'] OF CHAR;
digit_to_char :ARRAY [0C..11C] OF CHAR;
char_to_int :ARRAY ['0'..'9'] OF INTEGER;
mul_digits :ARRAY [0C..11C] OF ARRAY [0C..11C] OF INTEGER;
(* ----------------- *)
PROCEDURE IsAN( a :NumType ) :BOOLEAN;
(* is the large number A Number, has it been initialized ? *)
BEGIN
RETURN a.an = an_value;
END IsAN;
(* ----------------- *)
PROCEDURE xtoi( x :REAL; k :INTEGER ) :REAL;
(* raise a real number to an integer power *)
VAR
y :REAL;
i :INTEGER;
BEGIN
y := 1.0; (* initialize, and default for k=0 *)
IF k > 0 THEN
FOR i := 1 TO k DO
y := y * x;
END;
ELSE
IF k < 0 THEN
FOR i := 1 TO k DO
y := y / x;
END;
END;
END;
RETURN y;
END xtoi;
(* -------------------------------------------- *)
PROCEDURE SetDigits( d :INTEGER );
BEGIN
IF d < 2 THEN
max_digits := 2;
status := too_few_digits;
ELSIF d > max_precision THEN
max_digits := max_precision;
status := too_many_digits;
ELSE
max_digits := d;
status := ok;
END;
END SetDigits;
(* -------------------------------------------- *)
PROCEDURE Add( a, b :NumType; VAR c :NumType );
BEGIN
IF IsAN( a ) & IsAN( b ) THEN
IF a.sign = b.sign THEN
AddUtility( a, b, c );
IF a.sign = negative THEN
ChangeSign( c );
END;
ELSE
IF a.sign = positive THEN
SubtractUtility( a, b, c );
ELSE
SubtractUtility( b, a, c );
END;
END;
ELSE
c := zero;
status := nan;
END;
END Add;
(* -------------------------------------------- *)
PROCEDURE Subtract( a, b :NumType; VAR c :NumType );
BEGIN
IF IsAN( a ) & IsAN( b ) THEN
IF a.sign = b.sign THEN
IF a.sign = positive THEN
SubtractUtility( a, b, c );
ELSE
SubtractUtility( b, a, c );
END;
ELSE
AddUtility( a, b, c );
IF a.sign = negative THEN
ChangeSign( c );
END;
END;
ELSE
c := zero;
status := nan;
END;
END Subtract;
(* -------------------------------------------- *)
PROCEDURE Multiply( a, b :NumType; VAR c :NumType );
VAR
i, j, product, carry :INTEGER;
ex1 :NumType;
BEGIN
IF IsAN( a ) & IsAN( b ) THEN
IF ( Compare( a, zero ) = 0 ) OR ( Compare( b, zero ) = 0 ) THEN
(* multiplication by zero *)
c := zero;
status := ok;
ELSIF Compare( a, one ) = 0 THEN
(* multiplication by one *)
c := b;
status := ok;
ELSIF Compare( b, one ) = 0 THEN
(* multiplication by one *)
c := a;
status := ok;
ELSE (* regular multiplication *)
status := ok;
c := zero;
FOR i := max_digits - 1 TO 0 BY -1 DO
ex1 := zero;
ex1.exp := a.exp + b.exp - i;
carry := 0;
FOR j := max_digits - 1 TO 0 BY -1 DO
product := mul_digits[a.man[j],b.man[i]] + carry;
ex1.man[j] := int_to_digit[product MOD 10];
carry := product DIV 10;
END;
(* check for a final carry *)
IF carry > 0 THEN
ShiftRight( ex1 );
ex1.man[0] := int_to_digit[carry];
Times10( ex1 );
END;
(* add up what we've got so far *)
AddUtility( c, ex1, c );
END;
(* adjust product sign *)
IF a.sign # b.sign THEN
ChangeSign( c );
END;
END;
ELSE
c := zero;
status := nan;
END;
END Multiply;
(* -------------------------------------------- *)
PROCEDURE Divide( a, b :NumType; VAR c :NumType );
VAR
i, quotient :INTEGER;
ex1, ex2 :NumType;
BEGIN
IF IsAN( a ) & IsAN( b ) THEN
status := ok;
IF Compare( b, zero ) = 0 THEN
status := divide_by_zero;
ELSIF Compare( a, zero ) = 0 THEN
(* dividend is zero *)
c := zero;
ELSIF Compare( b, one ) = 0 THEN
(* divide by one *)
c := a;
ELSE (* regular division *)
c := zero;
c.exp := a.exp - b.exp;
ex1 := a; (* let ex1 = ABS( a ) / magnitude( a ) *)
Abs( ex1 );
ex1.exp := 0;
ex2 := b; (* let ex2 = ABS( b ) / magnitude( b ) *)
Abs( ex2 );
ex2.exp := 0;
FOR i := 0 TO max_digits-1 DO
quotient := 0;
WHILE Compare( ex1, ex2 ) = 1 DO (* while ex1 > ex2 *)
quotient := quotient + 1;
SubtractUtility( ex1, ex2, ex1 );
END;
c.man[i] := int_to_digit[quotient];
Divide10( ex2 );
END;
(* normalize quotient *)
Norm( c );
(* adjust quotient sign *)
IF a.sign # b.sign THEN
ChangeSign( c );
END;
END;
ELSE
c := zero;
status := nan;
END;
END Divide;
(* -------------------------------------------- *)
PROCEDURE Abs( VAR a :NumType );
BEGIN
a.sign := positive;
END Abs;
(* -------------------------------------------- *)
PROCEDURE ChangeSign( VAR a :NumType );
BEGIN
IF a.sign = positive THEN
a.sign := negative;
ELSE
a.sign := positive;
END;
END ChangeSign;
(* -------------------------------------------- *)
PROCEDURE Times10( VAR a :NumType );
BEGIN
IF IsAN( a ) THEN
a.exp := a.exp + 1;
IF a.exp > max_exp THEN
status := overflow;
ELSE
status := ok;
END;
ELSE
a := zero;
status := nan;
END;
END Times10;
(* -------------------------------------------- *)
PROCEDURE Divide10( VAR a :NumType );
BEGIN
IF IsAN( a ) THEN
a.exp := a.exp - 1;
IF a.exp < min_exp THEN
status := underflow;
ELSE
status := ok;
END;
ELSE
a := zero;
status := nan;
END;
END Divide10;
(* -------------------------------------------- *)
PROCEDURE Compare( a, b :NumType ) :INTEGER;
VAR
done :BOOLEAN;
i :INTEGER;
BEGIN
IF IsAN( a ) & IsAN( b ) THEN
IF a.sign # b.sign THEN
IF a.sign = positive THEN
RETURN 1; (* a is +, b is -, so a > b *)
ELSE
RETURN -1; (* a is -, b is +, so a < b *)
END;
ELSE
(* a and b have the same sign *)
IF a.exp # b.exp THEN
IF a.exp > b.exp THEN
IF a.sign = positive THEN
RETURN 1; (* a > b *)
ELSE
RETURN -1; (* a < b *)
END;
ELSE
(* a and b have same sign, and a exponent <= b exponent *)
IF a.sign = positive THEN
RETURN -1;
ELSE
RETURN 1;
END;
END;
ELSE
(* a and b have the same sign, and the same exponent *)
done := FALSE;
i := 0;
(* compare each digit till a difference is found, or the end *)
WHILE ( i <= max_digits - 1 ) & NOT done DO
IF a.man[i] # b.man[i] THEN
done := TRUE;
ELSE
i := i + 1;
END;
END;
IF i > max_digits - 1 THEN
(* end reached, so all digits match *)
RETURN 0;
ELSE
(* compare those last digits which are different *)
IF a.man[i] < b.man[i] THEN
IF a.sign = positive THEN
RETURN -1; (* a<b *)
ELSE
RETURN 1; (* a>b *)
END;
ELSE
IF a.sign = positive THEN
RETURN 1; (* a>b *)
ELSE
RETURN -1;
END;
END;
END;
END;
END;
ELSE
RETURN -1; (* this is questionable *)
END;
END Compare;
(* -------------------------------------------- *)
PROCEDURE Norm( VAR a :NumType );
VAR
value_zero :BOOLEAN;
i :INTEGER;
BEGIN
IF IsAN( a ) THEN
(* check for a zero *)
i := 0;
value_zero := TRUE;
WHILE ( i <= max_digits - 1 ) & value_zero DO
IF a.man[i] # int_to_digit[0] THEN
value_zero := FALSE;
END;
i := i + 1;
END;
IF value_zero THEN
(* normalize zero *)
a.sign := positive;
a.exp := 0;
ELSE
(* shift mantissa to left until most signifigant digit is non-zero,
increment exponent with each shift *)
WHILE a.man[0] = int_to_digit[0] DO
ShiftLeft( a );
Divide10( a )
END;
END;
ELSE
a := zero;
status := nan;
END;
END Norm;
(* -------------------------------------------- *)
PROCEDURE StringToNum( s :ARRAY OF CHAR; VAR a :NumType );
VAR
done, neg_exp :BOOLEAN;
i, len :CARDINAL;
j, exponent :INTEGER;
second_digit :BOOLEAN;
(* ---------------- *)
PROCEDURE Numeral( c :CHAR ) :BOOLEAN;
(* is the character a number *)
BEGIN
RETURN ( c >= '0' ) & ( c <= '9' );
END Numeral;
BEGIN
status := ok;
done := TRUE;
a := zero;
exponent := 0;
neg_exp := FALSE;
len := Length( s );
IF len > 0 THEN
i := 0;
(* chek for a negative mantissa *)
IF s[i] = '-' THEN
a.sign := negative;
i := i + 1;
END;
(* skip leading zeros *)
WHILE ( i < len ) & ( s[i] = '0' ) DO
i := i + 1;
END;
IF i < len THEN
(* more than just zeros in this string *)
j := 0;
second_digit := FALSE;
(* handle digits to the left of the decimal point *)
WHILE ( i < len ) & done DO
IF Numeral( s[i] ) THEN
a.man[j] := char_to_digit[s[i]];
IF second_digit THEN
Times10( a );
ELSE
second_digit := TRUE;
END;
j := j + 1;
i := i + 1;
ELSE
done := FALSE;
END;
END;
IF NOT done THEN
(* non-numeric was found, could it be something that we want *)
IF ( s[i] = '.' ) OR ( CAP(s[i]) = 'E' ) THEN
IF s[i] = '.' THEN
i := i + 1;
done := TRUE;
(* get all the digits to the right of the decimal point *)
WHILE ( i < len ) & done DO
IF Numeral( s[i] ) THEN
a.man[j] := char_to_digit[s[i]];
j := j + 1;
i := i + 1;
ELSE
done := FALSE;
END;
END;
END;
IF NOT done THEN
IF CAP( s[i] ) = 'E' THEN
i := i + 1;
neg_exp := s[i] = '-';
IF ( s[i] = '+' ) OR neg_exp THEN
i := i + 1; (* skip the sign *)
END;
(* handle the exponent *)
exponent := 0;
done := TRUE;
WHILE ( i < len ) & done DO
IF Numeral( s[i] ) THEN
exponent := exponent * 10 + char_to_int[s[i]];
ELSE
done := FALSE;
END;
END;
IF NOT done THEN
status := invalid_digit;
ELSE
IF neg_exp THEN
a.exp := a.exp - exponent;
ELSE
a.exp := a.exp + exponent;
END;
END;
ELSE
status := invalid_digit;
END;
END;
ELSE
status := invalid_digit;
END;
END;
END;
END;
END StringToNum;
(* -------------------------------------------- *)
PROCEDURE NumToString( a :NumType; VAR s :ARRAY OF CHAR );
VAR
i, last :INTEGER;
s_exp :ARRAY [0..exp_max_len] OF CHAR;
BEGIN
(* create a null string *)
s[0] := 0C;
IF IsAN( a ) THEN
(* the sign *)
IF a.sign = negative THEN
Append( '-', s );
END;
Append( digit_to_char[a.man[0]], s );
Append( '.', s );
(* find the last non-zero digit *)
last := max_digits - 1;
WHILE ( last >= 1 ) & ( a.man[last] = 0C ) DO
last := last - 1;
END;
FOR i := 1 TO last DO
Append( digit_to_char[a.man[i]], s );
END;
(* the exponent *)
Append( 'E', s );
IntToString( a.exp, 0, s_exp );
Append( s_exp, s );
ELSE
Append( '0', s );
status := nan;
END;
END NumToString;
(* -------------------------------------------- *)
PROCEDURE AddUtility( a, b :NumType; VAR c :NumType );
VAR
i, j, carry, digit, result :INTEGER;
ex1, ex2 :NumType;
BEGIN
IF a.exp > b.exp THEN
ex1 := a;
ex2 := b;
ELSE
ex1 := b;
ex2 := a;
END;
c := zero;
c.exp := ex1.exp;
carry := 0;
FOR i := max_digits-1 TO 0 BY -1 DO
j := i + ex2.exp - ex1.exp; (* j = index to ex2 *)
(* check that j falls within bounds *)
IF ( j >= 0 ) & ( j <= max_digits -1 ) THEN
digit := digit_to_int[ex2.man[j]]; (* get digit from ex2 *)
ELSE
digit := 0; (* j is outside, use zero *);
END;
(* perform addition with carry *)
result := digit_to_int[ex1.man[i]] + digit + carry;
(* check for next carry *)
IF result >= 10 THEN
result := result MOD 10;
carry := 1;
ELSE
carry := 0;
END;
c.man[i] := int_to_digit[result];
END;
(* handle final carry *)
IF carry > 0 THEN
ShiftRight( c ); (* shift mantissa to right *)
c.man[0] := 1C; (* make first digit = 1 *)
Times10( c ); (* update exponent *)
END;
(* set status *)
IF c.exp > max_exp THEN
status := overflow;
ELSE
status := ok;
END;
END AddUtility;
(* -------------------------------------------- *)
PROCEDURE SubtractUtility( a, b:NumType; VAR c :NumType );
VAR
i, j, borrow, digit, result :INTEGER;
ex1, ex2 :NumType;
positive_result :BOOLEAN;
BEGIN
Abs( a ); Abs( b );
IF a.exp > b.exp THEN
ex1 := a;
ex2 := b;
ELSE
ex1 := b;
ex2 := a;
END;
positive_result := Compare( ex1, ex2 ) = 1;
c := zero;
c.exp := ex1.exp;
borrow := 0;
FOR i := max_digits-1 TO 0 BY -1 DO
j := i + ex2.exp - ex1.exp; (* j is index into ex2 *)
(* check that j falls within bounds *)
IF ( j >= 0 ) & ( j <= max_digits -1 ) THEN
digit := digit_to_int[ex2.man[j]]; (* get digit from ex2 *)
ELSE
digit := 0; (* j outside array, user zero *)
END;
(* perform subtraction with borrowing *)
result := digit_to_int[ex1.man[i]];
IF positive_result THEN
result := result - digit;
ELSE
result := digit - result;
END;
result := result - borrow;
(* check for borrow *)
IF result < 0 THEN
result := result + 10;
borrow := 1;
ELSE
borrow := 0;
END;
c.man[i] := int_to_digit[result];
END;
(* normalize *)
Norm( c );
(* adjust sign *)
IF Compare( a, b ) = -1 THEN
ChangeSign( c );
END;
END SubtractUtility;
(* -------------------------------------------- *)
PROCEDURE ShiftRight( VAR a :NumType );
VAR
i :INTEGER;
BEGIN
FOR i := max_digits TO 1 BY -1 DO
a.man[i] := a.man[i-1];
END;
a.man[0] := int_to_digit[0]; (* put zero in most signifigant position *)
END ShiftRight;
(* -------------------------------------------- *)
PROCEDURE ShiftLeft( VAR a :NumType );
VAR
i :INTEGER;
BEGIN
FOR i := -1 TO max_digits - 2 DO
a.man[i] := a.man[i+1];
END;
a.man[max_digits-1] := int_to_digit[0]; (* put zero in least signifigant position *)
END ShiftLeft;
BEGIN
(* precompute conversions *)
FOR i := 0 TO 9 DO
digit_to_int[ CHR(i) ] := i;
int_to_digit[ i ] := CHR(i);
END;
i := 0;
FOR c := '0' TO '9' DO
char_to_digit[ c ] := CHR( i );
digit_to_char[ CHR(i) ] := c;
char_to_int[ c ] := i;
i := i + 1;
END;
(* precompute tables of adding and multiplying together character digits *)
FOR c := 0C TO 11C DO
FOR d := 0C TO 11C DO
mul_digits[c,d] := digit_to_int[c] * digit_to_int[d];
END;
END;
(* set default for number size *)
SetDigits( max_precision );
zero.sign := positive;
zero.exp := 0;
FOR i := -1 TO max_precision DO
zero.man[i] := int_to_digit[0];
END;
zero.an := an_value;
one := zero;
one.man[0] := int_to_digit[1];
END LargeNumbers.