IMPLEMENTATION MODULE SUPARITH; FROM MYINOUT IMPORT WriteLn, WriteString; FROM SUPINT IMPORT Base, LongBase, Exposant, SuperInteger, InfoSuperInteger, TwoExp, IncreaseLen, CopyBlock, Create, Affect, Exchange, NumBits; FROM SUPCONV IMPORT IntToSup, SupToInt, SupToLongReal, RoundLongRealToSup; FROM NUM IMPORT OneLR, TwoLR, t, Sqrt, Pow, Trunc, Lg, Round, LIToLR, IToLR, MaxLongReal, ZeroLR; CONST BaseMoinsUn = Base - 1; (* Valeur de la base moins un. *) VAR NbBitsMaxLR : INTEGER; (* nombre de bits necessaires pour *) (* representer le plus grand reel *) (* possible en representation entiere *) Tampon1, Tampon2, (* Utilisees pour memoriser la valeur *) (* d'un pointeur, elles ne sont pas *) (* crees. *) P, R, S, T, U, V, W, X, Y, Z, PP, RR, SS, TT, UU, VV, WW, XX, YY, ZZ : SuperInteger; (* SuperInteger crees servant a des *) (* calculs internes. *) hLR, jLR, chLR : LONGREAL; (* Util. dans DivideLR et DivideRacLR *) l : CARDINAL; (* utilise comme variable de boucle *) (* dans l'initialisation de hLR *) (******************************************************************************) PROCEDURE ShiftLeft ( A : SuperInteger; k : CARDINAL; B : SuperInteger ); BEGIN WITH B^ DO IF A^.Digit^[A^.SigPos] = 0 THEN (* A valait deja 0, donc B = 0. *) SigPos := 0; Digit^[0] := 0; Negative := FALSE ELSE Negative := A^.Negative; SigPos := A^.SigPos + k; Affect (A, R); (* Pour le cas ou A = B *) (* pour ne pas perdre la valeur *) IF TwoExp[Power] <= SigPos THEN IncreaseLen (B) END; CopyBlock (R^.Digit^, Digit^, SigPos, INTEGER (k) ) END END END ShiftLeft; (******************************************************************************) PROCEDURE ShiftRight ( A : SuperInteger; k : CARDINAL; B : SuperInteger ); BEGIN WITH B^ DO IF A^.SigPos >= k THEN (* On n'obtient pas zero. *) Negative := A^.Negative; SigPos := A^.SigPos - k; Affect (A, R); (* Pour le cas ou A = B *) (* pour ne pas perdre la valeur *) IF TwoExp[Power] <= SigPos THEN IncreaseLen (B) END; CopyBlock (R^.Digit^, Digit^, SigPos, - INTEGER (k) ) ELSE (* On obtient zero. *) SigPos := 0; Negative := FALSE; Digit^[0] := 0 END END END ShiftRight; (******************************************************************************) PROCEDURE CutLeft ( A : SuperInteger; k : CARDINAL; B : SuperInteger ); BEGIN WITH B^ DO IF A^.SigPos >= k THEN SigPos := k-1; WHILE (A^.Digit^[SigPos] = 0) AND (SigPos # 0) DO DEC (SigPos) END ELSE SigPos := A^.SigPos END; IF A <> B THEN Negative := A^.Negative; IF TwoExp[Power] <= SigPos THEN IncreaseLen (B) END; CopyBlock (A^.Digit^, Digit^, SigPos, 0) END END END CutLeft; (******************************************************************************) PROCEDURE Log2 ( A : SuperInteger ) : LONGREAL; VAR equal : BOOLEAN; BEGIN WITH A^ DO IntToSup (1, XX); ShiftLeft (XX, SigPos+1, XX); RETURN (Lg (DivideLR (A, XX, equal)) + IToLR (Exposant * (SigPos + 1))) END END Log2; (******************************************************************************) PROCEDURE Even ( A : SuperInteger ) : BOOLEAN; BEGIN RETURN NOT ( ODD ( A^.Digit^[0] ) ) END Even; (******************************************************************************) PROCEDURE Positive ( A : SuperInteger ) : BOOLEAN; BEGIN WITH A^ DO IF Negative OR ( Digit^[SigPos] = 0 ) THEN RETURN FALSE ELSE RETURN TRUE END END END Positive; (******************************************************************************) PROCEDURE Negative ( A : SuperInteger ) : BOOLEAN; BEGIN WITH A^ DO IF Negative AND ( Digit^[SigPos] # 0 ) THEN RETURN TRUE ELSE RETURN FALSE END END END Negative; (******************************************************************************) PROCEDURE Zero ( A : SuperInteger ) : BOOLEAN; BEGIN WITH A^ DO IF Digit^[SigPos] = 0 THEN RETURN TRUE END END; RETURN FALSE END Zero; (******************************************************************************) PROCEDURE PosDif ( VAR AA, AB : ARRAY OF LONGINT; I : CARDINAL ) : CARDINAL; (* Procedure non disponible a l'usager. Pour deux SuperInteger ayant le *) (* meme nb. de digits, retourne la premiere Position a partir de la gauche *) (* ou les digits sont differents. *) BEGIN WHILE ((AA [I] = AB [I]) AND (I > 0)) DO DEC (I) END; RETURN I END PosDif; (******************************************************************************) PROCEDURE AbsSizeOrd ( A, B : SuperInteger ) : INTEGER; VAR I : CARDINAL; BEGIN WITH A^ DO IF SigPos < B^.SigPos THEN RETURN 1 ELSIF SigPos > B^.SigPos THEN RETURN -1 END; I := PosDif(Digit^, B^.Digit^, SigPos); IF Digit^[I] < B^.Digit^[I] THEN RETURN 1 ELSIF Digit^[I] > B^.Digit^[I] THEN RETURN -1 END END; RETURN 0 END AbsSizeOrd; (******************************************************************************) PROCEDURE Equal ( A, B : SuperInteger ) : BOOLEAN; VAR I : CARDINAL; BEGIN WITH A^ DO IF Negative = B^.Negative THEN (* Meme signe *) IF SigPos <> B^.SigPos THEN (* Nombre different de digits *) RETURN FALSE ELSE (* Meme nombre de digits *) I := PosDif ( Digit^, B^.Digit^, SigPos ); IF Digit^[I] = B^.Digit^[I] THEN RETURN TRUE ELSE RETURN FALSE END END ELSE (* Signe different *) RETURN FALSE END END END Equal; (******************************************************************************) PROCEDURE Greater ( A, B : SuperInteger ) : BOOLEAN; VAR I : CARDINAL; BEGIN WITH A^ DO IF Negative = B^.Negative THEN (* Meme signe *) IF SigPos = B^.SigPos THEN (* A et B ont le meme nombre de digits*) I := PosDif ( Digit^, B^.Digit^, SigPos ); IF Negative THEN IF Digit^[I] >= B^.Digit^[I] THEN RETURN FALSE END ELSE IF Digit^[I] <= B^.Digit^[I] THEN RETURN FALSE END END ELSE IF (SigPos < B^.SigPos) AND (NOT Negative) THEN (* B a plus de digits que A et ils sont positifs. *) RETURN FALSE END; IF (SigPos > B^.SigPos) AND Negative THEN (* A a plus de digits que B et ils sont negatifs. *) RETURN FALSE END; END ELSE (* Signe different *) IF Negative THEN RETURN FALSE END END END; RETURN TRUE END Greater; (******************************************************************************) PROCEDURE Smaller ( A, B : SuperInteger ) : BOOLEAN; VAR I : CARDINAL; BEGIN WITH A^ DO IF Negative = B^.Negative THEN (* Meme signe *) IF SigPos = B^.SigPos THEN (* ont le meme nombre de digits *) I := PosDif ( Digit^, B^.Digit^, SigPos ); IF Negative THEN IF Digit^[I] > B^.Digit^[I] THEN RETURN TRUE END ELSE IF Digit^[I] < B^.Digit^[I] THEN RETURN TRUE END END ELSE IF (SigPos < B^.SigPos) AND (NOT Negative) THEN (* B a plus de digits que A et ils sont positifs. *) RETURN TRUE END; IF (SigPos > B^.SigPos) AND Negative THEN (* A a plus de digits que B et ils sont negatifs. *) RETURN TRUE END; END ELSE (* Signe different *) IF Negative THEN RETURN TRUE END END END; RETURN FALSE END Smaller; (******************************************************************************) PROCEDURE ChangeSign ( A, B : SuperInteger ); BEGIN IF A # B THEN Affect (A, B) END; B^.Negative := NOT A^.Negative; IF B^.Digit^[B^.SigPos] = 0 THEN B^.Negative := FALSE END; (* ajout de JT pour eviter -zero *) END ChangeSign; (******************************************************************************) PROCEDURE AbsVal ( A, B : SuperInteger ); BEGIN IF A # B THEN Affect ( A, B ) END; B^.Negative := FALSE END AbsVal; (******************************************************************************) (* Procedure non disponible a l'usager. Soustrait AMax de AMin ou AMax >= AMin*) (* et met le resultat dans Tampon. *) PROCEDURE InsideSub ( VAR AMax, AMin, ATampon : ARRAY OF LONGINT; MaxSigPos, MinSigPos : CARDINAL ); VAR I : CARDINAL; BEGIN WITH Tampon1^ DO FOR I := 0 TO MinSigPos DO ATampon[I] := AMax[I] - AMin[I]; END; FOR I := MinSigPos + 1 TO MaxSigPos DO ATampon[I] := AMax[I]; END; FOR I := 0 TO MaxSigPos DO IF ATampon[I] < 0 THEN DEC ( ATampon[I+1] ); INC ( ATampon[I], LongBase ) END END; WHILE ( SigPos # 0 ) AND ( ATampon[SigPos] = 0 ) DO DEC ( SigPos ) END; IF ATampon[SigPos] = 0 THEN Negative := FALSE END END END InsideSub; (******************************************************************************) (* Procedure non disponible a l'usager. Additionne AMax et AMin, AMax doit *) (* avoir un nombre de digit >= a AMin et met le resultat dans Tampon. *) PROCEDURE InsideAdd ( VAR AMax, AMin, ATampon : ARRAY OF LONGINT; MaxPos, MinPos : CARDINAL ); VAR I : CARDINAL; BEGIN WITH Tampon1^ DO ATampon[SigPos] := 0; (* Fait l'addition sans retenue *) FOR I := 0 TO MinPos DO ATampon[I] := AMax[I] + AMin[I] END; (* Copie le reste chiffre le plus long*) (* dans la reponse *) FOR I := MinPos + 1 TO MaxPos DO ATampon[I] := AMax[I] END; FOR I := 0 TO MaxPos DO IF ATampon[I] >= LongBase THEN INC ( ATampon[I+1] ); DEC ( ATampon[I], LongBase ) END END; IF ( SigPos # 0 ) AND ( ATampon[SigPos] = 0 ) THEN DEC ( SigPos ) END END END InsideAdd; (******************************************************************************) PROCEDURE Add ( A, B, C : SuperInteger ); VAR I : CARDINAL; Min, Max : SuperInteger; BEGIN IF ( A # C ) AND ( B # C ) THEN (* On utilise pas de Tampon *) Tampon1 := C ELSE (* On utilise un Tampon *) Tampon1 := X END; WITH Tampon1^ DO IF A^.Negative <> B^.Negative THEN(* signe different on soutrait *) IF A^.SigPos > B^.SigPos THEN (* Trouve plus long et plus court *) Max := A; Min := B ELSE IF A^.SigPos < B^.SigPos THEN Max := B; Min := A ELSE I := PosDif ( A^.Digit^, B^.Digit^, A^.SigPos ); IF A^.Digit^[I] > B^.Digit^[I] THEN Max := A; Min := B ELSE Max := B; Min := A END END END; (* Ajuste le signe, et reserve de *) Negative := Max^.Negative; (* l'espace memoire si necessaire *) SigPos := Max^.SigPos; IF TwoExp[Power] <= SigPos THEN IncreaseLen ( Tampon1 ) END; InsideSub ( Max^.Digit^, Min^.Digit^, Digit^, Max^.SigPos, Min^.SigPos ) ELSE (* Meme signe on additionne *) IF A^.SigPos > B^.SigPos THEN Max := A; Min := B ELSE Max := B; Min := A END; (* Verifie si la longueur est suffi- *) (* sante *) Negative := Max^.Negative; SigPos := Max^.SigPos + 1; IF TwoExp[Power] <= SigPos THEN IncreaseLen ( Tampon1 ) END; InsideAdd ( Max^.Digit^, Min^.Digit^, Digit^, Max^.SigPos, Min^.SigPos ) END END; IF Tampon1 # C THEN (* On a utilise un Tampon, alors on *) Affect ( X, C ) (* copy le Tampon dans C *) END; IF C^.Digit^[C^.SigPos] = 0 THEN C^.Negative := FALSE END; (* ajout de JT pour eviter -zero *) END Add; (******************************************************************************) PROCEDURE Subtract ( A, B, C : SuperInteger ); VAR I : CARDINAL; Max, Min : SuperInteger; BEGIN IF ( A # C ) AND ( B # C ) THEN Tampon1 := C ELSE Tampon1 := X END; WITH Tampon1^ DO IF A^.Negative = B^.Negative THEN (* Meme signe on soustrait *) IF A^.SigPos > B^.SigPos THEN (* Trouve plus long et plus court *) Max := A; Min := B ELSE IF A^.SigPos < B^.SigPos THEN Max := B; Min := A ELSE I := PosDif ( A^.Digit^, B^.Digit^, A^.SigPos ); IF A^.Digit^[I] > B^.Digit^[I] THEN Max := A; Min := B ELSE Max := B; Min := A END END END; IF Max = A THEN (* Ajuste le signe *) Negative := A^.Negative ELSE Negative := NOT A^.Negative END; SigPos := Max^.SigPos; (* Alloue l'espace memoire si neces- *) (* saire *) IF TwoExp[Power] <= SigPos THEN IncreaseLen ( Tampon1 ) END; InsideSub ( Max^.Digit^, Min^.Digit^, Digit^, Max^.SigPos, Min^.SigPos ) ELSE (* Signe different on additionne *) IF A^.SigPos > B^.SigPos THEN Max := A; Min := B ELSE Max := B; Min := A END; (* Verifie si la longueur est suffi- *) (* sante *) Negative := A^.Negative; SigPos := Max^.SigPos + 1; IF TwoExp[Power] <= SigPos THEN IncreaseLen ( Tampon1 ) END; InsideAdd ( Max^.Digit^, Min^.Digit^, Digit^, Max^.SigPos, Min^.SigPos ) END END; IF Tampon1 # C THEN (* On a utilise un Tampon alors on *) Affect ( X, C ) (* copy le Tampon dans C *) END; IF C^.Digit^[C^.SigPos] = 0 THEN C^.Negative := FALSE END; (* ajout de JT pour eviter -zero *) END Subtract; (******************************************************************************) (* + ajmal 18-12-96 *) PROCEDURE InsideMultiply ( VAR AA, AB, ATampon : ARRAY OF LONGINT; VAR A, B : SuperInteger ); VAR I, J : CARDINAL; T, K : LONGINT; BEGIN (* Algorithme, voir Knuth p. 253, 254 *) WITH A^ DO (* Initialise *) FOR I := 0 TO SigPos DO ATampon[I] := 0 END; (* Multiplie *) FOR J := 0 TO B^.SigPos DO K := 0; FOR I := 0 TO SigPos DO T := AB[J] * AA[I] + ATampon[I+J] + K; K := T DIV LongBase; ATampon [I+J] := T - K * LongBase; END; ATampon [SigPos+J+1] := K END END; WITH Tampon1^ DO (* Enleve les 0 a gauche *) WHILE ( SigPos # 0 ) AND ( ATampon [SigPos] = 0 ) DO DEC ( SigPos ) END END END InsideMultiply; PROCEDURE Multiply ( A, B, C : SuperInteger ); (* deplacee a l'exterieur ajmal 18 Dec 1996 PROCEDURE InsideMultiply (VAR AA, AB, ATampon : ARRAY OF LONGINT); VAR I, J : CARDINAL; T, K : LONGINT; BEGIN * Algorithme, voir Knuth p. 253, 254 * WITH A^ DO * Initialise * FOR I := 0 TO SigPos DO ATampon[I] := 0 END; * Multiplie * FOR J := 0 TO B^.SigPos DO K := 0; FOR I := 0 TO SigPos DO T := AB[J] * AA[I] + ATampon[I+J] + K; K := T DIV LongBase; ATampon [I+J] := T - K * LongBase; END; ATampon [SigPos+J+1] := K END END; WITH Tampon1^ DO * Enleve les 0 a gauche * WHILE ( SigPos # 0 ) AND ( ATampon [SigPos] = 0 ) DO DEC ( SigPos ) END END END InsideMultiply; end ajmal *) BEGIN IF ( A # C ) AND ( B # C ) THEN Tampon1 := C ELSE Tampon1 := X END; WITH Tampon1^ DO (* Alloue l'espace memoire suffisante *) SigPos := A^.SigPos + B^.SigPos + 1; IF TwoExp[Power] <= SigPos THEN IncreaseLen ( Tampon1 ) END; (* * ajmal 18-12-96 *) (* InsideMultiply ( A^.Digit^, B^.Digit^, Digit^ ); *) InsideMultiply ( A^.Digit^, B^.Digit^, Digit^, A, B ); IF A^.Negative <> B^.Negative THEN Negative := TRUE ELSE Negative := FALSE END; END; IF Tampon1 # C THEN Affect ( X, C ) END; IF C^.Digit^[C^.SigPos] = 0 THEN C^.Negative := FALSE END; (* ajout de JT pour eviter -zero *) END Multiply; (******************************************************************************) PROCEDURE MultiplyIntSup ( I : LONGINT; B, C : SuperInteger ); BEGIN IntToSup ( I, U ); Multiply ( U, B, C ) END MultiplyIntSup; (******************************************************************************) (* + ajmal 18-12-96 *) PROCEDURE OneDigitDivision ( VAR AA, AQ : ARRAY OF LONGINT; VAR r, K, b : LONGINT; VAR m, J : CARDINAL ); (* Utilisee lorsque le diviseur a 1 digit et donc que le reste aura un *) (* digit. Cette procedure divise le SuperInteger contenu dans AA par b *) (* et met le quotient dans AQ et le reste dans r. *) BEGIN r := 0; FOR J := m TO 0 BY -1 DO K := r * LongBase + AA[J]; AQ[J] := K DIV b; r := K - AQ[J] * b END END OneDigitDivision; PROCEDURE Multiplybyd ( VAR AA, AB : ARRAY OF LONGINT; SigPos : CARDINAL; VAR K, d, T : LONGINT; VAR I : CARDINAL ); (* Utilisee lorsque le diviseur a plusieur digits pour normaliser, voir *) (* algorithme Knuth. Cette procedure multiplie le superinteger contenu *) (* dans AA par d et met le resultat dans AB. *) BEGIN K := 0; FOR I := 0 TO SigPos DO T := d * AA[I] + K; K := T DIV LongBase; AB[I] := T - K * LongBase END; AB[SigPos+1] := K END Multiplybyd; PROCEDURE ManyDigitsDivision ( VAR AX, AY, AZ : ARRAY OF LONGINT; VAR SigPos : CARDINAL; VAR J, n, I : CARDINAL; VAR b, K, q, T, r, d : LONGINT ); (* Utilisee lorsque le diviseur a plusieur digits pour diviser, voir *) (* algorithme Knuth. Cette procedure divise le superinteger contenu dans*) (* AX par AY et met le resultat dans AZ. AX contiendra le reste. *) BEGIN FOR J := SigPos TO n+1 BY -1 DO b := AY[n]; K := AX[J] * LongBase + AX[J-1]; IF AX[J] = b THEN (* Evalue le digit du quotient *) q := BaseMoinsUn; ELSE q := K DIV b END; WHILE q * AY[n-1] > (K - q * b) * LongBase + AX [J-2] DO DEC (q) END; K := 0; (* Multiplie le digit par le diviseur *) FOR I := J-n-1 TO J DO (* et soustrait du dividende *) T := q * AY[I+n+1-J] + K; K := T DIV LongBase; AX[I] := AX[I] - ( T - K * LongBase ); IF AX[I] < 0 THEN AX[I] := AX[I] + LongBase; DEC ( AX[I+1] ) END END; IF AX[J+1] = -1 THEN (* Cette etape se produit avec une probabilite de 2/Base. *) (* Le digit q est trop grand de 1, on le decremente *) (* et on additionne le diviseur au au dividende. *) DEC (q); FOR I := J-n-1 TO J DO AX[I] := AY[I+n+1-J] + AX[I]; IF AX[I] >= LongBase THEN INC ( AX[I+1] ); AX[I] := AX[I] - LongBase END END END; AZ [J-n-1] := q END; SigPos := n; r := 0; (* Le quotient est ok, on denormalise *) (* le reste. *) FOR J := n TO 0 BY -1 DO K := r * LongBase + AX[J]; AX[J] := K DIV d; r := K - AX[J] * d END; (* Enleve les zeros a gauche du reste *) WHILE (SigPos # 0) AND (AX [SigPos] = 0) DO DEC (SigPos) END END ManyDigitsDivision; PROCEDURE Divide (A, B, Q, R : SuperInteger); VAR n, (* Nb digits-1 du diviseur et/ou reste*) m, (* Nb digits-1 du quotient *) I, J : CARDINAL; (* Indices de boucle *) d, (* Facteur de normalisation *) b, (* Digit du diviseur *) q, (* Digit du quotient *) r, (* Reste partiel *) T, K : LONGINT; (* Variables temporaires de calcul *) QNegative, (* Signe du quotient *) RNegative : BOOLEAN; (* Signe du reste *) (* On utilise des Tampons lorsque le diviseur a plusieurs digits *) (* Y contient le diviseur normalise *) (* X contient le dividende normalise si A=R *) (* deplacee a l'exterieur ajmal 18 Dec 1996 PROCEDURE OneDigitDivision ( VAR AA, AQ : ARRAY OF LONGINT ); * Utilisee lorsque le diviseur a 1 digit et donc que le reste aura un * * digit. Cette procedure divise le SuperInteger contenu dans AA par b * * et met le quotient dans AQ et le reste dans r. * BEGIN r := 0; FOR J := m TO 0 BY -1 DO K := r * LongBase + AA[J]; AQ[J] := K DIV b; r := K - AQ[J] * b END END OneDigitDivision; end ajmal *) (* deplacee a l'exterieur ajmal 18 Dec 1996 PROCEDURE Multiplybyd ( VAR AA, AB : ARRAY OF LONGINT; SigPos : CARDINAL ); (* Utilisee lorsque le diviseur a plusieur digits pour normaliser, voir *) (* algorithme Knuth. Cette procedure multiplie le superinteger contenu *) (* dans AA par d et met le resultat dans AB. *) BEGIN K := 0; FOR I := 0 TO SigPos DO T := d * AA[I] + K; K := T DIV LongBase; AB[I] := T - K * LongBase END; AB[SigPos+1] := K END Multiplybyd; end ajmal *) (* deplacee a l'exterieur ajmal 18 Dec 1996 PROCEDURE ManyDigitsDivision ( VAR AX, AY, AZ : ARRAY OF LONGINT; VAR SigPos : CARDINAL ); * Utilisee lorsque le diviseur a plusieur digits pour diviser, voir * * algorithme Knuth. Cette procedure divise le superinteger contenu dans* * AX par AY et met le resultat dans AZ. AX contiendra le reste. * BEGIN FOR J := SigPos TO n+1 BY -1 DO b := AY[n]; K := AX[J] * LongBase + AX[J-1]; IF AX[J] = b THEN * Evalue le digit du quotient * q := BaseMoinsUn; ELSE q := K DIV b END; WHILE q * AY[n-1] > (K - q * b) * LongBase + AX [J-2] DO DEC (q) END; K := 0; * Multiplie le digit par le diviseur * FOR I := J-n-1 TO J DO * et soustrait du dividende * T := q * AY[I+n+1-J] + K; K := T DIV LongBase; AX[I] := AX[I] - ( T - K * LongBase ); IF AX[I] < 0 THEN AX[I] := AX[I] + LongBase; DEC ( AX[I+1] ) END END; IF AX[J+1] = -1 THEN * Cette etape se produit avec une probabilite de 2/Base. * * Le digit q est trop grand de 1, on le decremente * * et on additionne le diviseur au au dividende. * DEC (q); FOR I := J-n-1 TO J DO AX[I] := AY[I+n+1-J] + AX[I]; IF AX[I] >= LongBase THEN INC ( AX[I+1] ); AX[I] := AX[I] - LongBase END END END; AZ [J-n-1] := q END; SigPos := n; r := 0; * Le quotient est ok, on denormalise * * le reste. * FOR J := n TO 0 BY -1 DO K := r * LongBase + AX[J]; AX[J] := K DIV d; r := K - AX[J] * d END; * Enleve les zeros a gauche du reste * WHILE (SigPos # 0) AND (AX [SigPos] = 0) DO DEC (SigPos) END END ManyDigitsDivision; end ajmal *) BEGIN IF A^.SigPos < B^.SigPos THEN (* Le diviseur a un nombre de digits *) (* plus grand que le dividende *) Affect ( A, R ); (* Met le dividende dans le reste *) WITH Q^ DO (* Met le quotient a 0 *) Negative := FALSE; SigPos := 0; Digit^[0] := 0 END ELSE (* Le diviseur a un nombre de digits *) (* <= a celui du dividende *) WITH B^ DO n := SigPos; (* Nb digits-1 diviseur et/ou reste *) m := A^.SigPos - SigPos; (* Nb digits du quotient *) RNegative := A^.Negative; QNegative := NOT (Negative = RNegative); END; IF n = 0 THEN (* Le diviseur a 1 digit *) b := B^.Digit^[0]; IF b = 0 THEN (* Le diviseur est 0 *) WriteLn; WriteString ('ERROR Divide, division by 0'); WriteLn; HALT ELSE (* Le diviseur a 1 digit et est *) (* different de zero *) WITH Q^ DO SigPos := m; (* Ajuste l'espace memoire du quotient*) IF TwoExp[Power] <= SigPos THEN IncreaseLen ( Q ) END; (* * ajmal 18-12-96 *) (* OneDigitDivision (A^.Digit^, Digit^) *) OneDigitDivision (A^.Digit^, Digit^, r, K, b, m, J ) END; WITH R^ DO Digit^[0] := r; SigPos := 0 END END ELSE (* Le diviseur a plus d'un digit *) (* Determine le facteur de nromalisa- *) (* tion. voir Knuth p.257 et 258 *) d := LongBase DIV ( B^.Digit^[n] + 1 ); WITH Y^ DO SigPos := B^.SigPos + 1; (* Y possede l'espace requis *) IF TwoExp[Power] <= SigPos THEN IncreaseLen ( Y ) END; (* Normalise le diviseur B et le met *) (* dans le Tampon Y *) (* * ajmal 18-12-96 *) (* Multiplybyd ( B^.Digit^, Digit^, n ) *) Multiplybyd ( B^.Digit^, Digit^, n, K, d, T, I ) END; IF A # R THEN Tampon1 := R ELSE Tampon1 := X END; WITH Tampon1^ DO (* Normalise le dividende A et le met *) (* dans le Tampon1 ( reste ) *) SigPos := A^.SigPos + 2; (* Buf possede l'espace requis *) IF TwoExp[Power] <= SigPos THEN IncreaseLen ( Tampon1 ) END; Digit^[SigPos] := 0; DEC ( SigPos ); (* Multiplie le dividende par d *) (* * ajmal 18-12-96 *) (* Multiplybyd ( A^.Digit^, Digit^, A^.SigPos ); *) Multiplybyd ( A^.Digit^, Digit^, A^.SigPos, K, d, T, I ); END; WITH Q^ DO SigPos := m; (* Ajuste l'espace memoire du quotient*) IF TwoExp[Power] <= SigPos THEN IncreaseLen ( Q ) END; (* * ajmal 18-12-96 *) (* ManyDigitsDivision (Tampon1^.Digit^, Y^.Digit^, Digit^, Tampon1^.SigPos); *) ManyDigitsDivision (Tampon1^.Digit^, Y^.Digit^, Digit^, Tampon1^.SigPos, J, n, I, b, K, q, T, r, d ); END; Affect (Tampon1, R); END; WITH R^ DO (* Etablit le signe du reste *) IF (SigPos = 0) AND (Digit^[SigPos] = 0) THEN Negative := FALSE ELSE Negative := RNegative END END; WITH Q^ DO (* Ote le zero a gauche du quotient *) IF ( SigPos # 0 ) AND ( Digit^[SigPos] = 0 ) THEN DEC ( SigPos ) END; (* Etablit le signe du quotient *) IF ( SigPos = 0 ) AND ( Digit^[0] = 0 ) THEN Negative := FALSE ELSE Negative := QNegative END END END; IF Q^.Digit^[Q^.SigPos] = 0 THEN Q^.Negative := FALSE END; (* ajout de JT pour eviter -zero *) IF R^.Digit^[R^.SigPos] = 0 THEN R^.Negative := FALSE END; (* ajout de JT pour eviter -zero *) END Divide; (******************************************************************************) PROCEDURE DivideSupInt ( A : SuperInteger; B : LONGINT; Q : SuperInteger; VAR R : LONGINT ); BEGIN IntToSup ( B, U ); Divide ( A, U, Q, U ); SupToInt ( U, R ) END DivideSupInt; (******************************************************************************) PROCEDURE SucLR (x : LONGREAL) : LONGREAL; (* Retourne le LONGREAL suivant immediatement x *) BEGIN IF x < ZeroLR THEN RETURN -PredLR(-x) END; (* IF x = ZeroLR THEN RETURN d END; d a arranger *) RETURN x + jLR*x; END SucLR; PROCEDURE PredLR (x : LONGREAL) : LONGREAL; (* Retourne le LONGREAL precedant immediatement x *) BEGIN IF x < ZeroLR THEN RETURN -SucLR(-x) END; (* IF x = ZeroLR THEN RETURN -d END; d a arranger *) RETURN chLR*x; END PredLR; (******************************************************************************) PROCEDURE DivideLR (S, T : SuperInteger; VAR equal : BOOLEAN) : LONGREAL; VAR Res : LONGREAL; cte : LONGREAL; Itemp : INTEGER; BEGIN equal := TRUE; IF Zero(S) THEN RETURN ZeroLR END; IF (Negative(S) AND NOT Negative(T)) THEN ChangeSign(S, S); Res := - DivideLR (S, T, equal); IF NOT equal THEN Res := SucLR(Res) END; ChangeSign(S, S); RETURN Res; END; IF (Negative(T) AND NOT Negative(S)) THEN ChangeSign(T, T); Res := - DivideLR (S, T, equal); IF NOT equal THEN Res := SucLR(Res) END; ChangeSign(T, T); RETURN Res; END; IntToSup(2, P); Itemp := VAL(INTEGER,t) - VAL(INTEGER, NumBits(S)) + VAL(INTEGER, NumBits(T)); IF Itemp <= 0 THEN Itemp := -Itemp; PowerSupInt(P, Itemp, P); Affect (T, R); Multiply (R, P, R); Divide (S, R, R, V); IF NOT Zero(V) THEN equal := FALSE END; SupToLongReal (P, cte); IF (NumBits(R) > t) THEN DivideSupInt (R, 2, R, Itemp); IF Itemp # 0 THEN equal := FALSE END; cte := cte * TwoLR; END; IF (VAL(INTEGER,NumBits(R)) + Round(Lg(cte))) > NbBitsMaxLR THEN (* resultat non representable en reel *) WriteLn; WriteString ("ERROR in DivideLR: Value too large for a LONGREAL"); WriteLn; HALT END; SupToLongReal(R, Res); Res := Res * cte; ELSE PowerSupInt(P, Itemp, P); Affect (S, R); Multiply (R, P, R); Divide (R, T, R, V); IF NOT Zero(V) THEN equal := FALSE END; SupToLongReal (P, cte); IF (NumBits(R) > t) THEN DivideSupInt (R, 2, R, Itemp); IF Itemp # 0 THEN equal := FALSE END; cte := cte / TwoLR; END; SupToLongReal(R, Res); Res := Res / cte; END; RETURN Res; END DivideLR; (******************************************************************************) PROCEDURE DivideRacLR (S, T : SuperInteger) : LONGREAL; VAR Res, ResAnc, ResSuiv, cte : LONGREAL; Itemp : INTEGER; equal : BOOLEAN; BEGIN IF Zero(S) THEN RETURN ZeroLR END; IF Negative(S) AND Negative(T) THEN ChangeSign (S, S); ChangeSign (T, T); Res := DivideRacLR (S, T); ChangeSign (S, S); ChangeSign (T, T); RETURN Res; END; Itemp := VAL(INTEGER,NumBits(S)) - VAL(INTEGER,NumBits(T)) - NbBitsMaxLR + 1; IF Itemp < 0 THEN Itemp := 0 END; IF ODD(Itemp) THEN Itemp := Itemp + 1 END; IntToSup (2, W); PowerSupInt (W, Itemp, W); Quotient (S, W, W); Res := Sqrt (DivideLR (W, T, equal)); Res := Res * Pow(2.0, LIToLR(Itemp)/2.0); equal := FALSE; Itemp := t - Trunc (Lg (Res) + 1.0); IF Itemp <= 0 THEN Itemp := 0 END; cte := Pow (2.0, LIToLR (Itemp)); RoundLongRealToSup (cte, W); Multiply (W, W, W); Multiply (S, W, W); ResSuiv := Res * cte; ResAnc := Res; RoundLongRealToSup (ResSuiv, P); Multiply (P, P, P); Multiply(T, P, P); IF Smaller(P, W) THEN REPEAT ResAnc := Res; Res := SucLR(Res); ResSuiv := Res * cte; RoundLongRealToSup (ResSuiv, P); Multiply (P, P, P); Multiply(T, P, P); UNTIL NOT Smaller(P, W); IF Equal(P, W) THEN equal := TRUE ELSE Res := ResAnc END; ELSIF Greater(P, W) THEN REPEAT Res := PredLR(Res); ResSuiv := Res * cte; RoundLongRealToSup (ResSuiv, P); Multiply (P, P, P); Multiply(T, P, P); UNTIL NOT Greater(P, W); IF Equal(P, W) THEN equal := TRUE END; ELSE equal := TRUE END; RETURN Res; END DivideRacLR; (******************************************************************************) PROCEDURE Quotient ( A, B, Q : SuperInteger ); BEGIN Divide ( A, B, Q, U ); END Quotient; (******************************************************************************) PROCEDURE Modulo ( A, B, R : SuperInteger ); BEGIN IF B^.Negative THEN (* deuxieme operande < 0 *) WriteLn; WriteString ('ERROR Modulo, Negative second operand'); WriteLn; HALT END; Divide ( A, B, U, R ); IF R^.Negative THEN Add ( B, R, R ) END END Modulo; (******************************************************************************) PROCEDURE PowerSupInt ( A : SuperInteger; I : LONGCARD; C : SuperInteger ); BEGIN IF I < 1 THEN (* L'exposant est 0 *) IntToSup ( 1, C ) ELSE (* L'exposant > 0 *) Affect ( A, Y ); IntToSup ( 1, Z ); WHILE I > 1 DO IF ODD ( I ) THEN Multiply ( Z, Y, Z ) END; Multiply ( Y, Y, Y ); I := I DIV 2 END; Multiply ( Y, Z, C ) END END PowerSupInt; (******************************************************************************) PROCEDURE GCD (A, B, C : SuperInteger); VAR I : CARDINAL; a, b, c, d, q, t, w, v : LONGINT; BEGIN IF ( A # C ) AND ( B # C ) THEN Tampon2 := C ELSE Tampon2 := W END; WITH A^ DO (* Trouve le plus grand et le met *) IF SigPos > B^.SigPos THEN (* dans W, le plus petit dans V *) AbsVal ( A, Tampon2 ); AbsVal ( B, V ) ELSE IF SigPos < B^.SigPos THEN AbsVal ( B, Tampon2 ); AbsVal ( A, V ) ELSE I := PosDif ( Digit^, B^.Digit^, SigPos ); IF Digit^[I] > B^.Digit^[I] THEN AbsVal ( A, Tampon2 ); AbsVal ( B, V ) ELSE IF Digit^[I] < B^.Digit^[I] THEN AbsVal ( B, Tampon2 ); AbsVal ( A, V ) ELSE AbsVal ( A, Tampon2 ); V^.SigPos := 0; V^.Digit^[0] := 0 END END END END END; WITH V^ DO (* Algorithme, voir Knuth p. 329, 330 *) WHILE Digit^[SigPos] # 0 DO IF Tampon2^.SigPos < 2 THEN (* On peut calculer le pgcd en simple *) SupToInt ( Tampon2, w ); (* precision *) SupToInt ( V, v ); WHILE v # 0 DO a := w MOD v; w := v; v := a END; SigPos := 0; Digit^[0] := 0; IntToSup ( w, Tampon2 ) ELSE (* On calcul le pgcd en multiple *) b := 0; (* precision *) IF Tampon2^.SigPos - SigPos < 2 THEN a := 1;c := 0;d := 1; w := Tampon2^.Digit^[Tampon2^.SigPos] * LongBase + Tampon2^.Digit^[Tampon2^.SigPos-1]; IF Tampon2^.SigPos = SigPos THEN v := Digit^[SigPos] * LongBase + Digit^[SigPos-1]; ELSE v := Digit^[SigPos] END; WHILE (( v + c ) # 0) AND (( v + d ) # 0) DO q := ( a + w ) DIV ( v + c ); IF q = ( b + w ) DIV ( v + d ) THEN t := a - q * c; a := c; c := t; t := b - q * d; b := d; d := t; t := w - q * v; w := v; v := t; ELSE v := -c END END END; IF b = 0 THEN Divide ( Tampon2, V, S, T ); Affect ( V, Tampon2 ); Affect ( T, V ) ELSE MultiplyIntSup ( a, Tampon2, T ); MultiplyIntSup ( b, V, Z ); MultiplyIntSup ( c, Tampon2, S ); MultiplyIntSup ( d, V, Y ); Add ( T, Z, Tampon2 ); Add ( S, Y, V) END END END END; IF Tampon2 # C THEN Affect ( W, C ) END; IF C^.Digit^[C^.SigPos] = 0 THEN C^.Negative := FALSE END; (* ajout de JT pour eviter -zero *) END GCD; (******************************************************************************) PROCEDURE Euclide (A, B : SuperInteger; VAR C, D, E, F, G : SuperInteger); BEGIN Affect (A, G); Affect (B, ZZ); IntToSup (1, C); IntToSup (0, D); IntToSup (0, E); IntToSup (1, F); IF Zero (A) THEN Exchange (G, ZZ); Exchange (C, E); Exchange (D, F); RETURN END; WHILE NOT Zero (ZZ) DO Exchange (G, ZZ); Exchange (C, E); Exchange (D, F); Quotient (ZZ, G, XX); ChangeSign (XX, XX); Multiply (XX, G, YY); Add (YY, ZZ, ZZ); Multiply (XX, C, YY); Add (YY, E, E); Multiply (XX, D, YY); Add (YY, F, F); END END Euclide; (******************************************************************************) PROCEDURE MultMod (A, B, M, C : SuperInteger); BEGIN Multiply (A, B, W); Modulo (W, M, C) END MultMod; (******************************************************************************) PROCEDURE PowerMod (A, B, M, C : SuperInteger); BEGIN IF B^.Negative THEN (* L'exposant est negatif, erreur *) WriteLn; WriteString ('ERROR PowerMod, Negative power'); WriteLn; HALT ELSE (* L'exposant est positif *) WITH S^ DO Affect ( B, S ); IntToSup ( 1, T ); IF Digit^[SigPos] = 0 THEN (* L'exposant est 0 *) Modulo ( T, M, C ) ELSE (* L'exposant > 0 *) IntToSup ( 2, R ); Affect ( A, V ); WHILE ( SigPos # 0 ) OR ( Digit^[SigPos] > 1 ) DO IF ODD ( Digit^[0] ) THEN MultMod ( T, V, M, T ) END; MultMod ( V, V, M, V ); Quotient ( S, R, S ) END; MultMod ( V, T, M, C ) END END END; IF C^.Digit^[C^.SigPos] = 0 THEN C^.Negative := FALSE END; (* ajout de JT pour eviter -zero *) END PowerMod; (******************************************************************************) PROCEDURE DivideMod ( A, B, M, C : SuperInteger ) : BOOLEAN; BEGIN Euclide (B, M, SS, TT, UU, VV, WW); Divide (A, WW, UU, VV); IF NOT Zero (VV) THEN RETURN FALSE ELSE Multiply (SS, UU, C); Modulo (C, M, C); RETURN TRUE END END DivideMod; (******************************************************************************) BEGIN (* Creation des tampons necessaires pour les operations mathematiques. *) Create (P, 6); Create (R, 2); Create (S, 6); Create (T, 6); Create (U, 6); Create (V, 6); Create (W, 6); Create (X, 6); Create (Y, 6); Create (Z, 6); Create (PP, 6); Create (RR, 6); Create (SS, 6); Create (TT, 6); Create (UU, 6); Create (VV, 6); Create (WW, 6); Create (XX, 6); Create (YY, 6); Create (ZZ, 6); NbBitsMaxLR := Round (Lg (MaxLongReal)); hLR := OneLR; FOR l := 1 TO t DO hLR := hLR / TwoLR END; jLR := hLR + hLR / TwoLR; chLR := OneLR - hLR; END SUPARITH.