Proper error messages instead of "foo" and "bar"
[sbcl.git] / contrib / sb-gmp / gmp.lisp
blobd054e29c6d23767d0b5c73d2cce1b30ac110ff05
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*
58 (in-package "SB-GMP")
60 (defvar *gmp-disabled* nil)
62 (defconstant +bignum-raw-area-offset+
63 (- (* sb-vm:bignum-digits-offset sb-vm:n-word-bytes)
64 sb-vm:other-pointer-lowtag))
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+))
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)))
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.")))
86 (defvar *gmp-features* nil)
87 (defvar *gmp-version* nil)
89 ;; We load only the library right now to avoid undefined alien
90 ;; style warnings
91 (%load-gmp)
94 ;;; types and initialization
96 (define-alien-type gmp-limb
97 #-(and win32 x86-64) unsigned-long
98 #+(and win32 x86-64) unsigned-long-long)
100 (deftype ui ()
101 #-(and win32 x86-64) '(unsigned-byte #.sb-vm:n-word-bits)
102 #+(and win32 x86-64) '(unsigned-byte 32))
104 (deftype si ()
105 #-(and win32 x86-64) '(signed-byte #.sb-vm:n-word-bits)
106 #+(and win32 x86-64) '(signed-byte 32))
108 (define-alien-type nil
109 (struct gmpint
110 (mp_alloc int)
111 (mp_size int)
112 (mp_d (* gmp-limb))))
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."
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))
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)
133 (the unsigned-byte (%normalize-bignum b count))))
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)))
144 ;;; conversion functions that also copy from GMP to SBCL bignum space
145 (declaim (inline gmp-z-to-bignum))
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))))
157 (declaim (inline blength bassert)
158 (ftype (function (integer) (values bignum-length &optional)) blength)
159 (ftype (function (integer) (values bignum &optional)) bassert))
161 (defun blength (a)
162 (declare (optimize (speed 3) (space 3) (safety 0)))
163 (etypecase a
164 (fixnum 1)
165 (t (%bignum-length a))))
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)))
173 ;;;; rationals
174 (define-alien-type nil
175 (struct gmprat
176 (mp_num (struct gmpint))
177 (mp_den (struct gmpint))))
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.
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))))
193 (define-alien-routine __gmpz_clear void
194 (x (* (struct gmpint))))
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))))
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))))
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))))
229 (define-twoarg-mpz-funs (__gmpz_sqrt
230 __gmpz_nextprime))
232 (define-threearg-mpz-funs (__gmpz_add
233 __gmpz_sub
234 __gmpz_mul
235 __gmpz_mod
236 __gmpz_gcd
237 __gmpz_lcm))
239 (define-fourarg-mpz-funs (__gmpz_cdiv_qr
240 __gmpz_fdiv_qr
241 __gmpz_tdiv_qr
242 __gmpz_powm))
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))
256 (define-alien-routine __gmpz_mul_2exp void
257 (r (* (struct gmpint)))
258 (b (* (struct gmpint)))
259 (e unsigned-long))
261 (define-alien-routine __gmpz_fdiv_q_2exp void
262 (r (* (struct gmpint)))
263 (b (* (struct gmpint)))
264 (e unsigned-long))
266 (define-alien-routine __gmpz_pow_ui void
267 (r (* (struct gmpint)))
268 (b (* (struct gmpint)))
269 (e unsigned-long))
271 (define-alien-routine __gmpz_probab_prime_p int
272 (n (* (struct gmpint)))
273 (reps int))
275 (define-alien-routine __gmpz_fac_ui void
276 (r (* (struct gmpint)))
277 (a unsigned-long))
279 (define-alien-routine __gmpz_2fac_ui void
280 (r (* (struct gmpint)))
281 (a unsigned-long))
283 (define-alien-routine __gmpz_mfac_uiui void
284 (r (* (struct gmpint)))
285 (n unsigned-long)
286 (m unsigned-long))
288 (define-alien-routine __gmpz_primorial_ui void
289 (r (* (struct gmpint)))
290 (n unsigned-long))
292 (define-alien-routine __gmpz_remove unsigned-long
293 (r (* (struct gmpint)))
294 (x (* (struct gmpint)))
295 (f (* (struct gmpint))))
297 (define-alien-routine __gmpz_bin_ui void
298 (r (* (struct gmpint)))
299 (n (* (struct gmpint)))
300 (k unsigned-long))
302 (define-alien-routine __gmpz_fib2_ui void
303 (r (* (struct gmpint)))
304 (a (* (struct gmpint)))
305 (b unsigned-long))
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))))
319 (define-threearg-mpq-funs (__gmpq_add
320 __gmpq_sub
321 __gmpq_mul
322 __gmpq_div))
325 ;;;; SBCL interface
327 ;;; utility macros for GMP mpz variable and result declaration and
328 ;;; incarnation of associated SBCL bignums
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-kernel::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))))))))
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))))))
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-kernel::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)))))))
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)))
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)))))
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)))))
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)))))
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))))
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))))))
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))))))
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))))))
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))))
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))))))
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))))
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)))))
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)))))
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)))))
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)))))
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)))
532 (defgmpfun mpz-nextprime (a)
533 (with-gmp-mpz-results (prime)
534 (with-mpz-vars ((a ga))
535 (__gmpz_nextprime (addr prime) (addr ga)))))
537 (defgmpfun mpz-fac (n)
538 (check-type n ui)
539 (with-gmp-mpz-results (fac)
540 (__gmpz_fac_ui (addr fac) n)))
542 (defgmpfun %mpz-2fac (n)
543 (check-type n ui)
544 (with-gmp-mpz-results (fac)
545 (__gmpz_2fac_ui (addr fac) n)))
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)))
553 (defgmpfun %mpz-primorial (n)
554 (check-type n ui)
555 (with-gmp-mpz-results (r)
556 (__gmpz_primorial_ui (addr r) n)))
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)))
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)))
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)))
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))))
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)))
607 ;;;; Random bignum (mpz) generation
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))))
617 (declaim (inline __gmp_randinit_mt
618 __gmp_randinit_lc_2exp
619 __gmp_randseed
620 __gmp_randseed_ui
621 __gmpz_urandomb
622 __gmpz_urandomm))
624 (define-alien-routine __gmp_randinit_mt void
625 (s (* (struct gmprandstate))))
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))
633 (define-alien-routine __gmp_randseed void
634 (s (* (struct gmprandstate)))
635 (sd (* (struct gmpint))))
637 (define-alien-routine __gmp_randseed_ui void
638 (s (* (struct gmprandstate)))
639 (sd unsigned-long))
641 (define-alien-routine __gmpz_urandomb void
642 (r (* (struct gmpint)))
643 (s (* (struct gmprandstate)))
644 (bcnt unsigned-long))
646 (define-alien-routine __gmpz_urandomm void
647 (r (* (struct gmpint)))
648 (s (* (struct gmprandstate)))
649 (n (* (struct gmpint))))
651 (defstruct (gmp-rstate (:constructor %make-gmp-rstate))
652 (ref (make-alien (struct gmprandstate))
653 :type (alien (* (struct gmprandstate))) :read-only t))
655 (declaim (sb-ext:maybe-inline make-gmp-rstate
656 make-gmp-rstate-lc
657 rand-seed
658 random-bitcount
659 random-int))
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))
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))
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))))
694 (error "SEED must be a positive integer")))))
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))))
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))))))
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)))
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))))))
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)))))))))))
807 (defmpqfun mpq-add __gmpq_add)
808 (defmpqfun mpq-sub __gmpq_sub)
809 (defmpqfun mpq-mul __gmpq_mul)
810 (defmpqfun mpq-div __gmpq_div)
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))
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))
841 *gmp-disabled*)
842 (orig-mul a b)
843 (mpz-mul a b)))
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))
852 *gmp-disabled*)
853 (orig-truncate a b)
854 (mpz-tdiv a b)))
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)))
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)))
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)))
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)))
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)))
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)))
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))
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)))
937 (mpz-pow base power))))))
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))
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))
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))
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))
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))
1006 (load-gmp)