contrib/sb-gmp/gmp.lisp

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