IMPLEMENTATION MODULE NUM; (******************************************************************************* Environnement : Solaris 2.3 Modula-2 ------------- Programmeurs : Pierre L'Ecuyer, Benoit Martin. ----------- *******************************************************************************) IMPORT lMathLib; FROM MYINOUT IMPORT Write, WriteLn, WriteString, WriteInt, ReadInt, ReadCard, WriteCard, ReadLn, ReadLongReal, WriteLongInt, Error; FROM ConLR IMPORT LRealToStr, StrToLReal; FROM ConNum IMPORT StrToSigned_32; FROM Strings IMPORT Length; FROM lMathLib IMPORT dpower, dsqrt, dexp, dln, dsin, dcos; VAR i : INTEGER; (* This should be a constant, but some compilers do not allow it. *) RcpLn2 : LONGREAL; (* = 1 / Ln (2). *) Cor : LONGREAL; Str : ARRAY [0..31] OF CHAR; c, Dot : CARDINAL; cc : CARDINAL; OK : BOOLEAN; PROCEDURE WriteFormatted (S : ARRAY OF CHAR; n : INTEGER); VAR i : CARDINAL; z : INTEGER; BEGIN z := Length(S); IF z >= n THEN z:= 0; ELSE z := n-z END; FOR i:= 1 TO z DO Write (' ') END; WriteString (S); END WriteFormatted; PROCEDURE S (x : LONGREAL) : LONGREAL; (* Si x est positif, calcule le plus petit LONGREAL strictement *) (* superieur a x. Verifie S(0) = 0 et S(-x) = -S(x). *) BEGIN RETURN x + jLR*x END S; PROCEDURE Pow (x, y : LONGREAL) : LONGREAL; BEGIN IF x >= 0.0 THEN RETURN dpower (x, y); ELSE IF ModLR (y, 2.0) = 0.0 THEN RETURN dpower (-x, y); ELSE RETURN -dpower (-x,y); END; END; END Pow; PROCEDURE Sqrt (x : LONGREAL) : LONGREAL; BEGIN RETURN dsqrt (x) END Sqrt; PROCEDURE Exp (x : LONGREAL) : LONGREAL; BEGIN RETURN dexp (x) END Exp; PROCEDURE Ln (x : LONGREAL) : LONGREAL; BEGIN RETURN dln (x) END Ln; PROCEDURE Lg (x : LONGREAL) : LONGREAL; BEGIN RETURN dln (x) * RcpLn2 END Lg; PROCEDURE Sin (x : LONGREAL) : LONGREAL; BEGIN RETURN dsin (x) END Sin; PROCEDURE Cos (x : LONGREAL) : LONGREAL; BEGIN RETURN dcos (x) END Cos; PROCEDURE IToLI (i : INTEGER ) : LONGINT; BEGIN RETURN VAL (LONGINT, i) END IToLI; PROCEDURE IToLR (i : INTEGER ) : LONGREAL; BEGIN RETURN LFLOAT (i) END IToLR; PROCEDURE LIToLR (i : LONGINT) : LONGREAL; BEGIN RETURN LFLOAT (i) END LIToLR; PROCEDURE LCToLR (c : LONGCARD) : LONGREAL; BEGIN RETURN LFLOAT (c) END LCToLR; PROCEDURE REALToLR (x : REAL) : LONGREAL; BEGIN RETURN LFLOAT (x) END REALToLR; PROCEDURE StringToLR (S : ARRAY OF CHAR) : LONGREAL; VAR x : LONGREAL; BEGIN StrToLReal (S, x, OK); RETURN x; END StringToLR; PROCEDURE LRToString (x : LONGREAL; i : INTEGER; S : ARRAY OF CHAR); BEGIN LRealToStr (x, i, S) END LRToString; PROCEDURE StringToLI (S : ARRAY OF CHAR) : LONGINT; VAR Ok : BOOLEAN; s : LONGINT; BEGIN StrToSigned_32 (S, 10, s, Ok); RETURN s END StringToLI; PROCEDURE ModLR (x, y : LONGREAL) : LONGREAL; VAR z : LONGREAL; BEGIN IF y <= 0.0 THEN Error ("ERROR ModLR in NUM: y < 0.0") END; z := ABS (x/y); IF z > MaxIntLongReal THEN Error ("ERROR ModLR in NUM: x/y not representable") END; IF x >= 0.0 THEN RETURN x - y * TruncLR (z) END; z := x + y * TruncLR (z); IF z < 0.0 THEN z := z + y END; RETURN z END ModLR; PROCEDURE TruncLR (x : LONGREAL) : LONGREAL; VAR s : INTEGER; Value : LONGREAL; BEGIN Value := ABS ( x ); IF Value <= MaxLongIntLR THEN RETURN LFLOAT (TRUNC (x)) END; IF Value > MaxIntLongReal THEN RETURN ( x ); END; (* MaxLongIntLR < Value <= MaxIntLongReal *) s := TRUNC (Ln (Value) * RcpLn2); (* we must ensure that 2**s =< Value =< 2**(s+1) , because of the *) (* precision problems caused by the Ln operation *) IF Value >= TwoExp [s+1] THEN INC (s) ELSE IF Value < TwoExp [s] THEN DEC (s) END END; (* LgMaxLongInt1 <= s <= LgMaxIntLongReal, 2**s =< Value =< 2**(s+1) *) Value := Value - TwoExp [s]; WHILE s > LgMaxLongInt1 DO DEC ( s ); IF Value >= TwoExp [s] THEN Value := Value - TwoExp [s] END END; (* we have s =< LgMaxLongInt1 and 0 =< Value < 2^LgMaxLongInt1 *) Value := Value - LFLOAT (TRUNC (Value)); (* Value = decimal part of ABS (x) *) IF x >= ZeroLR THEN RETURN (x - Value) ELSE RETURN (x + Value) END; END TruncLR; PROCEDURE RoundLR ( x : LONGREAL ) : LONGREAL; VAR Value : LONGREAL; BEGIN IF ABS (x) > MaxLongIntLR THEN Value := TruncLR (x); IF ABS (x - Value) >= Point5LR THEN IF x >= ZeroLR THEN RETURN Value+OneLR ELSE RETURN Value-OneLR END; ELSE (* ABS (x - TruncDble (x)) < 0.5 *) RETURN Value END ELSE (* ABS (x) <= MaxLongInt *) IF x > ZeroLR THEN RETURN LFLOAT (TRUNC (x + Point5LR)) ELSE RETURN LFLOAT (TRUNC (x - Point5LR)) END END END RoundLR; PROCEDURE Trunc (x : LONGREAL) : LONGINT; BEGIN IF ABS ( x ) <= MaxLongIntLR THEN RETURN TRUNC (x); ELSE WriteLR (x, 20, 2, 2); Error ('ERROR Trunc in NUM: Value too large for a LONGINT'); END END Trunc; PROCEDURE Round (x : LONGREAL) : LONGINT; BEGIN IF ABS ( x ) <= MaxLongIntLR THEN IF x > ZeroLR THEN RETURN TRUNC (x + Point5LR) ELSE RETURN TRUNC (x - Point5LR) END ELSE Error ('ERROR Round in NUM: Value too large for a LONGINT'); END END Round; PROCEDURE ReadLI ( VAR i : LONGINT ); BEGIN ReadInt (i) END ReadLI; PROCEDURE ReadLC ( VAR c : LONGCARD ); BEGIN ReadCard (c) END ReadLC; PROCEDURE ReadLR ( VAR x : LONGREAL ); BEGIN ReadLongReal (x); END ReadLR; PROCEDURE WriteLI (i : LONGINT; c : CARDINAL ); BEGIN WriteInt (i, c) END WriteLI; PROCEDURE WriteLC (i : LONGCARD; c : CARDINAL ); BEGIN WriteCard (i, c) END WriteLC; PROCEDURE WriteLR (x : LONGREAL; I, J, K : CARDINAL); VAR spaces : CARDINAL; BEGIN (*IF (J > 0) THEN Cor := 0.5 * Power (10.0, - FLOAT (J)); IF (x < 0.0) THEN Cor := - Cor END (* Note: on fait cette correction si LongRealToStr tronque. ??? *) (********** VERIFIER **************) (* Il n'y a pas de correction a faire sur les SUN *) ELSE Cor := 0.0 END; *) (***** patch parce qu'il y a un bug dans LRealToStr *******) IF ABS (x) > 1.0E50 THEN LRealToStr (x (*+ Cor*), -INTEGER(J), Str); ELSE LRealToStr (x (*+ Cor*), J, Str); END; IF (Length (Str) <= I) THEN IF x=0.0 THEN WriteFormatted (Str, I); RETURN; END; (* Verifier s'il y a suffisamment de chiffres significatifs. *) c := 0; Dot := 1; WHILE ((Str [c] < '1') OR (Str [c] > '9')) AND (c < Length (Str)) DO IF Str [c] = '.' THEN Dot := 0 END; INC (c) END; IF x = TruncLR(x) THEN Dot := 0; END; (* c'est un entier *) IF (Length (Str) - c - Dot >= K) THEN spaces := I - Length (Str); FOR c := 1 TO spaces DO Write (' '); END; FOR c := 0 TO Length(Str)-1 DO Write (Str [c]);END; RETURN END END; (* On doit utiliser la notation scientifique. *) LRealToStr (x, - VAL (INTEGER, J), Str); IF I<=Length(Str) THEN FOR c := 0 TO Length(Str)-1 DO Write(Str[c]) END; ELSE spaces := I - Length (Str); FOR c := 1 TO spaces DO Write(' ');END; FOR c := 0 TO Length(Str)-1 DO Write (Str [c]);END; END; END WriteLR; BEGIN ZeroLR := StringToLR ("0.0"); Point5LR := StringToLR ("0.5"); OneLR := StringToLR ("1.0"); TwoLR := StringToLR ("2.0"); ThreeLR := StringToLR ("3.0"); FourLR := StringToLR ("4.0"); RcpLn2 := OneLR / Ln (TwoLR); Pi := StringToLR ("3.1415926535897932"); e := StringToLR ("2.7182818284590452"); TwoExp [0] := OneLR; FOR i := 1 TO MaxTwoExp DO TwoExp [i] := TwoLR * TwoExp [i-1] END; MaxLongIntLR := StringToLR ("2147483647.0"); MaxIntLongReal := TwoExp [t] - OneLR; MaxLongReal := StringToLR ("8.9E306"); hLR := OneLR; FOR i := 1 TO t DO hLR := hLR/TwoLR END; jLR := hLR + hLR/TwoLR; chLR := OneLR - hLR; END NUM.