;;; -*- mode: lisp; syntax: common-lisp; package: gmp; base: 10 -*- ;; author: richard fateman, jan, 2006 ;; 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")) (defpackage :gmp (:use :cl)) (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 #|trying to clear the air about different versions of gmp 4.1.4 for different architectures. consider the computation of 20,000! by the program k4, below. how much does it matter which library we use? we tested this first bunch on a pentium 4. (we do the computation 10 times) (time (progn (dotimes (i 10)(k4 20000)) nil)) with generic code for any i86, written in c: 1,219 ms. pentium-generic/gmp.dll with some code specifically for pentium 0: 1,250 ms. p0/gmp.dll with some code specifically for pentium 3: 641 ms. p3/gmp.dll with some code specifically for pentium 4: 437 ms. p4/gmp.dll using a pentium 3 with generic code for any i86, written in c: 2,083 ms. pentium-generic/gmp.dll with some code specifically for pentium 0: 792 ms. p0/gmp.dll with some code specifically for pentium 3: 791 ms. p3/gmp.dll with some code specifically for pentium 4: p4/gmp.dll oh, using mpz_fac, directly from gmp 661 ms. |# ;;(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_mul_2exp "__gmpz_mul_2exp");; set target to op1*2^op2 ((target mpz (simple-array (signed-byte 32) (3)))(op1 mpz (simple-array (signed-byte 32) (3)))(op2 :long)) ;**** :returning :void :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_d "__gmpz_get_d") ((x mpz (simple-array (signed-byte 32) (3)))) :returning :double :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;; 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 (mpz_mul "__gmpz_mul") ((target mpz (simple-array (signed-byte 32) (3))) (op1 mpz (simple-array (signed-byte 32) (3))) (op2 mpz (simple-array (signed-byte 32) (3)))) :returning :void :arg-checking nil :call-direct t :strings-convert nil) (ff:def-foreign-call (mpz_add "__gmpz_add") ((target mpz (simple-array (signed-byte 32) (3)))(op1 mpz (simple-array (signed-byte 32) (3)))(op2 mpz (simple-array (signed-byte 32) (3)))) :returning :void :arg-checking nil :call-direct t) (ff:def-foreign-call (mpz_addmul "__gmpz_addmul");;rop <-rop+op1*op2 ((rop mpz (simple-array (signed-byte 32) (3)))(op1 mpz (simple-array (signed-byte 32) (3)))(op2 mpz (simple-array (signed-byte 32) (3)))) :returning :void :arg-checking nil :call-direct t) (ff:def-foreign-call;; convert fixnum y to mpz and store in x (mpz_set_si "__gmpz_set_si") ((x mpz (simple-array (signed-byte 32) (3))) (y :fixnum) ) :returning :void :arg-checking nil :call-direct t ) ;; sgn is not available because mpz_sgn is a macro (ff:def-foreign-call (mpz_neg "__gmpz_neg");; set targ to -op1 ((targ mpz (simple-array (signed-byte 32) (3)))(op1 mpz (simple-array (signed-byte 32) (3)))) :returning :void :arg-checking nil :call-direct t) ;; there are about 165 gmpz signed integer functions, ;; about 35 gmpq rational functions, ;; function entries etc should all be in the gmp.h file ) ;; if we want to use other parts of gmp, e.g. rationals or reals, we ;; could overload separately; e.g. another file for mpfr for ;; floating point versions. mpfr makes more sense than integers if we ;; really want to use sine/cosine/ etc or maybe even interval ;; versions. just gmpz has some functions of interest to number theory. ;; a gmpz object p represents an integer: p in z, the integers. ;; we set up the lisp system to keep track of all accessible gmpz ;; objects with its own garbage collector, but in an indirect fashion. ;; when a gmpz object is found, by the garbage collector, to be ;; inaccessible, the lisp gc can't just free it; it doesn't know how. ;; instead we signal lisp to call the gmp library program to clear its ;; own memory. this is done through the use of finalizations, which ;; are not required by the ansi standard but seem to be a feature in ;; many common lisp implementations. here we are using allegro ;; cl. other implementations which include finalizations equivalent to ;; allegro's include lispworks and gnu clisp. ;; we define gmpz, gmp integer structure as a structure with only one ;; entry to allow us to use defmethod. the one entry turns ;; out to be an array that is the gmpz datum. (defstruct gmpz z) ;;here z is the "guts" of the repr. a 3-word array whose data ;;are managed by gmp library. ;;how to print a gmpz number so that it is distinguishable from ;;a lisp number? set this format. the first choice makes 123{g}. ;;the second choice makes it indistinguishable from lisp. (defvar gmpzformat "~a\{g\}" ) ;; 123{g} ;(defvar gmpzformat "~a" ) ;;another possibility (defmethod print-object ((a gmpz) stream)(format stream gmpzformat (create_string_from_mpz a))) ;; storage deallocation for the gmpz objects requires some delicacy. the ;; guts of a gmz object is not garbage-collected automatically by ;; lisp because it is allocated in gmp-library controlled space. ;; but when the wrapper is about to be reclaimed because no one ;; is looking at it, the finalization method is run on the guts of it, ;; which tells the gmp library to deallocate the major part of the ;; number-- the sequence of "limbs" constituting the integer. this is ;; handled in alloc-gmpz below. (defun alloc-gmpz() ;; provides an empty gmp object for us to use, init to zero (let ((inside (make-array 3 :element-type '(signed-byte 32) :initial-element 0))) (mpz_init inside) ;;this next line sets up lisp to keep track of this object. ;;it sets up gc for what to do later, in case this object ;;is not accessible: the entries in the array "inside" will go away. ;;the array "inside" as well as the gmpz header will be gc'd later. (excl:schedule-finalization inside #'mpz_clear) (make-gmpz :z inside))) (defun alloc-gmpz2(n) ;; same as alloc-gmpz but initialize to n bits long (let ((inside (make-array 3 :element-type '(signed-byte 32) :initial-element 0))) (mpz_init2 inside n) (excl:schedule-finalization inside #'mpz_clear) (make-gmpz :z inside))) ;; note that mpz_clear is called only when a gc g has found that some ;; given gmpz object z is inaccessible from lisp. at the end of that ;; gc, the finalization program is run on z, so the guts are returned ;; to the gmp storage allocator. at the next gc after g, the object z ;; goes away as well. ;; create space for a number and put zero in it. (defun create_mpz_zero() (alloc-gmpz)) (defun copy_mpz (g) ; g is an inside.. ;new for simplegmp (let* ((r (alloc-gmpz)) (inside (gmpz-z r))) (mpz_add inside inside g) r)) (defmethod create_mpz_from_string((s string)) ;; s is a string like "123" (let* ((r (alloc-gmpz)) (inside (gmpz-z r))) (mpz_set_str inside s 10) ;; base 10 number conversion r)) (defmethod create_mpz_from_string((s integer)) ;; s isn't a string... (lisp2gmpz s)) (defun cs(x)(create_mpz_from_string x)) ;short name for above ;; given mpz number 123, return "123" (defmethod create_string_from_mpz((m gmpz)) (mpz_get_str 0 10 (gmpz-z m)) ) ;;; convert gmp number to a lisp number ;;; this is an inefficient hack. How often is it used though? (defmethod gmpz2lisp((x gmpz))(parse-integer (create_string_from_mpz x))) ;;; convert lisp string or number to gmpz. Bignums treated later. (defun into(x)(lisp2gmpz x)) (defmethod lisp2gmpz((x gmpz)) x) (defmethod lisp2gmpz((x string));; we hope, a string of decimal digits (let* ((a (alloc-gmpz)) (r (gmpz-z a))) (mpz_set_str r x 10) a)) (defmethod lisp2gmpz ((x fixnum)) (let* ((a (alloc-gmpz)) (r (gmpz-z a))) (mpz_set_si r x) a)) ;;; the program below is a slow hack, which sort of writes and then ;;; re-reads the number into a lisp number. Decomposing a lisp number ;;; word by word and clumping these together to make a GMPZ is ;;; done below in a more elaborate way using btoa. ;; we could convert floats which happen to be integers, too. e.g. (defmethod lisp2gmpz ((x float)) (lisp2gmpz (round x))) ;; round trip test, returns yes if successful ;;(defun rtgmpzlisp (x) (if (= x (gmpz2lisp(lisp2gmpz x))) 'yes 'no)) (defmethod create_double_from_mpz((m gmpz)) ;; creates a double precision version of mpz (mpz_get_d (gmpz-z m))) ;;(defun cd(m)(create_double_from_mpz m)) ; short version of above ;; rational division is not possible with ring of integers. ;;(defarithmetic / mpz_fdiv_qr) ; not yet linked up : returns quotient and remainder both ;; 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. (eval-when (compile load)(require :llstructs)) (defun nth-bigit (x n) (sys:memref x -14 (* 2 n) :unsigned-word)) (defun bignum-size(x) ;; in 16-bit quantities (excl::bm_size x)) (defmethod lisp2gmpz ((x integer)) ;; but not fixnum (let* ((out (create_mpz_zero)) ;the wrapper (r (gmpz-z out)) ;the inside (asize (bignum-size x)) (out2(create_mpz_zero)) (h (gmpz-z out2))) (do ((i (1- asize)(1- i))) ((< i 0) r) (declare (fixnum i)) (mpz_set_si h (nth-bigit x i)) (mpz_mul_2exp r r 16) (mpz_add r r h)) (if (< x 0) (setf (aref r 1)(- (aref r 1)))) out))