;; author: richard fateman, jan, 2006 ;;; extracted from the generic arith stuff. 7/2010 ;;; hacked to be faster. ACL specific. ;; this file is based in part on lisp/loadgmp2.cl ?2002? ;; which also has polynomial stuff in it. ;; working univariate polynomial stuff is in lisp/polymultx.cl though. ;(eval-when '(load) (require "ga")(provide "gmp")) ;(in-package :gmp) (eval-when (compile load) (declaim (optimize (speed 3) (safety 0) (space 0) (compilation-speed 0))) ) ;; if i learned about uffi i suppose it could be done in a uniform way ;; for any of the participating common lisp implementations. ;; load a working gmp library. ;;(load "h:/lisp/generic/p4/gmp.dll") ;; if this is a pentium 4 ;;(load "h:/lisp/generic/p3/gmp.dll") ;; if this is a pentium 3 ;;(load "h:/lisp/generic/p0/gmp.dll") ;; if this is a pentium ;;(load "gmp.dll") ;; one should work.. ;; the next section's syntax is particular for allegro cl. (eval-when (compile load eval) (ff:def-foreign-type mpz (:array :int)) (ff:def-foreign-call (mpz_init "__gmpz_init") ((x (cl::* :int))) :returning :int :arg-checking nil :call-direct t) (ff:def-foreign-call (mpz_init2 "__gmpz_init2") ((x (cl::* :int)) (len :int) ) ;; space for how many bits? :returning :int :arg-checking nil :call-direct t) (ff:def-foreign-call (mpz_fac "__gmpz_fac_ui") ;; factorial ((x (cl::* :int)) (n :int) ) ;; compute n! :returning :int :arg-checking nil :call-direct t) (ff:def-foreign-call (mpz_set_str "__gmpz_set_str") ;changes value of mpz x. ((x mpz (simple-array (signed-byte 32) (3))) (s (cl::* :char)) (base :int)) :strings-convert t :returning :int :arg-checking nil) ;; :call-direct t (ff:def-foreign-call (mpz_get_str "__gmpz_get_str") ((s :int) (base :int) (op mpz (simple-array (signed-byte 32) (3)))) ;;i think it is ok to let gmp allocate the space for this string. :returning ((cl::* :char) ) :arg-checking nil :call-direct t) (ff:def-foreign-call (mpz_get_d "__gmpz_get_d") ((x mpz (simple-array (signed-byte 32) (3)))) :returning :double :arg-checking nil :call-direct t) (ff:def-foreign-call;; remove gmp number from memory (mpz_clear "__gmpz_clear") ((x mpz (simple-array (signed-byte 32) (3)))) :returning :int :arg-checking nil :call-direct t) (ff:def-foreign-call;; compare (mpz_cmp "__gmpz_cmp") ((x mpz (simple-array (signed-byte 32) (3))) (y mpz (simple-array (signed-byte 32) (3))) ); returns pos if x>y, zero if x=y, neg if x0. If b=0, what do we do? (defmethod qr((a gmpz)(b gmpz)) (let ((r (create_mpz_zero)) (q (create_mpz_zero))) (mpz_fdiv_qr (gmpz-z q)(gmpz-z r)(gmpz-z a)(gmpz-z b)) (values r q)) ) ;;(defarithmetic expt) hm. can't do negative exponent. ;; not sure we want to do ;; (expt GMP GMP) (defun gmp_gcd (a b) (let((ans (alloc-gmpz))) (setf a (gmpz-z(into a))) (setf b (gmpz-z(into b))) ;; make both into gmpz numbers if not already (mpz_gcd (gmpz-z ans) a b) ans)) (defun gmp_gcdext (a b) (let((ss (alloc-gmpz)) (tt (alloc-gmpz)) (gg (alloc-gmpz))) (setf a (gmpz-z(into a))) (setf b (gmpz-z(into b))) ;; make both into gmpz numbers if not already (mpz_gcdext (gmpz-z gg) (gmpz-z ss)(gmpz-z tt) a b) (values gg ss tt))) ;returns g,s,t such that g=a*s+b*t and g>0 ;; programs that decompose an Allegro CL number into an ;; array of 16 bit quantities suitable for piling into GMP are ;; in the appendix to papers/polysbyGMP.tex. Copied here, too. #+ignore (defmethod lisp2gmpz ((x integer)) ;; but not fixnum (let* ((out (create_mpz_zero)) ;the wrapper (r (gmpz-z out)) ;the inside (a (btoa x));; bignum to array. (out2(create_mpz_zero)) (h (gmpz-z out2))) (do ((i (1- (length a))(1- i))) ((< i 0) r) (declare (fixnum i)) (mpz_set_si h (aref a i)) (mpz_mul_2exp r r 16) (mpz_add r r h)) (if (< x 0) (setf (aref r 1)(- (aref r 1)))) out)) #+ignore (defmethod lisp2gmpz ((x integer)) ;; but not fixnum (let* ((out (create_mpz_zero)) ;the wrapper (r (gmpz-z out)) ;the inside (a (btoa x));; bignum to array. ) (do ((i (1- (length a))(1- i))) ((< i 0) r) (declare (fixnum i)) (mpz_mul_2exp r r 16) (mpz_add_ui r r (aref a i))) (if (< x 0) (setf (aref r 1)(- (aref r 1)))) out)) (defmethod lisp2gmpz ((x integer)) ; x is lisp bignum, answer is gmp number (declare (optimize (speed 3)(safety 0))) (let* ((out (create_mpz_zero)) ;the wrapper (r (gmpz-z out))) ;the inside (if (fixnump x)(mpz_add_ui r r x) (do ((i (1-(excl::bm_size x))(1- i))) ((< i 0) r) (declare (fixnum i)) (mpz_mul_2exp r r 16) (mpz_add_ui r r (nth-bigitm x (+ i i))) )) (if (< x 0) (setf (aref r 1)(- (aref r 1)))) out)) ;;B to A: Bignum or fixnum to Array of bigits, 16 bits (defun btoa(x) (if (excl::fixnump x) (vector (abs x)) (let* ((size (excl::bm_size x)) ;; size of an ACL bignum in 16-bit quantities (ans (make-array size :element-type '(unsigned-byte 16)))) (declare (optimize (speed 3)(safety 0)(debug 0))(fixnum size)) (do ((i 0 (1+ i))) ((= i size) ans) (declare (fixnum i)) (setf (aref ans i)(nth-bigit x i)))))) (eval-when (compile load)(require :llstructs)) (defun nth-bigit (x n) (sys:memref x -14 (* 2 n) :unsigned-word)) ;; this is the place to add more interface program declarations. ;; we're not even using all the ones above, yet. ;; there are about 165 gmpz signed integer functions, ;; about 35 gmpq rational functions, ;; function entries etc should all be in the gmp.h file