Artifact bb94f111089ecc81ab0e81c0f5fe92b0a47d7cb8d76187cd671f9b1e71848479:
- File
psl-1983/3-1/util/bigbig.red
— part of check-in
[eb17ceb7f6]
at
2020-04-21 19:40:01
on branch master
— Add Reduce 3.0 to the historical section of the archive, and some more
files relating to version sof PSL from the early 1980s. Thanks are due to
Paul McJones and Nelson Beebe for these, as well as to all the original
authors.git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/historical@5328 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 26177) [annotate] [blame] [check-ins using] [more...]
- File
psl-1983/util/bigbig.red
— part of check-in
[eb17ceb7f6]
at
2020-04-21 19:40:01
on branch master
— Add Reduce 3.0 to the historical section of the archive, and some more
files relating to version sof PSL from the early 1980s. Thanks are due to
Paul McJones and Nelson Beebe for these, as well as to all the original
authors.git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/historical@5328 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 26177) [annotate] [blame] [check-ins using]
% BIGBIG.RED - Vector based BIGNUM package with INUM operations % M. L. Griss & B Morrison % 25 June 1982. % % Revision log: % 20 Dec: % MLG, changed TrimBigNUM to TrimBigNum1 in BhardDivide % 14 Dec: % Changed by MLG to put LOAD and IMPORTS in BUILD file % A. C . Norman - adjstments to many routines! % in particular corrections to BHardDivide (case D6 utterly wrong), % and adjustments to BExpt (for performance) and all logical % operators (for treatment of negative inputs); % 31 August 1982: % Copyright (C) 1982, A. C. Norman, B. Morrison, M. Griss % --------------------------------------------------------------- % ----------------------- % A bignum will be a VECTOR of Bigits: (digits in base BigBase): % [BIGPOS b1 ... bn] or [BIGNEG b1 ... bn]. BigZero is thus [BIGPOS] % All numbers are positive, with BIGNEG as 0 element to indicate negatives. Fluid '(BBase!* BBits!* LogicalBits!* WordHi!* WordLow!* Digit2Letter!* FloatHi!* FloatLow!* SysHi!* SysLo!* Carry!* OutputBase!*); % -------------------------------------------------------------------------- % -------------------------------------------------------------------------- % Support functions: % % U, V, V1, V2 for arguments are Bignums. Other arguments are usually % fix/i-nums. lisp procedure setbits x; % % This function sets the globals for big bignum package. % "x" should be total # of bits per word. <<BBits!*:=iquotient(isub1 x,2); % Total number of bits per word used. BBase!*:=TwoPower BBits!*; % "Beta", where n=A0 + A1*beta + A2*(beta^2)... WordHi!*:=BNum Isub1 BBase!*; % Highest value of Ai WordLow!*:=BMinus WordHi!*; % Lowest value of Ai LogicalBits!*:=ISub1 BBase!*; % Used in LAnd,Lor, etc. SysHi!*:=bsub1 btwopower isub1 x; % Largest representable Syslisp integer. SysLo!*:=BMinus BAdd1 SysHi!*; % Smallest representable Syslisp integer. BBase!*>>; lisp procedure BignumP (V); VectorP V and ((V[0] eq 'BIGPOS) or (V[0] eq 'BIGNEG)); lisp procedure NonBigNumError(V,L); StdError BldMsg(" Expect a BIGNUM in %r, given %p%n",L,V); lisp procedure BSize V; (BignumP V and UpbV V) or 0; lisp procedure GtPOS N; % Creates a positive Bignum with N "Bigits". Begin Scalar B; B:=MkVect N; IPutV(B,0,'BIGPOS); Return B; End; lisp procedure GtNeg N; % Creates a negative Bignum with N "Bigits". Begin Scalar B; B:=MkVect N; IPutV(B,0,'BIGNEG); Return B; End; lisp procedure TrimBigNum V3; % Truncate trailing 0. If Not BignumP V3 then NonBigNumError(V3,'TrimBigNum) else TrimBigNum1(V3,BSize V3); lisp procedure TrimBigNum1(V3,L3); % V3 is a bignum and L3 is the position in it of the highest % possible non-zero digit. Truncate V3 to remove leading zeros, % and if this leaves V3 totally zero make its sign positive; Begin While IGreaterP(L3,0) and IZeroP IGetV(V3,L3) do L3:=ISub1 L3; If IZerop Bsize TruncateVector(V3,L3) then IPutV(V3,0,'BIGPOS); return V3; end; lisp procedure big2sys U; if BLessP(U, SysLo!*) or BGreaterP(U, SysHi!*) then Error(99,list(U," is too large to be a Syslisp integer for BIG2SYS")) else begin scalar L,Sn,res,I; L:=BSize U; if IZeroP L then return 0; Sn:=BMinusP U; res:=IGetV(U,L); I:=ISub1 L; while not IZeroP I do <<res:=ITimes2(res, bbase!*); res:=IPlus2(res, IGetV(U,I)); I:=ISub1 I>>; if Sn then Res:=IMinus Res; return Res; end; lisp procedure TwoPower N; %fix/i-num 2**n 2**n; lisp procedure BTwoPower N; % gives 2**n; n is fix/i-num; result BigNum if not (fixp N or BignumP N) then NonIntegerError(N, 'BTwoPower) else begin scalar quot, rem, V; if bignump N then n:=big2sys n; quot:=Quotient(N,Bbits!*); rem:=Remainder(N,Bbits!*); V:=GtPOS(IAdd1 quot); IFor i:=1:quot do IPutV(v,i,0); IPutV(V,IAdd1 quot,twopower rem); return TrimBigNum1(V,IAdd1 quot); end; lisp procedure BZeroP V1; IZerop BSize V1 and not BMinusP V1; lisp procedure BOneP V1; Not BMinusP V1 and IOneP (BSize V1) and IOneP IGetV(V1,1); lisp procedure BAbs V1; if BMinusP V1 then BMinus V1 else V1; lisp procedure BMax(V1,V2); if BGreaterP(V2,V1) then V2 else V1; lisp procedure BMin(V1,V2); if BLessP(V2,V1) then V2 else V1; lisp procedure BExpt(V1,N); % V1 is Bignum, N is fix/i-num if not fixp N then NonIntegerError(N,'BEXPT) else if IZeroP N then int2B 1 else if IOneP N then V1 else if IMinusP N then BQuotient(int2B 1,BExpt(V1,IMinus N)) else begin scalar V2; V2 := BExpt(V1,IQuotient(N,2)); if IZeroP IRemainder(N,2) then return BTimes2(V2,V2) else return BTimes2(BTimes2(V2,V1),V2) end; % --------------------------------------- % Logical Operations % % All take Bignum arguments lisp procedure BLOr(V1,V2); % The main body of the OR code is only obeyed when both arguments % are positive, and so the result will be positive; if BMinusp V1 or BMinusp V2 then BLnot BLand(BLnot V1,BLnot V2) else begin scalar L1,L2,L3,V3; L1:=BSize V1; L2:=BSize V2; IF L2>L1 then <<L3:=L2; L2:=L1;L1:=L3; V3:=V2; V2:=V1;V1:=V3>>; V3:=GtPOS L1; IFor I:=1:L2 do IPutV(V3,I,ILor(IGetV(V1,I),IGetV(V2,I))); IFor I:=(IAdd1 L2):L1 do IPutV(V3,i,IGetV(V1,I)); Return V3 end; lisp procedure BLXor(V1,V2); % negative arguments are coped with using the identity % LXor(a,b) = LNot LXor(Lnot a,b) = LNor LXor(a,Lnot b); begin scalar L1,L2,L3,V3,S; if BMinusp V1 then << V1 := BLnot V1; S := t >>; if BMinusp V2 then << V2 := BLnot V2; S := not S >>; L1:=BSize V1; L2:=BSize V2; IF L2>L1 then <<L3:=L2; L2:=L1;L1:=L3; V3:=V2; V2:=V1;V1:=V3>>; V3:=GtPOS L1; IFor I:=1:L2 do IPutV(V3,I,ILXor(IGetV(V1,I),IGetV(V2,I))); IFor I:=(IAdd1 L2):L1 do IPutV(V3,i,IGetV(V1,I)); V1:=TrimBigNum1(V3,L1); if S then V1:=BLnot V1; return V1 end; % Not Used Currently: % % lisp Procedure BLDiff(V1,V2); % ***** STILL NEEDS ADJUSTING WRT -VE ARGS ***** % begin scalar V3,L1,L2; % L1:=BSize V1; % L2:=BSize V2; % V3:=GtPOS(max(L1,L2)); % IFor i:=1:min(L1,L2) do % IPutV(V3,i,ILAnd(IGetV(V1,i),ILXor(LogicalBits!*,IGetV(V2,i)))); % if IGreaterP(L1,L2) then IFor i:=(IAdd1 L2):L1 do IPutV(V3,i,IGetV(V1,i)); % if IGreaterP(L2,L1) then IFor i:=(IAdd1 L1):L2 do IPutV(V3,i,0); % return TrimBigNum1(V3,max(L1,L2)); % end; lisp procedure BLAnd(V1,V2); % If both args are -ve the result will be too. Otherwise result will % be positive; if BMinusp V1 and BMinusp V2 then BLnot BLor(BLnot V1,BLnot v2) else begin scalar L1,L2,L3,V3; L1:=BSize V1; L2:=BSize V2; L3:=Min(L1,L2); V3:=GtPOS L3; if BMinusp V1 then IFor I:=1:L3 do IPutV(V3,I,ILand(ILXor(Logicalbits!*,IGetV(V1,I)), IGetV(V2,I))) else if BMinusp V2 then IFor I:=1:L3 do IPutV(V3,I,ILand(IGetV(V1,I), ILXor(Logicalbits!*,IGetV(V2,I)))) else IFor I:=1:L3 do IPutV(V3,I,ILand(IGetV(V1,I),IGetV(V2,I))); return TrimBigNum1(V3,L3); End; lisp procedure BLNot(V1); BMinus BSmallAdd(V1,1); lisp procedure BLShift(V1,V2); % This seems a grimly inefficient way of doing things given that % the representation of big numbers uses a base that is a power of 2. % However it will do for now; if BMinusP V2 then BQuotient(V1, BTwoPower BMinus V2) else BTimes2(V1, BTwoPower V2); % ----------------------------------------- % Arithmetic Functions: % % U, V, V1, V2 are Bignum arguments. lisp procedure BMinus V1; % Negates V1. if BZeroP V1 then V1 else begin scalar L1,V2; L1:=BSize V1; if BMinusP V1 then V2 := GtPOS L1 else V2 := GtNEG L1; IFor I:=1:L1 do IPutV(V2,I,IGetV(V1,I)); return V2; end; % Returns V1 if V1 is strictly less than 0, NIL otherwise. % lisp procedure BMinusP V1; if (IGetV(V1,0) eq 'BIGNEG) then V1 else NIL; % To provide a conveninent ADD with CARRY. lisp procedure AddCarry A; begin scalar S; S:=IPlus2(A,Carry!*); if IGeq(S,BBase!*) then <<Carry!*:= 1; S:=IDifference(S,BBase!*)>> else Carry!*:=0; return S; end; lisp procedure BPlus2(V1,V2); begin scalar Sn1,Sn2; Sn1:=BMinusP V1; Sn2:=BMinusP V2; if Sn1 and Not Sn2 then return BDifference2(V2,BMinus V1,Nil); if Sn2 and Not Sn1 then return BDifference2(V1,BMinus V2,Nil); return BPlusA2(V1,V2,Sn1); end; lisp procedure BPlusA2(V1,V2,Sn1); % Plus with signs pre-checked and begin scalar L1,L2,L3,V3,temp; % identical. L1:=BSize V1; L2:=BSize V2; If IGreaterP(L2,L1) then <<L3:=L2; L2:=L1;L1:=L3; V3:=V2; V2:=V1;V1:=V3>>; L3:=IAdd1 L1; If Sn1 then V3:=GtNeg L3 else V3:=GtPOS L3; Carry!*:=0; IFor I:=1:L2 do <<temp:=IPlus2(IGetV(V1,I),IGetV(V2,I)); IPutV(V3,I,AddCarry temp)>>; temp:=IAdd1 L2; IFor I:=temp:L1 do IPutV(V3,I,AddCarry IGetV(V1,I)); IPutV(V3,L3,Carry!*); % Carry Out Return TrimBigNum1(V3,L3); end; lisp procedure BDifference(V1,V2); if BZeroP V2 then V1 else if BZeroP V1 then BMinus V2 else begin scalar Sn1,Sn2; Sn1:=BMinusP V1; Sn2:=BMinusP V2; if (Sn1 and Not Sn2) or (Sn2 and Not Sn1) then return BPlusA2(V1,BMinus V2,Sn1); return BDifference2(V1,V2,Sn1); end; lisp procedure SubCarry A; begin scalar S; S:=IDifference(A,Carry!*); if ILessP(S,0) then <<Carry!*:=1; S:=IPlus2(BBase!*,S)>> else Carry!*:=0; return S; end; Lisp procedure BDifference2(V1,V2,Sn1); % Signs pre-checked and identical. begin scalar i,L1,L2,L3,V3; L1:=BSize V1; L2:=BSize V2; if IGreaterP(L2,L1) then <<L3:=L1;L1:=L2;L2:=L3; V3:=V1;V1:=V2;V2:=V3; Sn1:=not Sn1>> else if L1 Eq L2 then <<i:=L1; while (IGetV(V2,i) Eq IGetV(V1,i) and IGreaterP(i,1)) do i:=ISub1 i; if IGreaterP(IGetV(V2,i),IGetV(V1,i)) then <<L3:=L1;L1:=L2;L2:=L3; V3:=V1;V1:=V2;V2:=V3;Sn1:=not Sn1>> >>; if Sn1 then V3:=GtNEG L1 else V3:=GtPOS L1; carry!*:=0; IFor I:=1:L2 do IPutV(V3,I,SubCarry IDifference(IGetV(V1,I),IGetV(V2,I))); IFor I:=(IAdd1 L2):L1 do IPutV(V3,I,SubCarry IGetV(V1,I)); return TrimBigNum1(V3,L1); end; lisp procedure BTimes2(V1,V2); begin scalar L1,L2,L3,Sn1,Sn2,V3; L1:=BSize V1; L2:=BSize V2; if IGreaterP(L2,L1) then <<V3:=V1; V1:=V2; V2:=V3; % If V1 is larger, will be fewer L3:=L1; L1:=L2; L2:=L3>>; % iterations of BDigitTimes2. L3:=IPlus2(L1,L2); Sn1:=BMinusP V1; Sn2:=BMinusP V2; If (Sn1 and Sn2) or not(Sn1 or Sn2) then V3:=GtPOS L3 else V3:=GtNEG L3; IFor I:=1:L3 do IPutV(V3,I,0); IFor I:=1:L2 do BDigitTimes2(V1,IGetV(V2,I),L1,I,V3); return TrimBigNum1(V3,L3); end; Lisp procedure BDigitTimes2(V1,V2,L1,I,V3); % V1 is a bignum, V2 a fixnum, L1=BSize L1, I=position of V2 in a bignum, % and V3 is bignum receiving result. I affects where in V3 the result of % a calculation goes; the relationship is that positions I:I+(L1-1) % of V3 receive the products of V2 and positions 1:L1 of V1. % V3 is changed as a side effect here. begin scalar J,carry,temp1,temp2; if zerop V2 then return V3 else << carry:=0; IFor H:=1:L1 do << temp1:=ITimes2(IGetV(V1,H),V2); temp2:=IPlus2(H,ISub1 I); J:=IPlus2(IPlus2(temp1,IGetV(V3,temp2)),carry); IPutV(V3,temp2,IRemainder(J,BBase!*)); carry:=IQuotient(J,BBase!*)>>; IPutV(V3,IPlus2(L1,I),carry)>>; % carry should be < BBase!* here return V3; end; Lisp procedure BSmallTimes2(V1,C); % V1 is a BigNum, C a fixnum. % Assume C positive, ignore sign(V1) % also assume V1 neq 0. if ZeroP C then return GtPOS 0 % Only used from BHardDivide, BReadAdd. else begin scalar J,carry,L1,L2,L3,V3; L1:=BSize V1; L2:=IPlus2(IQuotient(C,BBase!*),L1); L3:=IAdd1 L2; V3:=GtPOS L3; carry:=0; IFor H:=1:L1 do << J:=IPlus2(ITimes2(IGetV(V1,H),C),carry); IPutV(V3,H,IRemainder(J,BBase!*)); carry:=IQuotient(J,BBase!*)>>; IFor H:=(IAdd1 L1):L3 do << IPutV(V3,H,IRemainder(J:=carry,BBase!*)); carry:=IQuotient(J,BBase!*)>>; return TrimBigNum1(V3,L3); end; lisp procedure BQuotient(V1,V2); car BDivide(V1,V2); lisp procedure BRemainder(V1,V2); cdr BDivide(V1,V2); % BDivide returns a dotted pair, (Q . R). Q is the quotient and R is % the remainder. Both are bignums. R is of the same sign as V1. %; smacro procedure BSimpleQuotient(V1,L1,C,SnC); car BSimpleDivide(V1,L1,C,SnC); smacro procedure BSimpleRemainder(V1,L1,C,SnC); cdr BSimpleDivide(V1,L1,C,SnC); lisp procedure BDivide(V1,V2); begin scalar L1,L2,Q,R,V3; L2:=BSize V2; If IZerop L2 then error(99, "Attempt to divide by 0 in BDIVIDE"); L1:=BSize V1; If ILessP(L1,L2) or (L1 Eq L2 and ILessP(IGetV(V1,L1),IGetV(V2,L2))) % This also takes care of case then return (GtPOS 0 . V1); % when V1=0. if IOnep L2 then return BSimpleDivide(V1,L1,IGetV(V2,1),BMinusP V2); return BHardDivide(V1,L1,V2,L2); end; % C is a fixnum (inum?); V1 is a bignum and L1 is its length. % SnC is T if C (which is positive) should be considered negative. % Returns quotient . remainder; each is a bignum. % lisp procedure BSimpleDivide(V1,L1,C,SnC); begin scalar I,P,R,RR,Sn1,V2; Sn1:=BMinusP V1; if (Sn1 and SnC) or not(Sn1 or SnC) then V2:=GtPOS L1 else V2:=GtNEG L1; R:=0; I:=L1; While not IZeroP I do <<P:=IPlus2(ITimes2(R,BBase!*),IGetV(V1,I)); % Overflow. IPutV(V2,I,IQuotient(P, C)); R:=IRemainder(P, C); I:=ISub1 I>>; If Sn1 then RR:=GtNeg 1 else RR:=GtPOS 1; IPutV(RR,1,R); return (TrimBigNum1(V2,L1) . TrimBigNum1(RR,1)); end; lisp procedure BHardDivide(U,Lu,V,Lv); % This is an algorithm taken from Knuth. begin scalar U1,V1,A,D,LCV,LCV1,f,f2,J,K,Lq,carry,temp, LL,M,N,N1,P,Q,QBar,SnU,SnV,U2; N:=Lv; N1:=IAdd1 N; M:=IDifference(Lu,Lv); Lq:=IAdd1 M; % Deal with signs of inputs; SnU:=BMinusP U; SnV:=BMinusp V; % Note that these are not extra-boolean, i.e. % for positive numbers MBinusP returns nil, for % negative it returns its argument. Thus the % test (SnU=SnV) does not reliably compare the signs of % U and V; if SnU then if SnV then Q := GtPOS Lq else Q := GtNEG Lq else if SnV then Q := GtNEG Lq else Q := GtPOS Lq; U1 := GtPOS IAdd1 Lu; % U is ALWAYS stored as if one digit longer; % Compute a scale factor to normalize the long division; D:=IQuotient(BBase!*,IAdd1 IGetV(V,Lv)); % Now, at the same time, I remove the sign information from U and V % and scale them so that the leading coefficeint in V is fairly large; carry := 0; IFor i:=1:Lu do << temp := IPlus2(ITimes2(IGetV(U,I),D),carry); IPutV(U1,I,IRemainder(temp,BBase!*)); carry := IQuotient(temp,BBase!*) >>; Lu := IAdd1 Lu; IPutV(U1,Lu,carry); V1:=BSmallTimes2(V,D); % So far all variables contain safe values, % i.e. numbers < BBase!*; IPutV(V1,0,'BIGPOS); if ILessp(Lv,2) then NonBigNumError(V,'BHARDDIVIDE); % To be safe; LCV := IGetV(V1,Lv); LCV1 := IGetv(V1,ISub1 Lv); % Top two digits of the scaled V accessed once % here outside the main loop; % Now perform the main long division loop; IFor I:=0:M do << J:=IDifference(Lu,I); % J>K; working on U1[K:J] K:=IDifference(J,N1); % in this loop. A:=IGetV(U1,J); P := IPlus2(ITimes2(A,BBase!*),IGetv(U1,Isub1 J)); % N.B. P is up to 30 bits long. Take care! ; if A Eq LCV then QBar := ISub1 BBase!* else QBar := Iquotient(P,LCV); % approximate next digit; f:=ITimes2(QBar,LCV1); f2:=IPlus2(ITimes2(IDifference(P,ITimes2(QBar,LCV)),BBase!*), IGetV(U1,IDifference(J,2))); while IGreaterP(f,f2) do << % Correct most overshoots in Qbar; QBar:=ISub1 QBar; f:=IDifference(f,LCV1);; f2:=IPlus2(f2,ITimes2(LCV,BBase!*)) >>; carry := 0; % Ready to subtract QBar*V1 from U1; IFor L:=1:N do << temp := IPlus2( Idifference( IGetV(U1,IPlus2(K,L)), ITimes2(QBar,IGetV(V1,L))), carry); carry := IQuotient(temp,BBase!*); temp := IRemainder(temp,BBase!*); if IMinusp temp then << carry := ISub1 carry; temp := IPlus2(temp,BBase!*) >>; IPutV(U1,IPlus2(K,L),temp) >>; % Now propagate borrows up as far as they go; LL := IPlus2(K,N); while (not IZeroP carry) and ILessp(LL,J) do << LL := IAdd1 LL; temp := IPlus2(IGetV(U1,LL),carry); carry := IQuotient(temp,BBase!*); temp := IRemainder(temp,BBase!*); if IMinusP temp then << carry := ISub1 carry; temp := IPlus2(temp,BBase!*) >>; IPutV(U1,LL,temp) >>; if not IZerop carry then << % QBar was still wrong - correction step needed. % This should not happen very often; QBar := ISub1 QBar; % Add V1 back into U1; carry := 0; IFor L := 1:N do << carry := IPlus2( IPlus2(IGetV(U1,Iplus2(K,L)), IGetV(V1,L)), carry); IPutV(U1,IPlus2(K,L),IRemainder(carry,BBase!*)); carry := IQuotient(carry,BBase!*) >>; LL := IPlus2(K,N); while ILessp(LL,J) do << LL := IAdd1 LL; carry := IPlus2(IGetv(U1,LL),carry); IPutV(U1,LL,IRemainder(carry,BBase!*)); carry := IQuotient(carry,BBase!*) >> >>; IPutV(Q,IDifference(Lq,I),QBar) >>; % End of main loop; U1 := TrimBigNum1(U1,IDifference(Lu,M)); f := 0; f2 := 0; % Clean up potentially wild values; if not BZeroP U1 then << % Unnormalize the remainder by dividing by D if SnU then IPutV(U1,0,'BIGNEG); if not IOnep D then << Lu := BSize U1; carry := 0; IFor L:=Lu step -1 until 1 do << P := IPlus2(ITimes2(carry,BBase!*),IGetV(U1,L)); IPutv(U1,L,IQuotient(P,D)); carry := IRemainder(P,D) >>; P := 0; if not IZeroP carry then BHardBug("remainder when unscaling", U,V,TrimBigNum1(U1,Lu),TrimBigNum1(Q,Lq)); U1 := TrimBigNum1(U1,Lu) >> >>; Q := TrimBigNum1(Q,Lq); % In case leading digit happened to be zero; P := 0; % flush out a 30 bit number; % Here, for debugging purposes, I will try to validate the results I % have obtained by testing if Q*V+U1=U and 0<=U1<V. I Know this slows things % down, but I will remove it when my confidence has improved somewhat; % if not BZerop U1 then << % if (BMinusP U and not BMinusP U1) or % (BMinusP U1 and not BMinusP U) then % BHardBug("remainder has wrong sign",U,V,U1,Q) >>; % if not BAbs U1<BAbs V then BHardBug("remainder out of range",U,V,U1,Q) % else if not BZerop(BDifference(BPlus2(BTimes2(Q,V),U1),U)) then % BHardBug("quotient or remainder incorrect",U,V,U1,Q); return (Q . U1) end; lisp procedure BHardBug(msg,U,V,R,Q); % Because the inputs to BHardDivide are probably rather large, I am not % going to rely on BldMsg to display them; << Prin2T "***** Internal error in BHardDivide"; Prin2 "arg1="; Prin2T U; Prin2 "arg2="; Prin2T V; Prin2 "computed quotient="; Prin2T Q; Prin2 "computed remainder="; Prin2T R; StdError msg >>; lisp procedure BGreaterP(U,V); if BMinusP U then if BMinusP V then BUnsignedGreaterP(V,U) else nil else if BMinusP V then U else BUnsignedGreaterP(U,V); lisp procedure BLessp(U,V); if BMinusP U then if BMinusP V then BUnsignedGreaterP(U,V) else U else if BMinusP V then nil else BUnsignedGreaterP(V,U); lisp procedure BGeq(U,V); if BMinusP U then if BMinusP V then BUnsignedGeq(V,U) else nil else if BMinusP V then U else BUnsignedGeq(U,V); lisp procedure BLeq(U,V); if BMinusP U then if BMinusP V then BUnsignedGeq(U,V) else U else if BMinusP V then nil else BUnsignedGeq(V,U); lisp procedure BUnsignedGreaterP(U,V); % Compare magnitudes of two bignums; begin scalar Lu,Lv,I; Lu := BSize U; Lv := BSize V; if not (Lu eq Lv) then << if IGreaterP(Lu,Lv) then return U else return nil >>; while IGetV(U,Lv) eq IGetV(V,Lv) and IGreaterP(Lv,1) do Lv := ISub1 Lv; if IGreaterP(IGetV(U,Lv),IGetV(V,Lv)) then return U else return nil end; symbolic procedure BUnsignedGeq(U,V); % Compare magnitudes of two unsigned bignums; begin scalar Lu,Lv; Lu := BSize U; Lv := BSize V; if not (Lu eq Lv) then << if IGreaterP(Lu,Lv) then return U else return nil >>; while IGetV(U,Lv) eq IGetV(V,Lv) and IGreaterP(Lv,1) do Lv := ISub1 Lv; If IGreaterP(IGetV(V,Lv),IGetV(U,Lv)) then return nil else return U end; lisp procedure BAdd1 V; BSmallAdd(V,1); lisp procedure BSub1 U; BSmallDiff(U,1); % ------------------------------------------------ % Conversion to Float: lisp procedure FloatFromBigNum V; if BZeroP V then 0.0 else if BGreaterP(V, FloatHi!*) or BLessp(V, FloatLow!*) then Error(99,list("Argument, ",V," to FLOAT is too large")) else begin scalar L,Res,Sn,I; L:=BSize V; Sn:=BMinusP V; Res:=float IGetv(V,L); I:=ISub1 L; While not IZeroP I do << Res:=res*BBase!*; Res:=Res +IGetV(V,I); I:=ISub1 I>>; if Sn then Res:=minus res; return res; end; % ------------------------------------------------ % Input and Output: Digit2Letter!* := % Ascii values of digits and characters. '[48 49 50 51 52 53 54 55 56 57 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90]; % OutputBase!* is assumed to be positive and less than 37. lisp procedure BChannelPrin2(Channel,V); If not BignumP V then NonBigNumError(V, 'BPrin) %need? else begin scalar quot, rem, div, result, resultsign, myobase; myobase:=OutputBase!*; resultsign:=BMinusP V; div:=BSimpleDivide(V,Bsize V,OutputBase!*,nil); quot:=car div; rem:=cdr div; if Bzerop rem then rem:=0 else rem:=IGetV(rem,1); result:=rem . result; while Not BZeroP quot do <<div:=BSimpleDivide(quot,Bsize quot,OutputBase!*,nil); quot:=car div; rem:=cdr div; if Bzerop rem then rem:=0 else rem:=IGetV(rem,1); result:=rem . result>>; if resultsign then channelwritechar(Channel,char !-); if myobase neq 10 then <<ChannelWriteSysInteger(channel,myobase,10); ChannelWriteChar(Channel, char !#)>>; For each u in result do ChannelWriteChar(Channel, IGetV(digit2letter!*,u)); OutputBase!*:=myobase; return; end; lisp procedure BRead(s,radix,sn); % radix is < Bbase!* %s=string of digits, radix=base, sn=1 or -1 begin scalar sz, res, ch; sz:=size s; res:=GtPOS 1; ch:=indx(s,0); if IGeq(ch,char A) and ILeq(ch,char Z) then ch:=IPlus2(IDifference(ch,char A),10); if IGeq(ch,char 0) and ILeq(ch,char 9) then ch:=IDifference(ch,char 0); IPutV(res,1,ch); IFor i:=1:sz do <<ch:=indx(s,i); if IGeq(ch,char A) and ILeq(ch,char Z) then ch:=IDifference(ch,IDifference(char A,10)); if IGeq(ch,char 0) and ILeq(ch,char 9) then ch:=IDifference(ch,char 0); res:=BReadAdd(res, radix, ch)>>; if iminusp sn then res:=BMinus res; return res; end; lisp procedure BReadAdd(V, radix, ch); << V:=BSmallTimes2(V, radix); V:=BSmallAdd(V,ch)>>; lisp procedure BSmallAdd(V,C); %V big, C fix. if IZerop C then return V else if Bzerop V then return int2B C else if BMinusp V then BMinus BSmallDiff(BMinus V, C) else if IMinusP C then BSmallDiff(V, IMinus C) else begin scalar V1,L1; Carry!*:=C; L1:=BSize V; V1:=GtPOS(IAdd1 L1); IFor i:=1:L1 do IPutV(V1,i,addcarry IGetV(V,i)); if IOneP carry!* then IPutV(V1,IAdd1 L1,1) else return TrimBigNum1(V1,L1); return V1 end; lisp procedure BNum N; % temporary? Creates a Bignum of one digit, value N. begin scalar B; if IZerop n then return GtPOS 0 else if IMinusp N then <<b:=GtNEG 1; n:= IMinus n>> else b:=GtPos 1; IPutV(b,1,N); Return b; end; lisp procedure BSmallDiff(V,C); %V big, C fix if IZerop C then V else if BZeroP V then int2B IMinus C else if BMinusP V then BMinus BSmallAdd(BMinus V, C) else if IMinusP C then BSmallAdd(V, IMinus C) else begin scalar V1,L1; Carry!*:=C; L1:=BSize V; V1:=GtPOS L1; IFor i:=1:L1 do IPuTV(V1,i,subcarry IGetV(V,i)); if not IZeroP carry!* then StdError BldMsg(" BSmallDiff V<C %p %p%n",V,C); return TrimBigNum1(V1,L1); end; lisp procedure int2B n; % Temporary? Creates BigNum of value N. if not fixp n then NonIntegerError(n, 'int2B) else if ILessP(n,Bbase!*) then BNum n else begin scalar Str,ind,rad,Sn,r; Str:=bldmsg("%w",n); % like an "int2string" if indx(str,0)=char '!- then <<Sn:=-1; str:=sub(str,1,ISub1 (size str))>> else Sn:=1; IFor i:=0:size str do if indx(str,i)=char '!# then ind:=i; if ind then <<r:=sub(str,0,ISub1 ind); rad:=0; IFor i:=0:size r do rad:=IPlus2(ITimes2(rad,10),IDifference(indx(r,i),char 0)); str:=sub(str,IAdd1 ind,IDifference(size str,IAdd1 ind))>> else rad:=10; return Bread(str,rad,sn); end; %----------------------------------------------------- % "Fix" for Bignums lisp procedure bigfromfloat X; if fixp x or bigp x then x else begin scalar bigpart,floatpart,power,sign,thispart; if minusp X then <<sign:=-1; X:=minus X>> else sign:=1; bigpart:=bnum 0; while neq(X, 0) and neq(x,0.0) do << if X < bbase!* then << bigpart:=bplus2(bigpart, bnum fix x); X:=0 >> else <<floatpart:=x; power:=0; while floatpart>=bbase!* do % get high end of number. <<floatpart:=floatpart/bbase!*; power:=power + bbits!* >>; thispart:=btimes2(btwopower power, bnum fix floatpart); X:=X- floatfrombignum thispart; bigpart:=bplus2(bigpart, thispart) >> >>; if minusp sign then bigpart := bminus bigpart; return bigpart; end; if_system(VAX, <<setbits 32; FloatHi!*:=btimes2(bdifference(btwopower 67, btwopower 11), btwopower 60);% Largest representable float. FloatLow!*:=BMinus FloatHi!*>>); if_system(PDP10, <<setbits 36; FloatHi!*:=btimes2(bsub1 btwopower 62, btwopower 65); FloatLow!*:=BMinus FloatHi!*>>); % End of BIGBIG.RED ;