reduce-historical/psl-1983/util/bigbig.red

% 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 ;