;;; Copyright Igor Rivin, 2004, 2005, All rights reserved.

;;; The functions in this file compute the determinant of matrix by Gaussian Elimination. 
;;;They also compute the determinant of an integer matrix by reducing mod
;;;p for a sequence of largish primes p , and then chinese-remaindering
;;;the results. 
(load "utils.ss")
(load "numbers.ss")
(load "math.ss")
(define make-gaussdet 
  (lambda (mul plus minus inv zerop?)
    (let* ((lincomb
           (lambda (l1 l2 coeff)
             (map (lambda (x y) (plus x (mul coeff y))) l1 l2)))
           (lincombspec
            (lambda (l1 l2)
              (lincomb l1 l2 (car l1)))))
     
      (lambda (mat)
        (let gd-iter ((numrots 0) (thedet 1) (themat mat))
          (if (null? (car themat))
              thedet
              (let ((a11 (caar themat)))
                (cond ((null? (cdr themat)) (mul thedet a11))
                      ((zerop? a11)
                       (let ((len (length themat)))
                         (if (= numrots len)
                             0
                             (if (even? len)
                                 (gd-iter (+ numrots 1) (minus thedet) (simplerot themat))
                                 (gd-iter (+ numrots 1) thedet (simplerot themat))))))
                      (else (let* ((multiplier (minus (inv a11)))
                                   (newfirst (map (lambda (x) (mul x multiplier))
                                                  (car themat)))
                                   (newmat (map cdr (map (lambda(x) 
                                                           (lincombspec x newfirst))
                                                         (cdr themat)))))
                                                       (gd-iter 0 (mul thedet a11) newmat)))))))))))
  

(define hadamardbound
  (lambda (mat)
    (apply * (map ll2norm mat))))

(define justenough ;; the code assumes that the product of all the numbers in numlist exceeds bound.
  (lambda (numlist bound)
    (let j-iter ((theprod 1) (thenums '()) (therest numlist))
      (if (> theprod bound)
          thenums
          (let ((firstnum (car therest)))
            (j-iter (* theprod firstnum) (cons firstnum thenums) (cdr therest)))))))

(define moddet
  (lambda( p)
    (make-gaussdet (lambda (x y) (modtimes x y p))
                   (lambda (x y) (modplus x y p))
                   (lambda (x) (modminus x p))
                   (lambda (x) (modinv x p))
                   (lambda (x) (modzerop x p)))))

(define can-primes (reverse (vec-sieve (expt 2 15))))

(define cdet
  (lambda (mat primelist)
    (let* ((thebound (* 2 (hadamardbound mat)))
           (theprimes (justenough primelist thebound))
           (thefuncs (map moddet theprimes))
           (results (map (lambda (x) (x mat)) thefuncs))   
           (tmpres (chineseremainder (transpose2 (list results theprimes)))))
      (if (< (* 2 (car tmpres)) (cadr tmpres)) (car tmpres) (- (car tmpres) (cadr tmpres))))))
    
(define chinesedet
  (lambda (mat)
    (cdet mat can-primes)))