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