;; Using NTL from lisp

;;(load "c:/lisp/ntl/debug/ntl.dll")  ;; not really there.

;; mostly we are interested in NTL ZZX, polynomials in one
;; variable with coefficients that are arbitrary precision integers.
;; and also ZZXP,  coefficients mod p.

;; we may also need pairs for some purposes.

(eval-when (compile read eval load)(require :llstructs))
  ;;bignum or fixnum to array of unsigned-byte-32 bigits.


(defparameter buffer32  
    (make-array 1000 :element-type '(unsigned-byte 32) :initial-element 0
		:allocation :static-reclaimable))

(defparameter buffer16  
    (make-array 2000 :element-type '(unsigned-byte 16) :initial-element 0
		:allocation 
		:static-reclaimable))


(defun nth-bigit (x n) (sys:memref x -14 (* 2 n) :unsigned-word))
(defun bignum-size(x) ;; number of 16-bit quantities "bigits"
  (excl::bm_size x))

;; btoa16: input is a signed bignum.
;; result is 3 values. An array of 16bit words, the number of words,  a separate sign.
;; btoa16 works for any bignum or fixnum.  The elements in the array
;; represent the MAGNITUDE of the number, in 16 bit chunks.
;; the array a[0] is the least-significant 16 bits, etc.
;; atob16 converts the array, size,  sign, back to the bignum
;; btoa16 could be much much smaller, but we want to be efficient.

(defun btoa16(x)				
  ;;lisp bignum or fixnum to array of unsigned-byte-16 bigits.
  (declare (special buffer16))
  (let ((sign 1))
    (assert (numberp x))
    (if (< x 0)(setf sign -1 x (- x)))
    ;;fixnum could be up to 2^29-1 in this lisp.
    ;;I suppose this should be "portabilized" for other fixnum sizes.
    (if (fixnump x) (if (< x #.(expt 2 16))
			(values (make-array 1 :element-type '(unsigned-byte 16) 
					    :initial-element x)
				1
				sign)
		      (values (make-array 2 :element-type '(unsigned-byte 16) 
					  :initial-contents
					  (list 
					   (logand x #.(1-(expt 2 16)))
					   (ash x -16)))
			      2
			      sign))
      ;; we know it is really a lisp bignum.			
      (let ((size (bignum-size x))
	    ;; we re-use buffer16
	    ;; Later, grow buffer16 if needed.
	    (a16 buffer16))
	(declare (fixnum size))
	;; grow the buffer16 if needed
	(unless (< size (length buffer16))
	  (setf buffer16
		(make-array size :element-type '(unsigned-byte 16)
			    :initial-element 0 :allocation  :static-reclaimable))
	  (setf a16 buffer16))
	
	;; this loop does all the work
	(do ((i 0 (1+ i)));; 16 bits at a time
	    ((>= i size) (values a16 size sign))
	  (declare (fixnum i))
	  (setf (aref a16 i)
		(nth-bigit x i)))))))
  
 #+ignore		;; just for testing
(defun atob16(a size sign)	;; declared array of 16-bit numbers to bignum
  (assert (typep a '(simple-array (unsigned-byte 16))))
  (let ((num 0))
    (do ((i (1- size)  (1- i)))
	((< i  0) (if (= sign 1) num (- num))) ;tack on sign
      (setf num (+ (aref a i) (ash num 16))))))

;; ZZX is polynomial in one (unnamed) variable, e.g. x.
;; Make a ZZX from array of long coefficients (unsigned-byte 32)
;; Not completely general: longer coefficients not allowed here.

(ff:def-foreign-call 
 (LMP "LMakePoly")
 ((coefs (* :long)) (size :long))
 :returning  :foreign-address)

;; Make an array of coefs from a ZZX.

(ff:def-foreign-call
 (LUMP "LUnMakePoly")
 ((coefs (* :long))(poly  :foreign-address))
 :returning  :long)

(ff:def-foreign-call
 (LSetCoeffZ "LSetCoeffZ") ;works
 ((poly  :foreign-address) (power :long) (val :foreign-address))
 :returning  :void)

(ff:def-foreign-call ;works
 (LSetCoeffL "LSetCoeffL")
 ((poly  :foreign-address) (power :long) (val :long))
 :returning  :void)

(ff:def-foreign-call
 (LGetCoeffZ "LGetCoeffZ")
 ((poly  :foreign-address) (power :long))
 :returning  :foreign-address)


;; return ZZX product of two ZZX

(ff:def-foreign-call 
 (LMul2 "LMul2")
 ((poly1  :foreign-address)
  (poly2  :foreign-address))
 :returning  :foreign-address)

;; clobber ZZX target with the product of two ZZX

(ff:def-foreign-call 
 (LMul3 "LMul3")
 ((target :foreign-address)
  (poly1  :foreign-address)
  (poly2  :foreign-address))
 :returning  :void)

;; ZZ is arbitrary precision integer

;; return ZZ product of two ZZ
(ff:def-foreign-call 
    (LMulZZ2 "LMulZZ2")
    ((zz1  :foreign-address)
     (zz2  :foreign-address))
  :returning  :foreign-address)
;; clobber ZZ target with the product of two ZZ
(ff:def-foreign-call 
    (LMulZZ3 "LMulZZ3")
    ((target :foreign-address)
     (zz1  :foreign-address)
     (zz2  :foreign-address))
  :returning  :foreign-address)

;; similarly for addition.

(ff:def-foreign-call 
    (LAddZZ2 "LAddZZ2")
    ((zz1  :foreign-address)
     (zz2  :foreign-address))
  :returning  :foreign-address)

(ff:def-foreign-call 
    (LAddZZ3 "LAddZZ3")
    ((target :foreign-address)
     (zz1  :foreign-address)
     (zz2  :foreign-address))
  :returning  :foreign-address)

(ff:def-foreign-call 
    (LSizeZZX "LSizePoly")
    ((azz :foreign-address))
  :returning :long)

;; Convert a ZZ into an array of 32-bit unsigned bigits. Preserves the ZZ.
(ff:def-foreign-call 
    (LUMZS "LUnMakeZZSafe")
    ((bigits (* :long)) (azz :foreign-address))
  :returning :long)

;; Convert a long into a ZZ.
(ff:def-foreign-call (LMakeZZL "LMakeZZL")((z :long))
  :returning :foreign-address)

;; Convert a ZZ into a long.
(ff:def-foreign-call (LUnMakeZZL "LUnMakeZZL")((z :foreign-address))
  :returning :long)

(ff:def-foreign-call 
    (LNegateZZ "LNegateZZ") ;;z1 := -z2
    ((z1 :foreign-address)
     (z2 :foreign-address))
  :returning :void)


(ff:def-foreign-call  ;; z1 := z2 shifted left amt bits.
    (LeftShiftZZ3 "LeftShiftZZ3")
    ((z1 :foreign-address)
     (z2 :foreign-address)
     (amt :long))
  :returning :void)


(defun makentl16(arr size sign)
  ;; works on all inputs. Using 32-bit version fails on margins
  ;(assert (typep arr '(simple (unsigned-byte 16) (*))))
  (let ((sum (LMakeZZL 0)))
	
    (do ((i (1- size)(1- i)))
	((< i 0) nil)
      (declare (fixnum i))
      (LeftShiftZZ3 sum sum 16)
      (LAddZZ3 sum sum (LMakeZZL (aref arr i))))
    (unless (= sign 1)(LNegateZZ sum sum))
    sum))

#+ignore ;; round-trip test
(defun rt16(alispinteger)
 (NTLZZ2Lispint (Lispint2NTLZZ alispinteger)))

;; the next 2 functions' names are kind of long. Abbreviate them yourself.

(defun Lispint2NTLZZ(h)
  (assert (realp h) (h) "I can't convert ~s into an NTL ZZ" h)
  (multiple-value-bind (ar size sign)
      (btoa16 (round h)) ; works for floats,too.
    (makentl16 ar size sign)))

(defun NTLZZ2Lispint(p)			;p is ptr to NTL ZZ
    ;; in lisp, p looks like an integer, though it is a foreign address.
    (let ((size(LUMZS buffer32 p))
	  (sum 0))
    ;;  (format t "~%buffer= ~s" buffer32)
      (do ((i (1-(abs size)) (1- i)))
	  ((< i 0)(* (signum size) sum))
	(declare (fixnum i))
	(setf sum (+(aref buffer32 i) (* sum #.(expt 2 32)))))))

;; put together polynomials of arbitrary lisp bignum coefficients

(defun NTLZZ2LispintX(p)
  ;;p is ptr to NTL ZZ. Changes the NTL ZZ object to zero

  
    (let ((size(LUMZ buffer32 p))
	  (sum 0))
      (do ((i (1-(abs size)) (1- i)))
	  ((< i 0)(* (signum size) sum))
	(declare (fixnum i))
	(setf sum (+(aref buffer32 i) (* sum #.(expt 2 32)))))))

(defun lisp2NTLZZX(ar)
  (assert (typep ar 'array))
  (let ((poly (LMP (vector) 0)))	; an empty polynomial
    (do((i (1- (length ar)) (1- i)))
	((< i 0) poly)
      (LSetCoeffZ poly i (Lispint2NTLZZ (aref ar i) )))))

(defun NTLZZX2lisp(p)
  (let* ((size (LSizeZZX p))
	 (poly (make-array size)))
    
    (do((i (1- size) (1- i)))
	((< i 0) poly)
      (setf (aref poly i)
	(NTLZZ2Lispint(LGetCoeffZ p i))))))


;;;; factoring stuff

(ff:def-foreign-call
    LFactorPoly  ;;
    ((z1 :foreign-address) ;;  a ZZX
     (factors (* :long))
     (multiplicities  (* :long))
     (content :foreign-address))
   :returning :long)

;; array of multiplicities muls,
;; array of ZZXs factors

(defun pfact(apoly)
  (let* ((ntlpoly(lisp2NTLZZX apoly))
	 (maxdeg (length apoly))
	 (muls (make-array maxdeg :element-type '(signed-byte 32)
			 :initial-element 0
			 :allocation :static-reclaimable))
	 (facts (make-array maxdeg :element-type '(signed-byte 32)
			 :initial-element 0
			 :allocation :static-reclaimable))
	 (c (Lispint2NTLZZ 0)) ;the content
	 (thesize (LFactorPoly ntlpoly facts muls c)))

    (list 'content= (NTLZZ2Lispint c)
	  'factors=
	  (map 'list #'NTLZZX2lisp  ;; the list of factors
		     
		     (subseq facts 0 thesize))
	  'multiplicities=
	  (subseq muls 0 thesize);; the list of multiplicities
	  )))

#|
;; miscellaneous setup, testing.
(load "h:\\lisp\\ntldll\\debug\\ntldll.dll")
(defun setup() (load #P"ntldll\\debug\\ntldll.dll")(load "h:/lisp/loadntl.cl"))
(defun unset() (ff:unload-foreign-library "ntldll/debug/ntldll.dll"))
(defun unset() (ff:unload-foreign-library "thedll.dll"))
(defun setup() (load "thedll.dll")(load "h:/lisp/loadntl.cl"))

;; :ld f:/WinNTL-5.4/thedll/Debug/thedll.dll
;; these tests below all work fine 5/20/05
;; polynomial of longs
(setf a (make-array 10))
(dotimes (i 10 a)(setf (aref a i) (* (expt -1 i)i)))
(setf poly1 (lisp2NTLZZX a))
(dotimes (i 10 a)(setf (aref a i) (* 2 i)))
(setf poly2 (lisp2NTLZZX a))
(NTLZZX2lisp poly1)
(NTLZZX2lisp poly2)

;; demo how you can set and retrieve ZZ coefficients
;; in a polynomial
(LSetCoeffZ poly1 3 (Lispint2NTLZZ (expt 10 20)))
(NTLZZ2Lispint(LGetCoeffZ poly1 3))
;; show big numbers in a list
(setf v (make-array 30))
(dotimes (i 30)(setf (aref v i)(expt -4 i)))
(setf vvv (lisp2NTLZZX v))
(setf retv(NTLZZX2lisp vvv))
(every #'= v retv) ;; should be t
;; try squaring vvv:
(LMul3 vvv vvv vvv)
(NTLZZX2lisp vvv)

;;; try doing some factoring.
;; here is (2*x-4)*(3*x^2-1)^2  ...
;;                      5       4       3       2
;;(D2) 		    18 x  - 36 x  - 12 x  + 24 x  + 2 x - 4

(pfact (vector -4 2 24 -12 -36 18))






Further comments:
To do this "right" in lisp one should presumably declare
lisp objects which are of type NTL_ZZ, NTL_ZZX  etc, and
they should be suitably tagged in Common Lisp Object System (CLOS).
Each of these objects would then have suitable
conversion etc, so that one could use the same operator,
say mul to multiply whatever you had in mind.
Printing, parsing etc, could then be done more uniformly.
Calling the wrong function, e.g. LMulZZ3 instead of LMul3
will generally do something like segfault.

Drudge work to do, if we want to use more NTL resources,
just write the definitions following the patterns in LispWins.cpp.
For programs in NTL that return pairs or vectors or ...
some thought would have to be expended on how to represent
them in Lisp. Lists? Array?  Also how to convey the answers
between NTL and Lisp. e.g. Simplest primitives would probably
be: find the length L of a vector V, extract items
from V[i], i=0..(L-1).

To make polynomials whose coefficients are
NTL ZZ,  we create a zero polynomial, then by inserting
the high-order term it is filled out. We then continue to
use set/get coeff to fill in the rest.

Garbage collection could also be associated
with NTL objects, especially in CLOS; follow model used in
GMP. When the Lisp GC thinks there are no references to an NTL
object, make NTL erase it.

pfact sets up a call to polynomial factoring. Look at the
code in LispFact.cpp also.


|#


;; a swinnerton-dyer polynomial 
;;x^16 - 136*x^14 + 6476*x^12 - 141912*x^10 + 1513334*x^8 - 7453176*x^6 + 13950764*x^4 - 5596840*x^2 + 46225


;;(reverse (vector 1 0 -136 0 6476 0  -141912 0 1513334 0 -7453176 0 13950764		 0  -5596840 0  6225))