1%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2% 3% File: PU:NBIG30a.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: Experimental 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% 30-June-1993 (Herbert Melenk) 41% changed bigit length to full word size 42% 20-June-1989 (Winfried Neun) 43% changed bread so that long numbers can be read with less bignum overhead 44% 24-Jan-89 (Herbert Melenk) 45% enlarged bigit length up to word size - 2 46% 47%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 48 49%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 50% A bignum will be a VECTOR of Bigits: (digits in base BigBase): 51% [BIGPOS b1 ... bn] or [BIGNEG b1 ... bn]. BigZero is thus [BIGPOS] 52% All numbers are positive, with BIGNEG as 0 element to indicate negatives. 53% 54% The base is set to wordsize-2; for multiplication and division a 55% double length integer arithmetic is required. 56% 57% NOTE that bbase* is no longer an inum! So we use open coded functions 58% to access the values of bbase* and logicalbits* 59%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 60 61(compiletime (load muls32 fast-vector vfvect inum if-system double32)) 62 63(compiletime (if_system i386 (load arith387))) 64 65%------------------------- control of debug output ------------------------ 66 67(compiletime 68 (progn 69 70 (setq bigtest nil) 71 72 (dm bigtest(u) 73 (if bigtest 74 (cons 'progn (cdr u)) 75 '(nullum))) 76 77 (put 'nullum 'opencode '((*move (reg 1)(reg 1)))) 78 79)) 80 81%------------------------- arithmetic macros ----------------------------- 82 83(compiletime (progn 84 85 (ds +w(a b)(iplus2 a b)) % add words 86 (ds +x(a b)(addAndSetCarry a b)) % add extended (set carry*) 87 (ds +c(a b)(addWithCarry a b)) % add with carry and set carry* 88 (ds +x+c(a b) (addAndAddCarry a b)) % add a and b and add carry to 89 % *second-value* 90 91 (ds -c(a b)(subtractwithborrow a b)) 92 93 (ds >u(a b)(ugreaterp a b)) % unsigned greaterp 94 95 (remd '1+) 96 97 (ds 1+(x) (iadd1 x)) 98)) 99%------------------------- variables ------------------------------------- 100 101(fluid '( 102 betahi!* % Largest BetaNum as Inum 103 betalow!* % Smallest BetaNum as Inum 104 bigsyslow!* % Smallest SYSINT in BIG format 105 floatsyshi!* % Largest SYSINT in Float format 106 floatsyslow!* % Smallest SYSINT in Float format 107 floatbbase!* % bbase as a float 108 floatbbase-2* % 2.0 ** (bbits - 2) 109 bigfloathi!* % Largest Float in BIG format 110 bigfloatlow!* % Smallest Float in BIG format 111 staticbig!* % Warray for conversion of SYS to BIG 112 bone!* % A one 113 bzero!* % A zero 114 bbits!* % Number of Bits in BBASE!* 115 digit2letter!* carry!* outputbase!* 116 bigitsPerMantissa* % bigits needed for float conversion 117 *Karatsubabound* % number of bigits for Karatsuba start 118 float-shift* % for conversion bigit -> float 119 ieee-hidden-bit* 120 *second-value*)) 121 122(compiletime 123 (when (not bigtest) 124 (flag '(setbits trimbignum1 big2sysaux btwopower bexpt blor 125 blxor bland blnot blshift blrshift bllshift bminus 126 bplus2 bplusa2 bdifference bdifference2 127 btimes2 bdigittimes2 bsmalltimes2 bkaratsuba bwords 128 bwordsshiftleft bshiftandaddinto bdifference2inplace 129 bquotient bremainder bdivide-trivialtest 130 bsimpledivide bharddivide bhardbug 131 bgreaterp blessp bgeq bleq bunsignedgreaterp bunsignedgeq 132 badd1 bsub1 bdivide1000000ip bread breadadd bsmalladd 133 bnum bnumaux bsmalldiff biggcdn0 biggcdn1 bigit2float 134 ) 'internalfunction) 135)) 136 137(setq *Karatsubabound* 39) 138 139% -------------------------------------------------------------------------- 140% the bignum routines handle bit patterns which exceed the inf field. 141% If these items sleep in the stack when a garbage collection occurs, 142% the heap management will fail. Therefore we 143% - do all heap activities before the bigit operations start, 144% - clean up the stack when leaving a routine with bigit operations. 145% The macro "Cleanstack" does this job if included in a return expression: 146% it overwrites all stack positions and is transparent with respect to 147% the result lying in register 1. The macro code should be machine independent. 148% It will work correctly only if Nalloc* is set properly! 149% That is the case for 68020 e.g. with frames > 2 150 151(compiletime (load if-system)) 152 153(compiletime 154(if_system SPARC 155(progn 156 (defCmacro *CleanStack) 157 (de *CleanStack() 158 (prog (u) 159 (for (from i 1 nalloc* 1) 160 (do (setq u (cons `(*MOVE (reg 2) (FRAME ,i)) u)))) 161 (return (cons '(*MOVE (quote NIL) (reg 2)) u)))) 162 (put 'CleanStack1 'OpenCode 163 '((*CleanStack))) 164 (put 'CleanStack2 'OpenCode '((*MOVE (reg 1)(reg 1)))) % this is a dummy 165 (ds Cleanstack(r)(Cleanstack2 (Cleanstack1 r))) 166 ) 167 (progn 168 (defCmacro *CleanStack) 169 (de *CleanStack() 170 (prog (u) 171 (for (from i 1 nalloc* 1) 172 (do (setq u (cons `(*MOVE (reg 2) (FRAME ,i)) u)))) 173 (return (cons '(*MOVE (quote NIL) (reg 2)) u)))) 174 (put 'CleanStack1 'OpenCode 175 '((*CleanStack))) 176 (put 'CleanStack2 'OpenCode '((*MOVE (reg 1)(reg 1)))) % this is a dummy 177 (ds Cleanstack(r)(Cleanstack2 (Cleanstack1 r))) 178 ))) 179 180% -------------------------------------------------------------------------- 181 182% Support functions: 183% 184% U, V, V1, V2 for arguments are Bignums. Other arguments are usually 185% fix/i-nums. 186 187(de setbits (x) 188 % This function sets the globals for big bignum package. 189 % "x" should be total # of bits per word. 190 (prog (y) 191 (setq bbits!* x) % Total number of bits per word used. 192 % "Beta", where n=A0 + A1*beta + A2*(beta^2). 193 (setq floatbbase!* (float-expt 2.000 bbits!*)) 194 (setq floatbbase-2* (float-expt 2.000 (idifference bbits!* 2))) 195 % calculate the number of bigits needed for float conversion 196 (setq bigitsPerMantissa* 1) 197 (setq y floatbbase!*) 198 (while (not (eqn y (floatplus2 y floatbbase!*))) 199 (setq bigitsPerMantissa* (1+ bigitsPerMantissa*)) 200 (setq y (floattimes2 y floatbbase!*))) 201 (setq float-shift* (expt 2.0 (quotient bitsperword 2))) 202 % the general beta range is next integer below half word size 203 (setq betahi!* (isub1 (expt 2 (isub1 (quotient bitsperword 2))))) 204 (setq betalow!* (iminus betahi!*)) 205 (setq bone!* (bnum 1)) 206 (setq bzero!* (bnum 0)) 207)) 208 209% Smallest representable Syslisp integer. 210(de nonbignumerror (v l) 211 (stderror (bldmsg " Expect a BIGNUM in %r, given %p%n" l v))) 212 213(de bsize (v) 214 % Upper Limit of [BIGxxx a1 ... An] 215 (if (bigp v) (veclen (vecinf v)) 0)) 216 217(compiletime (ds bbsize (v) (veclen (vecinf v)))) 218 219(compiletime (ds bbminusp (v1) (eq (igetv v1 0) 'bigneg))) 220 221(de gtpos (n) 222 % Allocate [BIGPOS a1 ... an] 223 (prog nil 224 (setq n (gtwrds n)) 225 (iputv n 0 'bigpos) 226 (return (mkbign (vecinf n))))) 227 228(de gtneg (n) 229 % Allocate [BIGNEG a1 ... an] 230 (prog nil 231 (setq n (gtwrds n)) 232 (iputv n 0 'bigneg) 233 (return (mkbign (vecinf n))))) 234 235(de trimbignum (v3) 236 % truncate trailing 0 237 (if (not (bigp v3)) 238 (nonbignumerror v3 'trimbignum) 239 (trimbignum1 v3 (bbsize v3)))) 240 241(compiletime (ds bigaswrds (b) (mkwrds (inf b)))) 242 243(de trimbignum1 (b l3) 244 (prog (v3) 245 (setq v3 (bigaswrds b)) 246 (while (and (igreaterp l3 0) (izerop (igetv v3 l3))) 247 (setq l3 (isub1 l3))) 248 (if (izerop (upbw (truncatewords v3 l3))) 249 (return bzero!*) 250 (return b)))) 251 252%--------------------------- conversion --------------------------- 253 254(de big2sys (u) 255 % Convert a BIG to SYS, if in range 256 (prog(s b) 257 (setq b (bbsize u)) 258 (when (eq b 0) (return 0)) 259 (setq s (igetv u 1)) 260 (when (or (wgreaterp b 1) 261 (wlessp s 0)) 262 (continuableerror 99 "BIGNUM too large to convert to SYS" u)) 263 (return (if (bbminusp u) (wminus s) s)) 264 )) 265 266(de big2sysaux (u) (big2sys u)) 267 268(de sysx2int (u) 269 % convert a sysnumber plus *second-value* which contains the upper bits 270 % into a sysnumber, fixnumber or bignumber 271 % Here u and *second-value* are signed wordsize integers. 272 (prog (secundus tagu tags sign big) 273 (setq secundus *second-value*) % the overflow 274 275 (cond ((and (eq 0 secundus) (wgeq u 0)) (return (sys2int u))) 276 ((and (eq -1 secundus) (wlessp u 0)) (return (sys2int u))) 277 ((wlessp secundus 0) 278 (setq secundus (wnot secundus)) 279 (setq u (1+ (wnot u))) 280 (when (eq u 0) (setq secundus (1+ secundus))) % a carry 281 (setq sign t))) 282 283 % now all positive, must be hidden for gc ! 284 285 (setq tagu (tag u) u (inf u)) 286 (setq tags (tag secundus) secundus(inf secundus)) 287 (setq *second-value* 0) 288 289 (setq big (if sign (gtneg 2) (gtpos 2))) 290 (iputv big 2 (mkitem tags secundus)) 291 (iputv big 1 (mkitem tagu u)) 292 % upper word might be 0 293 (return (trimbignum1 big 2)))) 294 295(flag '(sysx2int) 'lose) 296%---------------------- aux arithmetic ----------------------------- 297 298(de twopower (n)(expt 2 n)) 299 300(de btwopower (n) 301 % gives 2**n; n is fix/i-num; result BigNum 302 (if (not (or (fixp n) (bigp n))) 303 (nonintegererror n 'btwopower) 304 (prog (quot rem v) 305 (when (bigp n) 306 (setq n (big2sys n))) 307 (setq quot (quotient n bbits!*)) 308 (setq rem (remainder n bbits!*)) 309 (setq v (gtpos (1+ quot))) 310 (vfor (from i 1 quot 1) (do (iputv v i 0))) 311 (iputv v (1+ quot) (twopower rem)) 312 (return (trimbignum1 v (1+ quot)))))) 313 314 315(compiletime 316 (progn 317 (ds bzerop (v1) (and (izerop (bbsize v1)) (not (bbminusp v1)))) 318 (ds bonep (v1) (and (not (bbminusp v1)) 319 (ionep (bbsize v1)) (ionep (igetv v1 1)))) 320 (ds babs (v1) (if (bbminusp v1) (bminus v1) v1)) 321 (ds bmax (v1 v2) (if (bgreaterp v2 v1) v2 v1)) 322 (ds bmin (v1 v2) (if (blessp v2 v1) v2 v1)) 323)) 324 325(de bexpt (v1 n) 326 % V1 is Bignum, N is fix/i-num 327 (cond ((not (fixp n)) (nonintegererror n 'bexpt)) 328 ((izerop n) bone!*) 329 ((ionep n) v1) 330 ((iminusp n) (bquotient bone!* (bexpt v1 (iminus n)))) 331 (t (prog (v2) 332 (setq v2 (bexpt v1 (iquotient n 2))) 333 (if (izerop (iremainder n 2)) 334 (return (btimes2 v2 v2)) 335 (return (btimes2 (btimes2 v2 v1) v2))))))) 336 337% --------------------------------------- 338% Logical Operations 339% 340% All take Bignum arguments 341(de blor (v1 v2) 342 % The main body of the OR code is only obeyed when both arguments 343 % are positive, and so the result will be positive; 344 (if (or (bbminusp v1) (bbminusp v2)) 345 (blnot (bland (blnot v1) (blnot v2))) 346 (prog (l1 l2 l3 v3) 347 (setq l1 (bbsize v1)) 348 (setq l2 (bbsize v2)) 349 (when (igreaterp l2 l1) (setq l3 l2) (setq l2 l1) (setq l1 l3) 350 (setq v3 v2) (setq v2 v1) (setq v1 v3)) 351 (setq v3 (gtpos l1)) 352 (vfor (from i 1 l2 1) 353 (do (iputv v3 i (ilor (igetv v1 i) (igetv v2 i))))) 354 (vfor (from i (1+ l2) l1 1) (do (iputv v3 i (igetv v1 i)))) 355 (return (cleanstack v3))))) 356 357(de blxor (v1 v2) 358 % negative arguments are coped with using the identity 359 % LXor(a,b) = LNot LXor(Lnot a,b) = LNor LXor(a,Lnot b); 360 (prog (l1 l2 l3 v3 s) 361 (when (bbminusp v1) (setq v1 (blnot v1)) (setq s t)) 362 (when (bbminusp v2) (setq v2 (blnot v2)) (setq s (not s))) 363 (setq l1 (bbsize v1)) 364 (setq l2 (bbsize v2)) 365 (when (igreaterp l2 l1) (setq l3 l2) (setq l2 l1) (setq l1 l3) 366 (setq v3 v2) (setq v2 v1) (setq v1 v3)) 367 (setq v3 (gtpos l1)) 368 (vfor (from i 1 l2 1) 369 (do (iputv v3 i (ilxor (igetv v1 i) (igetv v2 i))))) 370 (vfor (from i (1+ l2) l1 1) (do (iputv v3 i (igetv v1 i)))) 371 (setq v1 (trimbignum1 v3 l1)) 372 (when s (setq v1 (blnot v1))) 373 (return (cleanstack v1)))) 374 375(de bland (v1 v2) 376 % If both args are -ve the result will be too. Otherwise result will 377 % be positive; 378 (if (and (bbminusp v1) (bbminusp v2)) 379 (blnot (blor (blnot v1) (blnot v2))) 380 (prog (l1 l2 l3 v3) 381 (setq l1 (bbsize v1)) (setq l2 (bbsize v2)) (setq l3 (min l1 l2)) 382 (cond ((bbminusp v1) 383 % When one is negative, we have expand out to the 384 % size of the other one. Therefore, we use l2 as the 385 % size, not l3. When we exceed the size of the 386 % negative number, then we just use (logicalbits**). 387 (setq v3 (gtpos l2)) (setq l3 l2) 388 (vfor (from i 1 l2 1) 389 (do 390 (iputv v3 i 391 (iland (cond ((igreaterp i l1) (logicalbits**)) 392 (t (isub1 (igetv v1 i)))) 393 (igetv v2 i)))))) 394 ((bbminusp v2) 395 (setq v3 (gtpos l1)) (setq l3 l1) 396 (vfor (from i 1 l1 1) 397 (do 398 (iputv v3 i 399 (iland (igetv v1 i) 400 (cond ((igreaterp i l2)(logicalbits**)) 401 (t (isub1 (igetv v2 i))))))))) 402 403 (t (setq v3 (gtpos l3)) 404 (vfor (from i 1 l3 1) 405 (do 406 (iputv v3 i (iland (igetv v1 i) (igetv v2 i))))))) 407 (return(cleanstack (trimbignum1 v3 l3)))))) 408 409(de blnot (v1) (bminus (bsmalladd v1 1))) 410 411(de blshift(v1 n) 412 (if (iminusp n) (blrshift v1 (iminus n)) (bllshift v1 n))) 413 414% for both shifts we have to handle the case 415% shift amount = wordsize 416% separately because processors tend to handle a shift for 417% (wordsize) bits as nop. 418 419(de blrshift(v1 n) 420 % shift right n places 421 (prog (nw nb nr carry l1 v2 l2 x) 422 % split n into words and bits 423 (setq nw (iquotient n bbits*) nb (iminus (iremainder n bbits*)) 424 nr (+w bbits* nb)) 425 (setq l1 (bbsize v1)) 426 (setq l2 (idifference l1 nw)) 427 (when (ilessp l2 1) (return bzero*)) 428 (setq v2 (if (bbminusp v1)(gtneg l2)(gtpos l2))) 429 % for shifts we have to handle the case nb=0 430 % separately because processors tend to handle a shift for 431 % nr=(-wordsize) bits as nop. 432 (when (izerop nb) (go words)) 433 (setq carry 0) 434 (vfor (from i l2 1 -1) 435 (let j (+w i nw)) 436 (do 437 (progn 438 (setq x (igetv v1 j)) 439 (iputv v2 i(wor carry (wshift x nb))) 440 (setq carry (wshift x nr))))) 441 (go ret) 442 words 443 (vfor (from i l2 1 -1) 444 (let j (+w i nw)) 445 (do 446 (iputv v2 i (igetv v1 j)))) 447 ret 448 (return (cleanstack (trimbignum1 v2 l2))))) 449 450(de bllshift(v1 n) 451 % shift left n places 452 (prog (nw nb nr carry l1 v2 l2 x) 453 % split n into words and bits 454 (setq nw (iquotient n bbits*) nb (iremainder n bbits*) 455 nr (idifference nb bbits*)) % nr is <=0 !! 456 (setq l1 (bbsize v1)) 457 (setq l2 (+w l1 (1+ nw))) 458 (setq v2 (if (bbminusp v1)(gtneg l2)(gtpos l2))) 459 (setq carry 0) 460 % clear the trailing words 461 (vfor (from i 1 nw 1) (do (iputv v2 i 0))) 462 (when (izerop nb) (go words)) 463 (vfor (from i (1+ nw) (isub1 l2) 1) 464 (let j (idifference i nw)) 465 (do 466 (progn 467 (setq x (igetv v1 j)) 468 (iputv v2 i(wor carry (wshift x nb))) 469 (setq carry (wshift x nr))))) 470 (go ret) 471 words 472 (vfor (from i (1+ nw) (isub1 l2) 1) 473 (let j (idifference i nw)) 474 (do 475 (iputv v2 i(igetv v1 j)) )) 476 ret 477 (iputv v2 l2 carry) 478 (return (cleanstack (trimbignum1 v2 l2)))) ) 479 480% ----------------------------------------- 481% Arithmetic Functions: 482% U, V, V1, V2 are Bignum arguments. 483 484% --------------------------------- MINUS ------------------------------ 485 486(de bminus (v1) % Negates V1. 487 (if (bzerop v1) v1 488 (prog (l1 v2) 489 (setq l1 (bbsize v1)) 490 (if (bbminusp v1) 491 (setq v2 (gtpos l1)) 492 (setq v2 (gtneg l1))) 493 (vfor (from i 1 l1 1) (do (iputv v2 i (igetv v1 i)))) 494 (return (cleanstack v2))))) 495 496% --------------------------------- MINUSP ----------------------------- 497 498% Returns V1 if V1 is strictly less than 0, NIL otherwise. 499(de bminusp (v1) (if (eq (igetv v1 0) 'bigneg) v1 nil)) 500 501% ---------------------------------- PLUS ----------------------------- 502 503(de bplus2 (v1 v2) 504 (prog (sn1 sn2) 505 (setq sn1 (bbminusp v1)) (setq sn2 (bbminusp v2)) 506 (when (and sn1 (not sn2)) 507 (return (bdifference2 v2 (bminus v1) nil))) 508 (when (and sn2 (not sn1)) 509 (return (bdifference2 v1 (bminus v2) nil))) 510 (return (bplusa2 v1 v2 sn1)))) 511 512(compiletime 513 (put 'plusloop1 'opencode '( 514 (*alloc 5) 515 (*move (reg 1)(frame 1)) 516 (*move (reg 2)(frame 2)) 517 loop1 518 (*jumpwgreaterp (label done1)(reg 1)(reg 2)) 519 (*move (memory (reg 4) 0)(reg 1)) 520 (*move (memory (reg 5) 0)(reg 2)) 521 %(*call addwithcarry) 522 (*addwithcarry (reg 1)(reg 2)) 523 (*move (reg 1)(memory (reg 3) 0)) 524 (*wplus2 (reg 3) addressingunitsperitem) 525 (*wplus2 (reg 4) addressingunitsperitem) 526 (*wplus2 (reg 5) addressingunitsperitem) 527 (*move (frame 1)(reg 1)) 528 (*move (frame 2)(reg 2)) 529 (*wplus2 (reg 1) 1) 530 (*move (reg 1)(frame 1)) 531 (*jump (label loop1)) 532 done1 533 (*dealloc 5) 534))) 535 536(compiletime 537 (put 'plusloop2 'opencode '( 538 (*alloc 5) 539 (*move (reg 1)(frame 1)) 540 (*move (reg 2)(frame 2)) 541 loop2 542 (*jumpwgreaterp (label done2)(reg 1)(reg 2)) 543 (*move (memory (reg 4) 0)(reg 1)) 544 (*move 0 (reg 2)) 545 %(*call addwithcarry) 546 (*addwithcarry (reg 1)(reg 2)) 547 (*move (reg 1)(memory (reg 3) 0)) 548 (*wplus2 (reg 3) addressingunitsperitem) 549 (*wplus2 (reg 4) addressingunitsperitem) 550 (*move (frame 1)(reg 1)) 551 (*move (frame 2)(reg 2)) 552 (*wplus2 (reg 1) 1) 553 (*move (reg 1)(frame 1)) 554 (*jump (label loop2)) 555 done2 556 (*dealloc 5) 557))) 558 559 560(de bplusa2 (v1 v2 sn1) 561 % Plus with signs pre-checked and identical. 562 (prog (l1 l2 l3 v3 temp) 563 (setq l1 (bbsize v1)) (setq l2 (bbsize v2)) 564 (when (igreaterp l2 l1) (setq l3 l2) (setq l2 l1) (setq l1 l3) 565 (setq v3 v2) (setq v2 v1) (setq v1 v3)) 566 (setq l3 (1+ l1)) 567 (if sn1 (setq v3 (gtneg l3)) (setq v3 (gtpos l3))) 568 (setq carry!* 0) 569 %(vfor (from i 1 l2 1) 570 % (do (iputv v3 i (+c (igetv v1 i) (igetv v2 i)))) ) 571 (plusloop1 1 l2 572 (loc (igetv v3 1)) 573 (loc (igetv v1 1)) 574 (loc (igetv v2 1)) 575 ) 576 (setq temp (1+ l2)) 577 %(vfor (from i temp l1 1) (do (iputv v3 i (+c (igetv v1 i) 0)))) 578 (plusloop2 temp l1 579 (loc (igetv v3 temp)) 580 (loc (igetv v1 temp))) 581 (iputv v3 l3 carry!*) 582 (setq carry!* 0) % clear bigit 583 (return (cleanstack (trimbignum1 v3 l3))))) 584 585% ------------------------------ DIFFERENCE ---------------------------- 586 587(de bdifference (v1 v2) 588 (cond ((bzerop v2) v1) 589 ((bzerop v1) (bminus v2)) 590 (t (prog (sn1 sn2) 591 (setq sn1 (bbminusp v1)) 592 (setq sn2 (bbminusp v2)) 593 (when (or (and sn1 (not sn2)) (and sn2 (not sn1))) 594 (return (bplusa2 v1 (bminus v2) sn1))) 595 (return (bdifference2 v1 v2 sn1)))))) 596 597(compiletime 598 (put 'differenceloop0 'opencode '( 599 loop0 600 (*jumpeq (label done0)(reg 1) 0) 601 (*move (memory (reg 4) 0)(reg 3)) 602 (*move (memory (reg 5) 0)(reg 2)) 603 (*jumpnoteq (label done0) (reg 2)(reg 3)) 604 (*wdifference (reg 4) addressingunitsperitem) 605 (*wdifference (reg 5) addressingunitsperitem) 606 (*wdifference (reg 1) 1) 607 (*jump (label loop0)) 608 done0 609))) 610 611(compiletime 612 (put 'differenceloop1 'opencode '( 613 (*alloc 5) 614 (*move (reg 1)(frame 1)) 615 (*move (reg 2)(frame 2)) 616 loop1 617 (*jumpwgreaterp (label done1)(reg 1)(reg 2)) 618 (*move (memory (reg 4) 0)(reg 1)) 619 (*move (memory (reg 5) 0)(reg 2)) 620 %(*call subtractwithborrow) 621 (*subtractwithborrow (reg 1)(reg 2)) 622 (*move (reg 1)(memory (reg 3) 0)) 623 (*wplus2 (reg 3) addressingunitsperitem) 624 (*wplus2 (reg 4) addressingunitsperitem) 625 (*wplus2 (reg 5) addressingunitsperitem) 626 (*move (frame 1)(reg 1)) 627 (*move (frame 2)(reg 2)) 628 (*wplus2 (reg 1) 1) 629 (*move (reg 1)(frame 1)) 630 (*jump (label loop1)) 631 done1 632 (*dealloc 5) 633))) 634 635(compiletime 636 (put 'differenceloop2 'opencode '( 637 (*alloc 5) 638 (*move (reg 1)(frame 1)) 639 (*move (reg 2)(frame 2)) 640 loop2 641 (*jumpwgreaterp (label done2)(reg 1)(reg 2)) 642 (*move (memory (reg 4) 0)(reg 1)) 643 (*move 0 (reg 2)) 644 %(*call subtractwithborrow) 645 (*subtractwithborrow (reg 1)(reg 2)) 646 (*move (reg 1)(memory (reg 3) 0)) 647 (*wplus2 (reg 3) addressingunitsperitem) 648 (*wplus2 (reg 4) addressingunitsperitem) 649 (*move (frame 1)(reg 1)) 650 (*move (frame 2)(reg 2)) 651 (*wplus2 (reg 1) 1) 652 (*move (reg 1)(frame 1)) 653 (*jump (label loop2)) 654 done2 655 (*dealloc 5) 656))) 657 658 659(de bdifference2 (v1 v2 sn1) 660 % Signs pre-checked and identical. 661 (prog (i l1 l2 l3 v3) 662 (setq l1 (bbsize v1)) 663 (setq l2 (bbsize v2)) 664 (cond ((igreaterp l2 l1) (setq l3 l1) (setq l1 l2) (setq l2 l3) 665 (setq v3 v1)(setq v1 v2)(setq v2 v3)(setq sn1 (not sn1))) 666 ((eq l1 l2) 667 % (setq i l1) 668 % (while (and (eq (igetv v2 i) (igetv v1 i)) (igeq i 1)) 669 % (setq i (isub1 i))) 670 (setq i (differenceloop0 l1 0 0 671 (loc (igetv v2 l1)) 672 (loc (igetv v1 l2)))) 673 (when (eq i 0)(return (cleanstack bzero*))) 674 (when (>u (igetv v2 i) (igetv v1 i)) 675 (setq l3 l1) (setq l1 l2) (setq l2 l3) (setq v3 v1) 676 (setq v1 v2) (setq v2 v3) (setq sn1 (not sn1))))) 677 (if sn1 (setq v3 (gtneg l1)) (setq v3 (gtpos l1))) 678 (setq carry!* 0) 679 %(vfor (from i 1 l2 1) 680 % (do (iputv v3 i 681 % (-c (igetv v1 i) (igetv v2 i))))) 682 (differenceloop1 1 l2 683 (loc (igetv v3 1)) (loc (igetv v1 1)) (loc (igetv v2 1))) 684 %(vfor (from i (1+ l2) l1 1) 685 % (do (iputv v3 i (-c (igetv v1 i) 0)))) 686 (differenceloop2 (1+ l2) l1 687 (loc (igetv v3 1)) (loc (igetv v1 1))) 688 (setq carry!* 0) % clear bigit 689 (return (cleanstack (trimbignum1 v3 l1))))) 690 691% ------------------------------ TIMES --------------------------------- 692 693(compiletime 694 (put 'cleanloop 'opencode '( 695 clean-loop 696 (*jumpwgreaterp (label clean-done)(reg 1)(reg 2)) 697 (*move 0 (memory (reg 3) 0)) 698 (*wplus2 (reg 3) addressingunitsperitem) 699 (*wplus2 (reg 1) 1) 700 (*jump (label clean-loop)) 701 clean-done 702))) 703 704 705(de btimes2 (v1 v2) 706 (prog (l1 l2 l3 sn1 sn2 v3) 707 (setq l1 (bbsize v1)) (setq l2 (bbsize v2)) 708 (when (and(igreaterp l1 *Karatsubabound*) 709 (igreaterp l2 *Karatsubabound*)) 710 (return (bkaratsuba v1 l1 v2 l2))) 711 (when (igreaterp l2 l1) (setq v3 v1) (setq v1 v2) (setq v2 v3) 712 (setq l3 l1) (setq l1 l2) (setq l2 l3)) 713 % iterations of BDigitTimes2. 714 (setq l3 (+w l1 l2)) 715 (setq sn1 (bbminusp v1)) 716 (setq sn2 (bbminusp v2)) 717 (if (or (and sn1 sn2) (not (or sn1 sn2))) 718 (setq v3 (gtpos l3)) (setq v3 (gtneg l3))) 719 %(vfor (from i 1 l3 1) (do (iputv v3 i 0))) 720 (cleanloop 1 l3 (loc (igetv v3 1))) 721 (vfor (from i 1 l2 1) (do (bdigittimes2 v1 (igetv v2 i) l1 i v3))) 722 (return (cleanstack (trimbignum1 v3 l3))))) 723 724(de bdigittimes2 (v1 v2 l1 i v3) 725 % V1 is a bignum, V2 a fixnum, L1=BSize L1, I=position of V2 in a bignum, 726 % and V3 is bignum receiving result. I affects where in V3 the result of 727 % a calculation goes; the relationship is that positions I:I+(L1-1) 728 % of V3 receive the products of V2 and positions 1:L1 of V1. 729 % V3 is changed as a side effect here. 730 % 731 % upper bigit and carry are accumulated in *second-value* 732 733 (if (izerop v2) v3 734 (prog (j carry w temp1) 735 (setq carry 0) 736 (vfor (from h 1 l1 1) 737 (let temp2 (+w h (isub1 i))) 738 (do 739 (progn 740 (setq temp1 (wtimesdouble (igetv v1 h) v2)) 741 (setq w (+x+c temp1 carry)) 742 (iputv v3 temp2 (+x+c w (igetv v3 temp2))) 743 (setq carry *second-value*) 744 ))) 745 (iputv v3 (+w l1 i) carry) 746 % digit should be < BBase!* here 747 (setq *second-value* 0) % clear carry 748 (return (cleanstack v3))))) 749 750(de bsmalltimes2 (v1 c cc) 751 % V1 is a BigNum, C a fixnum. 752 % Assume C positive, ignore sign(V1) 753 % also assume V1 neq 0. 754 (cond ((and (izerop c)(izerop cc)) (return (gtpos 0))) 755 (t % Only used from BHardDivide, BReadAdd. 756 (prog (j l1 v3 carry) 757 (setq l1 (bbsize v1)) 758 (setq v3 (gtpos (1+ l1))) 759 (if (not(izerop cc)) % rejoin splitted word 760 (setq c (+w c (wshift cc (difference bitsperword 8))))) 761 (setq carry 0) (setq carry* 0) 762 (vfor (from h 1 l1 1) 763 (do 764 (progn 765 (setq j (+x (wtimesdouble (igetv v1 h) c) carry)) 766 (iputv v3 h j) 767 (setq carry (+w *second-value* carry*)) 768 ))) 769 (iputv v3 (1+ l1) carry) 770 (setq *second-value* 0) % clear carry 771 (return (cleanstack (trimbignum1 v3 (1+ l1)))))))) 772 773% Karatsuba method for multiplication of very large bignums 774% (e.g. from 350 decimal digits upwards) 775 776(de bkaratsuba(u lu v lv) 777 (prog (su sv n U0 V0 U1 V1 p1 p2 p3 uv) 778 % setq n to max(lu,lv)/2, rounded up 779 (setq n (if (igreaterp lu lv) lu lv)) 780 (setq n (if (izerop (wand n 1)) (wshift n -1)(1+(wshift n -1)))) 781 (setq U0 (bwords u lu 1 n) U1 (bwords u lu (1+ n) lu)) 782 (setq V0 (bwords v lv 1 n) V1 (bwords v lv (1+ n) lv)) 783 % the basic formula is 784 % uv := (2**2n + 2**n) U1 V1 + 785 % 2**n (U1 - U0)(V0 - V1) + 786 % (2**n + 1) U0 V0 787 (setq p1 (btimes2 U1 V1)) 788 (setq p3 (btimes2 U0 V0)) 789 % noo we can destory U0,..,V1 790 (setq p2 (btimes2 (bdifference2inPlace U1 U0 nil) 791 (bdifference2inPlace V0 V1 nil))) 792 (setq uv (bwordsshiftleft p1 (+w n n))) % the biggest part 793 (setq uv (bshiftAndAddInto p1 n uv)) 794 (setq uv (bshiftAndAddInto p2 n uv)) 795 (setq uv (bshiftAndAddInto p3 n uv)) 796 (setq uv (bshiftAndAddInto p3 0 uv)) 797 (when (not (eq (igetv u 0)(igetv v 0))) 798 (iputv uv 0 'bigneg)) 799 (return uv))) 800 801% auxiliary routines for Karatsuba 802 803(de bwords (u lu ilo ihi) 804 % select words ilo to ihi from bignum u with length lu 805 % result is a psitive bignum in any case 806 (prog (v lv x) 807 (when (igreaterp ihi lu)(setq ihi lu)) 808 (setq lv (1+(idifference ihi ilo))) 809 (when (ilessp lv 1) (return bzero*)) 810 (setq v (gtpos lv)) 811 (setq ilo(isub1 ilo)) 812 (vfor (from i 1 lv 1) 813 (let tmp (+w i ilo)) 814 (do (progn 815 (setq x (igetv u tmp)) 816 (iputv v i x)))) 817 (return (cleanstack (trimbignum1 v lv))))) 818 819 820(de bwordsshiftleft(u n) 821 (prog (lu v lv x) 822 (setq lu (bsize u)) 823 (setq lv (+w lu n)) 824 (setq v (if (bminusp u)(gtneg lv)(gtpos lv))) 825 (vfor (from i 1 n 1) (do (iputv v i 0))) 826 (vfor (from i (1+ n) lv 1) 827 (let tmp (idifference i n)) 828 (do (progn 829 (setq x (igetv u tmp)) 830 (iputv v i x)))) 831 (return (cleanstack (trimbignum1 v lv))))) 832 833 834(de bshiftAndAddInto (v2 n v1) 835 % v1 := v2*2**(n*bbits) + v1 836 % v1 is overwritten if it is large enough 837 (prog (l1 l2 l3 v3 temp) 838 (setq l1 (bbsize v1)) (setq l2 (bbsize v2)) 839 (when (and (bbminusp v2)(not(bbminusp v1))) 840 (setq v2 (bwordsshiftleft v2 n)) % v2 new 841 (iputv v2 0 'bigpos) % invert sign 842 (return (bdifference2inPlace v1 v2 nil))) 843 (when (and (bbminusp v1) (not (bminusp v2))) 844 (setq v2 (bwordsshiftleft v2 n)) % v2 new 845 (iputv v1 0 'bigpos) % destroys v1 here 846 (return (bdifference2inPlace v2 v1 nil))) 847 (when (igreaterp (+w l2 n) l1) 848 (return (bshiftAndAddInto v1 0 (bwordsshiftleft v2 n) ))) 849 % now signs are equal and v1 is long enough 850 (setq carry!* 0) 851 (setq l3(+w l2 n)) 852 (vfor (from i (1+ n) l3 1) 853 (let tmp (idifference i n)) 854 (do (progn (setq temp (+c (igetv v1 i) (igetv v2 tmp))) 855 (iputv v1 i temp)))) 856 (setq temp (1+ l3)) 857 (vfor (from i temp l1 1) (do (iputv v1 i (+c (igetv v1 i) 0)))) 858 % if there is a carry left we have to create a larger box 859 (when (not (izerop carry!*)) 860 (setq v3 (gtpos (1+ l1))) 861 (vfor (from i 1 l1 1) 862 (do (iputv v3 i (igetv v1 i)))) 863 (setq l1 (1+ l1)) 864 (iputv v3 l1 carry!*) 865 (setq v1 v3)) 866 (return (cleanstack v1)))) 867 868 869 870 871(de bdifference2inPlace (v1 v2 sn1) 872 % same as bdifference2, but v1 or v2 are overwritten by the result 873 % Signs pre-checked and identical. 874 (prog (i l1 l2 l3 v3) 875 (setq l1 (bbsize v1)) 876 (setq l2 (bbsize v2)) 877 (cond ((igreaterp l2 l1) (setq l3 l1) (setq l1 l2) (setq l2 l3) 878 (setq v3 v1) (setq v1 v2) (setq v2 v3) (setq sn1 (not sn1))) 879 ((eq l1 l2) (setq i l1) 880 (while (and (eq (igetv v2 i) (igetv v1 i)) 881 (igreaterp i 1)) 882 (setq i (isub1 i))) 883 (when (>u (igetv v2 i) (igetv v1 i)) 884 (setq l3 l1) (setq l1 l2) (setq l2 l3) (setq v3 v1) 885 (setq v1 v2) (setq v2 v3) (setq sn1 (not sn1))))) 886 (setq v3 v1) % overwriting! 887 (iputv v3 0 (if sn1 'bigneg 'bigpos)) 888 (setq carry!* 0) 889 (vfor (from i 1 l2 1) 890 (do (iputv v3 i (-c (igetv v1 i) (igetv v2 i))))) 891 (vfor (from i (1+ l2) l1 1) 892 (do (iputv v3 i (-c (igetv v1 i) 0)))) 893 (return (cleanstack (trimbignum1 v3 l1))))) 894 895 896 897 898% ------------------------ DIVISION ------------------------------------ 899 900(de bquotient (v1 v2) (car (bdivide v1 v2))) 901(de bremainder (v1 v2) (cdr (bdivide v1 v2))) 902 903% BDivide returns a dotted pair, (Q . R). Q is the quotient and R is 904% the remainder. Both are bignums. R is of the same sign as V1. 905 906(de bdivide (v1 v2) 907 (prog (l1 l2 q r v3) 908 (setq l2 (bbsize v2)) 909 (when (izerop l2) 910 (error 99 "Attempt to divide by 0 in BDIVIDE")) 911 (setq l1 (bbsize v1)) 912 (cond ((or (ilessp l1 l2) 913 (and (eq l1 l2) (bdivide-trivialtest v1 l1 v2 l2))) 914 (return (cons bzero* v1)))) 915 (when (ionep l2) 916 (return (bsimpledivide v1 l1 v2 (bbminusp v2)))) 917 (return (bharddivide v1 l1 v2 l2)))) 918 919(de bdivide-trivialtest(v1 l1 v2 l2) 920 (>u (igetv v2 l1) (igetv v1 l2))) 921 922% C is a fixnum (inum?); V1 is a bignum and L1 is its length. 923% SnC is T if C (which is positive) should be considered negative. 924% Returns quotient . remainder; each is a bignum. 925% 926(de bsimpledivide (v1 l1 c snc) 927 (prog (i p r rr sn1 v2 p) 928 (setq sn1 (bbminusp v1)) 929 (if (or (and sn1 snc) (not (or sn1 snc))) 930 (setq v2 (gtpos l1)) (setq v2 (gtneg l1))) 931 (if sn1 (setq rr (gtneg 1)) (setq rr (gtpos 1))) 932 (setq p (cons nil nil)) % do all heap requests before bigits come up 933 (when (bigp c)(setq c (igetv c 1))) % now no more gc can occur 934 (setq r 0) 935 (setq i l1) 936 (while (not (izerop i)) 937 (iputv v2 i (wquotientdouble r (igetv v1 i) c)) 938 (setq r *second-value*) 939 (setq i (isub1 i))) 940 (iputv rr 1 r) 941 (setq *second-value* 0) % clear carry 942 (rplaca p (trimbignum1 v2 l1)) (rplacd p (trimbignum1 rr 1)) % the pair 943 (return (cleanstack p)))) 944 945% boxes for double precision integer arithmetic 946(fluid '(box1 box2 box3 box4 box5 box6 box7)) 947(setq box1 (mkdouble)) (setq box2 (mkdouble)) (setq box3 (mkdouble)) 948(setq box4 (mkdouble)) (setq box5 (mkdouble)) (setq box6 (mkdouble)) 949(setq box7 (mkdouble)) 950 951(de bharddivide (u lu v lv) 952 % This is an algorithm taken from Knuth. 953 (prog (u1 v1 a d lcv lcv1 f f2 j k lq carry temp ll m n n1 p q qbar snu 954 snv u2 x result s w) 955 (setq n lv) 956 (setq n1 (1+ n)) 957 (setq m (idifference lu lv)) 958 (setq lq (1+ m)) 959 % Deal with signs of inputs; 960 (setq snu (bbminusp u)) 961 (setq snv (bbminusp v)) 962 % Note that these are not extra-boolean, i.e. 963 % for positive numbers MBinusP returns nil, for 964 % negative it returns its argument. Thus the 965 % test (SnU=SnV) does not reliably compare the signs of 966 % U and V; 967 (cond (snu (if snv (setq q (gtpos lq)) (setq q (gtneg lq)))) 968 (snv (setq q (gtneg lq))) 969 (t (setq q (gtpos lq)))) 970 (setq u1 (gtpos (1+ lu))) 971 (setq result (cons nil nil)) 972 % U is ALWAYS stored as if one digit longer; 973 % Compute a scale factor to normalize the long division; 974 (setq d (1+ (igetv v lv))) 975 (bigtest (print (list "leading digit" (sys2int d)))) 976 (setq d (if (weq d 0) 1 (wquotientdouble 1 0 d))) 977 (bigtest (print (list "scale factor " (sys2int d)))) 978 % Now, at the same time, I remove the sign information from U and V 979 % and scale them so that the leading coefficeint in V is fairly large; 980 % protect possible gc in bsmalltimes2 from bigits in the stack 981 (setq x (wshift d (minus(difference bitsperword 8))) 982 d (wshift (wshift d 8)-8) ) 983 (setq v1 (bsmalltimes2 v d x)) 984 (setq d (+w d (wshift x (difference bitsperword 8) ))) 985 986 (setq carry 0) 987 (vfor (from i 1 lu 1) 988 (do (progn (setq temp (+x (wtimesdouble (igetv u i) d) 989 carry)) 990 (iputv u1 i temp) 991 (setq carry (+w carry* *second-value*)) 992 ))) 993 (setq lu (1+ lu)) 994 (iputv u1 lu carry) 995 (bigtest (print (list "scaled u" (big2l u1)))) 996 (bigtest (print (list "scaled v" (big2l v1)))) 997 998 999 % So far all variables contain safe values, 1000 % i.e. numbers < BBase!*; 1001 (iputv v1 0 'bigpos) 1002 (when (ilessp lv 2) (nonbignumerror v 'bharddivide)) 1003 % To be safe; 1004 (setq lcv (igetv v1 lv)) 1005 (setq lcv1 (igetv v1 (isub1 lv))) 1006 % Top two digits of the scaled V accessed once 1007 % here outside the main loop; 1008 % Now perform the main long division loop; 1009 (ifor (from i 0 m 1) 1010 (do (prog () 1011 (setq j (idifference lu i)) 1012 (bigtest (print (list "===== loop i=" i))) 1013 1014 % J>K; working on U1[K:J] in this loop. 1015 (setq k (idifference j n1)) 1016 (setq a (igetv u1 j)) 1017 (setq p (filldouble a (igetv u1 (isub1 j)) box1)) 1018 (bigtest (print (list "p" (d2l p)))) 1019 (bigtest (print (list "lcv" (sys2int lcv)))) 1020 (bigtest (print (list "lcv1" (sys2int lcv1)))) 1021 (if (eq a lcv) 1022 (setq qbar -1) 1023 (setq qbar (quotientDouble2word p lcv))) 1024 (bigtest (print (list "qbar" (sys2int qbar)))) 1025 1026 adjust 1027 % approximate by next digit; 1028 % p = u(j)*b + u(j-1) 1029 % lcv = v(1) 1030 % lcv1 = v(2) 1031 % f = q*v(2) 1032 % f2 = [u(j)*b + u(j-1) - q*v(1)]*b + u(j-2) 1033 % |----------------------| = w 1034 (setq f (timesword2double qbar lcv1 box2)) 1035 (bigtest (print (list "f" (d2l f)))) 1036 (setq w (differenceDouble p 1037 (timesword2double qbar lcv box3) box4 )) 1038 (bigtest (print (list "w" (d2l w)))) 1039 (when (not (eq 0 (doublehigh w))) 1040 (go subtract)) % f2 > b^2 > f 1041 (setq f2 (filldouble 1042 (doublelow w) 1043 (igetv u1 (idifference j 2)) box3)) 1044 (bigtest (print (list "f2" (d2l f2)))) 1045 1046 % Correct most overshoots in Qbar; 1047 (when (doublegreaterp f f2) 1048 (setq qbar (isub1 qbar)) 1049 (bigtest(print(list "qbar adusted" (sys2int qbar)))) 1050 (go adjust)) 1051 1052 subtract 1053 % Ready to subtract QBar*V1 from U1; 1054 % in this loop "carry" stores locally the 1055 % borrow of the subtract part. 1056 (setq carry 0 s 0) 1057 (vfor (from l 1 n 1) 1058 (let aux (+w k l)) 1059 (do (progn 1060 (setq w (wtimesdouble (igetv v1 l) qbar)) 1061 (setq w (+x w s)) 1062 (setq s (+w *second-value* carry*)) 1063 (setq carry* carry) 1064 (iputv u1 aux (-c (igetv u1 aux) w)) 1065 (setq carry carry*) 1066 ))) 1067 1068 % propagate borrows up as far as they go; 1069 (setq ll (+w k n)) 1070 (while (and (or (not (izerop carry*)) 1071 (not (izerop s))) 1072 (ilessp ll j)) 1073 (setq ll (1+ ll)) 1074 (iputv u1 ll (-c (igetv u1 ll) s)) 1075 (setq s 0)) 1076 1077 1078 (unless (izerop carry*) 1079 % QBar was still wrong - correction step needed. 1080 % This should not happen very often; 1081 (bigtest (prin2t "2nd correction needed")) 1082 (setq qbar (isub1 qbar)) 1083 % Add V1 back into U1; 1084 (setq carry* 0) 1085 (vfor (from l 1 n 1) 1086 (let aux (+w k l)) 1087 (do (iputv u1 aux 1088 (+c (igetv u1 aux) (igetv v1 l)) )) )) 1089 (setq ll (+w k n)) 1090 (while (ilessp ll j) 1091 (setq ll (1+ ll)) 1092 (iputv u1 ll (+c (igetv u1 ll) 0)) ) 1093 (iputv q (idifference lq i) qbar)))) 1094 % End of main loop; 1095 (setq u1 (trimbignum1 u1 (idifference lu m))) 1096 (setq f 0) 1097 (setq f2 0) 1098 (setq q (trimbignum1 q lq)) 1099 % Clean up potentially wild values; 1100 (unless (bzerop u1) 1101 (when snu (iputv u1 0 'bigneg)) 1102 (unless (ionep d) 1103 (setq lu (bbsize u1)) 1104 (setq carry 0) 1105 (vfor (from l lu 1 -1) 1106 (do (progn 1107 (setq temp(wquotientdouble carry (igetv u1 l) d)) 1108 (iputv u1 l temp) 1109 (setq carry *second-value*)))) 1110 (setq u1(trimbignum1 u1 lu)) % mk 1111 (unless (izerop carry) 1112 (bhardbug "remainder when unscaling" u v v1 q u1)))) 1113 (setq *second-value* 0) % clear 1114 (rplaca result q)(rplacd result u1) 1115 (return (cleanstack result)))) 1116 1117(de bhardbug (msg u v v1 q r) 1118 % Because the inputs to BHardDivide are probably rather large, I am not 1119 % going to rely on BldMsg to display them; 1120 (progn (prin2t "***** Internal error in BHardDivide") 1121 (prin2 "arg1=") (prin2t u) 1122 (prin2 "arg2=") (prin2t v) 1123 (prin2 "scaled arg2=") (prin2t v1) 1124 (prin2 "computed quotient=") (prin2t q) 1125 (prin2 "computed remainder=") (prin2t r) 1126 (stderror msg))) 1127 1128% ------------------------------ comparisons ---------------------- 1129 1130(de bgreaterp (u v) 1131 (cond ((bbminusp u) (if (bbminusp v) (bunsignedgreaterp v u) nil)) 1132 ((bbminusp v) u) 1133 (t (bunsignedgreaterp u v)))) 1134 1135(de blessp (u v) 1136 (cond ((bbminusp u) (if (bbminusp v) (bunsignedgreaterp u v) u)) 1137 ((bbminusp v) nil) 1138 (t (bunsignedgreaterp v u)))) 1139 1140(de bgeq (u v) 1141 (cond ((bbminusp u) (if (bbminusp v) (bunsignedgeq v u) nil)) 1142 ((bbminusp v) u) 1143 (t (bunsignedgeq u v)))) 1144 1145(de bleq (u v) 1146 (cond ((bbminusp u) (if (bbminusp v) (bunsignedgeq u v) u)) 1147 ((bbminusp v) nil) 1148 (t (bunsignedgeq v u)))) 1149 1150(de bunsignedgreaterp (u v) 1151 % Compare magnitudes of two bignums; 1152 (prog (lu lv i) 1153 (setq lu (bbsize u)) (setq lv (bbsize v)) 1154 (unless (eq lu lv) (if (igreaterp lu lv) (return u) (return nil))) 1155 (while (and (eq (igetv u lv) (igetv v lv)) (igreaterp lv 1)) 1156 (setq lv (isub1 lv))) 1157 (if (>u (igetv u lv) (igetv v lv)) 1158 (return(cleanstack u)) (return (cleanstack nil))))) 1159 1160(de bunsignedgeq (u v) 1161 % Compare magnitudes of two unsigned bignums; 1162 (prog (lu lv) 1163 (setq lu (bbsize u)) (setq lv (bbsize v)) 1164 (unless (eq lu lv) 1165 (if (igreaterp lu lv) (return u) (return nil))) 1166 (while (and (eq (igetv u lv) (igetv v lv)) (igreaterp lv 1)) 1167 (setq lv (isub1 lv))) 1168 (if (>u (igetv v lv) (igetv u lv)) 1169 (return (cleanstack nil)) (return (cleanstack u))))) 1170 1171(de badd1 (v) (bsmalladd v 1)) 1172 1173(de bsub1 (u) (bsmalldiff u 1)) 1174 1175% ------------------------BIG<->FLOAT conversions ----------------- 1176 1177(de bigit2float(b) 1178 (prog(bh) 1179 (setq bh (wshift b (minus (quotient bitsperword 2)))) 1180 (setq b (wshift (wshift b (quotient bitsperword 2)) 1181 (minus (quotient bitsperword 2)))) 1182 (setq b (intfloat b)) 1183 (setq bh (floattimes2 (intfloat bh) float-shift*)) 1184 (return (floatplus2 b bh)) )) 1185 1186(de floatfrombignum (v) 1187 (cond ((bzerop v) 0.00000E+000) 1188 ((or (bgreaterp v bigfloathi!*) (blessp v bigfloatlow!*)) 1189 (error 99 (list "Argument, " v " to FLOAT is too large"))) 1190 (t (prog (res sn i j base) 1191 (setq i (bbsize v)) 1192 (setq j (idifference i bigitsPerMantissa*)) 1193 (when (ilessp j 1) (setq j 1)) 1194 (setq base (float-expt floatbbase!* (isub1 j))) 1195 (setq sn (bbminusp v)) 1196 (setq res (floattimes2 (bigit2float (igetv v j)) base)) 1197 (setq j(1+ j)) 1198 (while (not(igreaterp j i)) 1199 (setq base (floattimes2 base floatbbase!*)) 1200 (setq res (floatplus2 res 1201 (floattimes2 (bigit2float (igetv v j)) base))) 1202 (setq j(1+ j))) 1203 (when sn (setq res (minus res))) 1204 (return (cleanstack res)))))) 1205 1206(compiletime (setq system_list!* (cons bitsperword system_list!*))) 1207 1208(compiletime (if_system IEEE 1209 1210(progn 1211 1212(define-constant ieeeshift (difference 12 bitsperword)) 1213(define-constant signshift (difference 1 bitsperword)) 1214(define-constant ieeebias 1023) 1215(define-constant ieeemask 2047) 1216(ds floathiword(x) (floathighorder(inf x))) 1217(ds floatloword(x) (floatloworder(inf x))) 1218 1219(if_system 32 (progn 1220 1221(ds ieeezerop(u) 1222 (and (weq (floathiword u) 0) 1223 (weq (floatloword u) 0))) 1224 1225(ds 1226 ieeemant 1227 (f) 1228 ((lambda (lf) 1229 (lor 1230 (lshift 1231 (wor 1232 (wshift (wand (floathiword f) 16#fffff) 16#6) 1233 (wshift lf (minus 16#1a))) 1234 16#1a) 1235 (wand (lshift (minus 16#1) (minus 16#6)) lf))) 1236 (floatloword f))) 1237 1238) % if_system 32 1239 1240(progn % if_system 64 1241 1242(ds ieeezerop(u) (weq (floathiword u) 0)) 1243(ds ieeemant (f) (wand (floathiword f) 16#fffffffffffff)) 1244 1245)) % if_system 64 1246 1247(ds ieeeexpt(u) 1248 (wdifference(wand ieeemask 1249 (wshift(floathiword u)ieeeshift)) 1250 ieeebias)) 1251 1252(ds ieeesign (u) (wshift (floathiword u) signshift)) 1253 1254))) % compiletime if_system IEEE 1255 1256(if_system IEEE 1257 1258(de bigfromfloat (x) 1259 (cond 1260 ((or (fixp x) (bigp x)) x) 1261 ((ieeezerop x) bzero*) 1262 (t(prog (m e) 1263 (setq m (ieeemant x)) 1264 (setq e (ieeeexpt x)) 1265 (when (neq e (minus 16#3ff)) 1266 (setq m (lor m ieee-hidden-bit*))) 1267 (when (eq (ieeesign x) 1) 1268 (if (bigp m) 1269 (iputv m 0 'bigneg) 1270 (setq m (minus m)))) 1271 (return (lshift m (idifference e 52))))))) 1272 1273 1274% case: not IEEE 1275 1276(de bigfromfloat (x) 1277 (if (or (fixp x) (bigp x)) 1278 x 1279 (prog (bigpart floatpart power sign thispart) 1280 (if (minusp x) 1281 (progn (setq sign -1) 1282 (setq x (minus x))) 1283 (setq sign 1)) 1284 (setq bigpart bzero!*) 1285 (while (and (neq x 0) (neq x 0.00000E+000)) 1286 (if (lessp x floatbbase-2*) 1287 (progn (setq bigpart (bplus2 bigpart (bnumaux (int2sys (fix x))))) 1288 (setq x 0)) 1289 (progn (setq floatpart x) 1290 (setq power 0) 1291 (while (geq floatpart floatbbase-2*) 1292 % get high end of number. 1293 (progn 1294 (setq floatpart (quotient floatpart floatbbase-2*)) 1295 (setq power (+w power (difference bitsperword 2))))) 1296 (setq thispart 1297 (btimes2 (btwopower power) 1298 (bnumaux (int2sys (fix floatpart))))) 1299 (setq x (difference x (floatfrombignum thispart))) 1300 (setq bigpart (bplus2 bigpart thispart))))) 1301 (when (minusp sign) 1302 (setq bigpart (bminus bigpart))) 1303 (return bigpart)))) 1304 1305) % if_system IEEE 1306 1307 1308(de float-expt (v1 n) 1309 % V1 is float, N is positive i-num 1310 (cond ((izerop n) 1.0000) 1311 ((ionep n) v1) 1312 (t (prog (v2) 1313 (setq v2 (float-expt v1 (wshift n -1))) % quotient n 2 1314 (if (izerop (wand n 1)) % remainder n 2 1315 (return (floattimes2 v2 v2)) 1316 (return (floattimes2 (floattimes2 v2 v1) v2))))))) 1317 1318 1319 1320% -------------------------------- Input and Output ---------------- 1321 1322(setq digit2letter!* 1323 '[48 49 50 51 52 53 54 55 56 57 65 66 67 68 69 70 71 72 73 74 75 76 1324 77 78 79 80 81 82 83 84 85 86 87 88 89 90]) 1325 1326 1327(de bchannelprin2 (channel v1) 1328 ((lambda (v2 myobase i ob sn) 1329 (while (not (izerop i)) 1330 (progn (iputv v2 i (igetv v1 i)) (setq i (isub1 i)))) 1331 (setq sn (bbminusp v1)) 1332 1333 (cond (sn (channelwritechar channel (char !-)))) 1334 (cond 1335 ((neq myobase 10) 1336 (channelwritesysinteger channel myobase 10) 1337 (channelwritechar channel (char !#)))) 1338 (setq ob (if (eq ob 10) 1000000 1339 (progn 1340 (setq ob (itimes2 ob ob)) % ob should be < 6 bits 1341 (setq ob (itimes2 (itimes2 ob ob) ob))))) 1342 (bprin channel v2 ob) 1343 (setq outputbase!* myobase )) 1344 1345 (gtpos (bbsize v1)) 1346 outputbase* (bbsize v1) outputbase* 0)) 1347 1348 1349% divide v by 10**6, if the quotient is zero the recursion is terminated 1350% and the most siginifcant digit is printed (and so on down the stack) 1351% otherwise save the big digit (=remainder) and recurse 1352(de bprin (channel v ob) 1353 ((lambda (digit) 1354 (cond 1355 ((bzerop v) (channelwritesysinteger channel digit outputbase!*)) 1356 (t (bprin channel v ob) (printdigit channel digit)))) 1357 (bdivide1000000ip v ob))) 1358 1359% print a base 10**6 digit 1360(de printdigit (channel d) 1361 (prog (d1 d2 d3 d4 d5 d6 digl ob) 1362 1363 (setq digl digit2letter!* ob outputbase*) %local copy 1364 1365 (setq d (wdivide d ob)) (setq d1 *second-value*) 1366 (setq d (wdivide d ob)) (setq d2 *second-value*) 1367 (setq d (wdivide d ob)) (setq d3 *second-value*) 1368 (setq d (wdivide d ob)) (setq d4 *second-value*) 1369 (setq d (wdivide d ob)) (setq d5 *second-value*) 1370 (setq d (wdivide d ob)) (setq d6 *second-value*) 1371 1372 (channelwritechar channel (igetv digl d6)) 1373 (channelwritechar channel (igetv digl d5)) 1374 (channelwritechar channel (igetv digl d4)) 1375 (channelwritechar channel (igetv digl d3)) 1376 (channelwritechar channel (igetv digl d2)) 1377 (return (channelwritechar channel (igetv digl d1))))) 1378 1379% divide the bignum v1 (of length l1) by 100000, except the quotient is 1380% accumulated in the same place, the remainder, of course, ripples 1381% down to the bottom. Because the argument is modified there is no need 1382% to CONS up a result, but simply return the remainder. 1383 1384(de bdivide1000000ip (v1 c) 1385 (prog (i r l1) 1386 (setq l1 (bbsize v1)) 1387 (setq r 0) 1388 (setq i l1) 1389 (while (not (izerop i)) 1390 (iputv v1 i (wquotientdouble r (igetv v1 i) c)) 1391 (setq r *second-value*) 1392 (setq i (isub1 i))) 1393 (setq *second-value* 0) % clear carry 1394 (trimbignum1 v1 l1) 1395 (when (bbminusp v1) (setq r (iminus r))) 1396 (return (cleanstack r)))) 1397 1398 1399(de bread (s radix sn) 1400 % radix is < Bbase!* 1401 %s=string of digits, radix=base, sn=1 or -1 1402 (prog (sz accu res ch j million shortstring) 1403 (setq shortstring t) 1404 (setq j 2) 1405 (setq million 1) 1406 (setq sz (size s)) 1407 (setq res (gtpos 1)) 1408 (setq ch (indx s 0)) 1409 (when (and (igeq ch (char a)) (ileq ch (char z))) 1410 (setq ch (+w (idifference ch (char a)) 10))) 1411 (when (and (igeq ch (char 0)) (ileq ch (char 9))) 1412 (setq ch (idifference ch (char 0)))) 1413 (iputv res 1 0) 1414 (setq accu ch) 1415 (ifor (from i 1 sz 1) 1416 (do (progn (setq ch (indx s i)) 1417 (when (and (igeq ch (char a)) (ileq ch (char z))) 1418 (setq ch 1419 (idifference ch (idifference (char a) 10)))) 1420 (when (and (igeq ch (char 0)) (ileq ch (char 9))) 1421 (setq ch (idifference ch (char 0)))) 1422 (setq accu (+w (itimes2 radix accu) ch )) 1423 (if (eq j 6) % do breadadd in larger portions 1424 (progn 1425 (setq shortstring nil) 1426 (setq j 1) 1427 (setq res (breadadd res 1428 (itimes2 radix million) accu)) 1429 (setq million 1) 1430 (setq accu 0)) 1431 (setq j (1+ j )) 1432 (setq million (itimes2 million radix)))))) 1433 1434 (when shortstring (when (iminusp sn) (setq accu (iminus accu))) 1435 (return (bnum accu))) 1436 (setq res (breadadd res million accu)) 1437 (when (iminusp sn) (setq res (bminus res))) 1438 (return res))) 1439 1440(de breadadd (v radix ch) 1441 (setq v (bsmalltimes2 v radix 0)) 1442 (setq v (bsmalladd v ch))) 1443 1444(de bsmalladd (v c) 1445 %V big, C fix. 1446 (cond ((izerop c) (return v)) 1447 ((bzerop v) (return (int2big c))) 1448 ((bbminusp v) (bminus (bsmalldiff (bminus v) c))) 1449 ((iminusp c) (bsmalldiff v (iminus c))) 1450 (t (prog (v1 l1) 1451 (setq carry!* 0) 1452 (setq l1 (bbsize v)) 1453 (setq v1 (gtpos (1+ l1))) 1454 (vfor (from i 1 l1 1) 1455 (do 1456 (progn 1457 (iputv v1 i (+c c (igetv v i))) 1458 (setq c 0) 1459 ))) 1460 (if (ionep carry!*) 1461 (iputv v1 (1+ l1) 1) 1462 (return (cleanstack (trimbignum1 v1 l1)))) 1463 (return (cleanstack v1)))))) 1464 1465(de bnum (n) 1466 % Creates a Bignum of one BETA digit, value N. 1467 % N is POS or NEG 1468 (if (bigp n) n (bnumaux n))) 1469 1470(de bnumaux (n) 1471 % Creates a Bignum of one BIGIT value N. 1472 % N is POS or NEG 1473 (prog (b tagn sgn) 1474 (when (izerop n) (return (gtpos 0))) 1475 (setq sgn (iminusp n)) 1476 (setq tagn (tag n)) % bnumaux must be save at gc time! 1477 (setq n (inf n)) 1478 (cond (sgn (setq b (gtneg 1))) % e.g. here 1479 (t (setq b (gtpos 1)))) 1480 (setq n (mkitem tagn n)) 1481 (when sgn (setq n (iminus n))) 1482 (iputv b 1 n) 1483 (return b))) 1484 1485(de bsmalldiff (v c) 1486 %V big, C fix 1487 (cond ((izerop c) v) 1488 ((bzerop v) (int2big (iminus c))) 1489 ((bbminusp v) (bminus (bsmalladd (bminus v) c))) 1490 ((iminusp c) (bsmalladd v (iminus c))) 1491 (t (prog (v1 l1 cc) 1492 (setq cc c carry* 0) 1493 (setq l1 (bbsize v)) 1494 (setq v1 (gtpos l1)) 1495 (vfor (from i 1 l1 1) 1496 (do (progn 1497 (iputv v1 i (-c (igetv v i) cc)) 1498 (setq cc 0) 1499 ))) 1500 (unless (izerop carry!*) 1501 (stderror (bldmsg " BSmallDiff V<C %p %p%n" v c))) 1502 (setq carry!* 0) % clear bigit 1503 (return (cleanstack (trimbignum1 v1 l1))))))) 1504 1505(de int2big (n) 1506 % Creates BigNum of value N. 1507 % From any N, BETA,INUM,FIXNUM or BIGNUM 1508 (case (tag n) ((negint-tag posint-tag) (sys2big n)) 1509 ((fixnum-tag) (sys2big (fixval (fixinf n)))) 1510 ((bignum-tag) n) 1511 (nil (nonintegererror n 'int2big)))) 1512 1513% ------------------------- INSTALL --------------------------------- 1514 1515% Now Install Interfacing 1516 (prin2t '"SetupGlobals") 1517 (setbits bitsperword) 1518 (prin2t '" ... done") 1519 1520% ------------------------- input / output --------------------------- 1521 1522(loadtime (progn (setq staticbig!* (gtwarray 10)))) 1523 1524% Assume dont need more than 10 slots to represent a BigNum 1525% Version of SYSint 1526% -- Output--- 1527% MLG Change to interface to Recursive hooks, added for 1528% Prinlevel stuff 1529(copyd 'oldchannelprin1 'recursivechannelprin1) 1530(copyd 'oldchannelprin2 'recursivechannelprin2) 1531 1532(de recursivechannelprin1 (channel u level) 1533 (if (bigp u) 1534 (bchannelprin2 channel u) 1535 (oldchannelprin1 channel u level)) 1536 u) 1537 1538(de recursivechannelprin2 (channel u level) 1539 (if (bigp u) 1540 (bchannelprin2 channel u) 1541 (oldchannelprin2 channel u level)) 1542 u) 1543 1544% ------------------------ interface ------------------------------- 1545 1546 1547% simple conventions for integer types in nbig30: 1548 1549% 1) there is no legal bignum with less than two bigits. 1550% 2) there is no legal fixnum with more than 30 significant. 1551 1552% integers with less than 28 significant bits are sysints, and 1553% 28 .. 30 bits are fixnums. 1554 1555(compiletime 1556 (progn 1557 (ds checkifreallybigpair (uu) 1558 ((lambda (UUU) (rplaca UUU (checkifreallybig (car UUU))) 1559 (rplacd UUU (checkifreallybig (cdr UUU)))) uu)) 1560)) 1561 1562(de checkifreallybig (uu) 1563 % If BIGNUM result is in older FIXNUM or INUM range 1564 (let((s (bbsize uu))) 1565 (cond ((eq s 0) 0) 1566 ((or (wgreaterp s 1) % more than one word 1567 (wlessp (igetv uu 1) 0)) % 1st bit set 1568 uu) 1569 (t (sys2int (big2sysaux uu)))))) 1570 1571(de checkifreallybigpair (vv) 1572 % Used to process DIVIDE 1573 (cons (checkifreallybig (car vv)) (checkifreallybig (cdr vv)))) 1574 1575(de checkifreallybigornil (uu) 1576 % Used for EXTRA-boolean tests 1577 (if (null uu) 1578 nil 1579 (checkifreallybig uu))) 1580 1581 1582(de bigplus2 (u v) (checkifreallybig (bplus2 u v))) 1583(de bigdifference (u v) (checkifreallybig (bdifference u v))) 1584(de bigtimes2 (u v) (checkifreallybig (btimes2 u v))) 1585(de bigdivide (u v) (checkifreallybigpair (bdivide u v))) 1586(de bigquotient (u v) (checkifreallybig (bquotient u v))) 1587(de bigremainder (u v) (checkifreallybig (bremainder u v))) 1588(de bigland (u v) (checkifreallybig (bland u v))) 1589(de biglor (u v) (checkifreallybig (blor u v))) 1590(de biglxor (u v) (checkifreallybig (blxor u v))) 1591(de biglshift (u v) (checkifreallybig (blshift u v))) 1592 1593(de lshift (u v) 1594 (setq v (int2sys v)) % bigger numbers make no sense as shift amount 1595 (if (betap u) 1596 (cond ((wleq v (minus bitsperword)) 0) 1597 ((wlessp v 0) (sys2int (wshift u v))) 1598 ((wlessp v (iquotient bitsperword 2)) (sys2int (wshift u v))) 1599 (t (biglshift (sys2big u) v))) 1600 % Use int2big, not sys2big, since we might have fixnums. 1601 (biglshift (int2big u) v))) 1602 1603(copyd 'lsh 'lshift) 1604 1605(de biggreaterp (u v) (checkifreallybigornil (bgreaterp u v))) 1606(de biglessp (u v) (checkifreallybigornil (blessp u v))) 1607(de bigadd1 (u) (checkifreallybig (badd1 u))) 1608(de bigsub1 (u) (checkifreallybig (bsub1 u))) 1609(de biglnot (u) (checkifreallybig (blnot u))) 1610(de bigminus (u) (checkifreallybig (bminus u))) 1611(de bigminusp (u) (checkifreallybigornil (bminusp u))) 1612(de bigonep (u) (checkifreallybigornil (bonep u))) 1613(de bigzerop (u) (checkifreallybigornil (bzerop u))) 1614 1615% ---------------------------- GCDN ------------------------------- 1616% 1617% Lehmers algorithm for numerical gcd 1618% Knuth Vol 2 (2nd ed.) page 329 1619 1620(commentoutcode 1621 1622(de biggcdn (u v) 1623 % We are sure that both, u and v are "true" bignums with at least 1624 % two bigits. 1625 (prog(u^ v^ A B C D TT q ttt w lu lv) 1626 (when (bbminusp u)(setq u (bminus u))) 1627 (when (bbminusp v)(setq v (bminus v))) 1628 L1 (when (or (not(bigp u)) (not (bigp v))) 1629 (return (cleanstack (biggcdn1 u v)))) 1630 (setq lu (bbsize u) lv (bbsize v)) 1631 (cond ((igreaterp lu lv) (setq tt(remainder u v)) 1632 (setq u v v tt) 1633 (go L1)) 1634 ((ilessp lu lv) (setq v (remainder v u)) 1635 (go L1))) 1636 % now u and v have equal length (first bigit of same order) 1637 (setq u^ (igetv u lu) v^ (igetv v lv)) 1638 (setq A 1 B 0 C 0 D 1) 1639 L2 (when (or (izerop (+w v^ C)) (izerop (+w v^ D))) 1640 (go L4)) 1641 (setq q (iquotient(+w u^ A)(+w v^ C))) 1642 (when (not(weq q (iquotient(+w u^ B)(+w v^ D)))) 1643 (go L4)) 1644 L3 (setq TT (if (weq C 0) A (idifference A(itimes2 q C))) 1645 A C 1646 C TT 1647 TT (if (weq D 0) B (idifference B(itimes2 q D))) 1648 B D 1649 D TT 1650 TT (idifference u^ (itimes2 q v^)) 1651 u^ v^ 1652 v^ TT) 1653 (go L2) 1654 L4 (if (weq B 0) 1655 (if (weq 0 (setq TT (remainder u v))) (return (cleanstack v)) 1656 (setq u v v TT)) 1657 (progn 1658 (setq u^ (setq v^ (setq q 0))) 1659 (setq TT (cond((eq A 0) 0) 1660 ((eq A 1) u) 1661 (t (times A u)))) 1662 (setq TT (cond((eq B 0) TT) 1663 ((eq B 1) (plus TT v)) 1664 (t (plus TT (times B v))))) 1665 (setq w (cond((eq C 0) 0) 1666 ((eq C 1) u) 1667 (t (times C u)))) 1668 (setq w (cond ((eq D 0) w) 1669 ((eq D 1) (plus w v)) 1670 (t (plus w (times D v))))) 1671 (setq u Tt v w))) 1672 (go L1))) 1673 1674) 1675 1676(de biggcdn(u v) 1677 (when (bbminusp u)(setq u (bminus u))) 1678 (when (bbminusp v)(setq v (bminus v))) 1679 (cleanstack (biggcdn0 u v))) 1680 1681(de biggcdn0(u v) 1682 (cond ((eq v 0) u) 1683 ((or (not(bigp u)) (not (bigp v))) (biggcdn1 u v)) 1684 (t (biggcdn0 v (bigremainder u v))))) 1685 1686 1687(de biggcdn1(u v) % handles simple cases: at least one number not big 1688 % gcdn and wgcdn are defined in the arithmetic 1689 (cond((izerop v) u) 1690 ((and (intp u)(intp v))(wgcdn u v)) 1691 ((not(bigp u)) 1692 (cond((not (bigp v)) (gcdn u v)) 1693 (T(setq v (remainder v u)) 1694 (if (izerop v) u (gcdn u v) )))) 1695 (t %(not (bigp v)) 1696 (setq u (remainder u v)) 1697 (if (izerop u) v (gcdn v u v))))) 1698 1699% ---- Input ---- 1700(de makestringintolispinteger (s radix sn) 1701 (checkifreallybig (bread s radix sn))) 1702 1703(de int2sys (n) 1704 % Convert a random FIXed number to WORD Integer 1705 (case (tag n) ((posint-tag negint-tag) n) 1706 ((fixnum-tag) (fixval (fixinf n))) 1707 ((bignum-tag) (big2sysaux n)) 1708 (nil (nonnumber1error n 'int2sys)))) 1709 1710(de sys2big (n) 1711 % Convert a SYSint to a BIG 1712 % Must NOT use generic arith here 1713 % Careful that no GC if this BIGger than INUM 1714 (prog (sn a b) 1715 (when (weq n 0) (return bzero!*)) 1716 (when (wlessp n 0) (setq sn t)) 1717 % Careful handling of -N in case have largest NEG, not just 1718 % flip sign 1719 (when sn (if (izerop (wshift n 1)) (return bigsyslow!*) 1720 (setq n (iminus n)))) 1721 (setq a staticbig!*) % Grab the base 1722 (iputv a 1 n) 1723 (setq n 0) 1724 (setq b (if (izerop n) 1 2)) 1725 (setq b (if sn (gtneg b) (gtpos b))) 1726 (when (not(izerop n)) (iputv b 2 n)) 1727 (iputv b 1 (igetv a 1)) 1728 (return b))) 1729 1730% Coercion/Transfer Functions 1731(copyd 'oldfloatfix 'floatfix) 1732 1733(de floatfix (u) 1734 % Careful of sign and range 1735 (if (and (leq floatsyslow!* u) (leq u floatsyshi!*)) 1736 (oldfloatfix u) 1737 (bigfromfloat u))) 1738 1739(de betap (x) 1740 % test if NUMBER in reduced INUM range 1741 (if (intp x) (and (wleq x betahi!*) (wgeq x betalow!*)) nil)) 1742 1743(de betarangep (x) 1744 % Test if SYSINT in reduced INUM range 1745 (if (wleq x betahi!*) (wgeq x betalow!*) nil)) 1746 1747(de beta2p (x y) 1748 % Check for 2 argument arithmetic functions 1749 (when (betap x) (betap y))) 1750 1751% Also from .BUILD file 1752% Now install the important globals for this machine 1753 1754% Largest representable float. 1755 1756(if_system VAX 1757 (progn 1758 (setq BigFloatHi!* (btimes2 (bdifference (btwopower 67) (btwopower 11)) 1759 (btwopower 60))) 1760 (setq BigFloatLow!* (bminus BigFloatHi!*)))) 1761 1762(if_system IEEE 1763 (progn % IEEE double arithmetic 1764 (setq BigFloatHi!* (btimes2 (bsub1 (btwopower 53)) (btwopower 971))) 1765 (setq BigFloatLow!* (bminus BigFloatHi!*)))) 1766 1767 1768(prog (y syshi syslo) 1769 % Lowest value of Ai 1770 % here assume 2's complement 1771 (setq y (twopower (idifference bitsperword 2))) 1772 % eg, 36 bits, 2^35-1=2^34+2^34-1 1773 (setq syshi (plus y (difference y 1))) 1774 (setq y (minus y)) 1775 (setq syslo (plus y y)) 1776 (setq FloatSysHi!* (float SysHi)) 1777 (setq FloatSysLow!* (float SysLo)) 1778) 1779 1780(if_system IEEE 1781 (let ((one 1)) 1782 (setq ieee-hidden-bit* (lshift bone* 16#34)))) 1783 1784(setq bigsyslow!* (minus (bllshift bone!* bitsperword))) 1785 1786%------------------------ test material --------------------------- 1787 1788(bigtest 1789 1790(de l2big(l) 1791 % make bignum from list of numbers 1792 (prog(u n) 1793 (setq n (length(cdr l))) 1794 (setq u (if (igreaterp (car l) -1) 1795 (gtpos n) 1796 (gtneg n))) 1797 (setq l (reverse (cdr l))) 1798 (for (from i 1 n 1) 1799 (do (iputv u i (int2sys(nth l i))))) 1800 (return u))) 1801 1802(de big2l(u) 1803 %convert bignum in list 1804 (prog (l n) 1805 (setq l (list (if(bbminusp u) -1 0))) 1806 (setq n (bbsize u)) 1807 (for (from i 1 n 1) 1808 (do (setq l (cons (sys2int (igetv u i)) l)))) 1809 (return l) )) 1810 1811(de d2l(u) (list (sys2int (doublehigh u)) 1812 (sys2int (doublelow u)))) 1813 1814) 1815