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