Diff of /contrib/sb-gmp/gmp.lisp [000000] .. [1656e5] Maximize Restore

  Switch to unified view

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