1(defpackage "SB-GMP" 2 (:use "COMMON-LISP" "SB-ALIEN" "SB-BIGNUM") 3 ;; we need a few very internal symbols 4 (:import-from "SB-BIGNUM" 5 "%BIGNUM-0-OR-PLUSP" "%NORMALIZE-BIGNUM" 6 "NEGATE-BIGNUM-IN-PLACE") 7 (:export 8 ;; bignum integer operations 9 #:mpz-add 10 #:mpz-sub 11 #:mpz-mul 12 #:mpz-mod 13 #:mpz-mul-2exp ; shift left 14 #:mpz-cdiv 15 #:mpz-fdiv 16 #:mpz-fdiv-2exp ; shift right 17 #:mpz-tdiv 18 #:mpz-powm 19 #:mpz-pow 20 #:mpz-gcd 21 #:mpz-lcm 22 #:mpz-sqrt 23 #:mpz-probably-prime-p 24 #:mpz-nextprime 25 #:mpz-fac 26 ;; the following functions are GMP >= 5.1 only 27 #:mpz-2fac 28 #:mpz-mfac 29 #:mpz-primorial 30 ;; number theoretic functions 31 #:mpz-remove 32 #:mpz-bin 33 #:mpz-fib2 34 ;; random number generation 35 #:make-gmp-rstate 36 #:make-gmp-rstate-lc 37 #:rand-seed 38 #:random-bitcount 39 #:random-int 40 ;; ratio arithmetic 41 #:mpq-add 42 #:mpq-sub 43 #:mpq-mul 44 #:mpq-div 45 ;; (un)installer functions 46 ; these insert/remove the runtime patch in SBCL's runtime 47 #:install-gmp-funs 48 #:uninstall-gmp-funs 49 ; these also load/unload the shared library and setup/clear 50 ; hooks to handle core saves 51 #:load-gmp 52 #:unload-gmp 53 ;; special variables 54 #:*gmp-version* 55 #:*gmp-disabled* 56 )) 57 58(in-package "SB-GMP") 59 60(defvar *gmp-disabled* nil) 61 62(defconstant +bignum-raw-area-offset+ 63 (- (* sb-vm:bignum-digits-offset sb-vm:n-word-bytes) 64 sb-vm:other-pointer-lowtag)) 65 66(declaim (inline bignum-data-sap)) 67(defun bignum-data-sap (x) 68 (declare (type bignum x)) 69 (sb-sys:sap+ (sb-sys:int-sap (sb-kernel:get-lisp-obj-address x)) 70 +bignum-raw-area-offset+)) 71 72(defun try-load-shared-object (pathname) 73 (handler-case 74 (load-shared-object pathname :dont-save t) 75 (error (e) 76 (declare (ignore e)) 77 nil))) 78 79(defun %load-gmp () 80 (or (some #'try-load-shared-object 81 #-(or win32 darwin) '("libgmp.so" "libgmp.so.10" "libgmp.so.3") 82 #+darwin '("libgmp.dylib" "libgmp.10.dylib" "libgmp.3.dylib") 83 #+win32 '("libgmp.dll" "libgmp-10.dll" "libgmp-3.dll")) 84 (warn "GMP not loaded."))) 85 86(defvar *gmp-features* nil) 87(defvar *gmp-version* nil) 88 89;; We load only the library right now to avoid undefined alien 90;; style warnings 91(%load-gmp) 92 93 94;;; types and initialization 95 96(define-alien-type gmp-limb 97 #-(and win32 x86-64) unsigned-long 98 #+(and win32 x86-64) unsigned-long-long) 99 100(deftype ui () 101 #-(and win32 x86-64) '(unsigned-byte #.sb-vm:n-word-bits) 102 #+(and win32 x86-64) '(unsigned-byte 32)) 103 104(deftype si () 105 #-(and win32 x86-64) '(signed-byte #.sb-vm:n-word-bits) 106 #+(and win32 x86-64) '(signed-byte 32)) 107 108(define-alien-type nil 109 (struct gmpint 110 (mp_alloc int) 111 (mp_size int) 112 (mp_d (* gmp-limb)))) 113 114;; Section 3.6 "Memory Management" of the GMP manual states: "mpz_t 115;; and mpq_t variables never reduce their allocated space. Normally 116;; this is the best policy, since it avoids frequent 117;; reallocation. Applications that need to return memory to the heap 118;; at some particular point can use mpz_realloc2, or clear variables 119;; no longer needed." 120;; 121;; We can therefore allocate a bignum of sufficiant size and use the 122;; space for GMP computations without the need for memory transfer 123;; from C to Lisp space. 124(declaim (inline z-to-bignum z-to-bignum-neg)) 125 126(defun z-to-bignum (b count) 127 "Convert GMP integer in the buffer of a pre-allocated bignum." 128 (declare (optimize (speed 3) (space 3) (safety 0)) 129 (type bignum b) 130 (type bignum-length count)) 131 (if (zerop count) 132 0 133 (the unsigned-byte (%normalize-bignum b count)))) 134 135(defun z-to-bignum-neg (b count) 136 "Convert to twos complement int the buffer of a pre-allocated 137bignum." 138 (declare (optimize (speed 3) (space 3) (safety 0)) 139 (type bignum b) 140 (type bignum-length count)) 141 (negate-bignum-in-place b) 142 (the (integer * 0) (%normalize-bignum b count))) 143 144;;; conversion functions that also copy from GMP to SBCL bignum space 145(declaim (inline gmp-z-to-bignum)) 146 147(defun gmp-z-to-bignum (z b count) 148 "Convert and copy a positive GMP integer into the buffer of a 149pre-allocated bignum. The allocated bignum-length must be (1+ COUNT)." 150 (declare (optimize (speed 3) (space 3) (safety 0)) 151 (type (alien (* gmp-limb)) z) 152 (type bignum b) 153 (type bignum-length count)) 154 (dotimes (i count (%normalize-bignum b (1+ count))) 155 (%bignum-set b i (deref z i)))) 156 157(declaim (inline blength bassert) 158 (ftype (function (integer) (values bignum-length &optional)) blength) 159 (ftype (function (integer) (values bignum &optional)) bassert)) 160 161(defun blength (a) 162 (declare (optimize (speed 3) (space 3) (safety 0))) 163 (etypecase a 164 (fixnum 1) 165 (t (%bignum-length a)))) 166 167(defun bassert (a) 168 (declare (optimize (speed 3) (space 3) (safety 0))) 169 (etypecase a 170 (fixnum (make-small-bignum a)) 171 (t a))) 172 173;;;; rationals 174(define-alien-type nil 175 (struct gmprat 176 (mp_num (struct gmpint)) 177 (mp_den (struct gmpint)))) 178 179;;; memory initialization functions to support non-alloced results 180;;; since an upper bound cannot always correctly predetermined 181;;; (e.g. the memory required for the fib function exceed the number 182;;; of limbs that are be determined through the infamous Phi-relation 183;;; resulting in a memory access error. 184 185;; use these for non-prealloced bignum values, but only when 186;; ultimately necessary since copying back into bignum space a the end 187;; of the operation is about three times slower than the shared buffer 188;; approach. 189(declaim (inline __gmpz_init __gmpz_clear)) 190(define-alien-routine __gmpz_init void 191 (x (* (struct gmpint)))) 192 193(define-alien-routine __gmpz_clear void 194 (x (* (struct gmpint)))) 195 196 197;;; integer interface functions 198(defmacro define-twoarg-mpz-funs (funs) 199 (loop for i in funs collect `(define-alien-routine ,i void 200 (r (* (struct gmpint))) 201 (a (* (struct gmpint)))) 202 into defines 203 finally (return `(progn 204 (declaim (inline ,@funs)) 205 ,@defines)))) 206 207(defmacro define-threearg-mpz-funs (funs) 208 (loop for i in funs collect `(define-alien-routine ,i void 209 (r (* (struct gmpint))) 210 (a (* (struct gmpint))) 211 (b (* (struct gmpint)))) 212 into defines 213 finally (return `(progn 214 (declaim (inline ,@funs)) 215 ,@defines)))) 216 217(defmacro define-fourarg-mpz-funs (funs) 218 (loop for i in funs collect `(define-alien-routine ,i void 219 (r (* (struct gmpint))) 220 (a (* (struct gmpint))) 221 (b (* (struct gmpint))) 222 (c (* (struct gmpint)))) 223 into defines 224 finally (return `(progn 225 (declaim (inline ,@funs)) 226 ,@defines)))) 227 228 229(define-twoarg-mpz-funs (__gmpz_sqrt 230 __gmpz_nextprime)) 231 232(define-threearg-mpz-funs (__gmpz_add 233 __gmpz_sub 234 __gmpz_mul 235 __gmpz_mod 236 __gmpz_gcd 237 __gmpz_lcm)) 238 239(define-fourarg-mpz-funs (__gmpz_cdiv_qr 240 __gmpz_fdiv_qr 241 __gmpz_tdiv_qr 242 __gmpz_powm)) 243 244(declaim (inline __gmpz_mul_2exp 245 __gmpz_fdiv_q_2exp 246 __gmpz_pow_ui 247 __gmpz_probab_prime_p 248 __gmpz_fac_ui 249 __gmpz_2fac_ui 250 __gmpz_mfac_uiui 251 __gmpz_primorial_ui 252 __gmpz_remove 253 __gmpz_bin_ui 254 __gmpz_fib2_ui)) 255 256(define-alien-routine __gmpz_mul_2exp void 257 (r (* (struct gmpint))) 258 (b (* (struct gmpint))) 259 (e unsigned-long)) 260 261(define-alien-routine __gmpz_fdiv_q_2exp void 262 (r (* (struct gmpint))) 263 (b (* (struct gmpint))) 264 (e unsigned-long)) 265 266(define-alien-routine __gmpz_pow_ui void 267 (r (* (struct gmpint))) 268 (b (* (struct gmpint))) 269 (e unsigned-long)) 270 271(define-alien-routine __gmpz_probab_prime_p int 272 (n (* (struct gmpint))) 273 (reps int)) 274 275(define-alien-routine __gmpz_fac_ui void 276 (r (* (struct gmpint))) 277 (a unsigned-long)) 278 279(define-alien-routine __gmpz_2fac_ui void 280 (r (* (struct gmpint))) 281 (a unsigned-long)) 282 283(define-alien-routine __gmpz_mfac_uiui void 284 (r (* (struct gmpint))) 285 (n unsigned-long) 286 (m unsigned-long)) 287 288(define-alien-routine __gmpz_primorial_ui void 289 (r (* (struct gmpint))) 290 (n unsigned-long)) 291 292(define-alien-routine __gmpz_remove unsigned-long 293 (r (* (struct gmpint))) 294 (x (* (struct gmpint))) 295 (f (* (struct gmpint)))) 296 297(define-alien-routine __gmpz_bin_ui void 298 (r (* (struct gmpint))) 299 (n (* (struct gmpint))) 300 (k unsigned-long)) 301 302(define-alien-routine __gmpz_fib2_ui void 303 (r (* (struct gmpint))) 304 (a (* (struct gmpint))) 305 (b unsigned-long)) 306 307 308;; ratio functions 309(defmacro define-threearg-mpq-funs (funs) 310 (loop for i in funs collect `(define-alien-routine ,i void 311 (r (* (struct gmprat))) 312 (a (* (struct gmprat))) 313 (b (* (struct gmprat)))) 314 into defines 315 finally (return `(progn 316 (declaim (inline ,@funs)) 317 ,@defines)))) 318 319(define-threearg-mpq-funs (__gmpq_add 320 __gmpq_sub 321 __gmpq_mul 322 __gmpq_div)) 323 324 325;;;; SBCL interface 326 327;;; utility macros for GMP mpz variable and result declaration and 328;;; incarnation of associated SBCL bignums 329 330(defmacro with-mpz-results (pairs &body body) 331 (loop for (gres size) in pairs 332 for res = (gensym "RESULT") 333 collect `(when (> ,size sb-bignum::maximum-bignum-length) 334 (error "Size of result exceeds maxim bignum length")) into checks 335 collect `(,gres (struct gmpint)) into declares 336 collect `(,res (%allocate-bignum ,size)) 337 into resinits 338 collect `(setf (slot ,gres 'mp_alloc) (%bignum-length ,res) 339 (slot ,gres 'mp_size) 0 340 (slot ,gres 'mp_d) (bignum-data-sap ,res)) 341 into inits 342 collect `(if (minusp (slot ,gres 'mp_size)) ; check for negative result 343 (z-to-bignum-neg ,res ,size) 344 (z-to-bignum ,res ,size)) 345 into normlimbs 346 collect res into results 347 finally (return 348 `(progn 349 ,@checks 350 (let ,resinits 351 (sb-sys:with-pinned-objects ,results 352 (with-alien ,declares 353 ,@inits 354 ,@body 355 (values ,@normlimbs)))))))) 356 357(defmacro with-mpz-vars (pairs &body body) 358 (loop for (a ga) in pairs 359 for length = (gensym "LENGTH") 360 for plusp = (gensym "PLUSP") 361 for barg = (gensym "BARG") 362 for arg = (gensym "ARG") 363 collect `(,ga (struct gmpint)) into declares 364 collect `(,barg (bassert ,a)) into gmpinits 365 collect `(,plusp (%bignum-0-or-plusp ,barg (%bignum-length ,barg))) into gmpinits 366 collect `(,arg (if ,plusp ,barg (negate-bignum ,barg nil))) into gmpinits 367 collect `(,length (%bignum-length ,arg)) into gmpinits 368 collect arg into vars 369 collect `(setf (slot ,ga 'mp_alloc) ,length 370 (slot ,ga 'mp_size) 371 (progn ;; handle twos complements/ulong limbs mismatch 372 (when (zerop (%bignum-ref ,arg (1- ,length))) 373 (decf ,length)) 374 (if ,plusp ,length (- ,length))) 375 (slot ,ga 'mp_d) (bignum-data-sap ,arg)) 376 into gmpvarssetup 377 finally (return 378 `(with-alien ,declares 379 (let* ,gmpinits 380 (sb-sys:with-pinned-objects ,vars 381 ,@gmpvarssetup 382 ,@body)))))) 383 384(defmacro with-gmp-mpz-results (resultvars &body body) 385 (loop for gres in resultvars 386 for res = (gensym "RESULT") 387 for size = (gensym "SIZE") 388 collect size into sizes 389 collect `(,gres (struct gmpint)) into declares 390 collect `(__gmpz_init (addr ,gres)) into inits 391 collect `(,size (abs (slot ,gres 'mp_size))) 392 into resinits 393 collect `(when (> ,size (1- sb-bignum::maximum-bignum-length)) 394 (error "Size of result exceeds maxim bignum length")) into checks 395 collect `(,res (%allocate-bignum (1+ ,size))) 396 into resallocs 397 collect `(setf ,res (if (minusp (slot ,gres 'mp_size)) ; check for negative result 398 (- (gmp-z-to-bignum (slot ,gres 'mp_d) ,res ,size)) 399 (gmp-z-to-bignum (slot ,gres 'mp_d) ,res ,size))) 400 into copylimbs 401 collect `(__gmpz_clear (addr ,gres)) into clears 402 collect res into results 403 finally (return 404 `(with-alien ,declares 405 ,@inits 406 ,@body 407 (let ,resinits 408 (declare (type bignum-length ,@sizes)) 409 ,@checks 410 (let ,resallocs 411 ;; copy GMP limbs into result bignum 412 (sb-sys:with-pinned-objects ,results 413 ,@copylimbs) 414 ,@clears 415 (values ,@results))))))) 416 417;;; function definition and foreign function relationships 418(defmacro defgmpfun (name args &body body) 419 `(progn 420 (declaim (sb-ext:maybe-inline ,name)) 421 (defun ,name ,args 422 (declare (optimize (speed 3) (space 3) (safety 0)) 423 (type integer ,@args)) 424 ,@body))) 425 426 427;; SBCL/GMP functions 428(defgmpfun mpz-add (a b) 429 (with-mpz-results ((result (1+ (max (blength a) 430 (blength b))))) 431 (with-mpz-vars ((a ga) (b gb)) 432 (__gmpz_add (addr result) (addr ga) (addr gb))))) 433 434(defgmpfun mpz-sub (a b) 435 (with-mpz-results ((result (1+ (max (blength a) 436 (blength b))))) 437 (with-mpz-vars ((a ga) (b gb)) 438 (__gmpz_sub (addr result) (addr ga) (addr gb))))) 439 440(defgmpfun mpz-mul (a b) 441 (with-mpz-results ((result (+ (blength a) 442 (blength b)))) 443 (with-mpz-vars ((a ga) (b gb)) 444 (__gmpz_mul (addr result) (addr ga) (addr gb))))) 445 446(defgmpfun mpz-mul-2exp (a b) 447 (check-type b ui) 448 (with-mpz-results ((result (+ (1+ (blength a)) 449 (floor b sb-vm:n-word-bits)))) 450 (with-mpz-vars ((a ga)) 451 (__gmpz_mul_2exp (addr result) (addr ga) b)))) 452 453(defgmpfun mpz-mod (a b) 454 (with-mpz-results ((result (1+ (max (blength a) 455 (blength b))))) 456 (with-mpz-vars ((a ga) (b gb)) 457 (__gmpz_mod (addr result) (addr ga) (addr gb)) 458 (when (and (minusp (slot gb 'mp_size)) 459 (/= 0 (slot result 'mp_size))) 460 (__gmpz_add (addr result) (addr result) (addr gb)))))) 461 462(defgmpfun mpz-cdiv (n d) 463 (let ((size (1+ (max (blength n) 464 (blength d))))) 465 (with-mpz-results ((quot size) 466 (rem size)) 467 (with-mpz-vars ((n gn) (d gd)) 468 (__gmpz_cdiv_qr (addr quot) (addr rem) (addr gn) (addr gd)))))) 469 470(defgmpfun mpz-fdiv (n d) 471 (let ((size (1+ (max (blength n) 472 (blength d))))) 473 (with-mpz-results ((quot size) 474 (rem size)) 475 (with-mpz-vars ((n gn) (d gd)) 476 (__gmpz_fdiv_qr (addr quot) (addr rem) (addr gn) (addr gd)))))) 477 478(defgmpfun mpz-fdiv-2exp (a b) 479 (check-type b ui) 480 (with-mpz-results ((result (1+ (- (blength a) 481 (floor b sb-vm:n-word-bits))))) 482 (with-mpz-vars ((a ga)) 483 (__gmpz_fdiv_q_2exp (addr result) (addr ga) b)))) 484 485(defgmpfun mpz-tdiv (n d) 486 (let ((size (max (blength n) 487 (blength d)))) 488 (with-mpz-results ((quot size) 489 (rem size)) 490 (with-mpz-vars ((n gn) (d gd)) 491 (__gmpz_tdiv_qr (addr quot) (addr rem) (addr gn) (addr gd)))))) 492 493(defgmpfun mpz-pow (base exp) 494 (check-type exp (integer 0 #.most-positive-fixnum)) 495 (with-gmp-mpz-results (rop) 496 (with-mpz-vars ((base gbase)) 497 (__gmpz_pow_ui (addr rop) (addr gbase) exp)))) 498 499(defgmpfun mpz-powm (base exp mod) 500 (with-mpz-results ((rop (1+ (blength mod)))) 501 (with-mpz-vars ((base gbase) (exp gexp) (mod gmod)) 502 (__gmpz_powm (addr rop) (addr gbase) (addr gexp) (addr gmod))))) 503 504(defgmpfun mpz-gcd (a b) 505 (with-mpz-results ((result (min (blength a) 506 (blength b)))) 507 (with-mpz-vars ((a ga) (b gb)) 508 (__gmpz_gcd (addr result) (addr ga) (addr gb))))) 509 510(defgmpfun mpz-lcm (a b) 511 (with-mpz-results ((result (+ (blength a) 512 (blength b)))) 513 (with-mpz-vars ((a ga) (b gb)) 514 (__gmpz_lcm (addr result) (addr ga) (addr gb))))) 515 516(defgmpfun mpz-sqrt (a) 517 (with-mpz-results ((result (1+ (ceiling (blength a) 2)))) 518 (with-mpz-vars ((a ga)) 519 (__gmpz_sqrt (addr result) (addr ga))))) 520 521 522;;; Functions that use GMP-side allocated integers and copy the result 523;;; into a SBCL bignum at the end of the computation when the required 524;;; bignum length is known. 525(defun mpz-probably-prime-p (n &optional (reps 25)) 526 (declare (optimize (speed 3) (space 3) (safety 0)) 527 (type integer n reps)) 528 (check-type reps fixnum) 529 (with-mpz-vars ((n gn)) 530 (__gmpz_probab_prime_p (addr gn) reps))) 531 532(defgmpfun mpz-nextprime (a) 533 (with-gmp-mpz-results (prime) 534 (with-mpz-vars ((a ga)) 535 (__gmpz_nextprime (addr prime) (addr ga))))) 536 537(defgmpfun mpz-fac (n) 538 (check-type n ui) 539 (with-gmp-mpz-results (fac) 540 (__gmpz_fac_ui (addr fac) n))) 541 542(defgmpfun %mpz-2fac (n) 543 (check-type n ui) 544 (with-gmp-mpz-results (fac) 545 (__gmpz_2fac_ui (addr fac) n))) 546 547(defgmpfun %mpz-mfac (n m) 548 (check-type n ui) 549 (check-type m ui) 550 (with-gmp-mpz-results (fac) 551 (__gmpz_mfac_uiui (addr fac) n m))) 552 553(defgmpfun %mpz-primorial (n) 554 (check-type n ui) 555 (with-gmp-mpz-results (r) 556 (__gmpz_primorial_ui (addr r) n))) 557 558(defgmpfun mpz-remove-5.1 (n f) 559 (check-type f unsigned-byte 560 "handled by GMP prior to version 5.1") 561 (let* ((c 0) 562 (res (with-gmp-mpz-results (r) 563 (with-mpz-vars ((n gn) 564 (f gf)) 565 (setf c (__gmpz_remove (addr r) (addr gn) (addr gf))))))) 566 (values res c))) 567 568(defgmpfun mpz-remove (n f) 569 (let* ((c 0) 570 (res (with-gmp-mpz-results (r) 571 (with-mpz-vars ((n gn) 572 (f gf)) 573 (setf c (__gmpz_remove (addr r) (addr gn) (addr gf))))))) 574 (values res c))) 575 576(defun setup-5.1-stubs () 577 (macrolet ((stubify (name implementation &rest arguments) 578 `(setf (fdefinition ',name) 579 (if (member :sb-gmp-5.1 *gmp-features*) 580 (fdefinition ',implementation) 581 (lambda ,arguments 582 (declare (ignore ,@arguments)) 583 (error "~S is only available in GMP >= 5.1" 584 ',name)))))) 585 (stubify mpz-2fac %mpz-2fac n) 586 (stubify mpz-mfac %mpz-mfac n m) 587 (stubify mpz-primorial %mpz-primorial n)) 588 (unless (member :sb-gmp-5.1 *gmp-features*) 589 (setf (fdefinition 'mpz-remove) #'mpz-remove-5.1))) 590 591(defgmpfun mpz-bin (n k) 592 (check-type k ui) 593 (with-gmp-mpz-results (r) 594 (with-mpz-vars ((n gn)) 595 (__gmpz_bin_ui (addr r) (addr gn) k)))) 596 597(defgmpfun mpz-fib2 (n) 598 ;; (let ((size (1+ (ceiling (* n (log 1.618034 2)) 64))))) 599 ;; fibonacci number magnitude in bits is assymptotic to n(log_2 phi) 600 ;; This is correct for the result but appears not to be enough for GMP 601 ;; during computation (memory access error), so use GMP-side allocation. 602 (check-type n ui) 603 (with-gmp-mpz-results (fibn fibn-1) 604 (__gmpz_fib2_ui (addr fibn) (addr fibn-1) n))) 605 606 607;;;; Random bignum (mpz) generation 608 609;; we do not actually use the gestalt of the struct but need its size 610;; for allocation purposes 611(define-alien-type nil 612 (struct gmprandstate 613 (mp_seed (struct gmpint)) 614 (mp_alg int) 615 (mp_algdata (* t)))) 616 617(declaim (inline __gmp_randinit_mt 618 __gmp_randinit_lc_2exp 619 __gmp_randseed 620 __gmp_randseed_ui 621 __gmpz_urandomb 622 __gmpz_urandomm)) 623 624(define-alien-routine __gmp_randinit_mt void 625 (s (* (struct gmprandstate)))) 626 627(define-alien-routine __gmp_randinit_lc_2exp void 628 (s (* (struct gmprandstate))) 629 (a (* (struct gmpint))) 630 (c unsigned-long) 631 (exp unsigned-long)) 632 633(define-alien-routine __gmp_randseed void 634 (s (* (struct gmprandstate))) 635 (sd (* (struct gmpint)))) 636 637(define-alien-routine __gmp_randseed_ui void 638 (s (* (struct gmprandstate))) 639 (sd unsigned-long)) 640 641(define-alien-routine __gmpz_urandomb void 642 (r (* (struct gmpint))) 643 (s (* (struct gmprandstate))) 644 (bcnt unsigned-long)) 645 646(define-alien-routine __gmpz_urandomm void 647 (r (* (struct gmpint))) 648 (s (* (struct gmprandstate))) 649 (n (* (struct gmpint)))) 650 651(defstruct (gmp-rstate (:constructor %make-gmp-rstate)) 652 (ref (make-alien (struct gmprandstate)) 653 :type (alien (* (struct gmprandstate))) :read-only t)) 654 655(declaim (sb-ext:maybe-inline make-gmp-rstate 656 make-gmp-rstate-lc 657 rand-seed 658 random-bitcount 659 random-int)) 660 661(defun make-gmp-rstate () 662 "Instantiate a state for the GMP Mersenne-Twister random number generator." 663 (declare (optimize (speed 3) (space 3) (safety 0))) 664 (let* ((state (%make-gmp-rstate)) 665 (ref (gmp-rstate-ref state))) 666 (__gmp_randinit_mt ref) 667 (sb-ext:finalize state (lambda () (free-alien ref))) 668 state)) 669 670(defun make-gmp-rstate-lc (a c m2exp) 671 "Instantiate a state for the GMP linear congruential random number generator." 672 (declare (optimize (speed 3) (space 3))) 673 (check-type c ui) 674 (check-type m2exp ui) 675 (let* ((state (%make-gmp-rstate)) 676 (ref (gmp-rstate-ref state))) 677 (with-mpz-vars ((a ga)) 678 (__gmp_randinit_lc_2exp ref (addr ga) c m2exp)) 679 (sb-ext:finalize state (lambda () (free-alien ref))) 680 state)) 681 682(defun rand-seed (state seed) 683 "Initialize a random STATE with SEED." 684 (declare (optimize (speed 3) (space 3) (safety 0))) 685 (check-type state gmp-rstate) 686 (let ((ref (gmp-rstate-ref state))) 687 (cond 688 ((typep seed 'ui) 689 (__gmp_randseed_ui ref seed)) 690 ((typep seed '(integer 0 *)) 691 (with-mpz-vars ((seed gseed)) 692 (__gmp_randseed ref (addr gseed)))) 693 (t 694 (error "SEED must be a positive integer"))))) 695 696(defun random-bitcount (state bitcount) 697 "Return a random integer in the range 0..(2^bitcount - 1)." 698 (declare (optimize (speed 3) (space 3) (safety 0))) 699 (check-type state gmp-rstate) 700 (check-type bitcount ui) 701 (let ((ref (gmp-rstate-ref state))) 702 (with-mpz-results ((result (+ (ceiling bitcount sb-vm:n-word-bits) 2))) 703 (__gmpz_urandomb (addr result) ref bitcount)))) 704 705(defun random-int (state boundary) 706 "Return a random integer in the range 0..(boundary - 1)." 707 (declare (optimize (speed 3) (space 3))) 708 (check-type state gmp-rstate) 709 (let ((b (bassert boundary)) 710 (ref (gmp-rstate-ref state))) 711 (with-mpz-results ((result (1+ (%bignum-length b)))) 712 (with-mpz-vars ((b gb)) 713 (__gmpz_urandomm (addr result) ref (addr gb)))))) 714 715 716;;; Rational functions 717(declaim (inline %lsize)) 718(defun %lsize (minusp n) 719 (declare (optimize (speed 3) (space 3) (safety 0))) 720 "n must be a (potentially denormalized) bignum" 721 (let ((length (%bignum-length n))) 722 (when (zerop (%bignum-ref n (1- length))) 723 (decf length)) 724 (if minusp (- length) length))) 725 726(defmacro with-mpq-var (pair &body body) 727 (destructuring-bind (var mpqvar) pair 728 `(let* ((an (bassert (numerator ,var))) 729 (ad (bassert (denominator ,var))) 730 (asign (not (%bignum-0-or-plusp an (%bignum-length an))))) 731 (when asign 732 (setf an (negate-bignum an nil))) 733 (let* ((anlen (%lsize asign an)) 734 (adlen (%lsize NIL ad))) 735 (with-alien ((,mpqvar (struct gmprat))) 736 (sb-sys:with-pinned-objects (an ad) 737 (setf (slot (slot ,mpqvar 'mp_num) 'mp_size) anlen 738 (slot (slot ,mpqvar 'mp_num) 'mp_alloc) (abs anlen) 739 (slot (slot ,mpqvar 'mp_num) 'mp_d) 740 (bignum-data-sap an)) 741 (setf (slot (slot ,mpqvar 'mp_den) 'mp_size) adlen 742 (slot (slot ,mpqvar 'mp_den) 'mp_alloc) (abs adlen) 743 (slot (slot ,mpqvar 'mp_den) 'mp_d) 744 (bignum-data-sap ad)) 745 ,@body)))))) 746 747(defmacro defmpqfun (name gmpfun) 748 `(progn 749 (declaim (sb-ext:maybe-inline ,name)) 750 (defun ,name (a b) 751 (declare (optimize (speed 3) (space 3) (safety 0))) 752 (let ((size (+ (max (blength (numerator a)) 753 (blength (denominator a))) 754 (max (blength (numerator b)) 755 (blength (denominator b))) 756 3))) 757 (with-alien ((r (struct gmprat))) 758 (let ((num (%allocate-bignum size)) 759 (den (%allocate-bignum size))) 760 (sb-sys:with-pinned-objects (num den) 761 (setf (slot (slot r 'mp_num) 'mp_size) 0 762 (slot (slot r 'mp_num) 'mp_alloc) size 763 (slot (slot r 'mp_num) 'mp_d) (bignum-data-sap num)) 764 (setf (slot (slot r 'mp_den) 'mp_size) 0 765 (slot (slot r 'mp_den) 'mp_alloc) size 766 (slot (slot r 'mp_den) 'mp_d) (bignum-data-sap den)) 767 (let* ((an (bassert (numerator a))) 768 (ad (bassert (denominator a))) 769 (asign (not (%bignum-0-or-plusp an (%bignum-length an)))) 770 (bn (bassert (numerator b))) 771 (bd (bassert (denominator b))) 772 (bsign (not (%bignum-0-or-plusp bn (%bignum-length bn))))) 773 (when asign 774 (setf an (negate-bignum an nil))) 775 (when bsign 776 (setf bn (negate-bignum bn nil))) 777 (let* ((anlen (%lsize asign an)) 778 (adlen (%lsize NIL ad)) 779 (bnlen (%lsize bsign bn)) 780 (bdlen (%lsize NIL bd))) 781 (with-alien ((arga (struct gmprat)) 782 (argb (struct gmprat))) 783 (sb-sys:with-pinned-objects (an ad bn bd) 784 (setf (slot (slot arga 'mp_num) 'mp_size) anlen 785 (slot (slot arga 'mp_num) 'mp_alloc) (abs anlen) 786 (slot (slot arga 'mp_num) 'mp_d) 787 (bignum-data-sap an)) 788 (setf (slot (slot arga 'mp_den) 'mp_size) adlen 789 (slot (slot arga 'mp_den) 'mp_alloc) (abs adlen) 790 (slot (slot arga 'mp_den) 'mp_d) 791 (bignum-data-sap ad)) 792 (setf (slot (slot argb 'mp_num) 'mp_size) bnlen 793 (slot (slot argb 'mp_num) 'mp_alloc) (abs bnlen) 794 (slot (slot argb 'mp_num) 'mp_d) 795 (bignum-data-sap bn)) 796 (setf (slot (slot argb 'mp_den) 'mp_size) bdlen 797 (slot (slot argb 'mp_den) 'mp_alloc) (abs bdlen) 798 (slot (slot argb 'mp_den) 'mp_d) 799 (bignum-data-sap bd)) 800 (,gmpfun (addr r) (addr arga) (addr argb))))) 801 (locally (declare (optimize (speed 1))) 802 (sb-kernel::build-ratio (if (minusp (slot (slot r 'mp_num) 'mp_size)) 803 (z-to-bignum-neg num size) 804 (z-to-bignum num size)) 805 (z-to-bignum den size))))))))))) 806 807(defmpqfun mpq-add __gmpq_add) 808(defmpqfun mpq-sub __gmpq_sub) 809(defmpqfun mpq-mul __gmpq_mul) 810(defmpqfun mpq-div __gmpq_div) 811 812 813;;;; SBCL interface and integration installation 814(macrolet ((def (name original) 815 (let ((special (intern (format nil "*~A-FUNCTION*" name)))) 816 `(progn 817 (declaim (type function ,special) 818 (inline ,name)) 819 (defvar ,special (symbol-function ',original)) 820 (defun ,name (&rest args) 821 (apply (load-time-value ,special t) args)))))) 822 (def orig-mul multiply-bignums) 823 (def orig-truncate bignum-truncate) 824 (def orig-gcd bignum-gcd) 825 (def orig-lcm sb-kernel:two-arg-lcm) 826 (def orig-isqrt isqrt) 827 (def orig-two-arg-+ sb-kernel:two-arg-+) 828 (def orig-two-arg-- sb-kernel:two-arg--) 829 (def orig-two-arg-* sb-kernel:two-arg-*) 830 (def orig-two-arg-/ sb-kernel:two-arg-/) 831 (def orig-intexp sb-kernel::intexp)) 832 833;;; integers 834(defun gmp-mul (a b) 835 (declare (optimize (speed 3) (space 3)) 836 (type bignum a b) 837 (inline mpz-mul)) 838 (if (or (< (min (%bignum-length a) 839 (%bignum-length b)) 840 6) 841 *gmp-disabled*) 842 (orig-mul a b) 843 (mpz-mul a b))) 844 845(defun gmp-truncate (a b) 846 (declare (optimize (speed 3) (space 3)) 847 (type bignum a b) 848 (inline mpz-tdiv)) 849 (if (or (< (min (%bignum-length a) 850 (%bignum-length b)) 851 3) 852 *gmp-disabled*) 853 (orig-truncate a b) 854 (mpz-tdiv a b))) 855 856(defun gmp-lcm (a b) 857 (declare (optimize (speed 3) (space 3)) 858 (type integer a b) 859 (inline mpz-lcm)) 860 (if (or (and (typep a 'fixnum) 861 (typep b 'fixnum)) 862 *gmp-disabled*) 863 (orig-lcm a b) 864 (mpz-lcm a b))) 865 866(defun gmp-isqrt (n) 867 (declare (optimize (speed 3) (space 3)) 868 (type unsigned-byte n) 869 (inline mpz-sqrt)) 870 (if (or (typep n 'fixnum) 871 *gmp-disabled*) 872 (orig-isqrt n) 873 (mpz-sqrt n))) 874 875;;; rationals 876(defun gmp-two-arg-+ (x y) 877 (declare (optimize (speed 3) (space 3)) 878 (inline mpq-add)) 879 (if (and (or (typep x 'ratio) 880 (typep y 'ratio)) 881 (rationalp y) 882 (rationalp x) 883 (not *gmp-disabled*)) 884 (mpq-add x y) 885 (orig-two-arg-+ x y))) 886 887(defun gmp-two-arg-- (x y) 888 (declare (optimize (speed 3) (space 3)) 889 (inline mpq-sub)) 890 (if (and (or (typep x 'ratio) 891 (typep y 'ratio)) 892 (rationalp y) 893 (rationalp x) 894 (not *gmp-disabled*)) 895 (mpq-sub x y) 896 (orig-two-arg-- x y))) 897 898(defun gmp-two-arg-* (x y) 899 (declare (optimize (speed 3) (space 3)) 900 (inline mpq-mul)) 901 (if (and (or (typep x 'ratio) 902 (typep y 'ratio)) 903 (rationalp y) 904 (rationalp x) 905 (not *gmp-disabled*)) 906 (mpq-mul x y) 907 (orig-two-arg-* x y))) 908 909(defun gmp-two-arg-/ (x y) 910 (declare (optimize (speed 3) (space 3)) 911 (inline mpq-div)) 912 (if (and (rationalp x) 913 (rationalp y) 914 (not (eql y 0)) 915 (not *gmp-disabled*)) 916 (mpq-div x y) 917 (orig-two-arg-/ x y))) 918 919(defun gmp-intexp (base power) 920 (declare (inline mpz-mul-2exp mpz-pow)) 921 (check-type power (integer #.(1+ most-negative-fixnum) #.most-positive-fixnum)) 922 (cond 923 ((or (and (integerp base) 924 (< (abs power) 1000) 925 (< (blength base) 4)) 926 *gmp-disabled*) 927 (orig-intexp base power)) 928 (t 929 (cond ((minusp power) 930 (/ (the integer (gmp-intexp base (- power))))) 931 ((eql base 2) 932 (mpz-mul-2exp 1 power)) 933 ((typep base 'ratio) 934 (sb-kernel::%make-ratio (gmp-intexp (numerator base) power) 935 (gmp-intexp (denominator base) power))) 936 (t 937 (mpz-pow base power)))))) 938 939;;; installation 940(defun install-gmp-funs () 941 (sb-ext:without-package-locks 942 (macrolet ((def (destination source) 943 `(setf (fdefinition ',destination) 944 (fdefinition ',source)))) 945 (def sb-bignum:multiply-bignums gmp-mul) 946 (def sb-bignum:bignum-truncate gmp-truncate) 947 (def sb-bignum:bignum-gcd mpz-gcd) 948 (def sb-kernel:two-arg-lcm gmp-lcm) 949 (def sb-kernel:two-arg-+ gmp-two-arg-+) 950 (def sb-kernel:two-arg-- gmp-two-arg--) 951 (def sb-kernel:two-arg-* gmp-two-arg-*) 952 (def sb-kernel:two-arg-/ gmp-two-arg-/) 953 (def isqrt gmp-isqrt) 954 (def sb-kernel::intexp gmp-intexp))) 955 (values)) 956 957(defun uninstall-gmp-funs () 958 (sb-ext:without-package-locks 959 (macrolet ((def (destination source) 960 `(setf (fdefinition ',destination) 961 ,(intern (format nil "*~A-FUNCTION*" source))))) 962 (def multiply-bignums orig-mul) 963 (def bignum-truncate orig-truncate) 964 (def bignum-gcd orig-gcd) 965 (def sb-kernel:two-arg-lcm orig-lcm) 966 (def sb-kernel:two-arg-+ orig-two-arg-+) 967 (def sb-kernel:two-arg-- orig-two-arg--) 968 (def sb-kernel:two-arg-* orig-two-arg-*) 969 (def sb-kernel:two-arg-/ orig-two-arg-/) 970 (def isqrt orig-isqrt) 971 (def sb-kernel::intexp orig-intexp))) 972 (values)) 973 974(defun load-gmp (&key (persistently t)) 975 (setf *gmp-features* nil 976 *gmp-version* nil 977 *features* (set-difference *features* '(:sb-gmp :sb-gmp-5.0 :sb-gmp-5.1))) 978 (when persistently 979 (pushnew 'load-gmp sb-ext:*init-hooks*) 980 (pushnew 'uninstall-gmp-funs sb-ext:*save-hooks*)) 981 (let ((success (%load-gmp))) 982 (when success 983 (setf *gmp-version* (extern-alien "__gmp_version" c-string))) 984 (cond ((null *gmp-version*)) 985 ((string<= *gmp-version* "5.") 986 (warn "SB-GMP requires at least GMP version 5.0") 987 (setf success nil)) 988 (t 989 (pushnew :sb-gmp *gmp-features*) 990 (pushnew :sb-gmp-5.0 *gmp-features*) 991 (when (string>= *gmp-version* "5.1") 992 (pushnew :sb-gmp-5.1 *gmp-features*)) 993 (setf *features* (union *features* *gmp-features*)))) 994 (if success 995 (install-gmp-funs) 996 (uninstall-gmp-funs)) 997 (setup-5.1-stubs) 998 success)) 999 1000(defun unload-gmp () 1001 (setf sb-ext:*init-hooks* (remove 'load-gmp sb-ext:*init-hooks*)) 1002 (uninstall-gmp-funs) 1003 (setf sb-ext:*save-hooks* (remove 'uninstall-gmp-funs sb-ext:*save-hooks*)) 1004 (values)) 1005 1006(load-gmp) 1007