New contrib: SB-GMP
authorPaul Khuong <pvk@pvk.ca>
Fri, 28 Jun 2013 03:15:38 +0000 (23:15 -0400)
committerPaul Khuong <pvk@pvk.ca>
Fri, 28 Jun 2013 04:18:16 +0000 (00:18 -0400)
 This contrib was developed by Stephan Frank to replace some of our
 bignum and rational arithmetic code with calls to libgmp.  Simply
 loading the contrib will transparently accelerate arithmetic on
 large rationals when libgmp is available; if libgmp cannot be found,
 the contrib should change nothing.

 The contrib also wraps additional functions in GNU MP, so that they
 accept and return SBCL-native integers or ratios.  See GNU MP's manual
 for more information.

CREDITS
NEWS
contrib/sb-gmp/Makefile [new file with mode: 0644]
contrib/sb-gmp/README.md [new file with mode: 0644]
contrib/sb-gmp/gmp.lisp [new file with mode: 0644]
contrib/sb-gmp/sb-gmp.asd [new file with mode: 0644]
contrib/sb-gmp/tests.lisp [new file with mode: 0644]

diff --git a/CREDITS b/CREDITS
index 99a1b12..b2061fb 100644 (file)
--- a/CREDITS
+++ b/CREDITS
@@ -588,6 +588,10 @@ Andreas Fuchs:
   regressions of SBCL.  He assisted with the integration of the
   Unicode work.
 
+Stephan Frank:
+  He contributed the SB-GMP contrib to exploit libgmp in bignum and
+  ratio arithmetic.
+
 Nathan Froyd:
   He has fixed various bugs, and also done a lot of internal
   cleanup, not visible at the user level but important for
@@ -856,6 +860,7 @@ NS   Nikodemus Siivola
 PFD  Paul F. Dietz
 PRM  Pierre Mai
 PVE  Peter Van Eynde
+PK/PVK Paul-Virak Khuong
 PW   Paul Werkowski
 RAM  Robert MacLachlan
 RLT  Raymond Toy
diff --git a/NEWS b/NEWS
index 9f0015e..d5644be 100644 (file)
--- a/NEWS
+++ b/NEWS
@@ -1,5 +1,7 @@
 ;;;; -*- coding: utf-8; fill-column: 78 -*-
 changes relative to sbcl-1.1.8:
+  * new feature: the contrib SB-GMP links with libgmp at runtime to speed
+    up arithmetic on bignums and ratios.  (contributed by Stephan Frank)
   * enhancement: disassemble now annotates some previously missing static
     functions, like LENGTH.
   * enhancement: clean.sh now also cleans doc/internals.
diff --git a/contrib/sb-gmp/Makefile b/contrib/sb-gmp/Makefile
new file mode 100644 (file)
index 0000000..1b86400
--- /dev/null
@@ -0,0 +1,2 @@
+SYSTEM=sb-gmp
+include ../asdf-module.mk
diff --git a/contrib/sb-gmp/README.md b/contrib/sb-gmp/README.md
new file mode 100644 (file)
index 0000000..56ca695
--- /dev/null
@@ -0,0 +1,39 @@
+SB-GMP
+======
+
+GMP integration for SBCL
+
+This contrib enables the use of bignum computation routines from the
+GMP library for SBCL internal bignum calculations.
+
+Some functions can be transparently replaced within SBCL, namely:
+
+ - sb-bignum:multiply-bignum
+ - sb-bignum:bignum-truncate
+ - sb-bignum:bignum-gcd
+ - sb-kernel::two-arg-lcm
+ - cl:isqrt
+
+and for making use of the GMP rational arithmetic:
+
+ - sb-kernel::two-arg-+
+ - sb-kernel::two-arg--
+ - sb-kernel::two-arg-*
+ - sb-kernel::two-arg-/
+
+Most of the other SBCL bignum routines rely on these functions for the
+heavy computational lifting.
+
+However, SB-GMP also provides easy and transparent access to several
+other functions of the GMP library and may be able to replace other
+parts of the SBCL bignum machinery in the future. Refer to the GMP
+documentation for the respective function specification.
+
+The transparent SBCL integration is automatically activated by loading
+SB-GMP.  Some care is taken to behave correctly across core saves and
+loads; however, if this does not work correctly, `SB-GMP:UNLOAD-GMP`
+and `SB-GMP:LOAD-GMP` can be called to explicitly unload and reload
+libgmp and uninstall/reinstall the hooks.  The hooks can also be
+directly uninstalled and installed with `SB-GMP:UNINSTALL-GMP-FUNS`
+and `SB-GMP:INSTALL-GMP-FUNS`.  However, in most cases, binding
+`SB-GMP:*GMP-DISABLED*` to `T` should suffice.
diff --git a/contrib/sb-gmp/gmp.lisp b/contrib/sb-gmp/gmp.lisp
new file mode 100644 (file)
index 0000000..3dfefd6
--- /dev/null
@@ -0,0 +1,909 @@
+(defpackage "SB-GMP"
+  (:use "COMMON-LISP" "SB-ALIEN" "SB-BIGNUM")
+  ;; we need a few very internal symbols
+  (:import-from "SB-BIGNUM"
+                "%BIGNUM-0-OR-PLUSP" "%NORMALIZE-BIGNUM"
+                "NEGATE-BIGNUM-IN-PLACE")
+  (:export
+   ;; bignum integer operations
+   #:mpz-add
+   #:mpz-sub
+   #:mpz-mul
+   #:mpz-mod
+   #:mpz-cdiv
+   #:mpz-fdiv
+   #:mpz-tdiv
+   #:mpz-powm
+   #:mpz-pow
+   #:mpz-gcd
+   #:mpz-lcm
+   #:mpz-sqrt
+   #:mpz-probably-prime-p
+   #:mpz-nextprime
+   #:mpz-fac
+   ;; Following three are GMP >= 5.1 only
+   #:mpz-2fac
+   #:mpz-mfac
+   #:mpz-primorial
+   #:mpz-bin
+   #:mpz-fib2
+   ;; random number generation
+   #:make-gmp-rstate
+   #:make-gmp-rstate-lc
+   #:rand-seed
+   #:random-bitcount
+   #:random-int
+   ;; ratio arithmetic
+   #:mpq-add
+   #:mpq-sub
+   #:mpq-mul
+   #:mpq-div
+   ;; (un)installer functions
+   ; these insert/remove the runtime patch in SBCL's runtime
+   #:install-gmp-funs
+   #:uninstall-gmp-funs
+   ; these also load/unload the shared library and setup/clear
+   ; hooks to handle core saves
+   #:load-gmp
+   #:unload-gmp
+   ;; special variables
+   #:*gmp-version*
+   #:*gmp-disabled*
+   ))
+
+(in-package "SB-GMP")
+
+(defvar *gmp-disabled* nil)
+
+(defconstant +bignum-raw-area-offset+
+  (- (* sb-vm:bignum-digits-offset sb-vm:n-word-bytes)
+     sb-vm:other-pointer-lowtag))
+
+(declaim (inline bignum-data-sap))
+(defun bignum-data-sap (x)
+  (declare (type bignum x))
+  (sb-sys:sap+ (sb-sys:int-sap (sb-kernel:get-lisp-obj-address x))
+               +bignum-raw-area-offset+))
+
+(defun %load-gmp ()
+  (handler-case
+      (load-shared-object #-(or win32 darwin) "libgmp.so"
+                          #+darwin "libgmp.dylib"
+                          #+win32 "libgmp-10.dll"
+                          :dont-save t)
+    (error (e)
+      (warn "GMP not loaded (~a)" e)
+      (return-from %load-gmp nil)))
+  t)
+
+(defvar *gmp-features* nil)
+(defvar *gmp-version* nil)
+
+;; We load only the library right now to avoid undefined alien
+;; style warnings
+(%load-gmp)
+
+\f
+;;; types and initialization
+(define-alien-type nil
+    (struct gmpint
+            (mp_alloc int)
+            (mp_size int)
+            (mp_d (* unsigned-long))))
+
+;; Section 3.6 "Memory Management" of the GMP manual states: "mpz_t
+;; and mpq_t variables never reduce their allocated space. Normally
+;; this is the best policy, since it avoids frequent
+;; reallocation. Applications that need to return memory to the heap
+;; at some particular point can use mpz_realloc2, or clear variables
+;; no longer needed."
+;;
+;; We can therefore allocate a bignum of sufficiant size and use the
+;; space for GMP computations without the need for memory transfer
+;; from C to Lisp space.
+(declaim (inline z-to-bignum z-to-bignum-neg))
+
+(defun z-to-bignum (b count)
+  "Convert GMP integer in the buffer of a pre-allocated bignum."
+  (declare (optimize (speed 3) (space 3) (safety 0))
+           (type bignum-type b)
+           (type bignum-index count))
+  (if (zerop count)
+      0
+      (the unsigned-byte (%normalize-bignum b count))))
+
+(defun z-to-bignum-neg (b count)
+  "Convert to twos complement int the buffer of a pre-allocated
+bignum."
+  (declare (optimize (speed 3) (space 3) (safety 0))
+           (type bignum-type b)
+           (type bignum-index count))
+  (negate-bignum-in-place b)
+  (the (integer * 0) (%normalize-bignum b count)))
+
+;;; conversion functions that also copy from GMP to SBCL bignum space
+(declaim (inline gmp-z-to-bignum gmp-z-to-bignum-neg))
+
+(defun gmp-z-to-bignum (z b count)
+  "Convert and copy a positive GMP integer into the buffer of a
+pre-allocated bignum. The allocated bignum-length must be (1+ COUNT)."
+  (declare (optimize (speed 3) (space 3) (safety 0))
+           (type (alien (* unsigned-long)) z)
+           (type bignum-type b)
+           (type bignum-index count))
+  (dotimes (i count (%normalize-bignum b count))
+    (%bignum-set b i (deref z i))))
+
+(defun gmp-z-to-bignum-neg (z b count)
+  "Convert to twos complement and copy a negative GMP integer into the
+buffer of a pre-allocated bignum. The allocated bignum-length must
+be (1+ COUNT)."
+  (declare (optimize (speed 3) (space 3) (safety 0))
+           (type (alien (* unsigned-long)) z)
+           (type bignum-type b)
+           (type bignum-index count))
+  (let ((carry 0)
+        (add 1))
+    (declare (type (mod 2) carry add))
+    (dotimes (i count b)
+      (multiple-value-bind (value carry-tmp)
+          (%add-with-carry
+           (%lognot (deref z i)) add carry)
+        (%bignum-set b i value)
+        (setf carry carry-tmp
+              add 0)))))
+
+(declaim (inline blength bassert)
+         (ftype (function (integer) (values bignum-index &optional)) blength)
+         (ftype (function (integer) (values bignum &optional)) bassert))
+
+(defun blength (a)
+  (declare (optimize (speed 3) (space 3) (safety 0)))
+  (etypecase a
+    (fixnum 1)
+    (t (%bignum-length a))))
+
+(defun bassert (a)
+  (declare (optimize (speed 3) (space 3) (safety 0)))
+  (etypecase a
+    (fixnum (make-small-bignum a))
+    (t a)))
+
+;;;; rationals
+(define-alien-type nil
+    (struct gmprat
+            (mp_num (struct gmpint))
+            (mp_den (struct gmpint))))
+
+;;; memory initialization functions to support non-alloced results
+;;; since an upper bound cannot always correctly predetermined
+;;; (e.g. the memory required for the fib function exceed the number
+;;; of limbs that are be determined through the infamous Phi-relation
+;;; resulting in a memory access error.
+
+;; use these for non-prealloced bignum values, but only when
+;; ultimately necessary since copying back into bignum space a the end
+;; of the operation is about three times slower than the shared buffer
+;; approach.
+(declaim (inline __gmpz_init __gmpz_clear))
+(define-alien-routine __gmpz_init void
+  (x (* (struct gmpint))))
+
+(define-alien-routine __gmpz_clear void
+  (x (* (struct gmpint))))
+
+\f
+;;; integer interface functions
+(defmacro define-twoarg-mpz-funs (funs)
+  (loop for i in funs collect `(define-alien-routine ,i void
+                                 (r (* (struct gmpint)))
+                                 (a (* (struct gmpint))))
+          into defines
+        finally (return `(progn
+                           (declaim (inline ,@funs))
+                           ,@defines))))
+
+(defmacro define-threearg-mpz-funs (funs)
+  (loop for i in funs collect `(define-alien-routine ,i void
+                                 (r (* (struct gmpint)))
+                                 (a (* (struct gmpint)))
+                                 (b (* (struct gmpint))))
+          into defines
+        finally (return `(progn
+                           (declaim (inline ,@funs))
+                           ,@defines))))
+
+(defmacro define-fourarg-mpz-funs (funs)
+  (loop for i in funs collect `(define-alien-routine ,i void
+                                 (r (* (struct gmpint)))
+                                 (a (* (struct gmpint)))
+                                 (b (* (struct gmpint)))
+                                 (c (* (struct gmpint))))
+          into defines
+        finally (return `(progn
+                           (declaim (inline ,@funs))
+                           ,@defines))))
+
+\f
+(define-twoarg-mpz-funs (__gmpz_sqrt
+                         __gmpz_nextprime))
+
+(define-threearg-mpz-funs (__gmpz_add
+                           __gmpz_sub
+                           __gmpz_mul
+                           __gmpz_mod
+                           __gmpz_gcd
+                           __gmpz_lcm))
+
+(define-fourarg-mpz-funs (__gmpz_cdiv_qr
+                          __gmpz_fdiv_qr
+                          __gmpz_tdiv_qr
+                          __gmpz_powm))
+
+(declaim (inline __gmpz_pow_ui
+                 __gmpz_probab_prime_p
+                 __gmpz_fac_ui
+                 __gmpz_2fac_ui
+                 __gmpz_mfac_uiui
+                 __gmpz_primorial_ui
+                 __gmpz_bin_ui
+                 __gmpz_fib2_ui))
+
+(define-alien-routine __gmpz_pow_ui void
+  (r (* (struct gmpint)))
+  (b (* (struct gmpint)))
+  (e unsigned-long))
+
+(define-alien-routine __gmpz_probab_prime_p int
+  (n (* (struct gmpint)))
+  (reps int))
+
+(define-alien-routine __gmpz_fac_ui void
+  (r (* (struct gmpint)))
+  (a unsigned-long))
+
+(define-alien-routine __gmpz_2fac_ui void
+  (r (* (struct gmpint)))
+  (a unsigned-long))
+
+(define-alien-routine __gmpz_mfac_uiui void
+  (r (* (struct gmpint)))
+  (n unsigned-long)
+  (m unsigned-long))
+
+(define-alien-routine __gmpz_primorial_ui void
+  (r (* (struct gmpint)))
+  (n unsigned-long))
+
+(define-alien-routine __gmpz_bin_ui void
+  (r (* (struct gmpint)))
+  (n (* (struct gmpint)))
+  (k unsigned-long))
+
+(define-alien-routine __gmpz_fib2_ui void
+  (r (* (struct gmpint)))
+  (a (* (struct gmpint)))
+  (b unsigned-long))
+
+\f
+;; ratio functions
+(defmacro define-threearg-mpq-funs (funs)
+  (loop for i in funs collect `(define-alien-routine ,i void
+                                 (r (* (struct gmprat)))
+                                 (a (* (struct gmprat)))
+                                 (b (* (struct gmprat))))
+          into defines
+        finally (return `(progn
+                           (declaim (inline ,@funs))
+                           ,@defines))))
+
+(define-threearg-mpq-funs (__gmpq_add
+                           __gmpq_sub
+                           __gmpq_mul
+                           __gmpq_div))
+
+\f
+;;;; SBCL interface
+
+;;; utility macros for GMP mpz variable and result declaration and
+;;; incarnation of associated SBCL bignums
+(defmacro with-mpz-results (pairs &body body)
+  (loop for (gres size) in pairs
+        for res = (gensym "RESULT")
+        collect `(,gres (struct gmpint)) into declares
+        collect `(,res (%allocate-bignum ,size))
+          into resinits
+        collect `(setf (slot ,gres 'mp_alloc) (%bignum-length ,res)
+                       (slot ,gres 'mp_size) 0
+                       (slot ,gres 'mp_d) (bignum-data-sap ,res))
+          into inits
+        collect `(if (minusp (slot ,gres 'mp_size)) ; check for negative result
+                     (z-to-bignum-neg ,res ,size)
+                     (z-to-bignum ,res ,size))
+          into normlimbs
+        collect res into results
+        finally (return
+                  `(let ,resinits
+                     (sb-sys:with-pinned-objects ,results
+                       (with-alien ,declares
+                         ,@inits
+                         ,@body
+                         (values ,@normlimbs)))))))
+
+(defmacro with-mpz-vars (pairs &body body)
+  (loop for (a ga) in pairs
+        for length = (gensym "LENGTH")
+        for plusp = (gensym "PLUSP")
+        for barg = (gensym "BARG")
+        for arg = (gensym "ARG")
+        collect `(,ga (struct gmpint)) into declares
+        collect `(,barg (bassert ,a)) into gmpinits
+        collect `(,plusp (%bignum-0-or-plusp ,barg (%bignum-length ,barg))) into gmpinits
+        collect `(,arg (if ,plusp ,barg (negate-bignum ,barg nil))) into gmpinits
+        collect `(,length (%bignum-length ,arg)) into gmpinits
+        collect arg into vars
+        collect `(setf (slot ,ga 'mp_alloc) ,length
+                       (slot ,ga 'mp_size)
+                       (progn ;; handle twos complements/ulong limbs mismatch
+                         (when (zerop (%bignum-ref ,arg (1- ,length)))
+                           (decf ,length))
+                         (if ,plusp ,length (- ,length)))
+                       (slot ,ga 'mp_d) (bignum-data-sap ,arg))
+          into gmpvarssetup
+        finally (return
+                  `(with-alien ,declares
+                     (let* ,gmpinits
+                       (sb-sys:with-pinned-objects ,vars
+                         ,@gmpvarssetup
+                         ,@body))))))
+
+(defmacro with-gmp-mpz-results (resultvars &body body)
+  (loop for gres in resultvars
+        for res = (gensym "RESULT")
+        for size = (gensym "SIZE")
+        collect size into sizes
+        collect `(,gres (struct gmpint)) into declares
+        collect `(__gmpz_init (addr ,gres)) into inits
+        collect `(,size (1+ (abs (slot ,gres 'mp_size))))
+          into resinits
+        collect `(,res (%allocate-bignum ,size))
+          into resinits
+        collect `(setf ,res (if (minusp (slot ,gres 'mp_size)) ; check for negative result
+                                (gmp-z-to-bignum-neg (slot ,gres 'mp_d) ,res ,size)
+                                (gmp-z-to-bignum (slot ,gres 'mp_d) ,res ,size)))
+          into copylimbs
+        collect `(__gmpz_clear (addr ,gres)) into clears
+        collect res into results
+        finally (return
+                  `(with-alien ,declares
+                     ,@inits
+                     ,@body
+                     (let* ,resinits
+                       (declare (type bignum-index ,@sizes))
+                       ;; copy GMP limbs into result bignum
+                       (sb-sys:with-pinned-objects ,results
+                         ,@copylimbs)
+                       ,@clears
+                       (values ,@results))))))
+
+;;; function definition and foreign function relationships
+(defmacro defgmpfun (name args &body body)
+  `(progn
+     (declaim (sb-ext:maybe-inline ,name))
+     (defun ,name ,args
+       (declare (optimize (speed 3) (space 3) (safety 0))
+                (type integer ,@args))
+       ,@body)))
+
+\f
+;; SBCL/GMP functions
+(defgmpfun mpz-add (a b)
+  (with-mpz-results ((result (1+ (max (blength a)
+                                      (blength b)))))
+    (with-mpz-vars ((a ga) (b gb))
+      (__gmpz_add (addr result) (addr ga) (addr gb)))))
+
+(defgmpfun mpz-sub (a b)
+  (with-mpz-results ((result (1+ (max (blength a)
+                                      (blength b)))))
+    (with-mpz-vars ((a ga) (b gb))
+      (__gmpz_sub (addr result) (addr ga) (addr gb)))))
+
+(defgmpfun mpz-mul (a b)
+  (with-mpz-results ((result (+ (blength a)
+                                (blength b))))
+    (with-mpz-vars ((a ga) (b gb))
+      (__gmpz_mul (addr result) (addr ga) (addr gb)))))
+
+(defgmpfun mpz-mod (a b)
+  (with-mpz-results ((result (1+ (max (blength a)
+                                      (blength b)))))
+    (with-mpz-vars ((a ga) (b gb))
+      (__gmpz_mod (addr result) (addr ga) (addr gb))
+      (when (and (minusp (slot gb 'mp_size))
+                 (/= 0 (slot result 'mp_size)))
+        (__gmpz_add (addr result) (addr result) (addr gb))))))
+
+(defgmpfun mpz-cdiv (n d)
+  (let ((size (1+ (max (blength n)
+                       (blength d)))))
+    (with-mpz-results ((quot size)
+                       (rem size))
+      (with-mpz-vars ((n gn) (d gd))
+        (__gmpz_cdiv_qr (addr quot) (addr rem) (addr gn) (addr gd))))))
+
+(defgmpfun mpz-fdiv (n d)
+  (let ((size (1+ (max (blength n)
+                       (blength d)))))
+    (with-mpz-results ((quot size)
+                       (rem size))
+      (with-mpz-vars ((n gn) (d gd))
+        (__gmpz_fdiv_qr (addr quot) (addr rem) (addr gn) (addr gd))))))
+
+(defgmpfun mpz-tdiv (n d)
+  (let ((size (max (blength n)
+                   (blength d))))
+    (with-mpz-results ((quot size)
+                       (rem size))
+      (with-mpz-vars ((n gn) (d gd))
+        (__gmpz_tdiv_qr (addr quot) (addr rem) (addr gn) (addr gd))))))
+
+(defun mpz-pow (base exp)
+  (declare (optimize (speed 3) (space 3) (safety 0))
+           (type bignum-type base))
+  (check-type exp (unsigned-byte #.sb-vm:n-word-bits))
+  (with-gmp-mpz-results (rop)
+    (with-mpz-vars ((base gbase))
+      (__gmpz_pow_ui (addr rop) (addr gbase) exp))))
+
+(defgmpfun mpz-powm (base exp mod)
+  (with-mpz-results ((rop (1+ (blength mod))))
+    (with-mpz-vars ((base gbase) (exp gexp) (mod gmod))
+      (__gmpz_powm (addr rop) (addr gbase) (addr gexp) (addr gmod)))))
+
+(defgmpfun mpz-gcd (a b)
+  (with-mpz-results ((result (min (blength a)
+                                  (blength b))))
+    (with-mpz-vars ((a ga) (b gb))
+      (__gmpz_gcd (addr result) (addr ga) (addr gb)))))
+
+(defgmpfun mpz-lcm (a b)
+  (with-mpz-results ((result (+ (blength a)
+                                (blength b))))
+    (with-mpz-vars ((a ga) (b gb))
+      (__gmpz_lcm (addr result) (addr ga) (addr gb)))))
+
+(defgmpfun mpz-sqrt (a)
+  (with-mpz-results ((result (1+ (ceiling (blength a) 2))))
+    (with-mpz-vars ((a ga))
+      (__gmpz_sqrt (addr result) (addr ga)))))
+
+\f
+;;; Functions that use GMP-side allocated integers and copy the result
+;;; into a SBCL bignum at the end of the computation when the required
+;;; bignum length is known.
+(defun mpz-probably-prime-p (n &optional (reps 25))
+  (declare (optimize (speed 3) (space 3) (safety 0)))
+  (check-type reps fixnum)
+  (with-mpz-vars ((n gn))
+    (__gmpz_probab_prime_p (addr gn) reps)))
+
+(defun mpz-nextprime (a)
+  (declare (optimize (speed 3) (space 3) (safety 0)))
+  (with-gmp-mpz-results (prime)
+    (with-mpz-vars ((a ga))
+      (__gmpz_nextprime (addr prime) (addr ga)))))
+
+(defun mpz-fac (n)
+  (declare (optimize (speed 3) (space 3) (safety 0)))
+  (check-type n (unsigned-byte #.sb-vm:n-word-bits))
+  (with-gmp-mpz-results (fac)
+    (__gmpz_fac_ui (addr fac) n)))
+
+(defun %mpz-2fac (n)
+  (declare (optimize (speed 3) (space 3) (safety 0)))
+  (check-type n (unsigned-byte #.sb-vm:n-word-bits))
+  (with-gmp-mpz-results (fac)
+    (__gmpz_2fac_ui (addr fac) n)))
+
+(defun %mpz-mfac (n m)
+  (declare (optimize (speed 3) (space 3) (safety 0)))
+  (check-type n (unsigned-byte #.sb-vm:n-word-bits))
+  (check-type m (unsigned-byte #.sb-vm:n-word-bits))
+  (with-gmp-mpz-results (fac)
+    (__gmpz_mfac_uiui (addr fac) n m)))
+
+(defun %mpz-primorial (n)
+  (declare (optimize (speed 3) (space 3) (safety 0)))
+  (check-type n (unsigned-byte #.sb-vm:n-word-bits))
+  (with-gmp-mpz-results (r)
+    (__gmpz_primorial_ui (addr r) n)))
+
+(defun setup-5.1-stubs ()
+  (macrolet ((stubify (name implementation &rest arguments)
+               `(setf (fdefinition ',name)
+                      (if (member :sb-gmp-5.1 *gmp-features*)
+                          (fdefinition ',implementation)
+                          (lambda ,arguments
+                            (declare (ignore ,@arguments))
+                            (error "~S is only available in GMP >= 5.1"
+                                   ',name))))))
+    (stubify mpz-2fac %mpz-2fac n)
+    (stubify mpz-mfac %mpz-mfac n m)
+    (stubify mpz-primorial %mpz-primorial n)))
+
+(defun mpz-bin (n k)
+  (declare (optimize (speed 3) (space 3) (safety 0)))
+  (check-type k (unsigned-byte #.sb-vm:n-word-bits))
+  (with-gmp-mpz-results (r)
+    (with-mpz-vars ((n gn))
+      (__gmpz_bin_ui (addr r) (addr gn) k))))
+
+(defun mpz-fib2 (n)
+  (declare (optimize (speed 3) (space 3) (safety 0)))
+  ;; (let ((size (1+ (ceiling (* n (log 1.618034 2)) 64)))))
+  ;; fibonacci number magnitude in bits is assymptotic to n(log_2 phi)
+  ;; This is correct for the result but appears not to be enough for GMP
+  ;; during computation (memory access error), so use GMP-side allocation.
+  (check-type n (unsigned-byte #.sb-vm:n-word-bits))
+  (with-gmp-mpz-results (fibn fibn-1)
+    (__gmpz_fib2_ui (addr fibn) (addr fibn-1) n)))
+
+\f
+;;;; Random bignum (mpz) generation
+
+;; we do not actually use the gestalt of the struct but need its size
+;; for allocation purposes
+(define-alien-type nil
+    (struct gmprandstate
+            (mp_seed (struct gmpint))
+            (mp_alg int)
+            (mp_algdata (* t))))
+
+(declaim (inline __gmp_randinit_mt
+                 __gmp_randinit_lc_2exp
+                 __gmp_randseed
+                 __gmp_randseed_ui
+                 __gmpz_urandomb
+                 __gmpz_urandomm))
+
+(define-alien-routine __gmp_randinit_mt void
+  (s (* (struct gmprandstate))))
+
+(define-alien-routine __gmp_randinit_lc_2exp void
+  (s (* (struct gmprandstate)))
+  (a (* (struct gmpint)))
+  (c unsigned-long)
+  (exp unsigned-long))
+
+(define-alien-routine __gmp_randseed void
+  (s (* (struct gmprandstate)))
+  (sd (* (struct gmpint))))
+
+(define-alien-routine __gmp_randseed_ui void
+  (s (* (struct gmprandstate)))
+  (sd unsigned-long))
+
+(define-alien-routine __gmpz_urandomb void
+  (r (* (struct gmpint)))
+  (s (* (struct gmprandstate)))
+  (bcnt unsigned-long))
+
+(define-alien-routine __gmpz_urandomm void
+  (r (* (struct gmpint)))
+  (s (* (struct gmprandstate)))
+  (n (* (struct gmpint))))
+
+(defstruct (gmp-rstate (:constructor %make-gmp-rstate))
+  (ref (make-alien (struct gmprandstate))
+   :type (alien (* (struct gmprandstate))) :read-only t))
+
+(defun make-gmp-rstate ()
+  "Instantiate a state for the GMP Mersenne-Twister random number generator."
+  (declare (optimize (speed 3) (space 3)))
+  (let* ((state (%make-gmp-rstate))
+         (ref (gmp-rstate-ref state)))
+    (__gmp_randinit_mt ref)
+    (sb-ext:finalize state (lambda () (free-alien ref)))
+    state))
+
+(defun make-gmp-rstate-lc (a c m2exp)
+  "Instantiate a state for the GMP linear congruential random number generator."
+  (declare (optimize (speed 3) (space 3) (safety 0)))
+  (check-type c (unsigned-byte #.sb-vm:n-word-bits))
+  (check-type m2exp (unsigned-byte #.sb-vm:n-word-bits))
+  (let* ((state (%make-gmp-rstate))
+         (ref (gmp-rstate-ref state)))
+    (with-mpz-vars ((a ga))
+      (__gmp_randinit_lc_2exp ref (addr ga) c m2exp))
+    (sb-ext:finalize state (lambda () (free-alien ref)))
+    state))
+
+(defun rand-seed (state seed)
+  "Initialize a random STATE with SEED."
+  (declare (optimize (speed 3) (space 3) (safety 0)))
+  (check-type state gmp-rstate)
+  (let ((ref (gmp-rstate-ref state)))
+    (cond
+      ((typep seed '(unsigned-byte #.sb-vm:n-word-bits))
+       (__gmp_randseed_ui ref seed))
+      ((typep seed '(integer 0 *))
+       (with-mpz-vars ((seed gseed))
+         (__gmp_randseed ref (addr gseed))))
+      (t
+       (error "SEED must be a positive integer")))))
+
+(defun random-bitcount (state bitcount)
+  "Return a random integer in the range 0..(2^bitcount - 1)."
+  (declare (optimize (speed 3) (space 3) (safety 0)))
+  (check-type state gmp-rstate)
+  (check-type bitcount (unsigned-byte #.sb-vm:n-word-bits))
+  (let ((ref (gmp-rstate-ref state)))
+    (with-mpz-results ((result (+ (ceiling bitcount sb-vm:n-word-bits) 2)))
+      (__gmpz_urandomb (addr result) ref bitcount))))
+
+(defun random-int (state boundary)
+  "Return a random integer in the range 0..(boundary - 1)."
+  (declare (optimize (speed 3) (space 3) (safety 0)))
+  (check-type state gmp-rstate)
+  (let ((b (bassert boundary))
+        (ref (gmp-rstate-ref state)))
+    (with-mpz-results ((result (1+ (%bignum-length b))))
+      (with-mpz-vars ((b gb))
+        (__gmpz_urandomm (addr result) ref (addr gb))))))
+
+\f
+;;; Rational functions
+(declaim (inline %lsize))
+(defun %lsize (minusp n)
+  (declare (optimize (speed 3) (space 3) (safety 0)))
+  "n must be a (potentially denormalized) bignum"
+  (let ((length (%bignum-length n)))
+    (when (zerop (%bignum-ref n (1- length)))
+      (decf length))
+    (if minusp (- length) length)))
+
+(defmacro defmpqfun (name gmpfun)
+  `(progn
+     (declaim (sb-ext:maybe-inline ,name))
+     (defun ,name (a b)
+       (declare (optimize (speed 3) (space 3) (safety 0)))
+       (let ((size (+ (max (blength (numerator a))
+                           (blength (denominator a)))
+                      (max (blength (numerator b))
+                           (blength (denominator b)))
+                      3)))
+         (with-alien ((r (struct gmprat)))
+           (let ((num (%allocate-bignum size))
+                 (den (%allocate-bignum size)))
+             (sb-sys:with-pinned-objects (num den)
+               (setf (slot (slot r 'mp_num) 'mp_size) 0
+                     (slot (slot r 'mp_num) 'mp_alloc) size
+                     (slot (slot r 'mp_num) 'mp_d) (bignum-data-sap num))
+               (setf (slot (slot r 'mp_den) 'mp_size) 0
+                     (slot (slot r 'mp_den) 'mp_alloc) size
+                     (slot (slot r 'mp_den) 'mp_d) (bignum-data-sap den))
+               (let* ((an (bassert (numerator a)))
+                      (ad (bassert (denominator a)))
+                      (asign (not (%bignum-0-or-plusp an (%bignum-length an))))
+                      (bn (bassert (numerator b)))
+                      (bd (bassert (denominator b)))
+                      (bsign (not (%bignum-0-or-plusp bn (%bignum-length bn)))))
+                 (when asign
+                   (setf an (negate-bignum an nil)))
+                 (when bsign
+                   (setf bn (negate-bignum bn nil)))
+                 (let* ((anlen (%lsize asign an))
+                        (adlen (%lsize NIL ad))
+                        (bnlen (%lsize bsign bn))
+                        (bdlen (%lsize NIL bd)))
+                   (with-alien ((arga (struct gmprat))
+                                (argb (struct gmprat)))
+                     (sb-sys:with-pinned-objects (an ad bn bd)
+                       (setf (slot (slot arga 'mp_num) 'mp_size) anlen
+                             (slot (slot arga 'mp_num) 'mp_alloc) (abs anlen)
+                             (slot (slot arga 'mp_num) 'mp_d)
+                             (bignum-data-sap an))
+                       (setf (slot (slot arga 'mp_den) 'mp_size) adlen
+                             (slot (slot arga 'mp_den) 'mp_alloc) (abs adlen)
+                             (slot (slot arga 'mp_den) 'mp_d)
+                             (bignum-data-sap ad))
+                       (setf (slot (slot argb 'mp_num) 'mp_size) bnlen
+                             (slot (slot argb 'mp_num) 'mp_alloc) (abs bnlen)
+                             (slot (slot argb 'mp_num) 'mp_d)
+                             (bignum-data-sap bn))
+                       (setf (slot (slot argb 'mp_den) 'mp_size) bdlen
+                             (slot (slot argb 'mp_den) 'mp_alloc) (abs bdlen)
+                             (slot (slot argb 'mp_den) 'mp_d)
+                             (bignum-data-sap bd))
+                       (,gmpfun (addr r) (addr arga) (addr argb)))))
+                 (locally (declare (optimize (speed 1)))
+                   (sb-kernel::build-ratio (if (minusp (slot (slot r 'mp_num) 'mp_size))
+                                               (z-to-bignum-neg num size)
+                                               (z-to-bignum num size))
+                                           (z-to-bignum den size)))))))))))
+
+(defmpqfun mpq-add __gmpq_add)
+(defmpqfun mpq-sub __gmpq_sub)
+(defmpqfun mpq-mul __gmpq_mul)
+(defmpqfun mpq-div __gmpq_div)
+
+\f
+;;;; SBCL interface and integration installation
+(macrolet ((def (name original)
+             (let ((special (intern (format nil "*~A-FUNCTION*" name))))
+               `(progn
+                  (declaim (type function ,special)
+                           (inline ,name))
+                  (defvar ,special (symbol-function ',original))
+                  (defun ,name (&rest args)
+                    (apply (load-time-value ,special t) args))))))
+  (def orig-mul multiply-bignums)
+  (def orig-truncate bignum-truncate)
+  (def orig-gcd bignum-gcd)
+  (def orig-lcm sb-kernel:two-arg-lcm)
+  (def orig-isqrt isqrt)
+  (def orig-two-arg-+ sb-kernel:two-arg-+)
+  (def orig-two-arg-- sb-kernel:two-arg--)
+  (def orig-two-arg-* sb-kernel:two-arg-*)
+  (def orig-two-arg-/ sb-kernel:two-arg-/))
+
+;;; integers
+(defun gmp-mul (a b)
+  (declare (optimize (speed 3) (space 3))
+           (type bignum-type a b)
+           (inline mpz-mul))
+  (if (or (< (min (%bignum-length a)
+                  (%bignum-length b))
+             6)
+          *gmp-disabled*)
+      (orig-mul a b)
+      (mpz-mul a b)))
+
+(defun gmp-truncate (a b)
+  (declare (optimize (speed 3) (space 3))
+           (type bignum-type a b)
+           (inline mpz-tdiv))
+  (if (or (< (min (%bignum-length a)
+                  (%bignum-length b))
+             3)
+          *gmp-disabled*)
+      (orig-truncate a b)
+      (mpz-tdiv a b)))
+
+(defun gmp-lcm (a b)
+  (declare (optimize (speed 3) (space 3))
+           (type integer a b)
+           (inline mpz-lcm))
+  (if (or (and (typep a 'fixnum)
+               (typep b 'fixnum))
+          *gmp-disabled*)
+      (orig-lcm a b)
+      (mpz-lcm a b)))
+
+(defun gmp-isqrt (n)
+  (declare (optimize (speed 3) (space 3))
+           (type unsigned-byte n)
+           (inline mpz-sqrt))
+  (if (or (typep n 'fixnum)
+          *gmp-disabled*)
+      (orig-isqrt n)
+      (mpz-sqrt n)))
+
+;;; rationals
+(defun gmp-two-arg-+ (x y)
+  (declare (optimize (speed 3) (space 3))
+           (inline mpq-add))
+  (if (and (or (typep x 'ratio)
+               (typep y 'ratio))
+           (rationalp y)
+           (rationalp x)
+           (not *gmp-disabled*))
+      (mpq-add x y)
+      (orig-two-arg-+ x y)))
+
+(defun gmp-two-arg-- (x y)
+  (declare (optimize (speed 3) (space 3))
+           (inline mpq-sub))
+  (if (and (or (typep x 'ratio)
+               (typep y 'ratio))
+           (rationalp y)
+           (rationalp x)
+           (not *gmp-disabled*))
+      (mpq-sub x y)
+      (orig-two-arg-- x y)))
+
+(defun gmp-two-arg-* (x y)
+  (declare (optimize (speed 3) (space 3))
+           (inline mpq-mul))
+  (if (and (or (typep x 'ratio)
+               (typep y 'ratio))
+           (rationalp y)
+           (rationalp x)
+           (not *gmp-disabled*))
+      (mpq-mul x y)
+      (orig-two-arg-* x y)))
+
+(defun gmp-two-arg-/ (x y)
+  (declare (optimize (speed 3) (space 3))
+           (inline mpq-div))
+  (if (and (rationalp x)
+           (rationalp y)
+           (not (eql y 0))
+           (not *gmp-disabled*))
+      (mpq-div x y)
+      (orig-two-arg-/ x y)))
+
+;;; installation
+(defmacro with-package-locks-ignored (&body body)
+  `(handler-bind ((sb-ext:package-lock-violation
+                    (lambda (condition)
+                      (declare (ignore condition))
+                      (invoke-restart :ignore-all))))
+     ,@body))
+
+(defun install-gmp-funs ()
+  (with-package-locks-ignored
+      (macrolet ((def (destination source)
+                   `(setf (fdefinition ',destination)
+                          (fdefinition ',source))))
+        (def multiply-bignums gmp-mul)
+        (def bignum-truncate gmp-truncate)
+        (def bignum-gcd mpz-gcd)
+        (def sb-kernel:two-arg-lcm gmp-lcm)
+        (def sb-kernel:two-arg-+ gmp-two-arg-+)
+        (def sb-kernel:two-arg-- gmp-two-arg--)
+        (def sb-kernel:two-arg-* gmp-two-arg-*)
+        (def sb-kernel:two-arg-/ gmp-two-arg-/)
+        (def isqrt gmp-isqrt)))
+  (values))
+
+(defun uninstall-gmp-funs ()
+  (with-package-locks-ignored
+      (macrolet ((def (destination source)
+                   `(setf (fdefinition ',destination)
+                          ,(intern (format nil "*~A-FUNCTION*" source)))))
+        (def multiply-bignums orig-mul)
+        (def bignum-truncate orig-truncate)
+        (def bignum-gcd orig-gcd)
+        (def sb-kernel:two-arg-lcm orig-lcm)
+        (def sb-kernel:two-arg-+ orig-two-arg-+)
+        (def sb-kernel:two-arg-- orig-two-arg--)
+        (def sb-kernel:two-arg-* orig-two-arg-*)
+        (def sb-kernel:two-arg-/ orig-two-arg-/)
+        (def isqrt orig-isqrt)))
+  (values))
+
+(defun load-gmp (&key (persistently t))
+  (setf *gmp-features* nil
+        *gmp-version* nil
+        *features* (set-difference *features* '(:sb-gmp :sb-gmp-5.0 :sb-gmp-5.1)))
+  (when persistently
+    (pushnew 'load-gmp sb-ext:*init-hooks*)
+    (pushnew 'uninstall-gmp-funs sb-ext:*save-hooks*))
+  (let ((success (%load-gmp)))
+    (when success
+      (setf *gmp-version* (extern-alien "__gmp_version" c-string)))
+    (cond ((null *gmp-version*))
+          ((string<= *gmp-version* "5.")
+           (warn "SB-GMP requires at least GMP version 5.0")
+           (setf success nil))
+          (t
+           (pushnew :sb-gmp *gmp-features*)
+           (pushnew :sb-gmp-5.0 *gmp-features*)
+           (when (string>= *gmp-version* "5.1")
+             (pushnew :sb-gmp-5.1 *gmp-features*))
+           (setf *features* (union *features* *gmp-features*))))
+    (if success
+        (install-gmp-funs)
+        (uninstall-gmp-funs))
+    (setup-5.1-stubs)
+    success))
+
+(defun unload-gmp ()
+  (setf sb-ext:*init-hooks* (remove 'load-gmp sb-ext:*init-hooks*))
+  (uninstall-gmp-funs)
+  (setf sb-ext:*save-hooks* (remove 'uninstall-gmp-funs sb-ext:*save-hooks*))
+  (values))
+
+(load-gmp)
diff --git a/contrib/sb-gmp/sb-gmp.asd b/contrib/sb-gmp/sb-gmp.asd
new file mode 100644 (file)
index 0000000..8d595b0
--- /dev/null
@@ -0,0 +1,37 @@
+(defpackage #:sb-gmp-system (:use #:asdf #:cl))
+
+(in-package #:sb-gmp-system)
+
+(defsystem sb-gmp
+  :name "SB-GMP"
+  :version "0.1"
+  :description "bignum calculations for SBCL using the GMP library"
+  :serial t
+  :components ((:module sb-gmp
+                :pathname ""
+                :components ((:file "gmp")))))
+
+(defsystem sb-gmp-tests
+  :depends-on (sb-rt sb-gmp)
+  :components ((:file "tests")))
+
+(defmethod perform :after ((o load-op) (c (eql (find-system :sb-gmp))))
+  (provide 'sb-gmp))
+
+(defmethod perform ((o test-op) (c (eql (find-system :sb-gmp))))
+  (operate 'load-op 'sb-gmp)
+  (cond ((member :sb-gmp *features*)
+         (operate 'load-op 'sb-gmp-tests)
+         (operate 'test-op 'sb-gmp-tests))
+        (t
+         (warn "unable to test sb-gmp: libgmp unavailable"))))
+
+(defmethod perform ((o test-op) (c (eql (find-system :sb-gmp-tests))))
+  (multiple-value-bind (soft strict pending)
+      (funcall (intern "DO-TESTS" (find-package "SB-RT")))
+    (declare (ignorable pending))
+    (fresh-line)
+    (unless strict
+      (warn "ignoring expected failures in sb-gmp-tests"))
+    (unless soft
+      (error "sb-gmp-tests failed with unexpected failures"))))
diff --git a/contrib/sb-gmp/tests.lisp b/contrib/sb-gmp/tests.lisp
new file mode 100644 (file)
index 0000000..4555a1d
--- /dev/null
@@ -0,0 +1,197 @@
+(defpackage "SB-GMP-TESTS"
+  (:use "COMMON-LISP" "SB-GMP" "SB-RT"))
+
+(in-package "SB-GMP-TESTS")
+
+(defparameter *state* (make-gmp-rstate))
+(rand-seed *state* 1234)
+
+(defmacro defgenerator (name arguments &body body)
+  `(defun ,name ,arguments
+     (lambda () ,@body)))
+
+(defgenerator gen-mpz (&key (limbs 5) sign nonzero)
+  (let ((integer (random-bitcount *state*
+                                  (* limbs sb-vm:n-word-bits))))
+    (when (and nonzero (zerop integer))
+      (setf integer 1))
+    (ecase sign
+      ((+ nil) integer)
+      (- (- integer))
+      ((t random) (if (zerop (random 2))
+                      integer
+                      (- integer))))))
+
+(defun gen-mpq (&key (limbs 5) sign nonzero)
+  (let ((numerator (gen-mpz :limbs limbs :sign sign
+                            :nonzero nonzero))
+        (denominator (gen-mpz :limbs limbs :nonzero t)))
+    (lambda ()
+      (/ (funcall numerator) (funcall denominator)))))
+
+(defun maybe-apply (maybe-function &optional arguments)
+  (if (typep maybe-function '(or function symbol))
+      (apply maybe-function arguments)
+      maybe-function))
+
+(defun test-one-case (base tested &rest arguments)
+  (let* ((initial-hashes (mapcar #'sxhash arguments))
+         (base-values (let ((*gmp-disabled* t))
+                        (multiple-value-list (maybe-apply base arguments))))
+         (test-values (let ((*gmp-disabled* nil))
+                        (multiple-value-list (apply tested arguments))))
+         (final-hashes (mapcar #'sxhash arguments)))
+    (unless (and (= (length base-values) (length test-values))
+                 (every #'eql base-values test-values))
+      (error "Failed test: ~S returned ~S; expected ~S"
+             (cons tested arguments) test-values base-values))
+    (unless (every #'eql initial-hashes final-hashes)
+      (error "Failed test: ~S modified arguments ~{~A~^, ~} ~
+              (printed modified values)"
+             (cons tested arguments)
+             (loop for i upfrom 0
+                   for initial in initial-hashes
+                   for final in final-hashes
+                   unless (eql initial final)
+                     collect i))))
+  nil)
+
+;; Really just the most basic smoke test, otherwise
+;; build times ballon up a bit on slow machines.
+(defvar *iteration-count* 3)
+
+(defun test-n-cases (base tested &rest argument-generators)
+  (let ((*random-state* (sb-ext:seed-random-state 54321)))
+    (loop repeat *iteration-count* do
+      (apply 'test-one-case base tested
+             (mapcar #'maybe-apply argument-generators)))))
+
+(defmacro define-gmp-test ((name &key (repeat 1) limbs (gmp-seed 1234))
+                           &body body)
+  `(deftest ,name
+       (let ((*random-state* (sb-ext:seed-random-state 54321)))
+         (rand-seed *state* ,gmp-seed)
+         (handler-case
+             (dotimes (i ,repeat)
+               ;; try to get small failures first
+               (let ((limbs (case i
+                              (0 ,(subst `(lambda (x)
+                                            x 0)
+                                         'random
+                                         limbs))
+                              (1 ,(subst `(lambda (x)
+                                            (if (> x 1) 1 0))
+                                         'random
+                                         limbs))
+                              (t ,limbs))))
+                 (declare (ignorable limbs))
+                 ,@body))
+           (error (c)
+             (format t "~&~A~%" c)
+             nil)
+           (:no-error (&rest _) _ t)))
+     t))
+
+(define-gmp-test (mpz-add :repeat 7 :limbs (+ (random #xFFFFF) 2))
+  (test-n-cases '+ 'mpz-add
+                    (gen-mpz :limbs limbs :sign t)
+                    (gen-mpz :limbs limbs :sign t)))
+
+(define-gmp-test (mpz-sub :repeat 7 :limbs (+ (random #x1FFFF) 2))
+  (test-n-cases '- 'mpz-sub
+                (gen-mpz :limbs limbs :sign t)
+                (gen-mpz :limbs limbs :sign t)))
+
+(define-gmp-test (mpz-mul :repeat 7 :limbs (+ (random #x253F) 2))
+  (test-n-cases '* 'mpz-mul
+                (gen-mpz :limbs limbs :sign t)
+                (gen-mpz :limbs limbs :sign t)))
+
+(define-gmp-test (mpz-tdiv :repeat 7 :limbs (+ (random #x253F) 2))
+  (test-n-cases 'truncate 'mpz-tdiv
+                (gen-mpz :limbs limbs :sign t)
+                (gen-mpz :limbs limbs :sign t :nonzero t)))
+
+(define-gmp-test (mpz-fdiv :repeat 7 :limbs (+ (random #x253F) 2))
+  (test-n-cases 'floor 'mpz-fdiv
+                (gen-mpz :limbs limbs :sign t)
+                (gen-mpz :limbs limbs :sign t :nonzero t)))
+
+(define-gmp-test (mpz-cdiv :repeat 7 :limbs (+ (random #x253F) 2))
+  (test-n-cases 'ceiling 'mpz-cdiv
+                (gen-mpz :limbs limbs :sign t)
+                (gen-mpz :limbs limbs :sign t :nonzero t)))
+
+(define-gmp-test (mpz-gcd :repeat 7 :limbs (+ (random #x253F) 2))
+  (test-n-cases 'gcd 'mpz-gcd
+                (gen-mpz :limbs limbs :sign t)
+                (gen-mpz :limbs limbs :sign t)))
+
+(define-gmp-test (mpz-lcm :repeat 7 :limbs (+ (random #x253F) 2))
+  (test-n-cases 'lcm 'mpz-lcm
+                (gen-mpz :limbs limbs :sign t)
+                (gen-mpz :limbs limbs :sign t)))
+
+(define-gmp-test (isqrt :repeat 7 :limbs (+ (random #x253F) 2))
+  (test-n-cases 'isqrt 'mpz-sqrt (gen-mpz :limbs limbs)))
+
+(define-gmp-test (mpz-mod :repeat 7 :limbs (1+ (random #x253F)))
+  (test-n-cases 'mod 'mpz-mod
+                (gen-mpz :limbs limbs :sign t)
+                (gen-mpz :limbs limbs :sign t :nonzero t)))
+
+(define-gmp-test (mpz-powm :repeat 7 :limbs (1+ (random #x253F)))
+  (test-n-cases (lambda (base exponent mod)
+                  (let ((*gmp-disabled* nil)) ; atrociously slow otherwise
+                    (mod (expt base exponent) mod)))
+                'mpz-powm
+                (gen-mpz :limbs limbs :sign t)
+                (lambda ()
+                  (1+ (random 40)))
+                (gen-mpz :limbs (ceiling limbs 2) :nonzero t)))
+
+;; bugs that have been fixed
+(define-gmp-test (sign-conversion)
+  (test-one-case '+ 'mpz-add #x7FFFFFFFFFFFFFFF #x7FFFFFFFFFFFFFFF))
+(define-gmp-test (truncate-1)
+  (test-one-case 'truncate 'mpz-tdiv
+                 30951488519636377404900619671461408624764773310745985021994671444676860083493
+                 200662724990805535745252242839121922075))
+(define-gmp-test (truncate-2)
+  (test-one-case 'truncate 'mpz-tdiv
+                 320613729464106236061704728914573914390
+                 -285049280629101090500613812618405407883))
+
+(define-gmp-test (mpz-nextprime :repeat 7
+                                :gmp-seed 6234
+                                :limbs (1+ (random #x2F)))
+  (let ((a (gen-mpz :limbs limbs)))
+    (dotimes (i *iteration-count*)
+      (let* ((a (funcall a))
+             (p (mpz-nextprime a)))
+        (assert (>= p a))
+        (assert (plusp (mpz-probably-prime-p p)))))))
+
+(define-gmp-test (mpq-add :repeat 7 :limbs (1+ (random #x3FF))
+                          :gmp-seed 1235)
+  (test-n-cases '+ 'mpq-add
+                (gen-mpq :limbs limbs :sign t)
+                (gen-mpq :limbs limbs :sign t)))
+
+(define-gmp-test (mpq-sub :repeat 7 :limbs (1+ (random #x1FF))
+                          :gmp-seed 1235)
+  (test-n-cases '- 'mpq-sub
+                (gen-mpq :limbs limbs :sign t)
+                (gen-mpq :limbs limbs :sign t)))
+
+(define-gmp-test (mpq-mul :repeat 7 :limbs (1+ (random #x5FF))
+                          :gmp-seed 6235)
+  (test-n-cases '* 'mpq-mul
+                (gen-mpq :limbs limbs :sign t)
+                (gen-mpq :limbs limbs :sign t)))
+
+(define-gmp-test (mpq-div :repeat 7 :limbs (1+ (random #x3FF))
+                          :gmp-seed 7235)
+  (test-n-cases '/ 'mpq-div
+                (gen-mpq :limbs limbs :sign t)
+                (gen-mpq :limbs limbs :sign t)))