home *** CD-ROM | disk | FTP | other *** search
/ GEMini Atari / GEMini_Atari_CD-ROM_Walnut_Creek_December_1993.iso / files / language / examples / jacobi.lsp < prev    next >
Lisp/Scheme  |  1993-10-23  |  1KB  |  34 lines
  1. ;; Jacobi-Symbol-Berechnung
  2. ;; Bruno Haible Nov./Dez. 1987
  3.  
  4. (provide 'jacobi)
  5.  
  6. ; (jacobi a b) liefert für ganze Zahlen a,b mit ungeradem b>0 das Jacobi-Symbol
  7. ; (a / b). Es ist =0 genau dann, wenn a und b nicht teilerfremd sind.
  8. (defun jacobi (a b)
  9.   (setq a (mod a b)) ; erst 0<= a < b erzwingen.
  10.   (let ((v 1))
  11.     (loop ; (a/b) * v bleibt invariant
  12.        (cond ((= b 1) (return v))           ; Bei b=1 ist (a/b)=1.
  13.              ((zerop a) (return 0))         ; Bei b>1 und a=0 ist (a/b)=0.
  14.              ((> (* 2 a) b)                 ; Bei a > b/2 ist
  15.                 ; (a/b) = (-1/b) * ((b-a)/b) und (-1/b)=-1 falls b=3 mod 4.
  16.                 (setq a (- b a))
  17.                 (setq v (* v (ecase (mod b 4) ((1) 1) ((3) -1) ))) )
  18.              ((evenp a)                     ; Bei b>1 und a=2a' ist
  19.                 ; (a/b) = (2/b) * (a'/b) und (2/b)=-1 falls b=3,5 mod 8.
  20.                 (setq a (exquo a 2))
  21.                 (setq v (* v (ecase (mod b 8)
  22.                                ((1 7) 1) ((3 5) -1) ))) )
  23.              (t ; a und b ungerade, 0<a<b/2<b
  24.                 ; (a/b) = (-1)^((a-1)/4)((b-1)/4) * (b/a)
  25.                 (if (and (= (mod a 4) 3) (= (mod b 4) 3))
  26.                     (setq v (- v)))
  27.                 (rotatef a b)
  28.                 ; jetzt ist a > 2*b bekannt, setze a:=a mod b.
  29.                 (if (>= a (* 8 b))
  30.                     (setq a (mod a b))
  31.                     (do () ((< a b)) (setq a (- a b)) )
  32. ) ) )  )     )  )
  33.  
  34.