1%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2% 3% File: PU:NBIG0.SL 4% Description: Vector based BIGNUM package with INUM operations 5% Author: M. L. Griss & B Morrison 6% Created: 25 June 1982 7% Modified: 8% Mode: Lisp 9% Package: Utilities 10% Status: Open Source: BSD License 11% 12% (c) Copyright 1982, University of Utah 13% 14% Redistribution and use in source and binary forms, with or without 15% modification, are permitted provided that the following conditions are met: 16% 17% * Redistributions of source code must retain the relevant copyright 18% notice, this list of conditions and the following disclaimer. 19% * Redistributions in binary form must reproduce the above copyright 20% notice, this list of conditions and the following disclaimer in the 21% documentation and/or other materials provided with the distribution. 22% 23% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 24% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, 25% THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR 26% PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNERS OR 27% CONTRIBUTORS 28% BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 29% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 30% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 31% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 32% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 33% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 34% POSSIBILITY OF SUCH DAMAGE. 35% 36%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 37% 38% Revisions: 39% 40% 15-Nov-1988 (Tony Hearn) 41% Commented out bexpt, since expt in $pnk/easy-sl.sl supersedes it. 42% 25-Aug-1988 (C. Burdorf and T. Yamamoto) 43% Used MC68020 for maxfloatsize. Set up SPARC version for bigfloathi 44% and bigfloatlow initialization. 45% 22-Jan-88 (James Davenport for Will Galway) 46% floatfrombignum - don't recurse on negative numbers. 47% 05-Jan-88 (James Davenport for Mark Swanson) 48% bharddivide: various additional checks on dangerous values. 49% 04-Jan-88 (James Davenport for Mark Swanson) 50% bsimpledivide: set p to 0 before we might GC. 51% 09-Sep-87 (Leigh Stoller for Tony Hearn) 52% Added new printing routines for bignums. 53% 05-Sep-87 (Leigh Stoller) 54% Add German's reclaim fix. 55% 02-Sep-87 (Leigh Stoller for Russ Fish) 56% Use IntLShift instead of WShift on the VAX in redefinition of LShift. 57% 05-May-87 (Leigh Stoller) 58% Added Bob Kessler's fixes to bland lshift. 59% 23-Aug-84 09:39:05 (Brian Beach) 60% Made changes for use with micro-kernel: Removed (WCONST xxx) uses and 61% changed the two WFORs to IFORs. Removed (LOAD SYSLISP). 62% 26 Mar 1984 2146-PST (Nancy Kendzierski) 63% Commented out: (compiletime (if_system mc68000 (load c!-expr!-fix))) 64% since the file was removed when the new dreg compiler came into effect. 65% 05-Dec-83 17:50:15 (Nancy Kendzierski) 66% Added contents of .BUILD file. 67% 02-Dec-83 18:12:56 (Nancy Kendzierski) 68% Translated from Rlisp to Lisp. 69% 22-Jul-83 Nancy Kendzierski 70% Removed MC68000 *WTIMES2 fix, since it is now handled by hp-cmac.sl 71% [.BUILD file] 72% 10 March, 1983, MLG 73% LSH in Twopower replaced by 2**n 74% Fixed a bug in SYS2BIG that did not convert negative BIGNUMS correctly 75% 7 February 1983, MLG 76% Merged in NBIG1 (see its "revision history" below), plus clean-up. 77% Revision History of old NBIG1: 78% 28 Dec 1982, MLG: 79% Added BigZeroP and BigOneP for NArith 80% Changed Name to NBIG1.RED from BIGFACE 81% 22 Dec 1982, MLG: 82% Change way of converting from VECT to BIGN 83% Move Module dependency to .BUILD file 84% Changes for NEW-ARITH, involve name changes for MAKEFIXNUM 85% ISINUM, etc. 86% 21 December, 82: MLG 87% Change PRIN1 and PRIN2 hooks to refer to RecursiveChannelprinx 88% which changed in PK:PRINTERS.RED for prinlevel stuff 89% November: Variety of Bug Fixes by A. Norman 90% Use the BIGN tag for better Interface 91% 31 Dec 1982, MLG 92% Changed BNUM to check if arg ALREADY Big. Kludge 93% since new NARITH makes some things BIG earlier 94% since it calls the BIG funcs directly 95% 20 Dec 1982, MLG 96% Changed TrimBigNUM to TrimBigNum1 in BhardDivide 97% 14 Dec 1982, MLG 98% Changed to put LOAD and IMPORTS in BUILD file 99% 31 August 1982, A. C . Norman 100% Adjustments to many routines: in particular corrections to BHardDivide 101% (case D6 utterly wrong), and adjustments to BExpt (for performance) and 102% all logical operators (for treatment of negative inputs); 103% 104%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 105 106%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 107% A bignum will be a VECTOR of Bigits: (digits in base BigBase): 108% [BIGPOS b1 ... bn] or [BIGNEG b1 ... bn]. BigZero is thus [BIGPOS] 109% All numbers are positive, with BIGNEG as 0 element to indicate negatives. 110% BETA.RED - some values of BETA testing 111% On DEC-20, Important Ranges are: 112% -------------------------------- 113% POSBETA | 0 | n | 114% -------------------------------- 115% 19 17 bits 116% -------------------------------- 117% NEGBETA | -1 | | 118% -------------------------------- 119% 120% -------------------------------- 121% POSINT | 0 | 0 | | 122% -------------------------------- 123% 5 13 18 bits 124% -------------------------------- 125% NEGINT | -1 | -1 | | 126% -------------------------------- 127% Thus BETA: 2^17-1 -131072 ... 131071 128% INT 2^18-1 -262144 ... 262143 129% FIX 2^35-1 -34359738368 ... 34359738367 130% [Note that one bit used for sign in 36 bit word] 131%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 132 133(compiletime (load fast-vector inum if-system)) 134(compiletime (remprop 'addcarry 'opencode)) 135(compiletime (remprop 'subcarry 'opencode)) 136 137% For Bhard-Divide Bug; only for bignums, known to be buggy for other machines. 138% (compiletime (if_system mc68000 (load c!-expr!-fix))) 139 140(fluid '(bigbetahi!* % Largest BetaNum in BIG format 141 bigbetalow!* % Smallest BetaNum in BIG format 142 betahi!* % Largest BetaNum as Inum 143 betalow!* % Smallest BetaNum as Inum 144 syshi!* % Largest SYSINT in FixN format 145 syslow!* % Smallest SYSINT in FixN format 146 bigsyshi!* % Largest SYSINT in BIG format 147 bigsyslow!* % Smallest SYSINT in BIG format 148 floatsyshi!* % Largest SYSINT in Float format 149 floatsyslow!* % Smallest SYSINT in Float format 150 bbase!* % BETA, base of system 151 floatbbase!* % As a float 152 bigfloathi!* % Largest Float in BIG format 153 bigfloatlow!* % Smallest Float in BIG format 154 staticbig!* % Warray for conversion of SYS to BIG 155 bone!* % A one 156 bzero!* % A zero 157 bbits!* % Number of Bits in BBASE!* 158 logicalbits!* digit2letter!* carry!* outputbase!*)) 159 160% -------------------------------------------------------------------------- 161 162% -------------------------------------------------------------------------- 163 164% Support functions: 165% 166% U, V, V1, V2 for arguments are Bignums. Other arguments are usually 167% fix/i-nums. 168(ds putbig (b i val) % Access elements of a BIGNUM 169 (iputv b i val)) 170 171(ds getbig (b i) % Access elements of a BIGNUM 172 (igetv b i)) 173 174(de setbits (x) 175 % 176 % This function sets the globals for big bignum package. 177 % "x" should be total # of bits per word. 178 (prog (y) 179 (setq bbits!* (iquotient (isub1 x) 2)) 180 % Total number of bits per word used. 181 (setq bbase!* (twopower bbits!*)) 182 % "Beta", where n=A0 + A1*beta + A2*(beta^2). 183 (setq floatbbase!* (intfloat bbase!*)) 184 (setq logicalbits!* (isub1 bbase!*)) 185 % Used in LAnd,Lor, etc. 186 (setq betahi!* (isub1 bbase!*)) 187 (setq betalow!* (iminus bbase!*)) 188 (setq bone!* (bnum 1)) 189 (setq bzero!* (bnum 0)) 190 (setq bigbetahi!* (bnum betahi!*)) 191 % Highest value of Ai 192 (setq bigbetalow!* (bminus bigbetahi!*)) 193 % Lowest value of Ai 194 % here assume 2's complement 195 (setq y (twopower (idifference x 2))) 196 % eg, 36 bits, 2^35-1=2^34+2^34-1 197 (setq syshi!* (plus y (difference y 1))) 198 (setq y (minus y)) 199 (setq syslow!* (plus y y)) 200 (setq bigsyshi!* (bdifference (btwopower (isub1 x)) bone!*)) 201 % Largest representable Syslisp integer. 202 % Note that SYSPOS has leading 0, ie only x-1 active bits 203 (setq bigsyslow!* (bminus (bplus2 bone!* bigsyshi!*))))) 204 205% Smallest representable Syslisp integer. 206(de nonbignumerror (v l) 207 (stderror (bldmsg " Expect a BIGNUM in %r, given %p%n" l v))) 208 209(de bsize (v) 210 % Upper Limit of [BIGxxx a1 ... An] 211 (if (bigp v) 212 (veclen (vecinf v)) 213 0)) 214 215(de gtpos (n) 216 % Allocate [BIGPOS a1 ... an] 217 (prog nil 218 (setq n (mkvect n)) 219 (iputv n 0 'bigpos) 220 (return (mkbign (vecinf n))))) 221 222(de gtneg (n) 223 % Allocate [BIGNEG a1 ... an] 224 (prog nil 225 (setq n (mkvect n)) 226 (iputv n 0 'bigneg) 227 (return (mkbign (vecinf n))))) 228 229(de trimbignum (v3) 230 % truncate trailing 0 231 (if (not (bigp v3)) 232 (nonbignumerror v3 'trimbignum) 233 (trimbignum1 v3 (bsize v3)))) 234 235(de trimbignum1 (b l3) 236 (prog (v3) 237 (setq v3 (bigasvec b)) 238 (while (and (igreaterp l3 0) (izerop (igetv v3 l3))) 239 (setq l3 (isub1 l3))) 240 (if (izerop (upbv (truncatevector v3 l3))) 241 (return (gtpos 0)) 242 (return b)))) 243 244(de bigasvec (b) 245 % In order to see BIGITS 246 (mkvec (inf b))) 247 248(de vecasbig (v) 249 (mkbign (vecinf v))) 250 251(de big2sys (u) 252 % Convert a BIG to SYS, if in range 253 (if (or (blessp u bigsyslow!*) (bgreaterp u bigsyshi!*)) 254 (continuableerror 99 "BIGNUM too large to convert to SYS" u) 255 (big2sysaux u))) 256 257(de big2sysaux (u) 258 % Convert a BIGN that is in range to a SYSINT 259 (prog (l sn res) 260 (setq l (bsize u)) 261 (when (izerop l) 262 (return 0)) 263 (setq res (igetv u l)) 264 (setq l (isub1 l)) 265 (if (bminusp u) 266 (progn (setq res (minus res)) 267 (while (neq l 0) 268 (setq res (itimes2 res bbase!*)) 269 (setq res (idifference res (igetv u l))) 270 (setq l (isub1 l))) 271 nil) 272 (while (neq l 0) 273 (setq res (itimes2 res bbase!*)) 274 (setq res (iplus2 res (igetv u l))) 275 (setq l (isub1 l)))) 276 (return res))) 277 278(de twopower (n) 279 %fix/i-num 2**n 280 (expt 2 n)) 281 282(de btwopower (n) 283 % gives 2**n; n is fix/i-num; result BigNum 284 (if (not (or (fixp n) (bigp n))) 285 (nonintegererror n 'btwopower) 286 (prog (quot rem v) 287 (when (bigp n) 288 (setq n (big2sys n))) 289 (setq quot (quotient n bbits!*)) 290 (setq rem (remainder n bbits!*)) 291 (setq v (gtpos (iadd1 quot))) 292 (ifor (from i 1 quot 1) (do (iputv v i 0))) 293 (iputv v (iadd1 quot) (twopower rem)) 294 (return (trimbignum1 v (iadd1 quot)))))) 295 296(de bzerop (v1) 297 (and (izerop (bsize v1)) (not (bminusp v1)))) 298 299(de bonep (v1) 300 (and (not (bminusp v1)) (ionep (bsize v1)) (ionep (igetv v1 1)))) 301 302(de babs (v1) 303 (if (bminusp v1) 304 (bminus v1) 305 v1)) 306 307(de bmax (v1 v2) 308 (if (bgreaterp v2 v1) 309 v2 310 v1)) 311 312(de bmin (v1 v2) 313 (if (blessp v2 v1) 314 v2 315 v1)) 316 317% (de bexpt (v1 n) 318% % V1 is Bignum, N is fix/i-num 319% (cond ((not (fixp n)) (nonintegererror n 'bexpt)) 320% ((izerop n) bone!*) 321% ((ionep n) v1) 322% ((iminusp n) (bquotient bone!* (bexpt v1 (iminus n)))) 323% (t (prog (v2) 324% (setq v2 (bexpt v1 (iquotient n 2))) 325% (if (izerop (iremainder n 2)) 326% (return (btimes2 v2 v2)) 327% (return (btimes2 (btimes2 v2 v1) v2))))))) 328 329% --------------------------------------- 330% Logical Operations 331% 332% All take Bignum arguments 333(de blor (v1 v2) 334 % The main body of the OR code is only obeyed when both arguments 335 % are positive, and so the result will be positive; 336 (if (or (bminusp v1) (bminusp v2)) 337 (blnot (bland (blnot v1) (blnot v2))) 338 (prog (l1 l2 l3 v3) 339 (setq l1 (bsize v1)) 340 (setq l2 (bsize v2)) 341 (when (greaterp l2 l1) 342 (setq l3 l2) 343 (setq l2 l1) 344 (setq l1 l3) 345 (setq v3 v2) 346 (setq v2 v1) 347 (setq v1 v3)) 348 (setq v3 (gtpos l1)) 349 (ifor (from i 1 l2 1) 350 (do (iputv v3 i (ilor (igetv v1 i) (igetv v2 i))))) 351 (ifor (from i (iadd1 l2) l1 1) (do (iputv v3 i (igetv v1 i)))) 352 (return v3)))) 353 354(de blxor (v1 v2) 355 % negative arguments are coped with using the identity 356 % LXor(a,b) = LNot LXor(Lnot a,b) = LNor LXor(a,Lnot b); 357 (prog (l1 l2 l3 v3 s) 358 (when (bminusp v1) 359 (setq v1 (blnot v1)) 360 (setq s t)) 361 (when (bminusp v2) 362 (setq v2 (blnot v2)) 363 (setq s (not s))) 364 (setq l1 (bsize v1)) 365 (setq l2 (bsize v2)) 366 (when (greaterp l2 l1) 367 (setq l3 l2) 368 (setq l2 l1) 369 (setq l1 l3) 370 (setq v3 v2) 371 (setq v2 v1) 372 (setq v1 v3)) 373 (setq v3 (gtpos l1)) 374 (ifor (from i 1 l2 1) 375 (do (iputv v3 i (ilxor (igetv v1 i) (igetv v2 i))))) 376 (ifor (from i (iadd1 l2) l1 1) (do (iputv v3 i (igetv v1 i)))) 377 (setq v1 (trimbignum1 v3 l1)) 378 (when s 379 (setq v1 (blnot v1))) 380 (return v1))) 381 382% Not Used Currently: 383% 384% procedure BLDiff(V1,V2); 385% ***** STILL NEEDS ADJUSTING WRT -VE ARGS ***** 386% begin scalar V3,L1,L2; 387% L1:=BSize V1; 388% L2:=BSize V2; 389% V3:=GtPOS(max(L1,L2)); 390% IFor i:=1:min(L1,L2) do 391% IPutV(V3,i,ILAnd(IGetV(V1,i),ILXor(LogicalBits!*,IGetV(V2,i)))); 392% if IGreaterP(L1,L2) then IFor i:=(IAdd1 L2):L1 do IPutV(V3,i,IGetV(V1,i)); 393 394% if IGreaterP(L2,L1) then IFor i:=(IAdd1 L1):L2 do IPutV(V3,i,0); 395% return TrimBigNum1(V3,max(L1,L2)); 396% end; 397(de bland (v1 v2) 398 % If both args are -ve the result will be too. Otherwise result will 399 % be positive; 400 (if (and (bminusp v1) (bminusp v2)) 401 (blnot (blor (blnot v1) (blnot v2))) 402 (prog (l1 l2 l3 v3) 403 (setq l1 (bsize v1)) 404 (setq l2 (bsize v2)) 405 (setq l3 (min l1 l2)) 406 (cond ((bminusp v1) 407 % When one is negative, we have expand out to the 408 % size of the other one. Therefore, we use l2 as the 409 % size, not l3. When we exceed the size of the 410 % negative number, then we just use logicalbits*. 411 (setq v3 (gtpos l2)) 412 (setq l3 l2) 413 (ifor (from i 1 l2 1) 414 (do 415 (iputv v3 i 416 (iland (cond ((igreaterp i l1) logicalbits*) 417 (t (ilxor (isub1 (igetv v1 i)) 418 logicalbits*))) 419 (igetv v2 i)))))) 420 ((bminusp v2) 421 (setq v3 (gtpos l1)) 422 (setq l3 l1) 423 (ifor (from i 1 l1 1) 424 (do 425 (iputv v3 i 426 (iland (igetv v1 i) 427 (cond ((igreaterp i l2) logicalbits*) 428 (t (ilxor (isub1 (igetv v2 i)) 429 logicalbits*)))))))) 430 431 (t (setq v3 (gtpos l3)) 432 (ifor (from i 1 l3 1) 433 (do 434 (iputv v3 i (iland (igetv v1 i) (igetv v2 i))))))) 435 (return (trimbignum1 v3 l3))))) 436 437(de blnot (v1) 438 (bminus (bsmalladd v1 1))) 439 440(de blshift (v1 v2) 441 % This seems a grimly inefficient way of doing things given that 442 % the representation of big numbers uses a base that is a power of 2. 443 % However it will do for now; 444 (if (bminusp v2) 445 (bquotient v1 (btwopower (bminus v2))) 446 (btimes2 v1 (btwopower v2)))) 447 448% ----------------------------------------- 449% Arithmetic Functions: 450% 451% U, V, V1, V2 are Bignum arguments. 452(de bminus (v1) 453 % Negates V1. 454 (if (bzerop v1) 455 v1 456 (prog (l1 v2) 457 (setq l1 (bsize v1)) 458 (if (bminusp v1) 459 (setq v2 (gtpos l1)) 460 (setq v2 (gtneg l1))) 461 (ifor (from i 1 l1 1) (do (iputv v2 i (igetv v1 i)))) 462 (return v2)))) 463 464% Returns V1 if V1 is strictly less than 0, NIL otherwise. 465% 466(de bminusp (v1) 467 (if (eq (igetv v1 0) 'bigneg) 468 v1 469 nil)) 470 471% To provide a conveninent ADD with CARRY. 472(de addcarry (a) 473 (prog (s) 474 (setq s (iplus2 a carry!*)) 475 (if (igeq s bbase!*) 476 (progn (setq carry!* 1) 477 (setq s (idifference s bbase!*))) 478 (setq carry!* 0)) 479 (return s))) 480 481(de bplus2 (v1 v2) 482 (prog (sn1 sn2) 483 (setq sn1 (bminusp v1)) 484 (setq sn2 (bminusp v2)) 485 (when (and sn1 (not sn2)) 486 (return (bdifference2 v2 (bminus v1) nil))) 487 (when (and sn2 (not sn1)) 488 (return (bdifference2 v1 (bminus v2) nil))) 489 (return (bplusa2 v1 v2 sn1)))) 490 491(de bplusa2 (v1 v2 sn1) 492 % Plus with signs pre-checked and 493 (prog (l1 l2 l3 v3 temp) 494 % identical. 495 (setq l1 (bsize v1)) 496 (setq l2 (bsize v2)) 497 (when (igreaterp l2 l1) 498 (setq l3 l2) 499 (setq l2 l1) 500 (setq l1 l3) 501 (setq v3 v2) 502 (setq v2 v1) 503 (setq v1 v3)) 504 (setq l3 (iadd1 l1)) 505 (if sn1 506 (setq v3 (gtneg l3)) 507 (setq v3 (gtpos l3))) 508 (setq carry!* 0) 509 (ifor (from i 1 l2 1) 510 (do (progn (setq temp (iplus2 (igetv v1 i) (igetv v2 i))) 511 (iputv v3 i (addcarry temp))))) 512 (setq temp (iadd1 l2)) 513 (ifor (from i temp l1 1) (do (iputv v3 i (addcarry (igetv v1 i))))) 514 (iputv v3 l3 carry!*) 515 % Carry Out 516 (return (trimbignum1 v3 l3)))) 517 518(de bdifference (v1 v2) 519 (cond ((bzerop v2) v1) 520 ((bzerop v1) (bminus v2)) 521 (t (prog (sn1 sn2) 522 (setq sn1 (bminusp v1)) 523 (setq sn2 (bminusp v2)) 524 (when (or (and sn1 (not sn2)) (and sn2 (not sn1))) 525 (return (bplusa2 v1 (bminus v2) sn1))) 526 (return (bdifference2 v1 v2 sn1)))))) 527 528(de subcarry (a) 529 (prog (s) 530 (setq s (idifference a carry!*)) 531 (if (ilessp s 0) 532 (progn (setq carry!* 1) 533 (setq s (iplus2 bbase!* s))) 534 (setq carry!* 0)) 535 (return s))) 536 537(de bdifference2 (v1 v2 sn1) 538 % Signs pre-checked and identical. 539 (prog (i l1 l2 l3 v3) 540 (setq l1 (bsize v1)) 541 (setq l2 (bsize v2)) 542 (cond ((igreaterp l2 l1) (setq l3 l1) (setq l1 l2) (setq l2 l3) 543 (setq v3 v1) (setq v1 v2) (setq v2 v3) (setq sn1 (not sn1))) 544 ((eq l1 l2) (setq i l1) (while (and 545 (eq (igetv v2 i) 546 (igetv v1 i)) 547 (igreaterp i 1)) 548 (setq i (isub1 i))) 549 (when (igreaterp (igetv v2 i) (igetv v1 i)) 550 (setq l3 l1) 551 (setq l1 l2) 552 (setq l2 l3) 553 (setq v3 v1) 554 (setq v1 v2) 555 (setq v2 v3) 556 (setq sn1 (not sn1))))) 557 (if sn1 558 (setq v3 (gtneg l1)) 559 (setq v3 (gtpos l1))) 560 (setq carry!* 0) 561 (ifor (from i 1 l2 1) 562 (do (iputv v3 i 563 (subcarry (idifference (igetv v1 i) (igetv v2 i)))))) 564 (ifor (from i (iadd1 l2) l1 1) 565 (do (iputv v3 i (subcarry (igetv v1 i))))) 566 (return (trimbignum1 v3 l1)))) 567 568(de btimes2 (v1 v2) 569 (prog (l1 l2 l3 sn1 sn2 v3) 570 (setq l1 (bsize v1)) 571 (setq l2 (bsize v2)) 572 (when (igreaterp l2 l1) 573 (setq v3 v1) 574 (setq v1 v2) 575 (setq v2 v3) 576 % If V1 is larger, will be fewer 577 (setq l3 l1) 578 (setq l1 l2) 579 (setq l2 l3)) 580 % iterations of BDigitTimes2. 581 (setq l3 (iplus2 l1 l2)) 582 (setq sn1 (bminusp v1)) 583 (setq sn2 (bminusp v2)) 584 (if (or (and sn1 sn2) (not (or sn1 sn2))) 585 (setq v3 (gtpos l3)) 586 (setq v3 (gtneg l3))) 587 (ifor (from i 1 l3 1) (do (iputv v3 i 0))) 588 (ifor (from i 1 l2 1) (do (bdigittimes2 v1 (igetv v2 i) l1 i v3))) 589 (return (trimbignum1 v3 l3)))) 590 591(de bdigittimes2 (v1 v2 l1 i v3) 592 % V1 is a bignum, V2 a fixnum, L1=BSize L1, I=position of V2 in a bignum, 593 % and V3 is bignum receiving result. I affects where in V3 the result of 594 % a calculation goes; the relationship is that positions I:I+(L1-1) 595 % of V3 receive the products of V2 and positions 1:L1 of V1. 596 % V3 is changed as a side effect here. 597 (prog (j carry temp1 temp2) 598 (if (zerop v2) 599 (return v3) 600 (progn (setq carry 0) 601 (ifor (from h 1 l1 1) 602 (do 603 (progn (setq temp1 (itimes2 (igetv v1 h) v2)) 604 (setq temp2 (iplus2 h (isub1 i))) 605 (setq j 606 (iplus2 (iplus2 temp1 (igetv v3 temp2)) 607 carry)) 608 (iputv v3 temp2 (iremainder j bbase!*)) 609 (setq carry (iquotient j bbase!*))))) 610 (iputv v3 (iplus2 l1 i) carry))) 611 % carry should be < BBase!* here 612 (return v3))) 613 614(de bsmalltimes2 (v1 c) 615 % V1 is a BigNum, C a fixnum. 616 % Assume C positive, ignore sign(V1) 617 % also assume V1 neq 0. 618 (cond ((zerop c) (return (gtpos 0))) 619 (t % Only used from BHardDivide, BReadAdd. 620 (prog (j carry l1 l2 l3 v3) 621 (setq l1 (bsize v1)) 622 (setq l2 (iplus2 (iquotient c bbase!*) l1)) 623 (setq l3 (iadd1 l2)) 624 (setq v3 (gtpos l3)) 625 (setq carry 0) 626 (ifor (from h 1 l1 1) 627 (do 628 (progn (setq j 629 (iplus2 (itimes2 (igetv v1 h) c) carry)) 630 (iputv v3 h (iremainder j bbase!*)) 631 (setq carry (iquotient j bbase!*))))) 632 (ifor (from h (iadd1 l1) l3 1) 633 (do 634 (progn (iputv v3 h 635 (iremainder (setq j carry) bbase!*)) 636 (setq carry (iquotient j bbase!*))))) 637 (return (trimbignum1 v3 l3)))))) 638 639(de bquotient (v1 v2) 640 (car (bdivide v1 v2))) 641 642(de bremainder (v1 v2) 643 (cdr (bdivide v1 v2))) 644 645% BDivide returns a dotted pair, (Q . R). Q is the quotient and R is 646% the remainder. Both are bignums. R is of the same sign as V1. 647%; 648(ds bsimplequotient (v1 l1 c snc) (car (bsimpledivide v1 l1 c snc))) 649 650(ds bsimpleremainder (v1 l1 c snc) (cdr (bsimpledivide v1 l1 c snc))) 651 652(de bdivide (v1 v2) 653 (prog (l1 l2 q r v3) 654 (setq l2 (bsize v2)) 655 (when (izerop l2) 656 (error 99 "Attempt to divide by 0 in BDIVIDE")) 657 (setq l1 (bsize v1)) 658 (cond ((or (ilessp l1 l2) 659 (and (eq l1 l2) (ilessp (igetv v1 l1) (igetv v2 l2)))) 660 % This also takes care of case 661 (return (cons (gtpos 0) v1)))) 662 % when V1=0. 663 (when (ionep l2) 664 (return (bsimpledivide v1 l1 (igetv v2 1) (bminusp v2)))) 665 (return (bharddivide v1 l1 v2 l2)))) 666 667% C is a fixnum (inum?); V1 is a bignum and L1 is its length. 668% SnC is T if C (which is positive) should be considered negative. 669% Returns quotient . remainder; each is a bignum. 670% 671(de bsimpledivide (v1 l1 c snc) 672 (prog (i p r rr sn1 v2) 673 (setq sn1 (bminusp v1)) 674 (if (or (and sn1 snc) (not (or sn1 snc))) 675 (setq v2 (gtpos l1)) 676 (setq v2 (gtneg l1))) 677 (setq r 0) 678 (setq i l1) 679 (while (not (izerop i)) 680 (setq p (iplus2 (itimes2 r bbase!*) (igetv v1 i))) 681 % Overflow. 682 (iputv v2 i (iquotient p c)) 683 (setq r (iremainder p c)) 684 (setq i (isub1 i))) 685 (setq p 0) % In case we have to GC. JHD on advice from Mark Swanson 686 (if sn1 687 (setq rr (gtneg 1)) 688 (setq rr (gtpos 1))) 689 (iputv rr 1 r) 690 (return (cons (trimbignum1 v2 l1) (trimbignum1 rr 1))))) 691 692(de bharddivide (u lu v lv) 693 % This is an algorithm taken from Knuth. 694 (prog (u1 v1 a d lcv lcv1 f f2 j k lq carry temp ll m n n1 p q qbar snu 695 snv u2) 696 (setq n lv) 697 (setq n1 (iadd1 n)) 698 (setq m (idifference lu lv)) 699 (setq lq (iadd1 m)) 700 % Deal with signs of inputs; 701 (setq snu (bminusp u)) 702 (setq snv (bminusp v)) 703 % Note that these are not extra-boolean, i.e. 704 % for positive numbers MBinusP returns nil, for 705 % negative it returns its argument. Thus the 706 % test (SnU=SnV) does not reliably compare the signs of 707 % U and V; 708 (cond (snu (if snv 709 (setq q (gtpos lq)) 710 (setq q (gtneg lq)))) 711 (snv (setq q (gtneg lq))) 712 (t (setq q (gtpos lq)))) 713 (setq u1 (gtpos (iadd1 lu))) 714 % U is ALWAYS stored as if one digit longer; 715 % Compute a scale factor to normalize the long division; 716 (setq d (iquotient bbase!* (iadd1 (igetv v lv)))) 717 % Now, at the same time, I remove the sign information from U and V 718 % and scale them so that the leading coefficeint in V is fairly large; 719 720 (setq carry 0) 721 (ifor (from i 1 lu 1) 722 (do (progn (setq temp (iplus2 (itimes2 (igetv u i) d) carry)) 723 (iputv u1 i (iremainder temp bbase!*)) 724 (setq carry (iquotient temp bbase!*))))) 725 (setq lu (iadd1 lu)) 726 (iputv u1 lu carry) 727 % temp can have a >24 bit value, which is dangerous on 68000's 728 (setq temp 0) 729 (setq v1 (bsmalltimes2 v d)) 730 % So far all variables contain safe values, 731 % i.e. numbers < BBase!*; 732 (iputv v1 0 'bigpos) 733 (when (ilessp lv 2) 734 (nonbignumerror v 'bharddivide)) 735 % To be safe; 736 (setq lcv (igetv v1 lv)) 737 (setq lcv1 (igetv v1 (isub1 lv))) 738 % Top two digits of the scaled V accessed once 739 % here outside the main loop; 740 % Now perform the main long division loop; 741 (ifor (from i 0 m 1) 742 (do (progn (setq j (idifference lu i)) 743 % J>K; working on U1[K:J] 744 (setq k (idifference j n1)) 745 % in this loop. 746 (setq a (igetv u1 j)) 747 (setq p 748 (iplus2 (itimes2 a bbase!*) (igetv u1 (isub1 j)))) 749 % N.B. P is up to 30 bits long. Take care! ; 750 (if (eq a lcv) 751 (setq qbar (isub1 bbase!*)) 752 (setq qbar (iquotient p lcv))) 753 % approximate next digit; 754 (setq f (itimes2 qbar lcv1)) 755 (setq f2 756 (iplus2 757 (itimes2 (idifference p (itimes2 qbar lcv)) 758 bbase!*) 759 (igetv u1 (idifference j 2)))) 760 (while (igreaterp f f2) 761 % Correct most overshoots in Qbar; 762 (setq qbar (isub1 qbar)) 763 (setq f (idifference f lcv1)) 764 nil 765 (setq f2 (iplus2 f2 (itimes2 lcv bbase!*)))) 766 (setq carry 0) 767 % Ready to subtract QBar*V1 from U1; 768 (ifor (from l 1 n 1) 769 (do 770 (progn (setq temp 771 (iplus2 772 (idifference (igetv u1 (iplus2 k l)) 773 (itimes2 qbar (igetv v1 l))) 774 carry)) 775 (setq carry (iquotient temp bbase!*)) 776 (setq temp (iremainder temp bbase!*)) 777 (when (iminusp temp) 778 (setq carry (isub1 carry)) 779 (setq temp (iplus2 temp bbase!*))) 780 (iputv u1 (iplus2 k l) temp)))) 781 % Now propagate borrows up as far as they go; 782 (setq ll (iplus2 k n)) 783 (while (and (not (izerop carry)) (ilessp ll j)) 784 (setq ll (iadd1 ll)) 785 (setq temp (iplus2 (igetv u1 ll) carry)) 786 (setq carry (iquotient temp bbase!*)) 787 (setq temp (iremainder temp bbase!*)) 788 (when (iminusp temp) 789 (setq carry (isub1 carry)) 790 (setq temp (iplus2 temp bbase!*))) 791 (iputv u1 ll temp)) 792 (unless (izerop carry) 793 % QBar was still wrong - correction step needed. 794 % This should not happen very often; 795 (setq qbar (isub1 qbar)) 796 % Add V1 back into U1; 797 (setq carry 0) 798 (ifor (from l 1 n 1) 799 (do 800 (progn (setq carry 801 (iplus2 802 (iplus2 (igetv u1 (iplus2 k l)) 803 (igetv v1 l)) 804 carry)) 805 (iputv u1 (iplus2 k l) 806 (iremainder carry bbase!*)) 807 (setq carry (iquotient carry bbase!*))))) 808 (setq ll (iplus2 k n)) 809 (while (ilessp ll j) 810 (setq ll (iadd1 ll)) 811 (setq carry (iplus2 (igetv u1 ll) carry)) 812 (iputv u1 ll (iremainder carry bbase!*)) 813 (setq carry (iquotient carry bbase!*)))) 814 (iputv q (idifference lq i) qbar)))) 815 % End of main loop; 816 % zero out potentially dangerous values 817 (setq p 0) 818 (setq f 0) 819 (setq f2 0) 820 (setq u1 (trimbignum1 u1 (idifference lu m))) 821 % Clean up potentially wild values; 822 (unless (bzerop u1) 823 % Unnormalize the remainder by dividing by D 824 (when snu 825 (iputv u1 0 'bigneg)) 826 (unless (ionep d) 827 (setq lu (bsize u1)) 828 (setq carry 0) 829 (ifor (from l lu 1 -1) 830 (do (progn (setq p 831 (iplus2 (itimes2 carry bbase!*) (igetv u1 l))) 832 (iputv u1 l (iquotient p d)) 833 (setq carry (iremainder p d))))) 834 (setq p 0) 835 (unless (izerop carry) 836 (bhardbug "remainder when unscaling" u v (trimbignum1 u1 lu) 837 (trimbignum1 q lq))) 838 (setq u1 (trimbignum1 u1 lu)))) 839 (setq q (trimbignum1 q lq)) 840 % In case leading digit happened to be zero; 841 (setq p 0) 842 % flush out a 30 bit number; 843 % Here, for debugging purposes, I will try to validate the results I 844 845 % have obtained by testing if Q*V+U1=U and 0<=U1<V. I Know this slows things 846 847 % down, but I will remove it when my confidence has improved somewhat; 848 849 % if not BZerop U1 then << 850 % if (BMinusP U and not BMinusP U1) or 851 % (BMinusP U1 and not BMinusP U) then 852 % BHardBug("remainder has wrong sign",U,V,U1,Q) >>; 853 854 % if not BAbs U1<BAbs V then BHardBug("remainder out of range",U,V,U1,Q) 855 856 % else if not BZerop(BDifference(BPlus2(BTimes2(Q,V),U1),U)) then 857 858 % BHardBug("quotient or remainder incorrect",U,V,U1,Q); 859 (return (cons q u1)))) 860 861(de bhardbug (msg u v r q) 862 % Because the inputs to BHardDivide are probably rather large, I am not 863 % going to rely on BldMsg to display them; 864 (progn (prin2t "***** Internal error in BHardDivide") 865 (prin2 "arg1=") 866 (prin2t u) 867 (prin2 "arg2=") 868 (prin2t v) 869 (prin2 "computed quotient=") 870 (prin2t q) 871 (prin2 "computed remainder=") 872 (prin2t r) 873 (stderror msg))) 874 875(de bgreaterp (u v) 876 (cond ((bminusp u) (if (bminusp v) 877 (bunsignedgreaterp v u) 878 nil)) 879 ((bminusp v) u) 880 (t (bunsignedgreaterp u v)))) 881 882(de blessp (u v) 883 (cond ((bminusp u) (if (bminusp v) 884 (bunsignedgreaterp u v) 885 u)) 886 ((bminusp v) nil) 887 (t (bunsignedgreaterp v u)))) 888 889(de bgeq (u v) 890 (cond ((bminusp u) (if (bminusp v) 891 (bunsignedgeq v u) 892 nil)) 893 ((bminusp v) u) 894 (t (bunsignedgeq u v)))) 895 896(de bleq (u v) 897 (cond ((bminusp u) (if (bminusp v) 898 (bunsignedgeq u v) 899 u)) 900 ((bminusp v) nil) 901 (t (bunsignedgeq v u)))) 902 903(de bunsignedgreaterp (u v) 904 % Compare magnitudes of two bignums; 905 (prog (lu lv i) 906 (setq lu (bsize u)) 907 (setq lv (bsize v)) 908 (unless (eq lu lv) 909 (if (igreaterp lu lv) 910 (return u) 911 (return nil))) 912 (while (and (eq (igetv u lv) (igetv v lv)) (igreaterp lv 1)) 913 (setq lv (isub1 lv))) 914 (if (igreaterp (igetv u lv) (igetv v lv)) 915 (return u) 916 (return nil)))) 917 918(de bunsignedgeq (u v) 919 % Compare magnitudes of two unsigned bignums; 920 (prog (lu lv) 921 (setq lu (bsize u)) 922 (setq lv (bsize v)) 923 (unless (eq lu lv) 924 (if (igreaterp lu lv) 925 (return u) 926 (return nil))) 927 (while (and (eq (igetv u lv) (igetv v lv)) (igreaterp lv 1)) 928 (setq lv (isub1 lv))) 929 (if (igreaterp (igetv v lv) (igetv u lv)) 930 (return nil) 931 (return u)))) 932 933(de badd1 (v) 934 (bsmalladd v 1)) 935 936(de bsub1 (u) 937 (bsmalldiff u 1)) 938 939% ------------------------------------------------ 940% Conversion to Float: 941(de floatfrombignum (v) 942 (cond ((bzerop v) 0.00000E+000) 943 ((or (bgreaterp v bigfloathi!*) (blessp v bigfloatlow!*)) 944 (error 99 (list "Argument, " v " to FLOAT is too large"))) 945 (t (prog (l res sn i) 946 % Careful, do not want to call itself recursively 947 (setq l (bsize v)) 948 (setq sn (bminusp v)) 949 (setq res (intfloat (igetv v l))) 950 (setq i (isub1 l)) 951 (while (not (izerop i)) 952 (setq res (floattimes2 res floatbbase!*)) 953 (setq res (floatplus2 res (intfloat (igetv v i)))) 954 (setq i (isub1 i))) 955 (when sn 956 % used to be minus instead of floatdifference: 957 % that causes a recursion 958 (setq res (floatdifference 0.0 res))) 959 (return res))))) 960 961% ------------------------------------------------ 962% Input and Output: 963(setq digit2letter!* 964 '[48 49 50 51 52 53 54 55 56 57 65 66 67 68 69 70 71 72 73 74 75 76 965 77 78 79 80 81 82 83 84 85 86 87 88 89 90]) 966 967% a new routine for printing bignums 968% 969% bchannelprin2 is still the entry point of course, but not much else 970% 971% first make a copy of v1 (the number to be printed), output any 972% necessary base information then call bprin which actually does the work 973 974(de bchannelprin2 (channel v1) 975 ((lambda (v2 myobase i) 976 (while (not (izerop i)) 977 (progn (iputv v2 i (igetv v1 i)) (setq i (isub1 i)))) 978 (cond ((bminusp v1) (channelwritechar channel (char !-)))) 979 (cond 980 ((neq myobase 10) 981 (channelwritesysinteger channel myobase 10) 982 (channelwritechar channel (char !#)))) 983 (bprin channel v2) 984 (setq outputbase!* myobase)) 985 (gtpos (bsize v1)) %JAP 870802 986 outputbase!* 987 (bsize v1))) 988 989% divide v by 10**4, if the quotient is zero the recursion is terminated 990% and the most significant digit is printed (and so on down the stack) 991% otherwise save the big digit (=remainder) and recurse 992(de bprin (channel v) 993 ((lambda (digit) 994 (cond 995 ((bzerop v) (channelwritesysinteger channel digit outputbase!*)) 996 (t (bprin channel v) (printdigit channel digit)))) 997 (bdivide10000ip v))) 998 999% print a base 10**4 digit 1000(de printdigit (channel d) 1001 (prog (d1 d2 d3) 1002 (setq d1 (iremainder d 10)) (setq d (iquotient d 10)) 1003 (setq d2 (iremainder d 10)) (setq d (iquotient d 10)) 1004 (setq d3 (iremainder d 10)) (setq d (iquotient d 10)) 1005 (channelwritechar channel (igetv digit2letter!* (iremainder d 10))) 1006 (channelwritechar channel (igetv digit2letter!* d3)) 1007 (channelwritechar channel (igetv digit2letter!* d2)) 1008 (return (channelwritechar channel (igetv digit2letter!* d1))))) 1009 1010% divide the bignum v1 (of length l1) by 10000, except the quotient is 1011% accumulated in the same place, the remainder, of course, ripples 1012% down to the bottom. Because the argument is modified there is no need 1013% to CONS up a result, but simply return the remainder. 1014(de bdivide10000ip (v1) 1015 ((lambda (p r i) 1016 (while (not (izerop i)) 1017 (progn 1018 (setq p (iplus2 (itimes2 r bbase!*) (igetv v1 i))) 1019 (iputv v1 i (iquotient p 10000)) 1020 (setq r (iremainder p 10000)) 1021 (setq i (isub1 i)))) 1022 (trimbignum1 v1 (bsize v1)) 1023 (cond ((bminusp v1) (iminus r)) (t r))) 1024 nil 1025 0 1026 (bsize v1))) 1027 1028(de bread (s radix sn) 1029 % radix is < Bbase!* 1030 %s=string of digits, radix=base, sn=1 or -1 1031 (prog (sz res ch) 1032 (setq sz (size s)) 1033 (setq res (gtpos 1)) 1034 (setq ch (indx s 0)) 1035 (when (and (igeq ch (char a)) (ileq ch (char z))) 1036 (setq ch (iplus2 (idifference ch (char a)) 10))) 1037 (when (and (igeq ch (char 0)) (ileq ch (char 9))) 1038 (setq ch (idifference ch (char 0)))) 1039 (iputv res 1 ch) 1040 (ifor (from i 1 sz 1) 1041 (do (progn (setq ch (indx s i)) 1042 (when (and (igeq ch (char a)) (ileq ch (char z))) 1043 (setq ch 1044 (idifference ch (idifference (char a) 10)))) 1045 (when (and (igeq ch (char 0)) (ileq ch (char 9))) 1046 (setq ch (idifference ch (char 0)))) 1047 (setq res (breadadd res radix ch))))) 1048 (when (iminusp sn) 1049 (setq res (bminus res))) 1050 (return res))) 1051 1052(de breadadd (v radix ch) 1053 (setq v (bsmalltimes2 v radix)) 1054 (setq v (bsmalladd v ch))) 1055 1056(de bsmalladd (v c) 1057 %V big, C fix. 1058 (cond ((izerop c) (return v)) 1059 ((bzerop v) (return (int2big c))) 1060 ((bminusp v) (bminus (bsmalldiff (bminus v) c))) 1061 ((iminusp c) (bsmalldiff v (iminus c))) 1062 (t (prog (v1 l1) 1063 (setq carry!* c) 1064 (setq l1 (bsize v)) 1065 (setq v1 (gtpos (iadd1 l1))) 1066 (ifor (from i 1 l1 1) 1067 (do (iputv v1 i (addcarry (igetv v i))))) 1068 (if (ionep carry!*) 1069 (iputv v1 (iadd1 l1) 1) 1070 (return (trimbignum1 v1 l1))) 1071 (return v1))))) 1072 1073(de bnum (n) 1074 % Creates a Bignum of one BETA digit, value N. 1075 % N is POS or NEG 1076 (if (bigp n) 1077 n 1078 (bnumaux n))) 1079 1080(de bnumaux (n) 1081 % Creates a Bignum of one BIGIT value N. 1082 % N is POS or NEG 1083 (prog (b) 1084 (cond ((izerop n) (return (gtpos 0))) 1085 ((iminusp n) (setq b (gtneg 1)) (setq n (iminus n))) 1086 (t (setq b (gtpos 1)))) 1087 (iputv b 1 n) 1088 (return b))) 1089 1090(de bsmalldiff (v c) 1091 %V big, C fix 1092 (cond ((izerop c) v) 1093 ((bzerop v) (int2big (iminus c))) 1094 ((bminusp v) (bminus (bsmalladd (bminus v) c))) 1095 ((iminusp c) (bsmalladd v (iminus c))) 1096 (t (prog (v1 l1) 1097 (setq carry!* c) 1098 (setq l1 (bsize v)) 1099 (setq v1 (gtpos l1)) 1100 (ifor (from i 1 l1 1) 1101 (do (iputv v1 i (subcarry (igetv v i))))) 1102 (unless (izerop carry!*) 1103 (stderror (bldmsg " BSmallDiff V<C %p %p%n" v c))) 1104 (return (trimbignum1 v1 l1)))))) 1105 1106(de int2big (n) 1107 % Creates BigNum of value N. 1108 % From any N, BETA,INUM,FIXNUM or BIGNUM 1109 (case (tag n) ((negint-tag posint-tag) (sys2big n)) 1110 ((fixnum-tag) (sys2big (fixval (fixinf n)))) 1111 ((bignum-tag) n) 1112 (nil (nonintegererror n 'int2big)))) 1113 1114% Convert BIGNUMs to FLOAT 1115(de bigfromfloat (x) 1116 (if (or (fixp x) (bigp x)) 1117 x 1118 (prog (bigpart floatpart power sign thispart) 1119 (if (minusp x) 1120 (progn (setq sign -1) 1121 (setq x (minus x))) 1122 (setq sign 1)) 1123 (setq bigpart bzero!*) 1124 (while (and (neq x 0) (neq x 0.00000E+000)) 1125 (if (lessp x bbase!*) 1126 (progn (setq bigpart (bplus2 bigpart (bnum (fix x)))) 1127 (setq x 0)) 1128 (progn (setq floatpart x) 1129 (setq power 0) 1130 (while (geq floatpart bbase!*) 1131 % get high end of number. 1132 (progn (setq floatpart (quotient floatpart bbase!*)) 1133 (setq power (plus power bbits!*)))) 1134 (setq thispart 1135 (btimes2 (btwopower power) (bnum (fix floatpart)))) 1136 (setq x (difference x (floatfrombignum thispart))) 1137 (setq bigpart (bplus2 bigpart thispart))))) 1138 (when (minusp sign) 1139 (setq bigpart (bminus bigpart))) 1140 (return bigpart)))) 1141 1142% Now Install Interfacing 1143(de setupglobals () 1144 (prin2t '"SetupGlobals") 1145 (setbits bitsperword) 1146 (prin2t '" ... done") 1147 nil) 1148 1149(setupglobals) 1150 1151(loadtime 1152 (progn (setq staticbig!* (gtwarray 10)))) 1153 1154% Assume dont need more than 10 slots to represent a BigNum 1155% Version of SYSint 1156% -- Output--- 1157% MLG Change to interface to Recursive hooks, added for 1158% Prinlevel stuff 1159(copyd 'oldchannelprin1 'recursivechannelprin1) 1160 1161(copyd 'oldchannelprin2 'recursivechannelprin2) 1162 1163(de recursivechannelprin1 (channel u level) 1164 (if (bigp u) 1165 (bchannelprin2 channel u) 1166 (oldchannelprin1 channel u level)) 1167 u) 1168 1169(de recursivechannelprin2 (channel u level) 1170 (if (bigp u) 1171 (bchannelprin2 channel u) 1172 (oldchannelprin2 channel u level)) 1173 u) 1174 1175(de checkifreallybig (uu) 1176 % If BIGNUM result is in older FIXNUM or INUM range 1177 % Convert Back. 1178 %/ Need a faster test 1179 (if (or (blessp uu bigsyslow!*) (bgreaterp uu bigsyshi!*)) 1180 uu 1181 (sys2int (big2sysaux uu)))) 1182 1183(de checkifreallybigpair (vv) 1184 % Used to process DIVIDE 1185 (cons (checkifreallybig (car vv)) (checkifreallybig (cdr vv)))) 1186 1187(de checkifreallybigornil (uu) 1188 % Used for EXTRA-boolean tests 1189 (if (or (null uu) (blessp uu bigsyslow!*) (bgreaterp uu bigsyshi!*)) 1190 uu 1191 (sys2int (big2sysaux uu)))) 1192 1193(de bigplus2 (u v) 1194 (checkifreallybig (bplus2 u v))) 1195 1196(de bigdifference (u v) 1197 (checkifreallybig (bdifference u v))) 1198 1199(de bigtimes2 (u v) 1200 (checkifreallybig (btimes2 u v))) 1201 1202(de bigdivide (u v) 1203 (checkifreallybigpair (bdivide u v))) 1204 1205(de bigquotient (u v) 1206 (checkifreallybig (bquotient u v))) 1207 1208(de bigremainder (u v) 1209 (checkifreallybig (bremainder u v))) 1210 1211(de bigland (u v) 1212 (checkifreallybig (bland u v))) 1213 1214(de biglor (u v) 1215 (checkifreallybig (blor u v))) 1216 1217(de biglxor (u v) 1218 (checkifreallybig (blxor u v))) 1219 1220(de biglshift (u v) 1221 (checkifreallybig (blshift u v))) 1222 1223(de lshift (u v) 1224 (if (and (betap u) (betap v)) 1225 (cond ((wlessp v 0) 1226 (sys2int (if_system VAX (intlshift u v) (wshift u v)))) 1227 ((wlessp v bbits!*) 1228 (sys2int (if_system VAX (intlshift u v) (wshift u v)))) 1229 (t (biglshift (sys2big u) (sys2big v)))) 1230 % Use int2big, not sys2big, since we might have fixnums. 1231 (biglshift (int2big u) (int2big v)))) 1232 1233(copyd 'lsh 'lshift) 1234 1235(de biggreaterp (u v) 1236 (checkifreallybigornil (bgreaterp u v))) 1237 1238(de biglessp (u v) 1239 (checkifreallybigornil (blessp u v))) 1240 1241(de bigadd1 (u) 1242 (checkifreallybig (badd1 u))) 1243 1244(de bigsub1 (u) 1245 (checkifreallybig (bsub1 u))) 1246 1247(de biglnot (u) 1248 (checkifreallybig (blnot u))) 1249 1250(de bigminus (u) 1251 (checkifreallybig (bminus u))) 1252 1253(de bigminusp (u) 1254 (checkifreallybigornil (bminusp u))) 1255 1256(de bigonep (u) 1257 (checkifreallybigornil (bonep u))) 1258 1259(de bigzerop (u) 1260 (checkifreallybigornil (bzerop u))) 1261 1262% ---- Input ---- 1263(de makestringintolispinteger (s radix sn) 1264 (checkifreallybig (bread s radix sn))) 1265 1266(de int2sys (n) 1267 % Convert a random FIXed number to WORD Integer 1268 (case (tag n) ((posint-tag negint-tag) n) 1269 ((fixnum-tag) (fixval (fixinf n))) 1270 ((bignum-tag) (big2sysaux n)) 1271 (nil (nonnumber1error n 'int2sys)))) 1272 1273(de sys2big (n) 1274 % Convert a SYSint to a BIG 1275 % Must NOT use generic arith here 1276 % Careful that no GC if this BIGger than INUM 1277 (prog (sn a b) 1278 (when (weq n 0) 1279 (return (gtpos 0))) 1280 (setq a staticbig!*) 1281 % Grab the base 1282 (when (wlessp n 0) 1283 (setq sn t)) 1284 (setf (wgetv a 1) n) 1285 % Plant number 1286 (setq n 1) 1287 % now use N as counter 1288 % Careful handling of -N in case have largest NEG, not just 1289 % flip sign 1290 (if sn 1291 (progn (setq b (wminus bbase!*)) 1292 (while (wleq (wgetv a n) b) 1293 (setq n (wplus2 n 1)) 1294 (setf (wgetv a n) 1295 (wquotient (wgetv a (wdifference n 1)) bbase!*)) 1296 (setf (wgetv a (wdifference n 1)) 1297 (wdifference (wgetv a (wdifference n 1)) 1298 (wtimes2 (wgetv a n) bbase!*)))) 1299 (setq b (gtneg n)) 1300 (ifor (from i 1 n 1) 1301 (do (iputv b i (wminus (wgetv a i)))))) 1302 (progn (while (wgeq (wgetv a n) bbase!*) 1303 (setq n (wplus2 n 1)) 1304 (setf (wgetv a n) 1305 (wquotient (wgetv a (wdifference n 1)) bbase!*)) 1306 (setf (wgetv a (wdifference n 1)) 1307 (wdifference (wgetv a (wdifference n 1)) 1308 (wtimes2 (wgetv a n) bbase!*)))) 1309 (setq b (gtpos n)) 1310 (ifor (from i 1 n 1) (do (iputv b i (wgetv a i)))))) 1311 (return b))) 1312 1313% Coercion/Transfer Functions 1314(copyd 'oldfloatfix 'floatfix) 1315 1316(de floatfix (u) 1317 % Careful of sign and range 1318 (if (and (leq floatsyslow!* u) (leq u floatsyshi!*)) 1319 (oldfloatfix u) 1320 (bigfromfloat u))) 1321 1322(de betap (x) 1323 % test if NUMBER in reduced INUM range 1324 (if (intp x) 1325 (and (wleq x betahi!*) (wgeq x betalow!*)) 1326 nil)) 1327 1328(de betarangep (x) 1329 % Test if SYSINT in reduced INUM range 1330 (if (wleq x betahi!*) 1331 (wgeq x betalow!*) 1332 nil)) 1333 1334(de beta2p (x y) 1335 % Check for 2 argument arithmetic functions 1336 (when (betap x) 1337 (betap y))) 1338 1339% Also from .BUILD file 1340% Now install the important globals for this machine 1341 1342(if_system VAX 1343 (progn 1344 % Largest representable float. 1345 (setq BigFloatHi!* (btimes2 (bdifference (btwopower 67) (btwopower 11)) 1346 (btwopower 60))) 1347 (setq BigFloatLow!* (bminus BigFloatHi!*)))) 1348 1349(if_system SPARC 1350 (progn 1351 % HP9836 sizes, range 10^-308 .. 10 ^308 1352 % I guess: 10^308 = 2 ^1025 // 15.8 digits, IEEE double ~56 bits 1353 % Largest representable float. 1354 (setq BigFloatHi!* (btimes2 (bsub1 (btwopower 56)) 1355 (btwopower 961))) 1356 (setq BigFloatLow!* (bminus BigFloatHi!*)))) 1357 1358(if_system MC68000 1359 (progn 1360 % HP9836 sizes, range 10^-308 .. 10 ^308 1361 % I guess: 10^308 = 2 ^1025 // 15.8 digits, IEEE double ~56 bits 1362 % Largest representable float. 1363 (setq BigFloatHi!* (btimes2 (bsub1 (btwopower 56)) 1364 (btwopower 961))) 1365 (setq BigFloatLow!* (bminus BigFloatHi!*)))) 1366 1367(if_system MC88000 1368 (progn 1369 % HP9836 sizes, range 10^-308 .. 10 ^308 1370 % I guess: 10^308 = 2 ^1025 // 15.8 digits, IEEE double ~56 bits 1371 % Largest representable float. 1372 (setq BigFloatHi!* (btimes2 (bsub1 (btwopower 56)) 1373 (btwopower 961))) 1374 (setq BigFloatLow!* (bminus BigFloatHi!*)))) 1375 1376(if_system Convex 1377 (progn 1378 % HP9836 sizes, range 10^-308 .. 10 ^308 1379 % I guess: 10^308 = 2 ^1025 // 15.8 digits, IEEE double ~56 bits 1380 % Largest representable float. 1381 (setq BigFloatHi!* (btimes2 (bsub1 (btwopower 56)) 1382 (btwopower 961))) 1383 (setq BigFloatLow!* (bminus BigFloatHi!*)))) 1384 1385(if_system Mips 1386 (progn 1387 % HP9836 sizes, range 10^-308 .. 10 ^308 1388 % I guess: 10^308 = 2 ^1025 // 15.8 digits, IEEE double ~56 bits 1389 % Largest representable float. 1390 (setq BigFloatHi!* (btimes2 (bsub1 (btwopower 56)) 1391 (btwopower 961))) 1392 (setq BigFloatLow!* (bminus BigFloatHi!*)))) 1393 1394(if_system IBMRS 1395 (progn 1396 % HP9836 sizes, range 10^-308 .. 10 ^308 1397 % I guess: 10^308 = 2 ^1025 // 15.8 digits, IEEE double ~56 bits 1398 % Largest representable float. 1399 (setq BigFloatHi!* (btimes2 (bsub1 (btwopower 56)) 1400 (btwopower 961))) 1401 (setq BigFloatLow!* (bminus BigFloatHi!*)))) 1402 1403(if_system HP-RISC 1404 (progn 1405 % HP9836 sizes, range 10^-308 .. 10 ^308 1406 % I guess: 10^308 = 2 ^1025 // 15.8 digits, IEEE double ~56 bits 1407 % Largest representable float. 1408 (setq BigFloatHi!* (btimes2 (bsub1 (btwopower 56)) 1409 (btwopower 961))) 1410 (setq BigFloatLow!* (bminus BigFloatHi!*)))) 1411 1412(if_system PDP10 1413 (progn 1414 (setq BigFloatHi!* (btimes2 (bsub1 (btwopower 62)) (btwopower 65))) 1415 (setq BigFloatLow!* (bminus BigFloatHi!*)))) 1416 1417(setq FloatSysHi!* (float SysHi!*)) 1418(setq FloatSysLow!* (float SysLow!*)) 1419 1420% The following patch supplied by Herbert Melenk and Winfried Neun of 1421% the Konrad-Zuse Zentrum, Berlin, corrects a serious bug in the 1422% garbage collector that occurs when bignums are loaded. The patch 1423% must be used in compiled form. 1424 1425(copyd '!%!%reclaim '!%reclaim) 1426 1427(fluid '(gcvar1!* gcvar2!* arithargloc)) 1428 1429(de !%reclaim () 1430 (setq gcvar1!* (wgetv arithargloc 0)) 1431 (setq gcvar2!* (wgetv arithargloc 1)) 1432 (!%!%reclaim) 1433 (wputv arithargloc 0 gcvar1!*) 1434 (wputv arithargloc 1 gcvar2!*) 1435 nil) 1436 1437 1438