; gaussian integers, part 1

(define (gs a b) (cons a b))
(define (re x) (car x))
(define (im x) (cdr x))

(define (gauss a b)
  (when (not (integer? a))
    (error 'gauss "must be integer"))
  (when (not (integer? b))
    (error 'gauss "must be integer"))
  (gs a b))

(define (gauss-from-complex x)
  (gauss (real-part x) (imag-part x)))

(define (gauss-to-complex x)
  (make-rectangular (re x) (im x)))

(define (gauss-zero? x)
  (and (zero? (re x)) (zero? (im x))))

(define (gauss-unit? x)
  (or (and (= (abs (re x)) 1) (zero? (im x)))
      (and (zero? (re x)) (= (abs (im x)) 1))))

(define (gauss-conjugate x)
  (gs (re x) (- (im x))))

(define (gauss-norm x)
  (define (square x) (* x x))
  (+ (square (re x)) (square (im x))))

(define (gauss-eql? x y)
  (and (= (re x) (re y))
       (= (im x) (im y))))

(define (gauss-add . xs)
  (define (add x y)
    (gs (+ (re x) (re y))
        (+ (im x) (im y))))
  (let loop ((xs xs) (zs (gs 0 0)))
    (if (null? xs) zs
      (loop (cdr xs) (add (car xs) zs)))))

(define (gauss-negate x)
  (gs (- (re x)) (- (im x))))

(define (gauss-sub . xs)
  (define (sub x y)
    (gs (- (re x) (re y)) (- (im x) (im y))))
  (cond ((null? xs) (error 'gauss-sub "no operands"))
        ((null? (cdr xs)) (gauss-negate (car xs)))
        (else (let loop ((xs (cdr xs)) (zs (car xs)))
                (if (null? xs) zs
                  (loop (cdr xs) (sub zs (car xs))))))))

(define (gauss-mul . xs)
  (define (mul x y)
    (gs (- (* (re x) (re y))
           (* (im x) (im y)))
        (+ (* (re x) (im y))
           (* (im x) (re y)))))
  (let loop ((xs xs) (zs (gs 1 0)))
    (if (null? xs) zs
      (loop (cdr xs) (mul (car xs) zs)))))

(define (gauss-quotient num den)
  (let ((n (gauss-norm den))
        (r (+ (* (re num) (re den))
              (* (im num) (im den))))
        (i (- (* (re den) (im num))
              (* (re num) (im den)))))
    (gs (round (/ r n)) (round (/ i n)))))

(define (gauss-remainder num den quo)
  (gauss-sub num (gauss-mul den quo)))

; gaussian integers, part 2

(define (gauss-gcd x y)
  (if (gauss-zero? y) x
    (let* ((q (gauss-quotient x y))
           (r (gauss-remainder x y q)))
      (gauss-gcd y r))))

(define (gauss-divides? d n)
  (gauss-zero?
    (gauss-remainder n d
      (gauss-quotient n d))))

(define (gauss-coprime? x y)
  (gauss-unit? (gauss-gcd x y)))

(define (gauss-prime? x)
  (cond ((gauss-unit? x) #f)
        ((zero? (re x))
          (and (prime? (im x))
               (= (modulo (im x) 4) 3)))
        ((zero? (im x))
          (and (prime? (re x))
               (= (modulo (re x) 4) 3)))
        (else (prime? (gauss-norm x)))))

(define (gauss-factors x)
  (define (find-k p a)
    (if (= (expm a (/ (- p 1) 2) p) (- p 1))
        (expm a (/ (- p 1) 4) p)
        (find-k p (+ a 1))))
  (let loop ((g x) (qs (list))
             (ps (factors (gauss-norm x))))
    ;(display g) (display qs) (display ps) (newline)
    (cond ((null? ps)
            (if (gauss-eql? g (gs 1 0)) qs (cons g qs)))
          ((= (car ps) 2)
            (loop (gauss-quotient g (gs 1 1))
                  (cons (gs 1 1) qs) (cdr ps)))
          ((= (modulo (car ps) 4) 3)
            (let ((q (gs (car ps) 0)))
              (loop (gauss-quotient g q)
                    (cons q qs) (cddr ps))))
          (else
            (let* ((p (car ps))
                   (k (find-k p 2))
                   (q (gauss-gcd (gs p 0) (gs k 1)))
                   (z (gauss-quotient g q)))
              (when (not (gauss-zero? (gauss-remainder g q z)))
                (set! q (gauss-conjugate q))
                (set! z (gauss-quotient g q)))
              (loop z (cons q qs) (cdr ps)))))))

; library

(define (expm b e m)
  (define (m* x y) (modulo (* x y) m))
  (cond ((zero? e) 1)
        ((even? e) (expm (m* b b) (/ e 2) m))
        (else (m* b (expm (m* b b) (/ (- e 1) 2) m)))))

(define (last-pair xs) (if (null? (cdr xs)) xs (last-pair (cdr xs))))
(define (cycle . xs) (set-cdr! (last-pair xs) xs) xs)
(define wheel (cons 1 (cons 2 (cons 2 (cycle 4 2 4 2 4 6 2 6)))))

(define (prime? n)
  (let loop ((f 2) (ws wheel))
    (cond ((< n (* f f)) #t)
          ((zero? (modulo n f)) #f)
          (else (loop (+ f (car ws)) (cdr ws))))))

(define (factors n)
  (let loop ((n n) (f 2) (ws wheel) (fs (list)))
    (cond ((< n (* f f)) (if (= n 1) fs (reverse (cons n fs))))
          ((zero? (modulo n f)) (loop (/ n f) f ws (cons f fs)))
          (else (loop n (+ f (car ws)) (cdr ws) fs)))))

(display (gauss-factors (gauss 361 -1767))) (newline)
(display (apply gauss-mul (gauss-factors (gauss 361 -1767)))) (newline)

((-7 . -2) (19 . 0) (4 . 1) (2 . 1) (1 . 1))
(361 . -1767)