;; Jacobi-Symbol-Berechnung
;; Bruno Haible Nov./Dez. 1987

(provide 'jacobi)

; (jacobi a b) liefert f�r ganze Zahlen a,b mit ungeradem b>0 das Jacobi-Symbol
; (a / b). Es ist =0 genau dann, wenn a und b nicht teilerfremd sind.
(defun jacobi (a b)
  (setq a (mod a b)) ; erst 0<= a < b erzwingen.
  (let ((v 1))
    (loop ; (a/b) * v bleibt invariant
       (cond ((= b 1) (return v))           ; Bei b=1 ist (a/b)=1.
             ((zerop a) (return 0))         ; Bei b>1 und a=0 ist (a/b)=0.
             ((> (* 2 a) b)                 ; Bei a > b/2 ist
                ; (a/b) = (-1/b) * ((b-a)/b) und (-1/b)=-1 falls b=3 mod 4.
                (setq a (- b a))
                (setq v (* v (ecase (mod b 4) ((1) 1) ((3) -1) ))) )
             ((evenp a)                     ; Bei b>1 und a=2a' ist
                ; (a/b) = (2/b) * (a'/b) und (2/b)=-1 falls b=3,5 mod 8.
                (setq a (exquo a 2))
                (setq v (* v (ecase (mod b 8)
                               ((1 7) 1) ((3 5) -1) ))) )
             (t ; a und b ungerade, 0<a<b/2<b
                ; (a/b) = (-1)^((a-1)/4)((b-1)/4) * (b/a)
                (if (and (= (mod a 4) 3) (= (mod b 4) 3))
                    (setq v (- v)))
                (rotatef a b)
                ; jetzt ist a > 2*b bekannt, setze a:=a mod b.
                (if (>= a (* 8 b))
                    (setq a (mod a b))
                    (do () ((< a b)) (setq a (- a b)) )
) ) )  )     )  )