;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;;   LISP program for computing Fourier coefficients of Siegel Eisenstein
;;;   series using Katsurada's formula.
;;;   Copyright (C) 2001, 2016 Oliver D. King, with additions by
;;;   David S. Yuen and Cris Poor 
;;;
;;;   This program is free software: you can redistribute it and/or modify
;;;   it under the terms of the GNU General Public License as published by
;;;   the Free Software Foundation, either version 3 of the License, or
;;;   (at your option) any later version.
;;;
;;;   This program is distributed in the hope that it will be useful,
;;;   but WITHOUT ANY WARRANTY; without even the implied warranty of
;;;   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;;;   GNU General Public License for more details.
;;;
;;;   You should have received a copy of the GNU General Public License
;;;   along with this program.  If not, see <http://www.gnu.org/licenses/>.
;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;; This code implements methods for computing Fourier Coeffeicients of Siegel
;;; Eisentein series using the approach described in [Katsurada]. Some details
;;; of the implementation are discussed in [King]. The notation in the code
;;; generally follows [Katsurada], with references to other sources where relevant.
;;;
;;; References:
;;; [Katsurada] H. Katsurada, An explicit formula for Siegel series,  
;;;   American Journal of Mathematics, Vol. 121, No. 2, pp. 415-452, 1999
;;; [King] O.D. King, A mass formula for unimodular lattices with no roots,
;;;   Mathematics of Computation, Vol. 72, No. 242, pp. 839-863, 2003
;;; [Apostol] T.M. Apostol, Introduction to analytic number theory,
;;;   Springer-Verlag, New York, 1976
;;; [Serre] J.-P. Serre, A course in arithmatic, Springer-Verlag,
;;;   New York, 1973
;;; [SPLAG] J.H. Conway & N.J.A. Sloane,  Sphere packings, lattices and
;;;   groups (second edition), Springer-Verlag, New York, 1993
;;; [C&S] J.H. Conway & N.J.A. Sloane, Low-dimensional lattices. IV.
;;;   The mass formula, Proceedings of the Royal Society of London. Series A,
;;;   Mathematical and Physical Sciences, Vol. 419, No. 1857, pp. 259-286, 1988
;;; [Watson] G.L. Watson, Integral quadratic forms,
;;;      Cambridge University Press, 1960
;;;
;;;--------------------------- this is 80 characters wide ---------------------x

 
;;; (load "precomputed_misc.cl")  ; load list of primes up to 100000
;;; (load "precomputed_zetds.cl") ; load precomputed values of zeta_D(s)*d^(s-1)*sqrt(d)/pi^s
;;; Omit extension so preference is given to loading compiled code, if available.
;;; This behavior may be implementation-dependent, however; for clisp it appears
;;; that .fas (comiled) is given priority over .cl (source), unless it is older.

(load "precomputed_misc")    ; load list of primes up to 100000
(load "precomputed_zetads")  ; load precomputed values of zeta_D(s)*d^(s-1)*sqrt(d)/pi^s
 
;;; If rational input x = (p/q)^2, this finds p/q. (Some Common Lisp implementations
;;; will return rational value with ordinary sqrt, but this may be implementation dependent.)
;;; Note that numerator and denominator functions first put x in canonical reduced form.
(defun rationalsqrt (x) 
  (let* ((num (numerator x)) (den (denominator x)) 
         (facn (factorize num)) (facd (factorize den)) (ret 1))
    (dolist (term facn)
      (if term 
          (if (< 1 (car term)) 
                 (if (evenp (cadr term))
                     (setf ret (* ret (expt (car term) (/ (cadr term) 2)))) 
                   (print "sqrt not rational"))))) ; return nil?
    (dolist (term facd)
      (if term 
          (if (< 1 (car term)) 
                 (if (evenp (cadr term))
                     (setf ret (/ ret (expt (car term) (/ (cadr term) 2)))) 
                   (print "sqrt not rational")))))  ; return nil?
    (if (zerop num) (setf ret 0))
    ret))
  

;;; (defparameter +listofprimes+ '(2 3 5 7 11 13 17 19 23 29 31 37 41 43 47))

;;; Makes list of all primes less than or equal to n.
;;; Currently precomputed (for = 100000) and loaded from precomputed_misc
(defun makeprimelist (n)
  (let ((temp ()) (rootn (floor (sqrt n)))
        (sieve (make-array (+ 1 n) :initial-element 1)))
    (declare (special +listofprimes+))
    (setf (aref sieve 0) 0 (aref sieve 1) 0)
    (dotimes (i (+ 1 rootn))
      (if (= 1 (aref sieve i))
          (dotimes (j (- (floor n i) 1))
            (setf (aref sieve (* i (+ j 2))) 0))))
    (dotimes (i n) (if (= 1 (aref sieve i)) (setf temp (append temp (list i)))))
  (setf +listofprimes+ temp)))

;;; Prime factorization of non-negative integer n. 
;;; Return list of pairs ((p_1 m_1) (p_2 m_2) ...) in prime-power decomposition 
;;; n = p_1^m1 p2_m2 ..., or ((1 1)) for n = 1 or ((0 1)) for n = 0. 
;;; This is *not* an efficient general purpose algorithm, but was intended
;;; for factoring fairly small determinants of lattices.
(defun factorize (n)
  (let ((ret ()) (i 0))
    (declare (special +listofprimes+))
    (if (zerop n) (setf ret (list (list 0 1))) 
      (if (eq 1 n) (setf ret (list (list 1 1)))
        (do ((ls +listofprimes+ (cdr ls))) 
            ((or (null ls) (> (car ls) n)))
          (setf i (ord (car ls) n))
          (if (plusp i) 
              (setf ret (append ret (list (list (car ls) i))) 
                n (/ n (expt (car ls) i)))))))
    (if (> n 1) 
        (progn 
          (format t "~%going beyond primelist to factor ~d, may take awhile" n)
          (do ((j (car (last +listofprimes+)) (+ 2 j))) 
              ((= 1 n))
            (setf i (ord j n))
            (if (plusp i) 
                (setf ret (append ret (list (list j i))) 
                  n (/ n (expt j i)))))))   
    ret))


;; (defun primep n (member n +listofprimes+))

;;; Tests whether positive integer n is prime, returns t if yes, nil otherwise. 
;;; This is *not* an efficient general purpose primality testing algorithm.
(defun primep (n)
  (let ((facn (factorize n)))
    (and (= 1 (length facn)) (< 1 (car (car facn))) (= 1 (cadr (car facn))))))

;;; Factor b, compute product of Legendre symbols
(defun jacobisymbol (a b)
  (let ((ret 1) (facb (factorize b))) 
    (if (= b 1) (setf ret 1)
      (dolist (term facb)
        (if term 
            (setf ret (* ret (expt (legendresymbol a (car term)) (cadr term)))))))
    ret))

;;; Precomputed Bernoulli numbers for n=0 to 11
(defparameter +bernoulli+ '#(1 -1/2 1/6 0 -1/30 0 1/42 0 -1/30 0 5/66 0))

;;; Bernoulli number, [Apostol] p265. (In [Serre] VII.4.1 this is (-1)^(n+1) b_2n)
;;; This version avoids recomputing numbers
(defun bernoulli2 (n)
  (let ((temp) (m) (ret)) 
    (cond
     ((< n (length +bernoulli+)) (setf ret (aref +bernoulli+ n)))
     (t (progn 
          (setf temp (make-array (+ 1 n) :initial-element 0))
          (dotimes (j (length +bernoulli+)) 
            (setf (aref temp j) (aref +bernoulli+ j))) 
          (dotimes (j (- n -1 (length +bernoulli+)))
            (setf m (+ (length +bernoulli+) j))
            (dotimes (k m)
              (setf (aref temp m) 
                (- (aref temp m) (* (choose (+ 1 m) k) (aref temp k)))))
            (setf (aref temp m) (/ (aref temp m) (+ m 1))))
          (setf ret (aref temp n))
          (setf +bernoulli+ temp))))
    ret))


;;; Bernoulli number, [Apostol] p265. (In [Serre] VII.4.1 this is (-1)^(n+1) b_2n)
(defun bernoulli (n)
  (let ((temp) (m) (ret))
    (cond
     ((= n 0) (setf ret 1))
     ((= n 1) (setf ret -1/2))
     ((> n 1) (progn 
                (setf temp (make-array (+ 1 n) :initial-element 0))
                (setf (aref temp 0) 1 (aref temp 1) -1/2) 
                (dotimes (j (- n 1))
                  (setf m (+ 2 j))
                  (dotimes (k m)
                    (setf (aref temp m) 
                      (- (aref temp m) (* (choose (+ 1 m) k) (aref temp k)))))
                  (setf (aref temp m) (/ (aref temp m) (+ m 1))))
                (setf ret (aref temp n)))))
    ret))


;;; Bernoulli polynomials
(defun bernoullipoly (n x)
  (let ((ret 0)) 
    (dotimes (k (+ n 1)) 
      (setf ret (+ ret (* (choose n k) (bernoulli k) (expt x (- n k))))))
    ret))


;;; Bernoulli polynomials
(defun bernoullipoly2 (n x)
  (let ((ret 0) (k))
    (dotimes (j (+ n 1))
      (setf k (- n j))
      (setf ret (+ ret (* (choose n k) (bernoulli2 k) (expt x (- n k))))))
    ret))


;;; Minkowski Siegel mass of even unimodular lattices of dim n, [SPLAG] p409
(defun typeIImass (n)
  (let ((ret) (k (/ n 2)))
    (if (plusp (mod n 8)) (setf ret 0)
      (if (zerop n) (setf ret 1) 
        (progn
          (setf ret (/ (bernoulli2 k) n))
          (dotimes (j (- k 1)) 
            (setf ret (* ret (/ (bernoulli2 (* 2 (+ j 1))) (* 4 (+ j 1)))))))))
    (abs ret)))


;;; Minkowski Siegel mass of odd unimodular lattices of dim n, [SPLAG] p409
(defun typeImass (n)
  (let ((ret 1) (res (mod n 8)) (k (if (evenp n) (/ n 2) (/ (- n 1) 2))))
    (dotimes (j (- k 1)) (setf ret (* ret (bernoulli2 (* 2 (+ j 1))))))
    (if (oddp n) (setf ret (* ret (bernoulli2 (* 2 k)))))
    (cond ((= n 1) (setf ret 1/2))
          ((= res 0)  
           (progn 
             (setf ret (* ret (- 1 (expt 2 (- k))) (+ 1 (expt 2 (- 1 k))) (bernoulli2 k)))
             (setf ret (/ ret (* 2  (factorial k))))))
          ((or (= res 1) (= res 7))  
           (progn 
             (setf ret (* ret (+ (expt 2 k) 1)))
             (setf ret (/ ret (* (expt 2 (+ 1 (* 2 k))) (factorial k))))))
          ((or (= res 2) (= res 6)) 
           (progn 
             (setf ret (* ret (euler (- k 1))))
             (setf ret (/ ret (* (expt 2 (+ 1 (* 2 k))) (factorial (- k 1)))))))
          ((or (= res 3) (= res 5)) 
           (progn 
             (setf ret (* ret (- (expt 2 k) 1)))
             (setf ret (/ ret (* (expt 2 (+ 1 (* 2 k))) (factorial k))))))
          ((= res 4) 
           (progn 
             (setf ret (* ret (- 1 (expt 2 (- k))) (- 1 (expt 2 (- 1 k))) (bernoulli2 k)))
             (setf ret (/ ret (* 2  (factorial k)))))))
    (abs ret)))

;;; Precomputed absolute values of Euler numbers for indices 0,2,4,...,40; 
;;; (Odd indices have value zero).
;;; From Sloane's On-Line Encyclopedia of Integer Sequences (https://oeis.org/A000364)
;;; except for last 3, computed using euler2 function below. 
(defparameter +eulernumbers+  
  '#(1 1 5 61 1385 50521 2702765 199360981 19391512145 
     2404879675441 370371188237525 69348874393137901
     15514534163557086905 4087072509293123892361 
     1252259641403629865468285 441543893249023104553682821 
     177519391579539289436664789665 
     80723299235887898062168247453281
     41222060339517702122347079671259045
     23489580527043108252017828576198947741
     14851150718114980017877156781405826684425))


;;; Absolute value of Euler number, using precomputed value for n up to 40.
(defun euler (n) 
  (if (oddp n) 0 
    (if (< (/ n 2) (length +eulernumbers+)) (aref +eulernumbers+ (/ n 2))
      (abs (euler2 n)))))

;;; Euler numbers
(defun euler2 (n) 
  (let ((ret) (k) (vals))
    (if (oddp n) (setf ret 0)
      (progn (setf k (/ n 2))
        (setf vals (make-array (+ 1 k)))
        (dotimes (i (+ 1 k)) (setf ret (if (zerop i) -1 0))
          (dotimes (j i) 
            (setf ret (+ ret (* (choose (* 2 i) (* 2 j)) (aref vals j)))))
          (setf (aref vals i) (- ret)))
        (setf ret (aref vals k))))
    ret))


;;; Prints output of factorize in the format p_1^a_1 p_2^a_2 ... 
;;; to specified outstream (t for standard out)
(defun printfactors (facs outstrm)
  (dolist (term facs) 
    (if term 
        (cond ((= 0 (car term)) (format outstrm "0 "))
              ((= 1 (car term)) (format outstrm "1 "))
              ((= 1 (cadr term)) (format outstrm "~d " (car term)))
              (t (format outstrm "~d^~d " (car term) (cadr term))))))
  (values))


;;; Riemann zeta function for positive even n, omits the factor pi^n
(defun zeta (n)
  (* (expt -1 (+ 1 (/ n 2))) (expt 2 n) (bernoulli2 n) (/ 1 (* 2 (factorial n)))))


;;; Precomputed zeta(n), even n, 0 to 32, omits the factor pi^n; ignore odd n
(defparameter +zeta+
  '#(-1/2 0 1/6 0 1/90 0 1/945 0 1/9450 0 1/93555 0 691/638512875 0 
     2/18243225 0 3617/325641566250 0 43867/38979295480125 0 
     174611/1531329465290625 0 155366/13447856940643125 0 
     236364091/201919571963756521875 0
     1315862/11094481976030578125 0 
     6785560294/564653660170076273671875 0                  
     6892673020804/5660878804669082674070015625 0
     7709321041217/62490220571022341207266406250))

;;; Computes n! for non-negative integer n
(defun factorial (n)
  (declare (fixnum i n))
  (let ((ret 1)) 
    (dotimes (i n) (setf ret (* ret (+ 1 i))))
    ret))


;;; Binomial coefficient C(n,k)
(defun choose2 (n k)
  (declare (fixnum i n k))
  (let ((ret 1)) 
    (dotimes (i k) (setf ret (* ret (/ (- n i) (- k i)))))
    ret))


;;; Binomial coefficient C(n,k) using 1 rather than k divisions.
(defun choose (n k)
  (declare (fixnum i n k))
  (if (> k n) 0
    (let ((num 1) (denom 1)) 
      (dotimes (i k) (setf num (* num (- n i)) denom (* denom (- k i))))
      (/ num denom))))


;;; Dirichlet character, [C&S] p267
(defun dirichar (bigd m) (if (evenp m) 0 (jacobisymbol bigd m)))

;;; Modulus of Dirichlet character, C&S p267
(defun findmodulus (bigd)
  (let* ((done nil) (absd (abs bigd)) (vals (make-array (* 4 absd))) ;
         (ret (* 4 absd)))
    (dotimes (i (* 4 absd)) (setf (aref vals i) (dirichar bigd i))) 
    ;; (print vals)
    (do ((i 1 (+ 1 i)))
      ((or done (= i (* 4 absd))))
      (if (zerop (mod (* 4 absd) i))
          (do ((j (+ 1 i) (+ 1 j)))
              ((or (= j (* 4 absd))
                   (not (= (aref vals (mod j i)) (aref vals j))))
               (if (= j (* 4 absd)) (setf done t ret i))))))
    ret))


;;; Modulus and conductor of Dirichlet character (as list), [C&S] p267
(defun findconductor (bigd)
  (let* ((done nil) (modu (findmodulus bigd))  (ret modu)) ; 
    (do ((i 1 (+ 1 i)))
        ((or done (= i modu)))
      (if (zerop (mod modu i))
          (do ((j 1 (+ 1 j)))
              ((or (= j modu)
                   (not
                    (or
                     (> (gcd modu j) 1)
                     (plusp (mod (- j 1) i))
                     (= 1 (dirichar bigd j)))))
               (if (= j modu) (setf done t ret i))))))
    (list modu ret)))

;;; Primitive character, [C&S] p268
(defun primchar (bigd n)
  (let* ((mc (findconductor bigd)) (k (car mc)) (d (cadr mc)) (ret)) ;
    (if (< 1 (gcd d n)) (setf ret 0)
      (do ((q 1 (+ 1 q))) 
          ((and (primep (+ n (* q d))) (plusp (mod k (+ n (* q d)))))
           (setf ret (dirichar bigd (+ n (* q d)))))))
    ret))

;;; Primitive character, [C&S] p268
(defun primcharvec (bigd)
  (let* ((mc (findconductor bigd)) (k (car mc)) (d (cadr mc)) ; 
         (ret (make-array d :initial-element 8)))
    (dotimes (n d)
    (if (< 1 (gcd d n)) (setf (aref ret n) 0)
      (do ((q 1 (+ 1 q))) 
          ((and (primep (+ n (* q d))) (plusp (mod k (+ n (* q d)))))
           (setf (aref ret n) (dirichar bigd (+ n (* q d))))))))
    ret))


;;; Computes L-series zeta_D(s)*d^(s-1)*sqrt(d)/pi^s, with D=d*(-1)^s
;;; as in [C&S] table 6. Change to hash values after they are computed?
(defun zetad (smalld ss)
  (let* ((s (if (zerop ss) 1 ss)) (bigd (* smalld (expt -1 s)))
         (redd (/ bigd (expt 4 (ord 4 bigd))))
         (mc (findconductor redd)) (k (car mc)) (d (cadr mc))
         (charvec (primcharvec redd)) (val) (temp 0)) ; 
    (setf val (/ (- (expt 2 s)) 
                 (* (factorial s) (+ 1 (* (aref charvec (- d 1)) (expt -1 s))))))
    (if (not (= k d))
        (progn
          (dotimes (p k) 
            (if (and (primep p) (zerop (mod k p)))
                (setf val (* val (- 1 (/ (aref charvec (mod p d)) (expt p s)))))))
          (if (primep k)
              (setf val (* val (- 1 (/ (aref charvec (mod k d)) (expt k s))))))))
    (dotimes (r d) (setf temp (+ temp (* (aref charvec r) 
                                       (bernoullipoly2 s (/ r d))))))
    (setf val (* val temp))
    (setf val (* val (rationalsqrt (/ smalld d))))
    (if (< 1 (mod s 4)) (setf val (- val)))
    (setf val (* val (expt smalld (- s 1))))
    val))

;;; Product of F_p(B; X) polynomials for p | 2*det(B) evaluated at X = p^{-s},
;;; as in [Katsurada] Theorem 4.4. Here the input B is a [SPLAG]-style genus symbol,
;;; in the format described in the file precomputed_misc.cl.
;;; If half = t, the calculation is done for the half-integral lattice whose double
;;; has the given genus symbol. If interpolate = t, Lagrangian interpolation is used
;;; as in [King] p858.
(defun fpproduct (genus s interpolate half)
  (declare (fixnum s))
  (let* ((n (length genus)) ; (dim (+ (aref genus 1) (aref genus 2))) 
         (deter 1) (primelist '()) (wat) (fpbx) (fac) (ret 1) (X)) ; 
    (declare (fixnum n i p))
    (dotimes (i (- n 3)) 
      (setf deter (* deter (expt (car (aref genus (+ i 3))) 
                                 (abs (cadr (aref genus (+ i 3))))))))
    (setf fac (factorize deter))
    (dolist (term fac) 
      (if (and term (< 2 (car term))) 
          (setf primelist (append primelist (list (car term))))))
    (if (oddp (aref genus 2)) (setf deter (- deter)))
    (dolist (p primelist) 
      (setf X (expt p (- s)))
      (setf wat (genus2watsonp p genus))
      (setf fpbx (fpiter p wat X interpolate))
      (setf ret (* ret fpbx)))
    (setf X (expt 2 (- s))) 
    (if half (setf wat (halfwatson (genus2watsontwo genus)))
      (setf wat (genus2watsontwo genus)))
    (setf fpbx (f2iter wat X interpolate))
    (setf ret (* ret fpbx))
    ret))

;;; Computes c(n,k,B) from [Katsurada] p416. Here the input B is a [SPLAG]-style 
;;; genus symbol,in the format described in the file precomputed_misc.cl.
;;; If half = t, the calculation is done for the half-integral lattice whose double
;;; has the given genus symbol. If interpolate = t, Lagrangian interpolation is used
;;; as in [King] p858.
;;; NOTE: This assumes that B is positive definite with rank n, and k is even, with k > n. 
;;; Returns rational number, power of pi, and determinant of B (not half B).
;;; example: degree 6, weight 16
;;;     (cnk 6 16 #(2 6 0 (3  -1)  (4  -2  1  4)) t)
(defun cnk (n k genus half)
  (declare (fixnum n k))
  (let ((deter 1) (i) (powerofpi 0) (temp) (d) (facd) (s) (sq 1) (deter2) 
        (fpprod (if (plusp n) (fpproduct genus k (>= n k) t) 1)) ;; use interpolation when n>=k
        (ret (* (if (evenp (/ (* n k) 2)) 1 -1) 
                (expt 2 (* n (- k (/ (1- n) 2))))))) ; 
    (declare (fixnum i powerofpi i j s))
    (declare (special +zetadshash+))
    (dotimes (i (- (length genus) 3)) ; CHECK: should fixnum be declared in the macro?
      (setf deter (* deter (expt (car (aref genus (+ i 3))) 
                                 (abs (cadr (aref genus (+ i 3))))))))	
    (setf d deter) (setf deter2 deter)			 
    (if (oddp (aref genus 1)) (setf deter (- deter))) ; check
    ;; if n = 2*k and det not square, return 0 here? 
    (if half (setf deter (/ deter (expt 2 n))))
    ;; also require n < k here?
    (if (or (oddp k) (/= n (aref genus 1))) (setf ret nil) ; require even k and n = rank(genus).
      ;; else k even and n = rank(genus)
      (progn
         (setf ret (* ret fpprod))
         (setf ret (/ ret (zeta k)))
         (setf powerofpi (- powerofpi k))
         (dotimes (j (floor n 2))
           (setf ret (/ ret (zeta (* 2 (- k 1 j)))))
           (setf powerofpi (- powerofpi (* 2 (- k 1 j))))) 
         (dotimes (j n)
           (setf i (- (* 2 k) j))
           (setf powerofpi (+ powerofpi (floor i 2)))
           (setf ret (/ ret (gammahalf i))))
         (if (oddp n)
             (progn
               (setf ret (* ret (expt deter (* 1 (/ (- (* 2 k) n 1) 2)))))
               (if (= n (- (* 2 k) 1)) (setf ret (/ ret 2))))
           (progn        
             (setf facd (factorize d))
             (dolist (term facd)
               (if (and term (< 1 (car term)))
                   (setf sq (* sq (expt (car term) (floor (cadr term) 2))))))
             (setf d (/ d (* sq sq)))
             (setf s (- k (/ n 2)))
             (setf temp (* ;; (zetad d s)  
                           ;; use precomputed values if possible
                           (if (and (< s 17) (gethash d +zetadshash+))
                               (aref (gethash d +zetadshash+) s)
                             (zetad d s)) ;; CHECK: add this value to zetadshash?
                           (expt (/ deter d) (1- s))  ;; deter/d = sq*sq/2^n
                           (/ sq (if half (expt 2 (/ n 2)) 1))))
             (setf ret (* ret temp))
             (if (= n (* k 2)) (setf ret (- ret)))
             (setf temp (mod d 8))
             (if (plusp s)
                 (if (oddp (/ n 2))
                     (cond ((= temp 7) (setf ret (/ ret (- 1 (/ 1 (expt 2 s))))))
                           ((= temp 3) (setf ret (/ ret (+ 1 (/ 1 (expt 2 s)))))))
                   (cond ((= temp 1) (setf ret (/ ret (- 1 (/ 1 (expt 2 s))))))
                         ((= temp 5) (setf ret (/ ret (+ 1 (/ 1 (expt 2 s))))))))
               (if (/= 1 d) (setf ret 0)))))))
    ;; (format t "~d~5t~d~12t~d~30t" n  deter2  ret)  
    (values ret powerofpi deter2)))


;;; Solves mx=b for x, where m is a square Vandermonde matrix.
(defun gaussianelimination (m b)
  (let* ((n (length b)) (i) (j) (k) (temp))  
    (dotimes (i (- n 1)) 
      (dotimes (jj (- n i 1)) (setf j (+ jj i 1))
        (setf temp (/ (aref m j i) (aref m i i)))
        (dotimes (kk (- n i 1)) (setf k (+ kk i 1))
          (setf (aref m j k) (- (aref m j k) (* temp (aref m i k)))))
        (setf (aref b j) (- (aref b j) (* temp (aref b i))))))
    (dotimes (ii n) (setf i (- n ii 1))
      (dotimes (jj (- n i 1)) (setf j (+ jj i 1))
        (setf (aref b i) (- (aref b i) (* (aref b j) (aref m i j)))))
      (setf (aref b i) (/ (aref b i) (aref m i i))))
    b))

;;; Computes F_p(B; X) from [Katsurada] as symbolic polynomial for prime p. 
;;; Here wat is a p-adic diagonalization (or 2-adic decomposition) of B,
;;; as in [Watson], of dimension n. 
(defun fppoly (p wat n)
  (let* ((m (make-array (list n n))) (b (make-array n))) 
    (dotimes (i n) (dotimes (j n) (setf (aref m i j) (expt (+ 1 i) j))))
    (if (= p 2) (dotimes (i n) (setf (aref b i) (f2iter wat (+ i 1) t)))  
      (dotimes (i n)  (setf (aref b i) (fpiter p wat (+ i 1) t))))
    (gaussianelimination m b)))

;;; Computes F_p(B; X) from [Katsurada] as symbolic polynomials for all 
;;; relevant primes p. Here the input B is a [SPLAG]-style genus symbol,
;;; in the format described in the file precomputed_misc.cl.
;;; If half = t, the calculation is done for the half-integral lattice 
;;; whose double has the given genus symbol.
;;; This doesn't return anything; 
;;; it just prints out one polynomial per line in the format
;;; p #(c0 c1 c2 ... )
;;; where p is a prime dividing 2*det(B) and the ci 
;;; are coefficients in the polynomial c0 + c1*X + c2*X^2 + ...
(defun fppolys (genus half)
  (let* ((n (length genus))  (dim (+ (aref genus 1) (aref genus 2))) 
         (deter 1) (primelist '()) (wat) (fac))
    (declare (fixnum n dim i p)) ; check scope!
    (dotimes (i (- n 3)) 
      (setf deter (* deter (expt (car (aref genus (+ i 3))) 
                                 (abs (cadr (aref genus (+ i 3))))))))
    (setf fac (factorize deter))
    (dolist (term fac) 
      (if (and term (< 2 (car term))) 
          (setf primelist (append primelist (list (car term))))))
    (if (oddp (aref genus 1)) (setf deter (- deter))) ; not used? (aref genus 2)?
    (dolist (p primelist) 
      (setf wat (genus2watsonp p genus))
      (format t "~%~d  ~d" p (fppoly p wat (+ 7 dim))))  ; original + 1
    (if half (setf wat (halfwatson (genus2watsontwo genus)))
      (setf wat (genus2watsontwo genus)))
    (format t "~%~d  ~d" 2 (fppoly 2 wat (+ 7 dim)))     ; original + 1
    (values)))

;;; Similar to fppolys above, but does not print anything.
;;; Returns a list, each term of which is a pair (p #(c0 c1 ...)) 
;;; which can be printed in format {{p, c0+ c1*X + ...}, ...} 
;;; using function printpolys below
(defun fppolys2 (genus half)
  (let* ((n (length genus))  (dim (+ (aref genus 1) (aref genus 2))) 
         (deter 1) (primelist '()) (wat) (fac) (polylist '()))
    (declare (fixnum n dim i p)) ; check!
    (dotimes (i (- n 3)) 
      (setf deter (* deter (expt (car (aref genus (+ i 3))) 
                                 (abs (cadr (aref genus (+ i 3))))))))
    (setf fac (factorize deter))
    (dolist (term fac) 
      (if (and term (< 2 (car term))) 
          (setf primelist (append primelist (list (car term))))))
    (if (oddp (aref genus 1)) (setf deter (- deter))) ; not used? (aref genus 2)?
    (dolist (p primelist) 
      (setf wat (genus2watsonp p genus))
      (setf polylist (append polylist (list (list p (fppoly p wat (+ 7 dim)))))))  
    (if half (setf wat (halfwatson (genus2watsontwo genus)))
      (setf wat (genus2watsontwo genus)))
    (setf polylist (append polylist (list (list 2 (fppoly 2 wat (+ 7 dim))))))
    polylist))


;;; Evaluates the polynomial c0 + c1*x + c2*x^2 + ... at given x.
;;; Input argument coeffs is the vector of coefficients #(c0 c1 c2 ...)
(defun evalpoly (coeffs x)
  (let ((ret 0) (y 1))
    (dotimes (i (length coeffs)) 
      (if (not (zerop (aref coeffs i)))
          (setf ret (+ ret (* (aref coeffs i) y))))
      (setf y (* y x)))
    ret))


;;; Prints polynomial with coefs #(c0 c1 c2 ...) as c0 + c1*x + c2*x^2 ...
;;; to specified outstrm, with outstrm = t for standard out. 
;;; Terms with coefficient 0 are skipped, and - rather than + is used
;;; when the coefficient is negative.
(defun printpoly (outstrm coefs)
  (let ((leading t))
    (dotimes (i (length coefs))
      (if (not (zerop (aref coefs i)))
          (progn
            (if leading 
                (progn 
                  (setf leading nil)
                  (if (minusp (aref coefs i)) (format outstrm "-")))
              (if (<= 0 (aref coefs i)) 
                  (format outstrm "+") 
                (format outstrm "-")))
            (cond 
             ((= 0 i) (format outstrm "~d" (abs (aref coefs i))))
             ((= 1 i) (format outstrm "~dx" (abs (aref coefs i))))
             ((= 1 (abs (aref coefs i)))  (format outstrm "x^~d" i))
             (t (format outstrm "~dx^~d" (abs (aref coefs i)) i))))))
    (values)))


;;; Prints list of polynomials as output by function fppolys2 to outstrm, 
;;; in format {{p, c0 + c1*X + ...}, ...} 
;;; example: (printpolys t (fppolys2 #(2 6 0 (3  -1) (4  -2  1  4)) t))
(defun printpolys (outstrm polylist)
  (let ((firstp t))
   (format outstrm "{")
   (dolist (poly polylist)
      (if (not firstp) (format outstrm ","))
      (setf firstp nil)
      (format outstrm "{~d," (car poly))
      (printpoly outstrm (cadr poly))
      (format outstrm "}"))
   (format outstrm "}")
   (values)))


;;; Legendre symbol, p odd prime or p = 2 & a odd
(defun legendresymbol (a p) 
  (declare (fixnum p));; check size of a
  (cond ((> p 2)  
         (cond ((zerop (mod a p)) 0 )
               ((< 1 (mod (expt a (/ (1- p) 2)) p)) -1)
               (t 1)))
        ((= p 2) 
         (cond ((or (= 1 (mod a 8)) (= 7 (mod a 8))) 1)
               (t -1)))))

;;; Power of p dividing q, nonzero rational.
;;; Can get rid of one test since fractions are reduced.
(defun ord (p q)
  (declare (fixnum p)) 
  (let ((r 0) (c 0)) 
    (declare (fixnum r c))
    (setf c (numerator q))
    (if (zerop c) (progn (print "invalid: taking ord of 0") (setf c 1)))        ; return nil?
    (if (zerop (mod c p)) 
        (setf r (do ((i 0 (1+ i)) (j c (/ j p))) ((plusp (mod j p)) i))))
    (setf c (denominator q))
    (if (zerop (mod c p)) 
        (setf r (do ((i r (1- i)) (j c (/ j p))) ((plusp (mod j p)) i))))
    r))


;;; chi_p(a) from [Katsurada] p428, extended for rational a
(defun chip (p a) 
  (declare (fixnum p i))
  (let ((ret 0) (rn 0) (rd 0) (cn (numerator a)) (cd (denominator a)))
    (declare (fixnum ret rn rd cn cd))
    (if (zerop (mod cn p)) 
        (multiple-value-setq (rn cn) 
          (do ((i 0 (1+ i)) (j cn (/ j p))) 
              ((plusp (mod j p)) (values i j)))))
    (if (zerop (mod cd p)) 
        (multiple-value-setq (rd cd) 
          (do ((i 0 (1+ i)) (j cd (/ j p))) 
              ((plusp (mod j p)) (values i j)))))
    (setf rn (- rn rd)) 
    (setf cn (* cn cd))
    (cond ((and (> p 2) (evenp rn)) (setf ret (legendresymbol cn p)))
          ((and (= p 2) (evenp rn)) (cond ((= 1 (mod cn 8)) (setf ret 1)) 
                                          ((= 5 (mod cn 8)) (setf ret -1)))))
    ret))


;;; Hilbert symbol (a,b)_p, [Watson] p36
(defun hilbp (p a b) 
  (declare (fixnum p))
  (let ((pp (* p p)) 
        (r (* (numerator a) (denominator a)))                                   ; check denom prime to p?
        (s (* (numerator b) (denominator b)))
        (ret 0));; should end up 1 or -1 
    (declare (fixnum pp r s ret))
    (if (zerop (mod r pp)) 
        (setf r (do ((j r (/ j pp))) ((plusp (mod j pp)) j))))
    (if (zerop (mod s pp)) 
        (setf s (do ((j s (/ j pp))) ((plusp (mod j pp)) j))))
    (if (and (zerop (mod s p)) (zerop (mod r p))) (setq r (/ (* (- r) s) pp)))
    (if (zerop (mod r p)) (rotatef r s))
    (if (zerop (mod (* s r) p))
        (if (and (plusp (mod (* 2 r) p)) (zerop (mod s p))) 
            (setf ret (legendresymbol r p)) ;; else
          (if (and (= p 2) (zerop (mod s p)) 
                   (or (= 1 (mod r 8)) (= 1 (mod (+ r s) 8)))) 
              (setf ret 1) (setf ret -1))) ;; else p doesn't divide r*s
      (if (and (= 2 p) (= 1 (mod (* r s) 2)) ) 
          (if (or (= 1 (mod r 4)) (= 1 (mod s 4))) (setf ret 1) (setf ret -1))
        (if (and (> p 2) (plusp (mod (* r s ) p))) (setf ret 1))))              ; need last term?
    ret))   



;;; Vector of values delta_p as in [Katsurada] p430 for decomposition of B
(defun delbulk (p B) 
  (declare (fixnum p))
  (let* ((n (length B)) (invars (make-array n)) (i (- n 1)) (temp 1) (temp2))
    (declare (fixnum n i))
    (dolist (item (coerce (reverse B) 'list))
      (setf temp (* temp item))       
      ;; (setf temp2 (ord p (* temp (expt 2 (* 2 (floor (- n i) 2))))))         ; don't need 2^a for odd p 
      (setf temp2 (ord p temp))
      ;; (if (= 2 p) (setf temp2 (+ temp2 (* 2 (floor (- n i) 2)))))
      (if (evenp (- n i)) 
          (setf temp2 (* 2 (floor (+ temp2 1 (if (= 2 p) -1 0)) 2))))
      (setf (aref invars i) temp2) (setf i (1- i))) 
    invars))


;;; Vector of Hasse invariants h_p as in [Katsurada] p428 for decomposition of B
(defun hassebulk (p B)
  (declare (fixnum p))
  (let* ((n (length B)) (invars (make-array n)) (temp 1))
    (declare (fixnum n i j))
    (dotimes (i n)
      (dotimes (j (1+ i))
        (setf temp (* temp (hilbp p (aref B (- n i 1)) (aref B (- n j 1))))))
      (setf (aref invars (- n i 1)) temp))
    invars))


;;; Vector of values eta_p and xi_p as in [Katsurada] p428 for decomposition of B
(defun etaxibulk (p B)
  (declare (fixnum p))
  (let* ((n (length B)) (invars (make-array n)) (hassevec (hassebulk p B))
         (i (- n 1)) (temp 1))
    (declare (fixnum n i))
    (dolist (item (coerce (reverse B) 'list))
      (setf temp (* temp item)) (setf (aref invars i) temp) (setf i (1- i)))
    (dotimes (i n)
      (setf temp (aref invars i))
      (if (evenp (- n i)) (setf (aref invars i) 
                            (chip p (* (if (evenp (/ (- n i) 2)) 1 -1) temp)))
        (setf (aref invars i) 
          (* (aref hassevec i) 
             (hilbp p temp (* temp (if (evenp (/ (- n i 1) 2)) 1 -1)))
             (expt (hilbp p -1 -1) (/ (1- (* (- n i) (- n i))) 8))))))
    invars))


;;; Interpolates f(x) from array of vals f(1) f(p) f(p^2) ... f(p^deg) 
(defun interpolatep (p deg x vals)
  (declare (fixnum p deg))
  (let ((ppowers (make-array (1+ deg))) (temp 1) 
        (num 1) (ret 0) (ind -1) (diff))
    (declare (fixnum ind i j))
    (dotimes (i (1+ deg)) 
      (setf (aref ppowers i) temp)
      (setf diff (- x temp))
      (if (zerop diff) (setf ind i) (setf num (* num diff))) 
      (setf temp (* temp p)))
    (if (>= ind 0) (setf ret (aref vals ind))
      (dotimes (i (1+ deg))
        (setf temp 1)
        (dotimes (j (1+ deg))
          (if (/= i j) 
              (setf temp (* temp  (- (aref ppowers i) (aref ppowers j))))))
        (setf ret (+ ret (* (aref vals i) 
                            (/ num (* (- x (aref ppowers i)) temp)))))))
    ret))

;;; Computes F_p(B,Y) from [Katsurada] from p-adic diagonalization B of half-integral 
;;; matrix at specified Y. If interpolate=t, uses Lagrangian interpolation as in [King].
(defun fpiter (p B Y interpolate)
  (declare (fixnum p))
  (let* ((m (length B)) (etaxivec (etaxibulk p B)) (delvec (delbulk p B)) 
         (base (if interpolate 1 Y)) (c1) (c0) (X base) (i) (denom) 
         (del (aref delvec (1- m))) (deltilde) (xi) (xiprime) (etatilde) 
         (xitilde) (xitildeprime) (eta) (n) (sub1) (sub2) (sub3) (sub4) 
         (sub5) (sub6) (sub7) (degree (if interpolate (aref delvec 0) 0)) 
         (vals (make-array (+ 1 degree))) 
         (fvals (make-array (list m (+ m degree)))))
    (declare (fixnum i del deltilde xi xiprime etatilde xitilde 
                     xitildeprime eta n degree m))
    (dotimes (j (+ m degree)) 
      (if (plusp j) (setf X (* X p))) 
      (setf denom (- 1 (* p X)))
      (setf c1 (/ 1 denom ))  
      (setf c0 (/ (- (expt (* p X) (1+ del))) denom))
      (setf (aref fvals (1- m) j) (+ c1 c0)))
    (dotimes (k (1- m))
      (setf i (- m k 2)) 
      (setf n (+ k 2)) ;; check?
      (if (zerop (aref delvec i)) 
          (dotimes (j (+ i 1 degree)) (setf (aref fvals i j) 1))
        (if (evenp n) 
            (progn (setf del (aref delvec i))    
              (setf deltilde (aref delvec (1+ i)))
              (setf xi (aref etaxivec i))        
              (setf xiprime (+ 1 xi (- (* xi xi))))
              (setf etatilde (aref etaxivec (1+ i)))                            ; common subexpressions:
              (setf sub1 (expt p (/ n 2)))                                      ; p^(n/2)
              (setf sub2 (* sub1 sub1))                                         ; p^n
              (setf sub3 (* sub1 p xi))                                         ; p^(n/2 + 1) * xi
              (setf sub4 (* sub2 p))                                            ; p^(n+1)
              (setf sub5 (* sub1 xi))                                           ; p^(n/2) * xi
              (setf sub6 (+ del (- deltilde) (* xi xi)))                        ; del - deltilde + xi^2
              (setf sub7 (* xiprime etatilde (expt p (/ del 2))))               ; xiprime*etatilde*p^(del/2)
              (if (evenp xi) (setf sub7 (- sub7)))                              ; ... * (-1)^(xi+1)
              (setf X base)
              (dotimes (j (+ i 1 degree))
                (if (plusp j) (setf X (* p X))) 
                (setf denom (- 1 (* X X sub4)))   
                (setf c1 (/ (- 1 (* X sub5)) denom ))  
                (setf c0 (/ (* sub7 (- 1 (* X sub3)) 
                               (expt (* X sub1) sub6)) denom)) 
                (setf (aref fvals i j) (+ (* c1 (aref fvals (1+ i) (1+ j))) 
                                          (* c0 (aref fvals (1+ i) j))))))
          ;; else n odd
          (progn (setf del (aref delvec i)) 
            (setf deltilde (aref delvec (1+ i)))
            (setf xitilde (aref etaxivec (1+ i)))
            (setf xitildeprime (+ 1 xitilde (- (* xitilde xitilde))))
            (setf eta (aref etaxivec i))                                        ; common subexpressions:
            (setf sub1 (expt p (/ (1- n) 2)))                                   ; p^((n-1)/2)
            (setf sub2 (* sub1 p xitilde))                                      ; p^((n+1)/2) * xitilde
            (setf sub3 (- del deltilde -2 (* xitilde xitilde)))                 ; del - deltilde + 2 - xitilde^2
            (setf sub4 (* xitildeprime eta 
                          (expt p (/ (- (* del 2) deltilde -2) 2))))            ; (-1)^xitilde * ...
            (if (oddp xitilde) (setf sub4 (- sub4)))                            ; xitildeprime*eta*p^((2*del-deltilde+2)/2)
            (setf X base)
            (dotimes (j (+ i 1 degree))
              (if (plusp j) (setf X (* p X)))    
              (setf denom (- 1 (* X sub2)))     
              (setf c1 (/ 1 denom)) 
              (setf c0 (/ (* sub4 (expt (* X sub1) sub3)) denom))
              (setf (aref fvals i j) (+ (* c1 (aref fvals (1+ i) (1+ j))) 
                                        (* c0 (aref fvals (1+ i) j)))))))))
    (if interpolate 
        (progn
          (dotimes (i (1+ degree)) (setf (aref vals i) (aref fvals 0 i)))
          (values (interpolatep p degree Y vals) fvals))
      (values (aref fvals 0 0) fvals))))

;;; Computes F_2(B,Y) from [Katsurada] from 2-adic decomposition B of half-integral 
;;; matrix at specified Y. If interpolate=t, uses Lagrangian interpolation as in [King].
(defun f2iter (B Y interpolate)
  (let* ((m (length B)) (base (if interpolate 1 Y))
         (X) (temp) (foo) (goo) (n) (i) (temp2) (trm) 
         (fvals) (c1) (c0) (c10) (c11) (c20) (c21) (c212x) (c202x) 
         (denom) (denom2x) (sub0) (sub1) (sub2) (sub3) (sub4) (sub5) 
         (sub6) (sub7) (sub8) (sub9) (sub10) (degree) (vals)
         (xi) (xitilde) (xiprime) (xitildeprime) (xihat) (xihatprime) 
         (del) (deltilde) (delhat) (eta) (etatilde) (etahat) (sig) 
         (dimvec (make-array m))         (powervec (make-array m)) 
         (dvec (make-array m))           (delvec (make-array m)) 
         (detervec (make-array (1+ m)))  (auxdvec (make-array m)) 
         (auxdelvec (make-array m))      (etavec (make-array m)) 
         (xivec (make-array (1+ m)))     (hassevec (make-array (1+ m))))
    (setf temp 0)
    (dotimes (i m) 
      (if (< 1 (cadr (aref B (- m i 1)))) (setf temp (+ 2 temp)) 
        (setf temp (+ 1 temp)))
      (setf (aref dimvec (- m i 1)) temp))
    (setf temp 0)
    (setf (aref powervec 0) 0)
    (dotimes (i (1- m)) 
      (if (< 1 (cadr (aref B i))) (setf temp (+ 2 temp)) 
        (setf temp (+ 1 temp)))
      (setf (aref powervec (1+ i)) temp))
    (setf temp 1)
    (setf (aref detervec m) 1)
    (setf (aref hassevec m) 1) 
    (dotimes (j m)
      (setf i (- m j 1))
      (setf trm (aref B i))
      (setf temp2 (cadr trm))
      (cond ((= temp2 1) 
             (setf foo (* (expt 2 (car trm)) (caddr trm)))
             (setf (aref hassevec i) (* (aref hassevec (1+ i)) 
                                        (hilbp 2 foo temp) (hilbp 2 foo foo)))
             (setf temp (* temp foo)))
            ((= temp2 2) 
             (setf foo (* (expt 2 (car trm)) (caddr trm))) 
             (setf goo (* (expt 2 (car trm)) (cadddr trm)))
             (setf (aref hassevec i) (* (aref hassevec (1+ i))                  ; added mult below
                                        (hilbp 2 (* foo goo) temp) 
                                        (hilbp 2 foo foo)                         
                                        (hilbp 2 goo goo) (hilbp 2 foo goo)))
             (setf temp (* temp foo goo)))
            ((= temp2 3)  
             (setf foo (* (expt 4 (1- (car trm))) -1)) 
             (setf goo (expt 2 (1- (car trm))))
             (setf (aref hassevec i) (* (aref hassevec (1+ i)) 
                                        (hilbp 2 foo temp) (hilbp 2 goo -1)
                                        (hilbp 2 goo -1) -1))                   ; check this
             (setf temp (* temp foo)))
            ((= temp2 4)    
             (setf foo (* (expt 4 (1- (car trm))) 3)) 
             (setf goo (expt 2 (1- (car trm))))
             (setf (aref hassevec i) (* (aref hassevec (1+ i)) 
                                        (hilbp 2 foo temp) (hilbp 2 goo -1)
                                        (hilbp 2 goo 3)))                       ; check this 
             (setf temp (* temp foo))))
      (setf (aref detervec i) temp)
      ;; (setf temp2 (ord 2 (* (expt 2 (* 2 (floor (aref dimvec i) 2))) temp)))
      (setf temp2 (+ (* 2 (floor (aref dimvec i) 2)) (ord 2 temp)))
      (setf (aref dvec i) temp2)
      (setf (aref delvec i) 
        (if (evenp (aref dimvec i)) (* 2 (floor temp2 2)) temp2)))
    (setf (aref xivec m) 1)                                                     ; needed?
    (dotimes (i m) 
      (setf temp (aref dimvec i))
      (setf trm (aref B i))
      (if (< 1 (cadr trm))                                                      ; use stronger criterion?
          (progn 
            ;; (setf foo (ord 2 (* (expt 2 (* 2 (floor (1- temp) 2)))
            ;;                     (expt 2 (car trm)) (aref detervec (1+ i)))))
            (setf foo (+ (* 2 (floor (1- temp) 2)) (car trm) 
                         (ord 2 (aref detervec (1+ i)))))
            (setf (aref auxdvec i) foo)
            (setf (aref auxdelvec i) 
              (if (evenp temp) foo (* 2 (floor foo 2))))))
      (if (evenp temp) 
          (setf (aref xivec i) 
            (chip 2 (* (if (evenp (/ temp 2)) 1 -1) (aref detervec i))))))
    (dotimes (i m)
      (setf foo (aref dimvec i))
      (cond ((oddp foo)
             (setf goo (aref detervec i))
             (setf (aref etavec i) 
               (* (aref hassevec i) 
                  (hilbp 2 goo (if (evenp (/ (1- foo) 2)) goo (- goo))) 
                  (if (evenp (/ (1- (* foo foo)) 8)) 1 -1))))
            ((and (evenp foo) (< i (1- m))) 
             (setf goo (aref detervec (1+ i)))
             (cond ((and (= 2 (cadr (aref B i))) (evenp (aref dvec (1+ i))))
                    (setf (aref etavec i) 
                      (* (aref hassevec (1+ i))
                         (hilbp 2 (* goo (cadddr (aref B i)) ;; added 4/30
                                     (expt 2 (car (aref B i)))) 
                         (if (evenp (/ (- foo 2) 2)) goo (- goo))) 
                         (if (evenp (/ (1- (* (1- foo) (1- foo))) 8)) 1 -1))))
                   ((and (< 2 (cadr (aref B i))) 
                         (not (zerop (aref xivec (1+ i))))) 
                    (setf (aref etavec i) 
                      (* (aref hassevec (1+ i)) 
                         (hilbp 2 (expt 2 (car (aref B i))) 
                         (if (evenp (/ (- foo 2) 2)) goo (- goo))) 
                         (if (evenp (/ (1- (* (1- foo) (1- foo))) 8)) 1 -1))))
                   (t (setf (aref etavec i) 1))))
            (t (setf (aref etavec i) 1))))
    (setf temp (+ 3 (aref powervec (1- m))))
    (setf degree (if interpolate (aref delvec 0) 0))
    (setf fvals (make-array (list (1+ m) (+ temp degree))))
    (dotimes (j (+ temp degree)) (setf (aref fvals m j) 1))
    ;;(format t "dimvec    ~a ~%" dimvec)
    ;;(format t "powervec  ~a ~%" powervec)
    ;;(format t "delvec    ~a ~%" delvec)
    ;;(format t "auxdelvec ~a ~%" auxdelvec)
    ;;(format t "dvec      ~a ~%" dvec)
    ;;(format t "auxdvec   ~a ~%" auxdvec)
    ;;(format t "etavec    ~a ~%" etavec)
    ;;(format t "detervec  ~a ~%" detervec)
    ;;(format t "xivec     ~a ~%" xivec)
    ;;(format t "hassevec  ~a ~%" hassevec)
    (dotimes (k m)
      (setf i (- m k 1))
      (setf n (aref dimvec i))
      (setf trm (cadr (aref B i)))
      (if (and (< 2 trm) (zerop (car (aref B i))))
          (dotimes (j (+ 1 degree (aref powervec i))) (setf (aref fvals i j) 1))
        (cond ((and (evenp n) (< 1 trm))    
               (if (= 2 n)
                   (progn (setf xi (aref xivec i))
                     (setf xiprime (+ 1 xi (- (* xi xi))))
                     (setf xihat 1)
                     (setf xihatprime 1)
                     (setf deltilde (car (aref B i)))
                     (setf etatilde 1)
                     (setf del (aref delvec i))
                     (setf delhat 0)
                     (setf sig 0))
                 ;; else n>2
                 (progn (setf xi (aref xivec i))
                   (setf xiprime (+ 1 xi (- (* xi xi))))
                   (setf xihat (aref xivec (1+ i)))
                   (setf xihatprime (+ 1 xihat (- (* xihat xihat))))
                   (setf deltilde (aref auxdelvec i))
                   (setf etatilde (aref etavec i))
                   (setf del (aref delvec i))
                   (setf delhat (aref delvec (1+ i)))
                   (if (or (and (= 2 trm) (oddp (aref dvec i))) 
                           (and (< 2 trm) (zerop xihat))) 
                       (setf sig (/ (- (* 2 deltilde) del delhat -2) 2)) 
                     (setf sig 0))
                   ;; (format t "~% sig ~d " sig)
                   ))
               (setf sub0 (expt 2 (/ (- n 2) 2)))                               ; 2^((n-2)/2)
               (setf sub1 (* sub0 2))                                           ; 2^(n/2)
               (setf sub2 (* sub1 sub1))                                        ; 2^n
               (setf sub3 (* sub1 2 xi))                                        ; 2^(n/2 + 1) * xi
               (setf sub4 (* sub2 2))                                           ; 2^(n+1)
               (setf sub5 (* sub1 xi))                                          ; 2^(n/2) * xi
               (setf sub6 (+ del (- deltilde) (* xi xi) sig))                   ; del - deltilde + xi^2 + sig
               (setf sub7 (* xiprime etatilde (expt 2 (/ del 2))))              ; xiprime*etatilde*2^(del/2)
               (if (evenp xi) (setf sub7 (- sub7)))                             ; ... * (-1)^(xi+1)
               (setf sub8 (* sub1 xihat))                                       ; 2^(n/2) * xihat
               (setf sub9 (* xihatprime etatilde 
                             (expt 2 (/ (- (* 2 deltilde) 
                                           delhat -2 (* 2 sig)) 2))))
               (if (oddp xihat) (setf sub9 (- sub9)))
               (setf sub10 (- deltilde delhat -2 (* xihat xihat) sig))          ; deltilde-delhat+2-xihat^2-sig
               (setf X base)
               (dotimes (j (+ 1 degree (aref powervec i))) 
                 (if (plusp j) (setf X (* X 2)))  
                 (setf denom (- 1 (* X X sub4)))   
                 (setf c11 (/ (- 1 (* X sub5)) denom))   
                 (setf c10 (/ (* sub7 (- 1 (* X sub3)) 
                                 (expt (* X sub1) sub6)) denom)) 
                 (setf denom (- 1 (* X sub8)))
                 (setf denom2x (- 1 (* 2 X sub8)))    
                 (setf c21 (/ 1 denom))
                 (setf c20 (/ (* sub9 (expt (* X sub0) sub10)) denom))
                 (setf c212x (/ 1 denom2x))
                 (setf c202x (/ (* sub9 (expt (* 2 X sub0) sub10)) denom2x))
                 ;; (format t "~% ~d  ~d  ~d  ~d  ~d  ~d" 
                 ;;          c11 c10 c21 c20 c212x c202x)
                 (setf (aref fvals i j) 
                   (+ (* c11 c212x (aref fvals (+ i 1) (+ j 2))) 
                      (* (+ (* c11 c202x) (* c10 c21)) 
                         (aref fvals (+ i 1) (+ j 1))) 
                      (* c10 c20 (aref fvals (+ i 1) j))))))
              ((and (oddp n) (< 1 trm))   
               (setf xitilde (mod (1+ (aref auxdvec i)) 2))
               (if (= 2 (cadr (aref B i))) (setf xitilde 0)) ;; test 5/1!
               (setf xitildeprime 1)
               (setf etahat (aref etavec (1+ i)))
               (setf deltilde (aref auxdelvec i))
               (setf eta (aref etavec i))
               (setf del (aref delvec i))
               (setf delhat (aref delvec (1+ i)))
               (if (and (< 2 (cadr (aref B i))) (evenp (aref auxdvec i))) 
                   (setf sig 2) (setf sig 0))                                    ; use xitilde?
               (setf sub1 (expt 2 (/ (1- n) 2)))                                 ; 2^((n-1)/2)
               (setf sub2 (* sub1 2 xitilde))                                    ; 2^((n+1)/2) * xitilde
               (setf sub3 (- del deltilde -2 (* xitilde xitilde) (- sig)))       ; del - deltilde + 2 - xitilde^2 + sig
               (setf sub4 (* eta 
                             (expt 2 (/ (- (* del 2) deltilde -2 (- sig)) 2))))  ; (-1)^xitilde * ...
               (if (oddp xitilde) (setf sub4 (- sub4)))                          ; eta*2^((2*del -deltilde+2+sig)/2)        
               (setf sub5 (* sub1 sub1 2))                                       ; 2^n
               (setf sub6 (* sub1 xitilde))                                      ; 2^((n-1)/2) * xitilde
               (setf sub7 (* etahat (expt 2 (/ (- deltilde sig) 2))))            ; etahat*2^((deltilde-sig)/2)
               (if (evenp xitilde) (setf sub7 (- sub7)))                         ; ... * (-1)^(xitilde+1)
               (setf sub8 (- deltilde delhat (- (* xitilde xitilde)) sig))       ; deltilde-delhat+xitilde^2-sig
               (setf X base)
               (dotimes (j (+ 1 degree (aref powervec i))) 
                 (if (plusp j) (setf X (* X 2))) 
                 (setf denom (- 1 (* X sub2)))     
                 (setf c11 (/ 1 denom)) 
                 (setf c10 (/ (* sub4 (expt (* X sub1) sub3)) denom))
                 (setf denom (- 1 (* X X sub5)))
                 (setf denom2x (- 1 (* 4 X X sub5)))
                 (setf c21 (/ (- 1 (* X sub6)) denom))
                 (setf c20 (/ (* sub7 (- 1 (* X sub2)) 
                                 (expt (* X sub1) sub8)) denom))
                 (setf c212x (/ (- 1 (* 2 X sub6)) denom2x))
                 (setf c202x (/ (* sub7 (- 1 (* 2 X sub2)) 
                                   (expt (* 2 X sub1) sub8)) denom2x)) 
                 (setf (aref fvals i j) 
                   (+ (* c11 c212x (aref fvals (+ i 1) (+ j 2))) 
                      (* (+ (* c11 c202x) (* c10 c21)) 
                         (aref fvals (+ i 1) (+ j 1))) 
                      (* c10 c20 (aref fvals (+ i 1) j))))))
              ((and (evenp n) (= 1 trm))
               (setf del (aref delvec i)) 
               (setf deltilde (aref delvec (1+ i)))
               (setf xi (aref xivec i))
               (setf xiprime (+ 1 xi (- (* xi xi))))
               (setf etatilde (aref etavec (1+ i)))                             ; common subexpressions:
               (setf sub1 (expt 2 (/ n 2)))                                     ; 2^(n/2)
               (setf sub2 (* sub1 sub1))                                        ; 2^n
               (setf sub3 (* sub1 2 xi))                                        ; 2^(n/2 + 1) * xi
               (setf sub4 (* sub2 2))                                           ; 2^(n+1)
               (setf sub5 (* sub1 xi))                                          ; 2^(n/2) * xi
               (setf sub6 (+ del (- deltilde) (* xi xi)))                       ; del - deltilde + xi^2
               (setf sub7 (* xiprime etatilde (expt 2 (/ del 2))))              ; xiprime*etatilde*2^(del/2)
               (if (evenp xi) (setf sub7 (- sub7)))                             ; ... * (-1)^(xi+1)
               (setf X base)
               (dotimes (j (+ 1 degree (aref powervec i)))
                 (if (plusp j) (setf X (* X 2))) 
                 (setf denom (- 1 (* X X sub4)))  
                 (setf c1 (/ (- 1 (* X sub5)) denom ))  
                 (setf c0 (/ (* sub7 (- 1 (* X sub3)) 
                                (expt (* X sub1) sub6)) denom)) 
                 (setf (aref fvals i j) (+ (* c1 (aref fvals (1+ i) (1+ j))) 
                                           (* c0 (aref fvals (1+ i) j))))))
              ((and (oddp n) (= 1 trm))
               (if (= 1 n) 
                   (progn (setf del (aref delvec i)) 
                     (setf deltilde 0)
                     (setf xitilde 1)
                     (setf xitildeprime 1)
                     (setf eta 1)) 
                 ;; else n>1
                 (progn (setf del (aref delvec i)) 
                   (setf deltilde (aref delvec (1+ i)))
                   (setf xitilde (aref xivec (1+ i)))
                   (setf xitildeprime (+ 1 xitilde (- (* xitilde xitilde))))
                   (setf eta (aref etavec i))))                                 ; common subexpressions:
               (setf sub1 (expt 2 (/ (1- n) 2)))                                ; 2^((n-1)/2)
               (setf sub2 (* sub1 2 xitilde))                                   ; 2^((n+1)/2) * xitilde
               (setf sub3 (- del deltilde -2 (* xitilde xitilde)))              ; del - deltilde + 2 - xitilde^2
               (setf sub4 (* xitildeprime eta 
                             (expt 2 (/ (- (* del 2) deltilde -2) 2))))         ; (-1)^xitilde* ...
               (if (oddp xitilde) (setf sub4 (- sub4)))                         ; xitildeprime*eta*2^((2*del-deltilde+2)/2)
               (setf X base)
               (dotimes (j (+ 1 degree (aref powervec i)))
                 (if (plusp j) (setf X (* X 2))) 
                 (setf denom (- 1 (* X sub2)))     
                 (setf c1 (/ 1 denom)) 
                 (setf c0 (/ (* sub4 (expt (* X sub1) sub3)) denom))
                 (setf (aref fvals i j) (+ (* c1 (aref fvals (1+ i) (1+ j))) 
                                           (* c0 (aref fvals (1+ i) j)))))))))
    (if interpolate 
        (progn
          (setf vals (make-array (+ 1 degree)))
          (dotimes (i (1+ degree)) (setf (aref vals i) (aref fvals 0 i))) 
          (values (interpolatep 2 degree Y vals) fvals))
      (values (aref fvals 0 0) fvals))))


;;; Computes a p-adic diagonalization (as in [Watson]) from [SPLAG]-style genus symbol, odd prime
;;; The format of genus symbol is described in the file precomputed_misc.cl.
;;; If half = t, the calculation is done for the half-integral lattice whose double
;;; has the given genus symbol. If interpolate = t, Lagrangian interpolation is used
;;; as in [King] p858.
;;; examples: (genus2watsonp 11 #(2 10 0 (11 -1))) ; a10, prime 11
;;;           (genus2watsonp 5 #(2 10 0 (11 -1)))  ; a10, prime 5
(defun genus2watsonp (p genus) 
  (declare (fixnum p))
  (let* ((m (length genus)) (n (+ (aref genus 1) (aref genus 2)))               ; change to list?
         (diag (make-array n :initial-element 1)) (q) (nres 1) (epsn) 
         (epsprod 1) (dex 0) (determt 1) (j))
    (declare (fixnum m n nres epsn epsprod dex i j q))
    ;; (if (> 3 p) (print "use for odd primes p only")) 
    (do ((i 1 (1+ i))) ((> 1 (legendresymbol i p)) (setf nres i)))              ; nres is quadratic nonresidue
    (dotimes (i (- m 3)) ;; use dolist?
      (setf j (- m i 1)) 
      (setf q (car (aref genus j)))
      (setf epsn (cadr (aref genus j)))
      (if (zerop (mod q p)) 
          (progn (if (minusp epsn) (setq epsprod (- epsprod)))
            (dotimes (k (1- (abs epsn))) 
              (setf (aref diag dex) q) (setf dex (1+ dex)))
            (if (plusp epsn) (setf (aref diag dex) q) 
              (setf (aref diag dex) (* nres q)))
            (setf dex (1+ dex)))
        ;; else p doesn't divide q
        (setf determt (* determt (expt q (abs epsn))))))
    (if (oddp (aref genus 2)) (setf determt (- determt)))
    (if (/= (legendresymbol determt p) epsprod) 
        (if (<= dex (1- n)) (setf (aref diag (1- n)) nres) 
          (print "invalid genus"))) ; return nil?
    diag))

;;; Computes a 2-adic decomposition (as in [Watson]) from [SPLAG]-style genus symbol.
;;; The format of genus symbol is described in the file precomputed_misc.cl.
;;; (Check trace of part not divisible by 2 in type I lattice? Use exhanced genus symbol?)
;;; examples: (genus2watsontwo #(2 10 0 (11 -1))); a10
;;;           (genus2watsontwo #(2 11 0 (3 -1) (4 1 1 1))); a11
(defun genus2watsontwo (genus)
  (let* ((m (length genus)) (n (+ (aref genus 1) (aref genus 2))) 
         (diag '()) (antisqct) (pexcess 0) (dim) (q) (epsn) (ep) (i)
         (epsprod 1) (determt 1) (typ) (trc) (proto) (hhct) (yyct) (fac)
         (temp) (protoct) (a) (b) (c) (disc) (primelist '()))
    (declare (fixnum m n antisqct pexcess dim q epsn ep i epsprod
                     typ trc hhct yyct protoct a b c p k j))
    (dotimes (i (- m 3)) 
      (if (evenp (car (aref genus (+ i 3)))) 
          (setf n (- n (abs (cadr (aref genus (+ i 3)))))) 
        (setf determt (* determt (expt (car (aref genus (+ i 3))) 
                                       (abs (cadr (aref genus (+ i 3)))))))))
    (if (= 1 (aref genus 0))
        (progn   
          (setf fac (factorize determt))
          (dolist (term fac) 
            (if (and term (< 2 (car term))) 
                (setf primelist (append primelist (list (car term))))))
          (dolist (p primelist) 
            (setf temp 0) 
            (setf antisqct 0)
            (dotimes (i (- m 3))
              (setf q (car (aref genus (+ i 3))))
              (setf epsn (cadr (aref genus (+ i 3))))
              (if (zerop (mod q p)) 
                  (progn (setf temp (+ temp (* (abs epsn) (1- q)))) 
                    (if (and (minusp epsn) (oddp (ord p q))) 
                        (setf antisqct (1+ antisqct))))))
            (setf pexcess (+ pexcess (mod (+ temp (* 4 antisqct)) 8))))
          (setf pexcess (+ pexcess (aref genus 1) (- (aref genus 2))))))
    (dotimes (k (- m 2))
      (setf i (- m k 1))
      (if (= 2 i) (setf q 2) (setf q (car (aref genus i))))
      (if (evenp q) 
          (progn 
            (if (= 2 i)
                (progn (setf q 1) (setf dim n) (setf typ (aref genus 0)) 
                  (if (oddp (aref genus 2)) (setf determt (- determt)))
                  (setf temp (mod determt 8))
                  (if (or (and (or (= 1 temp) (= 7 temp)) (=  1 epsprod)) 
                          (and (or (= 3 temp) (= 5 temp)) (= -1 epsprod))) 
                      (setf ep 1) 
                    (setf ep -1))
                  (setf epsn (* ep dim)))
              ;; else i >= 3 
              (progn (setf epsn (cadr (aref genus i))) 
                (setf dim (abs epsn)) 
                (setf typ (caddr (aref genus i))) 
                (if (plusp epsn) (setf ep 1) 
                  (setf ep -1 epsprod (- epsprod))))) 
            (setf q (ord 2 q)) 
            (if (and (minusp epsn) (oddp q)) (setf pexcess (- pexcess 4)))
            (if (= 2 typ) 
                (progn (if (oddp dim) 
                           (print "shouldnt have type II comp of odd dim")); return nil?
                  (dotimes (j (1- (/ dim 2))) 
                    (setf diag (append diag (list (list (1+ q) 3)))))
                  (if (>= dim 2) 
                      (if (plusp ep) 
                          (setf diag (append diag (list (list (1+ q) 3)))) 
                        (setf diag (append diag (list (list (1+ q) 4)))))))
              ;; else typ=1
              (progn 
                (if (= i 2) (setf trc pexcess) 
                  (progn (setf trc (cadddr (aref genus i))) 
                    (setf pexcess (- pexcess trc))))
                (setf trc (mod trc 8))  
                (if (> trc 4) (setf trc (- trc 8)))
                (if (oddp (+ trc dim)) (print "incompatible trc and dim")); return nil?
                (cond ((= 1 dim) 
                       ;;(print "dim1")
                       (if (or (and (plusp  ep) (= 1 (abs trc)))
                               (and (minusp ep) (= 3 (abs trc))))
                           (setf diag (append diag (list (list q 1 trc))))  
                         (format t "incompat 1 dim e ~d and trc ~d" ep trc))); return nil?
                      ((= 2 dim)
                       ;;(print "dim2")
                       (cond ((and (plusp ep) (= 2 (abs trc))) 
                              (setf diag 
                                (append diag (list (list q 2 (/ trc 2) 
                                                         (/ trc 2))))))
                             ((and (plusp ep) (zerop trc)) 
                              (setf diag 
                                (append diag (list (list q 2 1 -1)))))
                             ((and (minusp ep) (= 2 (abs trc))) 
                              (setf diag 
                                (append diag 
                                        (list (list q 2 (- (/ trc 2)) 
                                                    (* trc 3 1/2))))))
                             ((and (minusp ep) (= 4 trc)) 
                              (setf diag (append diag (list (list q 2 1 3)))))
                             (t (print "incompat e and trc for dim 1 comp")))); return nil?
                      ((< 2 dim) 
                       (setf proto (make-array dim))
                       (dotimes (j dim) 
                         (setf (aref proto j) (if (evenp (- n j)) -1 1)))
                       (if (oddp dim) 
                           (if (plusp ep) 
                               (cond
                                ((= trc  1) 
                                 (setf (aref proto 0)  1) 
                                 (setf (aref proto 1)  1) 
                                 (setf (aref proto 2) -1))
                                ((= trc -1) 
                                 (setf (aref proto 0)  3) 
                                 (setf (aref proto 1)  3) 
                                 (setf (aref proto 2)  1))
                                ((= trc  3) 
                                 (setf (aref proto 0)  1) 
                                 (setf (aref proto 1)  1) 
                                 (setf (aref proto 2)  1))
                                ((= trc -3) 
                                 (setf (aref proto 0)  3) 
                                 (setf (aref proto 1)  3) 
                                 (setf (aref proto 2) -1))
                                (t (print "trc should be odd"))) ; return nil?
                             ;; else ep<0
                             (cond 
                              ((= trc  1) 
                               (setf (aref proto 0)  3) 
                               (setf (aref proto 1) -1) 
                               (setf (aref proto 2) -1))
                              ((= trc -1) 
                               (setf (aref proto 0) -3) 
                               (setf (aref proto 1)  1) 
                               (setf (aref proto 2)  1))
                              ((= trc  3) 
                               (setf (aref proto 0) -3) 
                               (setf (aref proto 1) -1) 
                               (setf (aref proto 2) -1))
                              ((= trc -3) 
                               (setf (aref proto 0) -3) 
                               (setf (aref proto 1)  1) 
                               (setf (aref proto 2) -1))
                              (t (print "trc should be odd")))) ; return nil?
                         ;; else dim even
                         (if (plusp ep) 
                             (cond
                              ((= trc  0) 
                               (setf (aref proto 0)  1) 
                               (setf (aref proto 1)  1) 
                               (setf (aref proto 2) -1) 
                               (setf (aref proto 3) -1))
                              ((= trc  2) 
                               (setf (aref proto 0)  1) 
                               (setf (aref proto 1)  1) 
                               (setf (aref proto 2)  1) 
                               (setf (aref proto 3) -1))
                              ((= trc -2) 
                               (setf (aref proto 0)  1) 
                               (setf (aref proto 1) -1) 
                               (setf (aref proto 2) -1) 
                               (setf (aref proto 3) -1))
                              ((= trc  4) 
                               (setf (aref proto 0)  1) 
                               (setf (aref proto 1)  1) 
                               (setf (aref proto 2)  1) 
                               (setf (aref proto 3)  1))
                              (t (print "trc should be even"))) ; return nil?
                           ;; else ep<0
                           (cond 
                            ((= trc  0) 
                             (setf (aref proto 0)  3) 
                             (setf (aref proto 1) -1) 
                             (setf (aref proto 2) -1) 
                             (setf (aref proto 3) -1))
                            ((= trc  2) 
                             (setf (aref proto 0)  3) 
                             (setf (aref proto 1)  1) 
                             (setf (aref proto 2) -1) 
                             (setf (aref proto 3) -1))
                            ((= trc -2) 
                             (setf (aref proto 0) -3) 
                             (setf (aref proto 1) -1) 
                             (setf (aref proto 2)  1) 
                             (setf (aref proto 3)  1))
                            ((= trc  4) 
                             (setf (aref proto 0)  3) 
                             (setf (aref proto 1)  1) 
                             (setf (aref proto 2)  1) 
                             (setf (aref proto 3) -1))
                            (t (print "trc should be even"))))) ; return nil?
                       (setf protoct dim) (setf hhct 0) (setf yyct 0)
                       (dotimes
                           (j (if (evenp dim) (1- (/ dim 2)) (/ (1- dim) 2)))
                         (setf a (aref proto (- protoct 1)))
                         (setf b (aref proto (- protoct 2)))
                         (setf c (aref proto (- protoct 3)))
                         (setf disc (* a b c))
                         (if (= 1 (hilbp 2 (- (* a b)) (- (* b c))))
                             (progn 
                               (setf (aref proto (- protoct 3)) 
                                 (mod (- disc) 8)) 
                               (setf hhct (1+ hhct)))
                           (progn 
                             (setf (aref proto (- protoct 3)) 
                               (mod (* 3 disc) 8))
                             (setf yyct (1+ yyct))))
                         (setf protoct (- protoct 2)))
                       (if (evenp yyct) (setf hhct (+ hhct yyct)) 
                         (progn (setf hhct (+ hhct yyct -1)) 
                           (setf diag (append diag (list (list (1+ q) 4))))))
                       (dotimes (j hhct) 
                         (setf diag (append diag (list (list (1+ q) 3)))))
                       (if (= 1 protoct) 
                           (setf diag 
                             (append diag (list (list q 1 (aref proto 0)))))
                         (setf diag 
                           (append diag (list (list q 2 (aref proto 0) 
                                                    (aref proto 1)))))))))))))
    (coerce diag 'vector)))




;;; Computes 2-adic decomposition of half-integral B from 2-adic decomposition of 2*B.
;;; Used for shifting from integral quadratic forms to half-integral quadratic forms. 
;;  examples: (genus2watsontwo #(2 11 0 (3 -1) (4 1 1 1))); a11
;;           (halfwatson (genus2watsontwo #(2 11 0 (3 -1) (4 1 1 1)))); a11
(defun halfwatson (B)
  (let* ((n (length B)) (ret (make-array n)) (temp))
    (declare (fixnum n i temp))
    (dotimes (i n)
      (setf (aref ret i) (copy-seq (aref B i)))
      (setf temp (car (aref ret i)))
      (if (plusp temp) (setf (car (aref ret i)) (1- temp)) 
        (print "wouldnt be half integral"))) ; return nil?
    ret))


;;; Sort vector representation of genus by prime powers, small to large.
;;; Can be used to sort arguments to mergesortedgenera below
(defun sortgenus (gen) 
  (concatenate 'vector (subseq gen 0 3) 
    (sort (subseq gen 3) #'(lambda (x y) (< (car x) (car y))))))


;;; Computes genus of direct sum from genus of each summands.
;;; Ugly. Re-written for presorted genus symbols below.
(defun mergegenera (gen1 gen2) 
  (let* ((m (length gen1)) (n (length gen2)) 
         (temp1) (temp2) (i) (j) (leftovers) (ret (make-array m)) 
         (gencop (make-array n :initial-element 1)))
    (declare (fixnum m n i j ii jj))
    (setf (aref ret 0) (min (aref gen1 0) (aref gen2 0)))
    (setf (aref ret 1) (+ (aref gen1 1) (aref gen2 1)))
    (setf (aref ret 2) (+ (aref gen1 2) (aref gen2 2)))
    (dotimes (ii (- m 3))
      (setf i (+ ii 3))
      (setf (aref ret i) (copy-seq (aref gen1 i)))
      (setf temp1 (aref ret i))
      (dotimes (jj (- n 3))  
        (setf j (+ jj 3))
        (setf temp2 (aref gen2 j))
        (if (and (= 1 (aref gencop j)) (= (car temp1) (car temp2)))
            (progn (setf (aref gencop j) 0)
              (setf (cadr (aref ret i)) 
                (* (if (plusp (cadr temp1)) 1 -1) 
                   (if (plusp (cadr temp2)) 1 -1)
                   (+ (abs (cadr temp1)) (abs (cadr temp2)))))
              (if (evenp (car temp1)) 
                  (cond ((and (= 1 (caddr temp1)) (= 1 (caddr temp2))) 
                         (setf (cadddr (aref ret i)) 
                           (mod (+ (cadddr temp1) (cadddr temp2)) 8)))
                        ((and (= 2 (caddr temp1)) (= 1 (caddr temp2))) 
                         (setf (aref ret i) 
                           (list (car temp1) (cadr (aref ret i)) 
                                 1 (cadddr temp2))))))))))
    (setf leftovers (make-array (- (reduce #'+ gencop) 3)))
    (setf j 0)
    (dotimes (i (- n 3)) 
      (if (= 1 (aref gencop (+ i 3))) 
          (progn (setf (aref leftovers j) (copy-seq (aref gen2 (+ i 3)))) 
            (setf j (1+ j)))))
    (concatenate 'vector ret leftovers)))

;;; Computes genus of direct sum from genus of each summands.
;;; Assumes each genus symbol is sorted and has no duplicates.
(defun mergesortedgenera (gen1 gen2) 
  (let* ((m (length gen1)) (n (length gen2)) (temp1) (temp2) (i 3) (j 3) 
         (epsn) (ret (make-array (+ n m -3))))
    (declare (fixnum m n i j epsn k))
    (setf (aref ret 0) (min (aref gen1 0) (aref gen2 0)))
    (setf (aref ret 1) (+ (aref gen1 1) (aref gen2 1)))
    (setf (aref ret 2) (+ (aref gen1 2) (aref gen2 2)))
    (if (> m 3) (setf temp1 (aref gen1 3)))
    (if (> n 3) (setf temp2 (aref gen2 3)))
    (do ((k 3 (1+ k))) ((and (= i m) (= j n)) (setf ret (subseq ret 0 k)))
      (cond ((= i m)  (setf (aref ret k) (copy-seq temp2)) (setf j (1+ j)) 
             (if (< j n) (setf temp2 (aref gen2 j))))
            ((= j n)  (setf (aref ret k) (copy-seq temp1)) (setf i (1+ i)) 
             (if (< i m) (setf temp1 (aref gen1 i))))
            ((> (car temp1) (car temp2))  
             (setf (aref ret k) (copy-seq temp2)) (setf j (1+ j)) 
             (if (< j n) (setf temp2 (aref gen2 j))))   
            ((> (car temp2) (car temp1))  
             (setf (aref ret k) (copy-seq temp1)) (setf i (1+ i)) 
             (if (< i m) (setf temp1 (aref gen1 i))))
            (t (setf epsn (* (if (plusp (cadr temp1)) 1 -1) 
                             (if (plusp (cadr temp2)) 1 -1)
                             (+ (abs (cadr temp1)) (abs (cadr temp2)))))
               (cond ((oddp (car temp1)) (setf (aref ret k)  
                                           (list (car temp1) epsn)))
                     ((and (= 1 (caddr temp1)) (= 1 (caddr temp2))) 
                      (setf (aref ret k) 
                        (list (car temp1) epsn 1 
                              (mod (+ (cadddr temp1) (cadddr temp2)) 8))))
                     ((and (= 2 (caddr temp1)) (= 1 (caddr temp2))) 
                      (setf (aref ret k) 
                        (list (car temp1) epsn 1 (cadddr temp2))))
                     ((and (= 1 (caddr temp1)) (= 2 (caddr temp2))) 
                      (setf (aref ret k) 
                        (list (car temp1) epsn 1 (cadddr temp1))))
                     (t (setf (aref ret k) (list (car temp1) epsn 2))))
               (setf i (1+ i)) (setf j (1+ j)) 
               (if (< i m) (setf temp1 (aref gen1 i))) 
               (if (< j n) (setf temp2 (aref gen2 j))))))
    ret))


;;; Computes genus of direct sum of root lattices from list of components.
;;; Currently limited to dimensions up to 40 for each component, as it uses
;;; precomputed genus of each component.
;;; example:
;;;     (multiple-value-setq (alist dlist elist) (string2roots "a1^8 a2 e8^2"))
;;;     (roots2genus alist dlist elist)
;;; same as:
;;;     (roots2genus '((1 8) (2 1)) '() '((8 2)))
;;; both return:
;;;     #(2 26 0 (2 8 1 0) (3 -1))
(defun roots2genus (rA rD rE)
  (let ((ret '#(2 0 0)))
    (declare (fixnum i))
    (declare (special +genA+ +genD+ +genE+))
    (dolist (trm rA) 
      (dotimes (i (cadr trm)) 
        (setf ret (mergesortedgenera ret (aref +genA+ (car trm))))))
    (dolist (trm rD) 
      (dotimes (i (cadr trm)) 
        (setf ret (mergesortedgenera ret (aref +genD+ (car trm))))))
    (dolist (trm rE) 
      (dotimes (i (cadr trm)) 
        (setf ret (mergesortedgenera ret (aref +genE+ (car trm))))))
    ret))

;;; Prints lists of components of root lattices in concise form to standard out.
;;; c.f. string2roots below.
;;; example: (printroots '((1 10) (10 2)) '((4 2)) '((8 1)))
(defun printroots (rA rD rE)
    (dolist (trm rA)  
      (if (= 1 (cadr trm)) (format t "a~d " (car trm))
        (format t "a~d^~d " (car trm) (cadr trm))))
    (dolist (trm rD)  
      (if (= 1 (cadr trm)) (format t "d~d " (car trm))
        (format t "d~d^~d " (car trm) (cadr trm))))
    (dolist (trm rE)  
      (if (= 1 (cadr trm)) (format t "e~d " (car trm))
        (format t "e~d^~d " (car trm) (cadr trm)))))


;;; Computes gamma(n/2) for n pos int. Omits the term pi^(1/2) when n odd.
(defun gammahalf (n)
  (declare (fixnum n i))
  (let ((ret 1))
    (declare (rational ret))
    (if (evenp n) (dotimes (i (1- (/ n 2))) (setf ret (* ret (1+ i))))
      (progn (dotimes (i (/ (1- n) 2)) (setf ret (* ret (- n i 1))))
        (setf ret (/ ret (expt 2 (1- n))))))
    ret))


;;; Not currently used
(defun gammap (p dim deter X)
  (declare (fixnum p dim));; deter?
  (let* ((n dim) (xi) (ret (- 1 X)) (XX (* X X)) (temp))
    (declare (fixnum n xi))
    ;;(format t "~%gammap ~a ~a ~a ~a" p dim deter X)
    (if (evenp n)
        (progn (setf xi (chip p (if (evenp (/ n 2)) deter (- deter))))
          (setf temp (/ 1 (- 1 (* xi X (expt p (/ n 2)))))) ; check nonzero?
          (setf ret (* ret temp)) ;; put in loop? do n/2 times? 
          (dotimes (i (/ n 2)) 
            (setf ret (* ret (- 1 (* (expt p (* 2 (1+ i))) XX))))))
      ;; else n odd
      (dotimes (i (/ (1- n) 2)) 
        (setf ret (* ret (- 1 (* (expt p (* 2 (1+ i))) XX))))))
    ret))

;;; Theta series of E_8 lattice, eg [SPLAG] chapter 4.8.1
;;; Returns 0 unless m is an even integer.
(defun e8coeff (m)
  (declare (fixnum m))
  (let ((j) (n (/ m 2)) (ret 0))
  (declare (fixnum n i j))
  (if (and (integerp n) (evenp m))
     (if (= 0 m) (setf ret 1) 
        (progn 
           (dotimes (i n) (setf j (1+ i))
              (if (zerop (mod n j)) (setf ret (+ ret (expt j 3)))))
           (setf ret (* 240 ret))))
     ;; else n odd
     (setf ret 0))
  ret))



;;; Converts string describing direct sum of root lattices to lists of components.
;;; Returns separate values for a, d, e components, in that order
;;; examples: (string2roots "a1^20 a10^2 d4^2 e8")
;;;   (multiple-value-setq (alist dlist elist) (string2roots "a1^20 a10^2 d4^2 e8"))
;;; Each component should have the form tn or tn^m, where t is type (a, d or e)
;;; n is the dimension, and m is the multiplicity (if ^m is omitted, taken to be 1). 
(defun string2roots (str)
  (declare (string str))
  (let ((n (length str)) (alist ()) (dlist ()) (elist ())
        (typ) (dim) (exp) (dex 0))
    (declare (fixnum n dim exp dex))
    (do ()
        ((= dex n))
      (if (char= (char str dex) #\space) (setf dex (1+ dex)))
      (setf exp 1)
      (setf typ (char str dex)) ; should be a, d or e
      (setf dex (1+ dex))
      (setf dim (digit-char-p (char str dex) 10)); first digit in dimension
      (setf dex (1+ dex))
      (do () 
      	  ((or (= dex n) (not (digit-char-p (char str dex) 10))))
              (progn (setf dim (+ (* 10 dim) (digit-char-p (char str dex) 10)))
                (setf dex (1+ dex))))         
      (if (and (< dex n) (char= (char str dex) #\^))
         (progn (setf dex (1+ dex))
            (setf exp 0)
            (do () 
      	       ((or (= dex n) (not (digit-char-p (char str dex) 10))))
                  (progn (setf exp (+ (* 10 exp) (digit-char-p (char str dex) 10)))
                      (setf dex (1+ dex))))))         
      (cond ((char= typ #\a) 
             (setf alist (append alist (list (list dim exp)))))
            ((char= typ #\d) 
             (setf dlist (append dlist (list (list dim exp)))))
            ((char= typ #\e) 
             (setf elist (append elist (list (list dim exp)))))))  
    (values alist dlist elist)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;