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

(load "utils.ss")

(define modinv 
  (lambda (x p)
    (cadr (egcd x p))))

(define modtimes
  (lambda (x y p)
    (remainder (* x y) p)))

(define modplus
  (lambda (x y p)
    (remainder (+ x y) p)))

(define modminus
  (lambda (x p)
    (remainder (- x) p)))

(define modzerop
  (lambda (x p)
    (= 0 (remainder x p))))

;; Clever but inefficient sieve.

(define make-cycler
  (lambda (n)
    (let ((thestate 0))
      (lambda (x) 
        (if (eq? x 'step)
          (begin
            (set! thestate (+ thestate 1))
            (if (= thestate n)
                (set! thestate 0))
            thestate)
          n)))))


(define make-sieve
  (lambda ()
    (let ((curno 2)
          (curcycles (list (make-cycler 2))))
      (lambda ()
        (set! curno (+ curno 1))
        (let ((isprime (apply * (map (lambda (x) (x 'step)) curcycles))))
          (if (= 0 isprime)
              (list curno #f)
              (begin
                (set! curcycles (cons (make-cycler curno) curcycles))
                (list curno #t))
              ))))))

;; Sieve of Eratosthenes.

(define vec-sieve
  (lambda (n)
    (let ((thenums (make-vector n #t))
          (sqrtn (sqrt n)))
      (vector-set! thenums 0 #f)
      (vector-set! thenums 1 #f)
      (do ((i 2 (+ i 1))) 
        ((> i sqrtn))
        (if (vector-ref thenums i)
            (do ((j (+ i i) (+ j i)))
              ((>= j n))
              (vector-set! thenums j #f))))
      (select2 identity (range n) (vector->list thenums)))))


    
    
(define chineseremainder2
  (lambda (l1 l2)
    (let ((r1 (car l1))
          (p1 (cadr l1))
          (r2 (car l2))
          (p2 (cadr l2)))
      (let ((e_gcd (egcd p1 p2)))
        (let ((q1 (cadr e_gcd))
              (q2 (caddr e_gcd))
              (theprod (* p1 p2)) )
          (list (remainder (+ (* r2 p1 q1) (* r1 p2 q2)) theprod) theprod))))))
              
(define chineseremainder 
  (lambda (l)
    (cond ((null? l) '(0 1))
          ((null? (cdr l)) (car l))
          (else (fold chineseremainder2 l '(0 1))))))



(define egcd-core
  (lambda (m n) ;;assumes that m n are positive and n > m
    (let egcdaux ((m m) (n n) (themat (mat-identity 2)))
      (if (= m 0)
          (cons n (reverse (car themat)))
          (let* ((q (quotient n m))
                 (r (remainder n m))
                 (newmat (list (list 0 1) (list 1 (- 0 q)))))
            (egcdaux r m (mat-mult newmat themat)))))))

; computes the extended gcd of m and n.

(define egcd
  (lambda (m n)
    (cond ((< m 0)
           (let ((r1 (egcd (- 0 m) n)))
             (list (car r1) (- (cadr r1)) (caddr r1))))
          ((< n 0)
           (let ((r1 (egcd m (- 0 n))))
             (list (car r1) (cadr r1) (- (caddr r1)))))
          ((< n m)
           (let ((r1 (egcd n m)))
             (list (car r1) (caddr r1) (cadr r1))))
          (else (egcd-core m n)))))