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-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 \f
  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
 137 bignum."
 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
 149 pre-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 \f
 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 \f
 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 \f
 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 \f
 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 \f
 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 \f
 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 \f
 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 \f
 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 \f
 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)