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   (declare (optimize (speed 3) (space 3) (safety 0))
 453            (type bignum-type base))
 454   (check-type exp (unsigned-byte #.sb-vm:n-word-bits))
 455   (with-gmp-mpz-results (rop)
 456     (with-mpz-vars ((base gbase))
 457       (__gmpz_pow_ui (addr rop) (addr gbase) exp))))
 458
 459 (defgmpfun mpz-powm (base exp mod)
 460   (with-mpz-results ((rop (1+ (blength mod))))
 461     (with-mpz-vars ((base gbase) (exp gexp) (mod gmod))
 462       (__gmpz_powm (addr rop) (addr gbase) (addr gexp) (addr gmod)))))
 463
 464 (defgmpfun mpz-gcd (a b)
 465   (with-mpz-results ((result (min (blength a)
 466                                   (blength b))))
 467     (with-mpz-vars ((a ga) (b gb))
 468       (__gmpz_gcd (addr result) (addr ga) (addr gb)))))
 469
 470 (defgmpfun mpz-lcm (a b)
 471   (with-mpz-results ((result (+ (blength a)
 472                                 (blength b))))
 473     (with-mpz-vars ((a ga) (b gb))
 474       (__gmpz_lcm (addr result) (addr ga) (addr gb)))))
 475
 476 (defgmpfun mpz-sqrt (a)
 477   (with-mpz-results ((result (1+ (ceiling (blength a) 2))))
 478     (with-mpz-vars ((a ga))
 479       (__gmpz_sqrt (addr result) (addr ga)))))
 480
 481 \f
 482 ;;; Functions that use GMP-side allocated integers and copy the result
 483 ;;; into a SBCL bignum at the end of the computation when the required
 484 ;;; bignum length is known.
 485 (defun mpz-probably-prime-p (n &optional (reps 25))
 486   (declare (optimize (speed 3) (space 3) (safety 0)))
 487   (check-type reps fixnum)
 488   (with-mpz-vars ((n gn))
 489     (__gmpz_probab_prime_p (addr gn) reps)))
 490
 491 (defun mpz-nextprime (a)
 492   (declare (optimize (speed 3) (space 3) (safety 0)))
 493   (with-gmp-mpz-results (prime)
 494     (with-mpz-vars ((a ga))
 495       (__gmpz_nextprime (addr prime) (addr ga)))))
 496
 497 (defun mpz-fac (n)
 498   (declare (optimize (speed 3) (space 3) (safety 0)))
 499   (check-type n (unsigned-byte #.sb-vm:n-word-bits))
 500   (with-gmp-mpz-results (fac)
 501     (__gmpz_fac_ui (addr fac) n)))
 502
 503 (defun %mpz-2fac (n)
 504   (declare (optimize (speed 3) (space 3) (safety 0)))
 505   (check-type n (unsigned-byte #.sb-vm:n-word-bits))
 506   (with-gmp-mpz-results (fac)
 507     (__gmpz_2fac_ui (addr fac) n)))
 508
 509 (defun %mpz-mfac (n m)
 510   (declare (optimize (speed 3) (space 3) (safety 0)))
 511   (check-type n (unsigned-byte #.sb-vm:n-word-bits))
 512   (check-type m (unsigned-byte #.sb-vm:n-word-bits))
 513   (with-gmp-mpz-results (fac)
 514     (__gmpz_mfac_uiui (addr fac) n m)))
 515
 516 (defun %mpz-primorial (n)
 517   (declare (optimize (speed 3) (space 3) (safety 0)))
 518   (check-type n (unsigned-byte #.sb-vm:n-word-bits))
 519   (with-gmp-mpz-results (r)
 520     (__gmpz_primorial_ui (addr r) n)))
 521
 522 (defun setup-5.1-stubs ()
 523   (macrolet ((stubify (name implementation &rest arguments)
 524                `(setf (fdefinition ',name)
 525                       (if (member :sb-gmp-5.1 *gmp-features*)
 526                           (fdefinition ',implementation)
 527                           (lambda ,arguments
 528                             (declare (ignore ,@arguments))
 529                             (error "~S is only available in GMP >= 5.1"
 530                                    ',name))))))
 531     (stubify mpz-2fac %mpz-2fac n)
 532     (stubify mpz-mfac %mpz-mfac n m)
 533     (stubify mpz-primorial %mpz-primorial n)))
 534
 535 (defun mpz-bin (n k)
 536   (declare (optimize (speed 3) (space 3) (safety 0)))
 537   (check-type k (unsigned-byte #.sb-vm:n-word-bits))
 538   (with-gmp-mpz-results (r)
 539     (with-mpz-vars ((n gn))
 540       (__gmpz_bin_ui (addr r) (addr gn) k))))
 541
 542 (defun mpz-fib2 (n)
 543   (declare (optimize (speed 3) (space 3) (safety 0)))
 544   ;; (let ((size (1+ (ceiling (* n (log 1.618034 2)) 64)))))
 545   ;; fibonacci number magnitude in bits is assymptotic to n(log_2 phi)
 546   ;; This is correct for the result but appears not to be enough for GMP
 547   ;; during computation (memory access error), so use GMP-side allocation.
 548   (check-type n (unsigned-byte #.sb-vm:n-word-bits))
 549   (with-gmp-mpz-results (fibn fibn-1)
 550     (__gmpz_fib2_ui (addr fibn) (addr fibn-1) n)))
 551
 552 \f
 553 ;;;; Random bignum (mpz) generation
 554
 555 ;; we do not actually use the gestalt of the struct but need its size
 556 ;; for allocation purposes
 557 (define-alien-type nil
 558     (struct gmprandstate
 559             (mp_seed (struct gmpint))
 560             (mp_alg int)
 561             (mp_algdata (* t))))
 562
 563 (declaim (inline __gmp_randinit_mt
 564                  __gmp_randinit_lc_2exp
 565                  __gmp_randseed
 566                  __gmp_randseed_ui
 567                  __gmpz_urandomb
 568                  __gmpz_urandomm))
 569
 570 (define-alien-routine __gmp_randinit_mt void
 571   (s (* (struct gmprandstate))))
 572
 573 (define-alien-routine __gmp_randinit_lc_2exp void
 574   (s (* (struct gmprandstate)))
 575   (a (* (struct gmpint)))
 576   (c unsigned-long)
 577   (exp unsigned-long))
 578
 579 (define-alien-routine __gmp_randseed void
 580   (s (* (struct gmprandstate)))
 581   (sd (* (struct gmpint))))
 582
 583 (define-alien-routine __gmp_randseed_ui void
 584   (s (* (struct gmprandstate)))
 585   (sd unsigned-long))
 586
 587 (define-alien-routine __gmpz_urandomb void
 588   (r (* (struct gmpint)))
 589   (s (* (struct gmprandstate)))
 590   (bcnt unsigned-long))
 591
 592 (define-alien-routine __gmpz_urandomm void
 593   (r (* (struct gmpint)))
 594   (s (* (struct gmprandstate)))
 595   (n (* (struct gmpint))))
 596
 597 (defstruct (gmp-rstate (:constructor %make-gmp-rstate))
 598   (ref (make-alien (struct gmprandstate))
 599    :type (alien (* (struct gmprandstate))) :read-only t))
 600
 601 (defun make-gmp-rstate ()
 602   "Instantiate a state for the GMP Mersenne-Twister random number generator."
 603   (declare (optimize (speed 3) (space 3)))
 604   (let* ((state (%make-gmp-rstate))
 605          (ref (gmp-rstate-ref state)))
 606     (__gmp_randinit_mt ref)
 607     (sb-ext:finalize state (lambda () (free-alien ref)))
 608     state))
 609
 610 (defun make-gmp-rstate-lc (a c m2exp)
 611   "Instantiate a state for the GMP linear congruential random number generator."
 612   (declare (optimize (speed 3) (space 3) (safety 0)))
 613   (check-type c (unsigned-byte #.sb-vm:n-word-bits))
 614   (check-type m2exp (unsigned-byte #.sb-vm:n-word-bits))
 615   (let* ((state (%make-gmp-rstate))
 616          (ref (gmp-rstate-ref state)))
 617     (with-mpz-vars ((a ga))
 618       (__gmp_randinit_lc_2exp ref (addr ga) c m2exp))
 619     (sb-ext:finalize state (lambda () (free-alien ref)))
 620     state))
 621
 622 (defun rand-seed (state seed)
 623   "Initialize a random STATE with SEED."
 624   (declare (optimize (speed 3) (space 3) (safety 0)))
 625   (check-type state gmp-rstate)
 626   (let ((ref (gmp-rstate-ref state)))
 627     (cond
 628       ((typep seed '(unsigned-byte #.sb-vm:n-word-bits))
 629        (__gmp_randseed_ui ref seed))
 630       ((typep seed '(integer 0 *))
 631        (with-mpz-vars ((seed gseed))
 632          (__gmp_randseed ref (addr gseed))))
 633       (t
 634        (error "SEED must be a positive integer")))))
 635
 636 (defun random-bitcount (state bitcount)
 637   "Return a random integer in the range 0..(2^bitcount - 1)."
 638   (declare (optimize (speed 3) (space 3) (safety 0)))
 639   (check-type state gmp-rstate)
 640   (check-type bitcount (unsigned-byte #.sb-vm:n-word-bits))
 641   (let ((ref (gmp-rstate-ref state)))
 642     (with-mpz-results ((result (+ (ceiling bitcount sb-vm:n-word-bits) 2)))
 643       (__gmpz_urandomb (addr result) ref bitcount))))
 644
 645 (defun random-int (state boundary)
 646   "Return a random integer in the range 0..(boundary - 1)."
 647   (declare (optimize (speed 3) (space 3) (safety 0)))
 648   (check-type state gmp-rstate)
 649   (let ((b (bassert boundary))
 650         (ref (gmp-rstate-ref state)))
 651     (with-mpz-results ((result (1+ (%bignum-length b))))
 652       (with-mpz-vars ((b gb))
 653         (__gmpz_urandomm (addr result) ref (addr gb))))))
 654
 655 \f
 656 ;;; Rational functions
 657 (declaim (inline %lsize))
 658 (defun %lsize (minusp n)
 659   (declare (optimize (speed 3) (space 3) (safety 0)))
 660   "n must be a (potentially denormalized) bignum"
 661   (let ((length (%bignum-length n)))
 662     (when (zerop (%bignum-ref n (1- length)))
 663       (decf length))
 664     (if minusp (- length) length)))
 665
 666 (defmacro defmpqfun (name gmpfun)
 667   `(progn
 668      (declaim (sb-ext:maybe-inline ,name))
 669      (defun ,name (a b)
 670        (declare (optimize (speed 3) (space 3) (safety 0)))
 671        (let ((size (+ (max (blength (numerator a))
 672                            (blength (denominator a)))
 673                       (max (blength (numerator b))
 674                            (blength (denominator b)))
 675                       3)))
 676          (with-alien ((r (struct gmprat)))
 677            (let ((num (%allocate-bignum size))
 678                  (den (%allocate-bignum size)))
 679              (sb-sys:with-pinned-objects (num den)
 680                (setf (slot (slot r 'mp_num) 'mp_size) 0
 681                      (slot (slot r 'mp_num) 'mp_alloc) size
 682                      (slot (slot r 'mp_num) 'mp_d) (bignum-data-sap num))
 683                (setf (slot (slot r 'mp_den) 'mp_size) 0
 684                      (slot (slot r 'mp_den) 'mp_alloc) size
 685                      (slot (slot r 'mp_den) 'mp_d) (bignum-data-sap den))
 686                (let* ((an (bassert (numerator a)))
 687                       (ad (bassert (denominator a)))
 688                       (asign (not (%bignum-0-or-plusp an (%bignum-length an))))
 689                       (bn (bassert (numerator b)))
 690                       (bd (bassert (denominator b)))
 691                       (bsign (not (%bignum-0-or-plusp bn (%bignum-length bn)))))
 692                  (when asign
 693                    (setf an (negate-bignum an nil)))
 694                  (when bsign
 695                    (setf bn (negate-bignum bn nil)))
 696                  (let* ((anlen (%lsize asign an))
 697                         (adlen (%lsize NIL ad))
 698                         (bnlen (%lsize bsign bn))
 699                         (bdlen (%lsize NIL bd)))
 700                    (with-alien ((arga (struct gmprat))
 701                                 (argb (struct gmprat)))
 702                      (sb-sys:with-pinned-objects (an ad bn bd)
 703                        (setf (slot (slot arga 'mp_num) 'mp_size) anlen
 704                              (slot (slot arga 'mp_num) 'mp_alloc) (abs anlen)
 705                              (slot (slot arga 'mp_num) 'mp_d)
 706                              (bignum-data-sap an))
 707                        (setf (slot (slot arga 'mp_den) 'mp_size) adlen
 708                              (slot (slot arga 'mp_den) 'mp_alloc) (abs adlen)
 709                              (slot (slot arga 'mp_den) 'mp_d)
 710                              (bignum-data-sap ad))
 711                        (setf (slot (slot argb 'mp_num) 'mp_size) bnlen
 712                              (slot (slot argb 'mp_num) 'mp_alloc) (abs bnlen)
 713                              (slot (slot argb 'mp_num) 'mp_d)
 714                              (bignum-data-sap bn))
 715                        (setf (slot (slot argb 'mp_den) 'mp_size) bdlen
 716                              (slot (slot argb 'mp_den) 'mp_alloc) (abs bdlen)
 717                              (slot (slot argb 'mp_den) 'mp_d)
 718                              (bignum-data-sap bd))
 719                        (,gmpfun (addr r) (addr arga) (addr argb)))))
 720                  (locally (declare (optimize (speed 1)))
 721                    (sb-kernel::build-ratio (if (minusp (slot (slot r 'mp_num) 'mp_size))
 722                                                (z-to-bignum-neg num size)
 723                                                (z-to-bignum num size))
 724                                            (z-to-bignum den size)))))))))))
 725
 726 (defmpqfun mpq-add __gmpq_add)
 727 (defmpqfun mpq-sub __gmpq_sub)
 728 (defmpqfun mpq-mul __gmpq_mul)
 729 (defmpqfun mpq-div __gmpq_div)
 730
 731 \f
 732 ;;;; SBCL interface and integration installation
 733 (macrolet ((def (name original)
 734              (let ((special (intern (format nil "*~A-FUNCTION*" name))))
 735                `(progn
 736                   (declaim (type function ,special)
 737                            (inline ,name))
 738                   (defvar ,special (symbol-function ',original))
 739                   (defun ,name (&rest args)
 740                     (apply (load-time-value ,special t) args))))))
 741   (def orig-mul multiply-bignums)
 742   (def orig-truncate bignum-truncate)
 743   (def orig-gcd bignum-gcd)
 744   (def orig-lcm sb-kernel:two-arg-lcm)
 745   (def orig-isqrt isqrt)
 746   (def orig-two-arg-+ sb-kernel:two-arg-+)
 747   (def orig-two-arg-- sb-kernel:two-arg--)
 748   (def orig-two-arg-* sb-kernel:two-arg-*)
 749   (def orig-two-arg-/ sb-kernel:two-arg-/))
 750
 751 ;;; integers
 752 (defun gmp-mul (a b)
 753   (declare (optimize (speed 3) (space 3))
 754            (type bignum-type a b)
 755            (inline mpz-mul))
 756   (if (or (< (min (%bignum-length a)
 757                   (%bignum-length b))
 758              6)
 759           *gmp-disabled*)
 760       (orig-mul a b)
 761       (mpz-mul a b)))
 762
 763 (defun gmp-truncate (a b)
 764   (declare (optimize (speed 3) (space 3))
 765            (type bignum-type a b)
 766            (inline mpz-tdiv))
 767   (if (or (< (min (%bignum-length a)
 768                   (%bignum-length b))
 769              3)
 770           *gmp-disabled*)
 771       (orig-truncate a b)
 772       (mpz-tdiv a b)))
 773
 774 (defun gmp-lcm (a b)
 775   (declare (optimize (speed 3) (space 3))
 776            (type integer a b)
 777            (inline mpz-lcm))
 778   (if (or (and (typep a 'fixnum)
 779                (typep b 'fixnum))
 780           *gmp-disabled*)
 781       (orig-lcm a b)
 782       (mpz-lcm a b)))
 783
 784 (defun gmp-isqrt (n)
 785   (declare (optimize (speed 3) (space 3))
 786            (type unsigned-byte n)
 787            (inline mpz-sqrt))
 788   (if (or (typep n 'fixnum)
 789           *gmp-disabled*)
 790       (orig-isqrt n)
 791       (mpz-sqrt n)))
 792
 793 ;;; rationals
 794 (defun gmp-two-arg-+ (x y)
 795   (declare (optimize (speed 3) (space 3))
 796            (inline mpq-add))
 797   (if (and (or (typep x 'ratio)
 798                (typep y 'ratio))
 799            (rationalp y)
 800            (rationalp x)
 801            (not *gmp-disabled*))
 802       (mpq-add x y)
 803       (orig-two-arg-+ x y)))
 804
 805 (defun gmp-two-arg-- (x y)
 806   (declare (optimize (speed 3) (space 3))
 807            (inline mpq-sub))
 808   (if (and (or (typep x 'ratio)
 809                (typep y 'ratio))
 810            (rationalp y)
 811            (rationalp x)
 812            (not *gmp-disabled*))
 813       (mpq-sub x y)
 814       (orig-two-arg-- x y)))
 815
 816 (defun gmp-two-arg-* (x y)
 817   (declare (optimize (speed 3) (space 3))
 818            (inline mpq-mul))
 819   (if (and (or (typep x 'ratio)
 820                (typep y 'ratio))
 821            (rationalp y)
 822            (rationalp x)
 823            (not *gmp-disabled*))
 824       (mpq-mul x y)
 825       (orig-two-arg-* x y)))
 826
 827 (defun gmp-two-arg-/ (x y)
 828   (declare (optimize (speed 3) (space 3))
 829            (inline mpq-div))
 830   (if (and (rationalp x)
 831            (rationalp y)
 832            (not (eql y 0))
 833            (not *gmp-disabled*))
 834       (mpq-div x y)
 835       (orig-two-arg-/ x y)))
 836
 837 ;;; installation
 838 (defmacro with-package-locks-ignored (&body body)
 839   `(handler-bind ((sb-ext:package-lock-violation
 840                     (lambda (condition)
 841                       (declare (ignore condition))
 842                       (invoke-restart :ignore-all))))
 843      ,@body))
 844
 845 (defun install-gmp-funs ()
 846   (with-package-locks-ignored
 847       (macrolet ((def (destination source)
 848                    `(setf (fdefinition ',destination)
 849                           (fdefinition ',source))))
 850         (def multiply-bignums gmp-mul)
 851         (def bignum-truncate gmp-truncate)
 852         (def bignum-gcd mpz-gcd)
 853         (def sb-kernel:two-arg-lcm gmp-lcm)
 854         (def sb-kernel:two-arg-+ gmp-two-arg-+)
 855         (def sb-kernel:two-arg-- gmp-two-arg--)
 856         (def sb-kernel:two-arg-* gmp-two-arg-*)
 857         (def sb-kernel:two-arg-/ gmp-two-arg-/)
 858         (def isqrt gmp-isqrt)))
 859   (values))
 860
 861 (defun uninstall-gmp-funs ()
 862   (with-package-locks-ignored
 863       (macrolet ((def (destination source)
 864                    `(setf (fdefinition ',destination)
 865                           ,(intern (format nil "*~A-FUNCTION*" source)))))
 866         (def multiply-bignums orig-mul)
 867         (def bignum-truncate orig-truncate)
 868         (def bignum-gcd orig-gcd)
 869         (def sb-kernel:two-arg-lcm orig-lcm)
 870         (def sb-kernel:two-arg-+ orig-two-arg-+)
 871         (def sb-kernel:two-arg-- orig-two-arg--)
 872         (def sb-kernel:two-arg-* orig-two-arg-*)
 873         (def sb-kernel:two-arg-/ orig-two-arg-/)
 874         (def isqrt orig-isqrt)))
 875   (values))
 876
 877 (defun load-gmp (&key (persistently t))
 878   (setf *gmp-features* nil
 879         *gmp-version* nil
 880         *features* (set-difference *features* '(:sb-gmp :sb-gmp-5.0 :sb-gmp-5.1)))
 881   (when persistently
 882     (pushnew 'load-gmp sb-ext:*init-hooks*)
 883     (pushnew 'uninstall-gmp-funs sb-ext:*save-hooks*))
 884   (let ((success (%load-gmp)))
 885     (when success
 886       (setf *gmp-version* (extern-alien "__gmp_version" c-string)))
 887     (cond ((null *gmp-version*))
 888           ((string<= *gmp-version* "5.")
 889            (warn "SB-GMP requires at least GMP version 5.0")
 890            (setf success nil))
 891           (t
 892            (pushnew :sb-gmp *gmp-features*)
 893            (pushnew :sb-gmp-5.0 *gmp-features*)
 894            (when (string>= *gmp-version* "5.1")
 895              (pushnew :sb-gmp-5.1 *gmp-features*))
 896            (setf *features* (union *features* *gmp-features*))))
 897     (if success
 898         (install-gmp-funs)
 899         (uninstall-gmp-funs))
 900     (setup-5.1-stubs)
 901     success))
 902
 903 (defun unload-gmp ()
 904   (setf sb-ext:*init-hooks* (remove 'load-gmp sb-ext:*init-hooks*))
 905   (uninstall-gmp-funs)
 906   (setf sb-ext:*save-hooks* (remove 'uninstall-gmp-funs sb-ext:*save-hooks*))
 907   (values))
 908
 909 (load-gmp)