% NBIG0.RED - Vector based BIGNUM package with INUM operations
% M. L. Griss & B Morrison, 25 June 1982.
% Copyright (C) 1982, A. C. Norman, B. Morrison, M. Griss
%
% Revision log:
% 10 March, 1983, MLG
% LSH in Twopower replaced by 2**n
% Fixed a bug in SYS2BIG that did not convert negative BIGNUMS correctly
% 7 February 1983, MLG
% Merged in NBIG1 (see its "revision history" below), plus clean-up.
% Revision History of old NBIG1:
% 28 Dec 1982, MLG:
% Added BigZeroP and BigOneP for NArith
% Changed Name to NBIG1.RED from BIGFACE
% 22 Dec 1982, MLG:
% Change way of converting from VECT to BIGN
% Move Module dependency to .BUILD file
% Changes for NEW-ARITH, involve name changes for MAKEFIXNUM
% ISINUM, etc.
% 21 December, 82: MLG
% Change PRIN1 and PRIN2 hooks to refer to RecursiveChannelprinx
% which changed in PK:PRINTERS.RED for prinlevel stuff
% November: Variety of Bug Fixes by A. Norman
% Use the BIGN tag for better Interface
%
% 31 Dec 1982, MLG
% Changed BNUM to check if arg ALREADY Big. Kludge
% since new NARITH makes some things BIG earlier
% since it calls the BIG funcs directly
% 20 Dec 1982, MLG
% Changed TrimBigNUM to TrimBigNum1 in BhardDivide
%
% 14 Dec 1982, MLG
% Changed to put LOAD and IMPORTS in BUILD file
%
% 31 August 1982, A. C . Norman
% Adjustments 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);
% ---------------------------------------------------------------
% -----------------------
% 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.
% BETA.RED - some values of BETA testing
% On DEC-20, Important Ranges are:
% --------------------------------
% POSBETA | 0 | n |
% --------------------------------
% 19 17 bits
% --------------------------------
% NEGBETA | -1 | |
% --------------------------------
%
% --------------------------------
% POSINT | 0 | 0 | |
% --------------------------------
% 5 13 18 bits
% --------------------------------
% NEGINT | -1 | -1 | |
% --------------------------------
% Thus BETA: 2^17-1 -131072 ... 131071
% INT 2^18-1 -262144 ... 262143
% FIX 2^35-1 -34359738368 ... 34359738367
% [Note that one bit used for sign in 36 bit word]
fluid '(BigBetaHi!* % Largest BetaNum in BIG format
BigBetaLow!* % Smallest BetaNum in BIG format
BetaHi!* % Largest BetaNum as Inum
BetaLow!* % Smallest BetaNum as Inum
SysHi!* % Largest SYSINT in FixN format
SysLow!* % Smallest SYSINT in FixN format
BigSysHi!* % Largest SYSINT in BIG format
BigSysLow!* % Smallest SYSINT in BIG format
FloatSysHi!* % Largest SYSINT in Float format
FloatSysLow!* % Smallest SYSINT in Float format
BBase!* % BETA, base of system
FloatBbase!* % As a float
BigFloatHi!* % Largest Float in BIG format
BigFloatLow!* % Smallest Float in BIG format
StaticBig!* % Warray for conversion of SYS to BIG
Bone!* % A one
Bzero!* % A zero
BBits!* % Number of Bits in BBASE!*
LogicalBits!*
Digit2Letter!*
Carry!*
OutputBase!*
);
% --------------------------------------------------------------------------
% --------------------------------------------------------------------------
% Support functions:
%
% U, V, V1, V2 for arguments are Bignums. Other arguments are usually
% fix/i-nums.
smacro procedure PutBig(b,i,val);
% Access elements of a BIGNUM
IputV(b,i,val);
smacro procedure GetBig(b,i);
% Access elements of a BIGNUM
IgetV(B,i);
procedure setbits x;
%
% This function sets the globals for big bignum package.
% "x" should be total # of bits per word.
Begin scalar y;
BBits!*:=iquotient(isub1 x,2); % Total number of bits per word used.
BBase!*:=TwoPower BBits!*; % "Beta", where n=A0 + A1*beta + A2*(beta^2).
FloatBbase!* := IntFloat Bbase!*;
LogicalBits!*:=ISub1 BBase!*; % Used in LAnd,Lor, etc.
BetaHi!*:=isub1 Bbase!*;
BetaLow!* :=Iminus Bbase!*;
Bone!* := Bnum 1;
Bzero!* := Bnum 0;
BigBetaHi!*:=BNum BetaHi!*; % Highest value of Ai
BigBetaLow!*:=BMinus BigBetaHi!*; % Lowest value of Ai
% here assume 2's complement
y:=TwoPower idifference (x,2); % eg, 36 bits, 2^35-1=2^34+2^34-1
SysHi!* :=y+(y-1);
y:=-y;
Syslow!* :=y+y;
BigSysHi!*:=bdifference(btwopower isub1 x,
Bone!*); % Largest representable Syslisp integer.
% Note that SYSPOS has leading 0, ie only x-1 active bits
BigSysLow!*:=BMinus BPlus2(Bone!*, BigSysHi!*);
% Smallest representable Syslisp integer.
end;
procedure NonBigNumError(V,L);
StdError BldMsg(" Expect a BIGNUM in %r, given %p%n",L,V);
procedure BSize V;
% Upper Limit of [BIGxxx a1 ... An]
If BigP V then VecLen VecInf V else 0;
procedure GtPOS N;
% Allocate [BIGPOS a1 ... an]
Begin
N:=MkVect N;
IPutV(N,0,'BIGPOS);
Return MkBigN Vecinf N;
End;
procedure GtNeg N;
% Allocate [BIGNEG a1 ... an]
Begin
N:=MkVect N;
IPutV(N,0,'BIGNEG);
Return MkBigN VecInf N;
End;
procedure TrimBigNum V3;
% truncate trailing 0
If Not BigP V3 then NonBigNumError(V3,'TrimBigNum)
else TrimBigNum1(V3,BSize V3);
procedure TrimBigNum1(B,L3);
Begin scalar v3;
V3:=BigAsVec B;
While IGreaterP(L3,0) and IZeroP IGetV(V3,L3) do L3:=ISub1 L3;
If IZerop UpBv TruncateVector(V3,L3) then return GtPOS 0
else return B;
end;
procedure BigAsVec B;
% In order to see BIGITS
MkVec Inf B;
procedure VecAsBig V;
MkBigN VecInf V;
Procedure BIG2Sys U;
% Convert a BIG to SYS, if in range
If Blessp(U,BigSysLow!*) or Bgreaterp(U,BigSysHi!*) then
ContinuableError(99,"BIGNUM too large to convert to SYS", U)
else Big2SysAux U;
procedure Big2SysAux U;
% Convert a BIGN that is in range to a SYSINT
begin scalar L,Sn,res;
L:=BSize U;
if IZeroP L then return 0;
res:=IGetV(U,L);
L:=ISub1 L;
If BMinusP U then
<<res:=-res;
while L neq 0 do <<res:=ITimes2(res, Bbase!*);
res:=IDifference(res, IGetV(U,L));
L:=ISub1 L>>;
>>
else
while L neq 0 do <<res:=ITimes2(res, Bbase!*);
res:=IPlus2(res, IGetV(U,L));
L:=ISub1 L>>;
return Res;
end;
procedure TwoPower N; %fix/i-num 2**n
2**n;
procedure BTwoPower N; % gives 2**n; n is fix/i-num; result BigNum
if not (fixp N or BigP N) then NonIntegerError(N, 'BTwoPower)
else begin scalar quot, rem, V;
if BigP 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;
procedure BZeroP V1;
IZerop BSize V1 and not BMinusP V1;
procedure BOneP V1;
Not BMinusP V1 and IOneP (BSize V1) and IOneP IGetV(V1,1);
procedure BAbs V1;
if BMinusP V1 then BMinus V1 else V1;
procedure BMax(V1,V2);
if BGreaterP(V2,V1) then V2 else V1;
procedure BMin(V1,V2);
if BLessP(V2,V1) then V2 else V1;
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 Bone!*
else if IOneP N then V1
else if IMinusP N then BQuotient(Bone!*,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
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;
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:
%
% 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;
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;
procedure BLNot(V1);
BMinus BSmallAdd(V1,1);
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.
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.
%
procedure BMinusP V1;
if (IGetV(V1,0) eq 'BIGNEG) then V1 else NIL;
% To provide a conveninent ADD with CARRY.
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;
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;
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;
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;
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;
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;
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;
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;
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;
procedure BQuotient(V1,V2);
car BDivide(V1,V2);
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);
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.
%
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;
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;
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 >>;
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);
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);
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);
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);
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;
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;
procedure BAdd1 V;
BSmallAdd(V, 1);
procedure BSub1 U;
BSmallDiff(U, 1);
% ------------------------------------------------
% Conversion to Float:
procedure FloatFromBigNum V;
if BZeroP V then 0.0
else if BGreaterP(V, BigFloatHi!*) or BLessp(V, BigFloatLow!*)
then Error(99,list("Argument, ",V," to FLOAT is too large"))
else begin scalar L,Res,Sn,I;
% Careful, do not want to call itself recursively
L:=BSize V;
Sn:=BMinusP V;
Res:=IntFloat IGetv(V,L);
I:=ISub1 L;
While not IZeroP I do << Res:=FloatTimes2(res,FloatBBase!*);
Res:=FloatPlus2(Res, IntFloat 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.
procedure BChannelPrin2(Channel,V);
If not BigP 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;
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;
procedure BReadAdd(V, radix, ch);
<< V:=BSmallTimes2(V, radix);
V:=BSmallAdd(V,ch)>>;
procedure BSmallAdd(V,C); %V big, C fix.
if IZerop C then return V
else if Bzerop V then return int2Big 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;
procedure BNum N;
% Creates a Bignum of one BETA digit, value N.
% N is POS or NEG
IF BIGP N then N else BnumAux N;
procedure BNumAux N;
% Creates a Bignum of one BIGIT value N.
% N is POS or NEG
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;
procedure BSmallDiff(V,C); %V big, C fix
if IZerop C then V
else if BZeroP V then int2Big 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;
on syslisp;
syslsp procedure int2Big n;
% Creates BigNum of value N.
% From any N, BETA,INUM,FIXNUM or BIGNUM
case tag n of
NEGINT,POSINT: sys2Big n;
FIXN: sys2Big fixval fixinf n;
BIGN: N;
default: NonIntegerError(n, 'int2Big);
End;
off syslisp;
% Convert BIGNUMs to FLOAT
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:=bzero!*;
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;
% Now Install Interfacing
on syslisp;
syslsp procedure SetUpGlobals;
<< Prin2t '"SetupGlobals";
SetBits BitsPerWord;
Prin2T '" ... done";>>;
off syslisp;
SetupGlobals();
LoadTime <<
StaticBig!*:=GtWarray 10>>;
% Assume dont need more than 10 slots to represent a BigNum
% Version of SYSint
% -- Output---
% MLG Change to interface to Recursive hooks, added for
% Prinlevel stuff
CopyD('OldChannelPrin1,'RecursiveChannelPrin1);
CopyD('OldChannelPrin2,'RecursiveChannelPrin2);
Procedure RecursiveChannelPrin1(Channel,U,Level);
<<if BigP U then BChannelPrin2(Channel,U)
else OldChannelPrin1(Channel, U,Level);U>>;
Procedure RecursiveChannelPrin2(Channel,U,level);
<<If BigP U then BChannelPrin2(Channel, U)
else OldChannelPrin2(Channel, U,level);U>>;
procedure checkifreallybig UU;
% If BIGNUM result is in older FIXNUM or INUM range
% Convert Back.
%/ Need a faster test
if BLessP(UU, BigSysLow!*) or BGreaterp(UU,BigSysHi!*) then UU
else Sys2Int Big2SysAux UU;
procedure checkifreallybigpair VV;
% Used to process DIVIDE
checkifreallybig car VV . checkifreallybig cdr VV;
procedure checkifreallybigornil UU;
% Used for EXTRA-boolean tests
if Null UU or BLessp(UU, BigSysLow!*) or BGreaterP(UU,BigSysHi!*) then UU
else Sys2Int Big2SysAux UU;
procedure BigPlus2(U,V);
CheckIfReallyBig BPlus2(U,V);
procedure BigDifference(U,V);
CheckIfReallyBig BDifference(U,V);
procedure BigTimes2(U,V);
CheckIfReallyBig BTimes2(U,V);
procedure BigDivide(U,V);
CheckIfReallyBigPair BDivide(U,V);
procedure BigQuotient(U,V);
CheckIfReallyBig BQuotient(U,V);
procedure BigRemainder(U,V);
CheckIfReallyBig BRemainder(U,V);
procedure BigLAnd(U,V);
CheckIfReallyBig BLand(U,V);
procedure BigLOr(U,V);
CheckIfReallyBig BLOr(U,V);
procedure BigLXOr(U,V);
CheckIfReallyBig BLXor(U,V);
procedure BigLShift(U,V);
CheckIfReallyBig BLShift(U,V);
on syslisp;
procedure Lshift(U,V);
If BetaP U and BetaP V
then (if V<0 then Sys2Int Wshift(U,V)
else if V< LispVar (BBits!* ) then Sys2Int Wshift(U,V)
else BigLshift(Sys2Big U, Sys2Big V) )
else BigLshift(Sys2Big U, Sys2Big V) ;
off syslisp;
Copyd('LSH,'Lshift);
procedure BigGreaterP(U,V);
CheckIfReallyBigOrNil BGreaterP(U,V);
procedure BigLessP(U,V);
CheckIfReallyBigOrNil BLessP(U,V);
procedure BigAdd1 U;
CheckIfReallyBig BAdd1 U;
procedure BigSub1 U;
CheckIfReallyBig BSub1 U;
procedure BigLNot U;
CheckIfReallyBig BLNot U;
procedure BigMinus U;
CheckIfReallyBig BMinus U;
procedure BigMinusP U;
CheckIfReallyBigOrNil BMinusP U;
procedure BigOneP U;
CheckIfReallyBigOrNil BOneP U;
procedure BigZeroP U;
CheckIfReallyBigOrNil BZeroP U;
% ---- Input ----
procedure MakeStringIntoLispInteger(S,Radix,Sn);
CheckIfReallyBig BRead(S,Radix,Sn);
on syslisp;
procedure Int2Sys N;
% Convert a random FIXed number to WORD Integer
case tag(N) of
POSINT,NEGINT: N;
FIXN: FixVal FixInf N;
BIGN: Big2SysAux N;
default: NonNumber1Error(N,'Int2SYS);
End;
syslsp procedure Sys2Big N;
% Convert a SYSint to a BIG
% Must NOT use generic arith here
% Careful that no GC if this BIGger than INUM
Begin scalar Sn, A, B;
If N=0 then return GtPos 0;
A:= LispVar StaticBig!*; % Grab the base
If N<0 then sn:=T;
A[1]:=N; % Plant number
N:=1; % now use N as counter
% Careful handling of -N in case have largest NEG, not just
% flip sign
If Sn then <<B:=-Bbase!*;
While A[n]<=B do
<<N:=N+1;
A[n]:=A[n-1]/Bbase!*;
A[n-1]:=A[n-1]-a[n]*Bbase!*>>;
B:=GtNeg N;
For i:=1:N do Iputv(B,i,-A[i])>>
else << While A[n]>=Bbase!* do
<<N:=N+1;
A[n]:=A[n-1]/Bbase!*;
A[n-1]:=A[n-1]-a[n]*Bbase!*>>;
B:= GtPos N;
For i:=1:N do IputV(B,i,A[i])>>;
Return B;
End;
off syslisp;
% Coercion/Transfer Functions
copyd('oldFloatFix,'FloatFix);
procedure FloatFix U;
% Careful of sign and range
If FloatSysLow!* <= U and U <= FloatSysHi!* then Oldfloatfix U
else bigfromfloat U;
on syslisp;
procedure BetaP x;
% test if NUMBER in reduced INUM range
If Intp x then (x <= Lispvar(betaHi!*)) and (x >= LispVar(betaLow!*))
else NIL;
procedure BetaRangeP x;
% Test if SYSINT in reduced INUM range
if (x <= Lispvar(betaHi!*)) then (x>=LispVar(betaLow!*)) else NIL;
procedure Beta2P(x,y);
% Check for 2 argument arithmetic functions
if BetaP x then BetaP y;
off syslisp;
End;
end;