; Bruno Haible 10.11.1987-12.12.1987, 25.2.1990
; Primfaktorzerleger f�r ganze Zahlen

(provide 'primzahl)
;(in-package 'primzahl)
;(export '(primteiler pfzv pdivisors divisors pfz prim-von-bis
;          *Kleinheit* prime probprime notprime factored *pfz-table* 
;          isprobprime isprime verffz setpfz))
(require 'avl) ; AVL-B„ume


(setq *print-circle* t) ; wegen der zyklischen Listen in fz-other


; (intlist a b) ergibt (a a+1 ... b-1 b), wenn a und b integers sind.
(defun intlist (a b)
  (do ((l nil (cons i l))
       (i b (1- i)))
      ((< i a) l)
) )


; exakter Quotient von Integers, schneller als / :
#-CLISP
(defun exquo (a b)
  (multiple-value-bind (q r) (floor a b)
    (unless (zerop r) (error "Quotient ~S/~S nicht exakt." a b))
    q
) )


(require '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.


; (pfz n) liefert die Primfaktorzerlegung von n>0, die Primfaktoren von n
; in aufsteigender Reihenfolge, evtl. mehrfach vertreten.

#|
; vorl„ufige Version des Primfaktorzerlegers
(defun pfz (n)
  (reverse
    (let ((l nil))         ; l = Liste der bereits gefundenen Faktoren
      (do () ((oddp n))
        (setq n (exquo n 2))
        (push 2 l))
      (let ((p 3))         ; p = Test-teiler (ungerade)
        (loop
          (if (= n 1) (return l))
          (setq s (isqrt n))
          (do () ((or (> p s) (zerop (mod n p))))
                 (incf p 2))
          (if (> p s) (return (cons n l)))
          (push p l)
          (setq n (exquo n p))
) ) ) ) )
|#


; (Miller-Rabin-test N) testet, ob eine gegebene Zahl n>1 wohl eine
; Primzahl ist.
; Falls das Ergebnis NIL ist, ist N garantiert zusammengesetzt,
; Eventuell kommt ein nichttrivialer Teiler von N als zweiter Wert heraus.
; Falls das Ergebnis T ist, ist N mit hoher Wahrscheinlichkeit prim.
; (Fehlerwahrscheinlichkeit < 1E-30 )
; Ist n gerade und >2, so ist das Ergebnis NIL.
(defun miller-rabin-test (n)
  ; Vorausschleife, die bereits 87% aller Zahlen faktorisiert
  (dolist (p '(2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67))
    (if (zerop (mod n p))
        (return-from miller-rabin-test (if (= n p) t (values nil p)))))
  ; erst jetzt geht's richtig los:
  (let (blist n1 s t1)
    (if (< n 2) (return-from miller-rabin-test nil))
    (setq blist (cond ((< n 2000) '(2))
                      ((< n 1300000) '(2 3))
                      ((< n 25000000) '(2 3 5))
                      ((< n 3200000000) '(2 3 5 7))
                      (t '(2 3 5 7 11 13 17 19 23 29
                           31 37 41 43 47 53 59 61 67 71
                           73 79 83 89 97 101 103 107 109 113
                           127 131 137 139 149 151 157 163 167 173
                           179 181 191 193 197 199 211 223 227 229))))
    ; mit 50 Primzahlen: P(falsch) < (1/4)^50 < 1E-30
    (setq n1 (- n 1))
    (setq s n1  t1 0)
    (do () ((oddp s))
           (setq s (exquo s 2)) (incf t1))
    ; n-1 = s * 2^t1 mit ungeradem s zerlegt.
    (dolist (b blist)
      (if (zerop (mod n b))
        (if (= n b) (return-from miller-rabin-test t)
                    (return-from miller-rabin-test (values nil b))))
      (do ((c (exptmod (mod b n) s n) (sqrmod c n))
           (c1 n1 c)
           (j 0 (1+ j)))
          ((or (= j t1) (= c 1))
           (cond ((not (= c 1)) ; sicher j=t1, aber c/=1
                  (return-from miller-rabin-test nil))
                 ((not (= c1 n1)) ; sicher c=1, c1/=n-1, also j>0
                  ; ggt(c1-1,n) * ggt(c1+1,n) = ggt(c1^2-1,n) = ggt(c-1,n) = n
                  (return-from miller-rabin-test (values nil (gcd n (- c1 1)))))
      )   )) ; Eine Iteration �ber b beendet
    )   ; blist abgearbeitet
    t   ; n wahrscheinlich prim
) )

; Hilfsfunktion f�r Miller-Rabin-Test:
; Quadrieren modulo n von x mit 0<=x<n
(defun sqrmod (x n)
  (mod (* x x) n))

; Hilfsfunktion f�r Miller-Rabin-Test:
; Potenzieren modulo n von a mit 0<=a<n und b>0
(defun exptmod (a b n)
  (let ((c 1))
    (loop    ; a^b*c mod n bleibt invariant
      (if (oddp b) (setq c (mod (* a c) n)))
      (if (= b 1) (return c))
      (setq a (mod (* a a) n) b (floor b 2))
) ) )


; verifizierter Primzahltest f�r kleine Zahlen n>0
(defun klprimtest (n)
  (cond ((= n 1) nil)
        ((evenp n) (= n 2))
        (t (let ((s (isqrt n)))
             (do ((i 3 (+ i 2)))
                 ((> i s) t)
               (if (zerop (mod n i)) (return nil))
) )     )  ) )


; Primfaktorzerlegung kleiner Zahlen n>0
(defun klpfz (n)
  (let ((l nil))
    (do () ((oddp n)) (setq n (ash n -1)) (push 2 l)) ; Zweier heraus
    (block Rest-prim
      (do ((p 3))
          ((= n 1))
        (loop ; in dieser Schleife wird ein Teiler gesucht
          (cond ((> (* p p) n) (return-from Rest-prim (push n l)))
                ((zerop (mod n p)) (push p l) (setq n (exquo n p)) (return))
                (t (incf p 2))
    ) ) ) )
    (nreverse l)
) )


; Datenstruktur: In *pfz-table* steht eine Tabelle von Primfaktorzerlegungs-
; ergebnissen. Zu jeder Zahl (sie ist der Schl�ssel) steht drin:
; Ob sie Primzahl ist, eine eventuelle Faktorzerlegung bzw. evtl. eine
; Primitivwurzel.

(defparameter *Kleinheit* 10000
  "Nur Zahlen >= *Kleinheit* werden in die PFZ-Tabelle aufgenommen.")
; *Kleinheit* sollte nur ge„ndert werden, wenn gleichzeitig
; (setq *pfz-table* nil) ausgef�hrt wird!
(assert (> *Kleinheit* 1))

(deftype fz-class () "Das Ergebnis einer Faktorzerlegung, ganz grob"
  '(member              ; ein Symbol:
            nil        ; (nur kurz nach der Initialisierung)
            prime      ; n ist prim
            probprime  ; n ist wahrscheinlich prim
            notprime   ; n ist garantiert zusammengesetzt, Faktoren unbekannt
            factored   ; n ist zusammengesetzt, schon faktorisiert
)  )

(defstruct fz     "Eine Faktorisierung als Ganzes"
  (num 2   :type integer :read-only t)  ; die zu faktorisierende Zahl, >1
  (cl  nil :type fz-class) ; Klassifikation
  (other nil) ; weitere Information:
     ; bei cl = prime : NIL oder eine Primitivwurzel mod n
     ; bei cl = probprime : Eine zyklische Liste von wartenden Funktionen.
     ;                      Jede von ihnen liefert NIL, wenn sie wieder aufge-
     ;                      rufen werden will, und (NIL . factor), falls sie
     ;                      die Zahl als zusammengesetzt nachgewiesen hat, evtl.
     ;                      einen nichttrivialen Faktor gefunden hat, und
     ;                      (T . Prw), falls sie die Zahl als prim nachgewiesen
     ;                      hat, mit evtl. einer Primitivwurzel Prw.
     ;                      (Eventuell zu Beginn die leere Liste.)
     ; bei cl = notprime : Eine zyklische Liste von wartenden Funktionen
     ;                     (Closures), die mitten in der Suche nach einem
     ;                     Faktor von n stecken. Jede von ihnen liefert NIL,
     ;                     wenn sie wieder aufgerufen werden will, und einen
     ;                     nichttrivialen Faktor, wenn sie einen gefunden hat.
     ;                     (Eventuell zu Beginn die leere Liste.)
     ; bei cl = factored : Ein Konstrukt vom Typ fzl, das die Faktorenzerlegung
     ;                     von n widerspiegelt.
)

(defstruct fzl "Faktorenzerlegung, Liste"
  (allprimes nil  :type symbol) ; Flag, ob alle Faktoren von factorlist
                                ; bekanntermažen Primzahlen sind. 
  (factorlist nil :type list) ; Eine Liste (f1 ... fk) von Faktoren von n,
        ; mit 1 < f1 <= ... <= fk und n = f1 * ... * fk. Die fj sind dabei
        ; (Pointer auf) Konstrukte vom Typ fz, bzw. bei fj<*Kleinheit* ist fj
        ; die Zahl selbst und Primzahl.
)


; (fz-key f) ergibt zu einem Element einer Faktorzerlegungsliste die
; wirkliche Zahl.
(defun fz-key (g)
  (if (integerp g) g (fz-num g))
)


; (fz-isprime f) stellt fest, ob ein Element einer Faktorzerlegungsliste
; prim ist.
(defun fz-isprime (f)
  (or (integerp f) (eq (fz-cl f) 'prime))
)


; (mergefz l1 l2) mischt zwei Faktorzerlegungen l1 und l2 von n1 bzw. n2,
; so daž eine Faktorzerlegung von n1*n2 herauskommt. Vorsicht: Destruktiv!
(defun mergefz (l1 l2)
  (merge 'list l1 l2 #'< :key #'fz-key)
)


; (fzl-recalc h) berechnet neu, ob alle Faktoren von h (h ist vom Typ fzl)
; prim sind, tr„gt diese Information in h ein und liefert sie als Ergebnis.
(defun fzl-recalc (h)
  (setf (fzl-allprimes h)
    (reduce #'(lambda (allp f) (and allp (fz-isprime f)))
            (fzl-factorlist h)
            :initial-value t
) ) )
  

(defun fz-comp (x y) (< (fz-num x) (fz-num y)))

(defun fz-eq-test (x y) (= (fz-num x) (fz-num y)))

(defvar *pfz-table* (avl:seq-to-avl nil #'fz-comp)
  "Die globale Tabelle, die alle Primfaktorzerlegungsergebnisse enth„lt")
; Die in *pfz-table* abgespeicherten Werte sind alle vom Typ fz.
; Die Ordnungsrelation ist durch die Zahl selbst (fz-num *) gegeben.


; Ausgeben der globalen Tabelle:
(defun pt (&optional (output-stream *standard-output*))
  (pprint (avl:avl-to-seq *pfz-table*) output-stream)
)


(defun pfz-table-lookup (n)
; ergibt die in *pfz-table* stehende Faktorzerlegung von n
; 1. Version (langsamer):
; (avl:member (make-fz :num n) *pfz-table* #'fz-comp
;   #'(lambda (x y) (and (fz-eq-test x y) y)) ))
; 2. Version (schneller, aber weniger portabel):
  (do ((tr *pfz-table*))
      ((or (null tr) (= n (fz-num (avl::node-value tr))))
       (cond (tr (avl::node-value tr))))
    (cond ((< n (fz-num (avl::node-value tr))) (setq tr (avl::node-left tr)))
          ((> n (fz-num (avl::node-value tr))) (setq tr (avl::node-right tr)))
          (t (error "Unvergleichbarkeit zweier Zahlen!"))
) ) )


; (pfz-table-insert b) f�gt zu *pfz-table* einen neuen Eintrag b (vom Typ fz)
; hinzu. Ergebnis ist b.
(defun pfz-table-insert (b)
  (setq *pfz-table* (avl:insert b *pfz-table* #'fz-comp #'fz-eq-test))
  b
)


; (pfz-table-lookup-insert n) sucht in *pfz-table* nach einem Eintrag zur Zahl
; n. Wird kein Eintrag gefunden, so wird ein neuer in die Tabelle eingef�gt.
; Das Ergebnis ist stets der entsprechende Eintrag in der Tabelle *pfz-table*.
(defun pfz-table-lookup-insert (n)
  (assert (>= n *Kleinheit*))
  (let ((b (pfz-table-lookup n)))
    (unless b (pfz-table-insert (setq b (make-fz :num n))))
    b
) )


; (recall-factors n) ergibt eine Liste aller Faktoren der nat�rlichen Zahl n>1,
; soweit sie bekannt sind, in einer sortierten Liste, im selben Format wie bei
; fzl. Der zweite Wert zeigt an, ob alles sicher Primzahlen sind.
(defun recall-factors (n)
  (if (< n *Kleinheit*)
      (values (klpfz n) t)
      (let ((b (pfz-table-lookup n)))
        (if (null b)
            ; nicht gefunden
            (progn
              (mr n)
              (values (list (pfz-table-lookup n)) nil)
            )
            ; gefunden
            (case (fz-cl b)
              ((prime) (values (list b) t))
              ((probprime notprime) (values (list b) nil))
              ((factored)
               (setq b (fz-other b))
               (values (fzl-factorlist b) (fzl-allprimes b))
              )
              (t (error "Falsch aufgebauter fz-Struct!"))
) )   ) )   )


; (build-fzl l) ergibt die Faktorzerlegung vom Typ fzl von n, wenn bekannt ist,
; daž das Produkt aller Elemente von l (alles nat�rliche Zahlen >1) = n ist.
(defun build-fzl (l)
  (let ((l1 nil) (allp t))
    (dolist (f l) (multiple-value-bind (f1 f2) (recall-factors f)
                    (setq l1 (mergefz l1 (copy-list f1)))
                    (setq allp (and allp f2))
    )             )
    (make-fzl :allprimes allp :factorlist l1)
) )


; F�hrt f�r eine Zahl >= *Kleinheit*, �ber die noch nichts bekannt ist,
; den Miller-Rabin-Test durch und speichert das Ergebnis in *pfz-table* ab.
; Bei Ergebnis t (also n wohl prim) ist n=2 oder n>2 ungerade.
(defun mr (n)
  (multiple-value-bind (pr fct) (miller-rabin-test n)
    (cond (pr ; wohl prim
              (pfz-table-insert (make-fz :num n :cl 'probprime))
              t)
          ; sonst sicher nicht prim
          (fct ; Faktor gefunden
              (pfz-table-insert
                (make-fz :num n :cl 'factored
                         :other (build-fzl (list fct (exquo n fct)))))
              nil)
          (t ; kein Faktor gefunden
              (pfz-table-insert (make-fz :num n :cl 'notprime))
              nil)
) ) )


; (isprobprime n) stellt fest, ob die nat�rliche Zahl n wahrscheinlich eine
; Primzahl ist.
; Bei geraden Zahlen >2 liefert dies stets NIL.
(defun isprobprime (n)
  (if (< n *Kleinheit*)
      (klprimtest n)
      (let ((b (pfz-table-lookup n)))
        (if b (member (fz-cl b) '(prime probprime)) ; bekannt
              ; sonst muž gerechnet werden:
              (mr n)
) )   ) )


; (isprime n) stellt fest, ob die nat�rliche Zahl n eine Primzahl ist.
(defun isprime (n)
  (if (< n *Kleinheit*)
      (klprimtest n)
      (let ((b (pfz-table-lookup n)))
        (if (and b (member (fz-cl b) '(prime))) ; bekannt
            (return-from isprime t))
        (if (and b (member (fz-cl b) '(notprime factored))) ; bekannt
            (return-from isprime nil))
        (assert (or (null b) (eq (fz-cl b) 'probprime)))
        (if (and (null b) (null (mr n))) (return-from isprime nil))
        ; jetzt ist bekannt, daž n wohl prim ist, steht in *pfz-table*.
        (if (null b) (setq b (pfz-table-lookup n)))
        (assert (and b (member (fz-cl b) '(prime probprime))))
        (if (eq (fz-cl b) 'prime) (return-from isprime t))
        (verffz b) ; Primzahlverdacht erh„rten
        (assert (not (eq (fz-cl b) 'probprime)))
        (eq (fz-cl b) 'prime)
) )   )


(defparameter nearly-infinite 1000000000000000000
  "Schier unendliche Zahl von Iterationen")


; Sei (fz-cl b) = 'probprime und (fz-other b) /= NIL.
; (tryprimeprove b) ruft die wartenden Funktionen in (fz-other b) der Reihe
; nach auf, bis ein Ergebnis kommt. Maximal k Iterationen (k fehlt->grenzenlos).
; Das Ergebnis ist die Anzahl der verbrauchten Iterationen.
(defun tryprimeprove (b &optional (k nearly-infinite) &aux (it 0) erg)
  (loop
     (when (>= it k) (return it))
     (incf it)
     (when (setq erg (funcall (first (fz-other b))))
       ; Ungewižheit beendet.
       ; erg = (NIL) oder (NIL . factor) oder (T) oder (T . Prw)
       (cond ((car erg) (setf (fz-cl b) 'prime)
                        (setf (fz-other b) (cdr erg)) )
             ((null (cdr erg)) (setf (fz-cl b) 'notprime)
                               (setf (fz-other b) nil) )
             (t (setq erg (cdr erg))
                (setf (fz-other b) nil)
                (setf (fz-cl b) 'factored)
                (setf (fz-other b)
                      (build-fzl (list erg (exquo (fz-num b) erg))))
       )     )
       (return it)
     )
     (pop (fz-other b)) ; Funktionen-Schlange rotieren
) )


; 1. Primzahlbeweiser: Bei jedem Aufruf wird eine gewisse Anzahl von
; Teilern durchprobiert, maximal bis zur Wurzel.
; [Allgemein bekannt.]
(defun tryprimeprover1 (n)
  (let ((p 1)
        (s (isqrt n)))
    #'(lambda ()
        (block tryprimeproving1
          (dotimes (i (min 100 (ceiling (- s p) 2)))
            (incf p 2)
            (if (zerop (mod n p))
                (return-from tryprimeproving1 (cons nil p)))
          )
          (if (>= p s) (cons t nil) nil)
) )   ) )


; 2. Primzahlbeweiser: Aus einer hinreichend guten Faktorzerlegung von n-1
; wird eine Zahl der Ordnung n-1 mod n bestimmt, also ist n prim.
; [Donald E. Knuth, The art of computer programming, Band 2, Seminumerical
; algorithms, Abschnitt 4.5.4, Aufgabe 26, Seite 397.]
(defun tryprimeprover2 (n)
  (let ((enough nil)
        f r plist
        (n1 (1- n)))
    (flet ((factored-enough (b) ; stellt fest, ob n-1 weit genug faktorisiert
             (setq enough                                                 ; ist.
               (and (= (fz-num b) n1)
                    (eq (fz-cl b) 'factored)
                    (progn
                      ; f sammelt die primen Faktoren p von n,
                      ; r sammelt die nichtprimen Faktoren von n,
                      ; plist = Menge aller primen Faktoren von n.
                      (setq f 1)
                      (setq r 1)
                      (setq plist nil)
                      (dolist (p (fzl-factorlist (fz-other b)))
                        (if (fz-isprime p)
                            (progn
                              (setq f (* f (fz-key p)))
                              (setq plist (adjoin (fz-key p) plist))
                            )
                            (setq r (* r (fz-key p)))
                      ) )
                      ; alles ok, wenn r<=f+1
                      (<= r (1+ f))
           ) ) )    )
           (tpp2-rest () ; f�hrt den Primzahlbeweis aus, wenn f,r,plist da.
             ; Bei r=1 k”nnte man sogar eine Primitivwurzel errechnen.
             ; [D. E. Knuth, Band 2, Abschnitt 4.5.4, Aufgabe 10, Seite 395.]
             (do ((x 2 (1+ x)))
                 ((null plist) (cons t nil))
               (unless (= (exptmod x n1 n) 1)
                 (format t
                   "~%Bei ~D hatte der Miller-Rabin-Test falsch geraten!~%"
                   n)
                 (return (cons nil nil))
               )
               (setq plist
                 (remove-if #'(lambda (p)
                                (= (gcd (1- (exptmod x (exquo n1 p) n)) n) 1) )
                            plist
          )) ) ) )
      (if (< n1 *Kleinheit*)
          #'(lambda ()
              (setq f n1)
              (setq r 1)
              (setq plist (primteiler n1))
              (tpp2-rest))
          (let ((b (pfz-table-lookup n1)))
            (when (null b)
                  (isprobprime n1) ; n1 in *pfz-table* eintragen
                  (setq b (pfz-table-lookup n1)) )
            #'(lambda ()
                (verffz b 1 :return-test #'factored-enough)
                ; Anmerkung: Immer nur eine Zeitscheibe ist wohl ung�nstig,
                ; weil dann stets der kleinste m”gliche Faktor von n1
                ; verfeinert wird.
                (when (or enough (factored-enough b))
                      (tpp2-rest)
              ) )
) ) ) )   )


; Sei (fz-cl b) = 'notprime und (fz-other b) /= NIL.
; (tryfactor b) ruft die wartenden Funktionen in (fz-other b) der Reihe nach
; auf, bis ein Ergebnis kommt. Maximal k Iterationen (k fehlt -> grenzenlos).
; Ergebnis ist die Anzahl der verbrauchten Iterationen.
(defun tryfactor (b &optional (k nearly-infinite) &aux (it 0) erg)
  (loop
    (when (>= it k) (return it))
    (incf it)
    (when (setq erg (funcall (first (fz-other b))))
      (setf (fz-other b) nil)
      (setf (fz-cl b) 'factored)
      (setf (fz-other b) (build-fzl (list erg (exquo (fz-num b) erg))))
      (return it))
    (pop (fz-other b)) ; Funktionen-Schlange rotieren
) )


; Alle Faktorisierungsalgorithmen bekommen eine ungerade Zahl n>1, die
; garantiert zusammengesetzt ist, und liefern eine parameterlose Funktion ab,
; die nach wiederholtem Aufruf einen Faktor von n liefert (und NIL, solange
; sie noch keinen Faktor gefunden hat).


; erster Faktorisierungsalgorithmus: Teilersuche bis zur Wurzel
; [Allgemein bekannt.]
(defun tryfactor1 (n)
  (let ((p 1)
        (s (isqrt n)))
    #'(lambda ()
        (block tryfactoring1
          (dotimes (i (min 100 (ceiling (- s p) 2)))
            (incf p 2)
            (if (zerop (mod n p)) (return-from tryfactoring1 p))
          )
          (if (>= p s) (error "Nicht-Primzahl ~D ohne Teiler." n))
          nil
) )   ) )


; zweiter Faktorisierungsalgorithmus: Pollardsche Iteration, in der Version
; von R.P. Brent.
; [Richard P. Brent, An improved Monte Carlo factorization algorithm,
;  BIT 20 (1980), 176-184, Seite 182].
(defun tryfactor2 (n &aux (m 100))
  (let ((state 1) ; Zustand: 1 bei der Initialisierung,
                  ;          2 bei den Iterationen ins Leere,
                  ;          3 bei den Iterationen in Bl”cken.
        (s 1) ; Iterationsparameter
        x     ; relative Anfangszahl
        y     ; Iterationswert
        r     ; wird immer wieder verdoppelt
        k     ; Anzahl der bisher ausgef�hrten Iterationen (0<=k<=r)
        q     ; bisher aufgelaufenes Produkt
       ) ; m = Blockl„nge
    (flet ((brentiter (x)
             (mod (+ (* x x) s) n)
          ))
      #'(lambda ()
          (let ((it 100)) ; Noch verbleibende Iterationen
            (loop
              (if (<= it 0) (return nil))
              (case state
                ((1) (setq y 0)
                     (setq x y)
                     (setq r 1)
                     (setq k 0)
                     (setq q 1)
                     (setq state 2)
                )
                ((2) (let ((i (min it (- r k))))
                       (dotimes (i1 i) (setq y (brentiter y)))
                       (incf k i)
                       (decf it i))
                     (when (= k r)
                       (setq k 0)
                       (setq state 3)
                )    )
                ((3) (let ((ys y)
                           (i (min it (- r k) m))
                           g)
                       (dotimes (i1 i)
                         (setq y (brentiter y))
                         (setq q (mod (* q (abs (- x y))) n))
                       )
                       (incf k i)
                       (decf it i)
                       (setq g (gcd q n))
                       (cond ((> g 1) ; Habe wohl etwas gefunden
                         (if (= g n) ; wirklich?
                           (loop ; nachiterieren
                             (setq ys (brentiter ys))
                             (decf it)
                             (setq g (gcd (- x ys) n))
                             (if (> g 1) (return))
                         ) )
                         (if (< g n) ; ja!
                           (return g))
                         ; nein ...
                         (incf s) ; ganz neuer Versuch
                         (setq state 1)
                        )
                        ((= k r) ; state 3 zu Ende
                         (setq x y)
                         (setq r (* 2 r))
                         (setq k 0)
                         (setq state 2)
                )    ) ))
) ) )   ) ) ) )


; dritter Faktorisierungsalgorithmus: Fermats Methode.
; [Donald E. Knuth, The art of computer programming, Band 2, Seminumerical
;  algorithms, Abschnitt 4.5.4, Algorithmus C, Seite 371.]
(defun tryfactor3 (n)
  (let* ((s (isqrt n)) 
         (x (1+ (* 2 s)))
         (y 1)
         (r (- (* s s) n)) )
    #'(lambda ()
        (dotimes (it 500)
          (cond ((> r 0) (decf r y) (incf y 2))
                ((< r 0) (incf r x) (incf x 2))
                (t (assert (zerop r)) (return (exquo (- x y) 2)))
) )   ) ) )

#|
; vierter Faktorisierungsalgorithmus: Siebmethode.
; [Donald E. Knuth, The art of computer programming, Band 2, Seminumerical
;  algorithms, Abschnitt 4.5.4, Algorithmus D, Seite 373.]
(defun tryfactor4 (n)
  (let* ((x (isqrt n))
         (mlist '#(3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61))
         (r (length mlist))
|#

; f�nfter Faktorisierungsalgorithmus: Kettenbruch.
; Auch Algorithmus von Kraitchik - Morrison - Brillhart genannt.
; [Donald E. Knuth, The art of computer programming, Band 2, Seminumerical
;  algorithms, Abschnitt 4.5.4, Algorithmus E, Seiten 380-384.]
; [Carl Pomerance: The quadratic sieve factoring algorithm, in:
;  Lecture Notes in Computer Science, Band 209, Advances in Cryptology,
;  Seiten 169-173.]
; [Marvin C. Wunderlich: A running time analysis of Brillhart's continued
;  fraction factoring method, in:
;  Lecture Notes in Mathematics, Band 751, Seiten 328-342.]
(defun tryfactor5 (n)
  (let* ((m (floor (exp (* 0.25 (sqrt (* (log n) (log (log n))))))))
         ; man kann den Faktor 0.25 auch vergr”žern, wenn der Platz reicht.
         (factorbase (make-array (list m) :element-type 'integer))
                       ; wird ein Array von m Primzahlen sein
         (eqnlist nil) ; wird eine Liste von Gleichungen sein, wobei die Array-
                       ; Komponenten in Gaužscher Treppennormalform �ber
                       ; GF(2) sind.
        )
    ; (assert (>= m 1))
    (flet ((base-factor (v) ; stellt fest, ob abs(v) ganz �ber der Faktorbasis
             ; zerf„llt. Wenn ja, also v = (-1)^e0 * ... * pm^em * v1^2, kommt
             ; ( #(e0 ... em) . v1) als Ergebnis.
             (block base-factoring
               (let ((a (make-array (list (1+ m)) :initial-element 0)))
                 (assert (not (zerop v)))
                 (when (minusp v) ; Ziehe Faktor (-1) heraus
                   (setq v (- v))
                   (setf (aref a 0) 1)
                 )
                 (dotimes (i m) ; Ziehe kleine Faktoren heraus
                   (let ((p (aref factorbase i)))
                     (do ()
                       ((cond ((zerop (mod v p))
                               (setq v (exquo v p))
                               (incf (aref a (1+ i)))
                               nil)
                              (t)
                 ) ) ) ))
                 (unless (= v 1) ; Ist der Rest ein Quadrat?
                   (let ((v2 (isqrt v)))
                     (if (= (* v2 v2) v)
                         (setq v v2)
                         (return-from base-factoring nil) ; nein -> wertlos
                 ) ) )
                 (cons a v)
           ) ) )
           (add-eqn (eqn) ; f�gt eine neue Gleichung eqn = (#(e0..em) x . v)
                          ; die x^2 = (-1)^e0 ... pm^em v^2 mod n bedeutet,
                          ; zum Wissen in eqnlist hinzu.
             ; Wird ein nichttrivialer Faktor von n gefunden, erfolgt ein THROW.
             (block add-eqn
               (do ((l1 nil)
                    (l2 eqnlist)
                    h ha a)
                 ; stets (append (reverse l1) l2) = eqnlist
                 ((endp l2) ; Liste eqnlist zu Ende untersucht.
                  (if ; Sind in eqn alle Exponenten gerade?
                      (some #'oddp (car eqn)) ; Nein: muž eqn behalten
                      (setq eqnlist (nreconc l1 (list eqn)))
                      ; Ja! Ist die Gleichung x^2 = y^2 mod n trivial ?
                      (let ((x (second eqn))
                            (y (let ((v (cddr eqn)))
                                 (dotimes (i m)
                                   (setq v (* v
                                     (expt (aref factorbase i)
                                           (exquo (aref (car eqn) (1+ i)) 2)
                                 ) )))
                                 v
                           ))  )
                        (assert (zerop (mod (- (* x x) (* y y)) n)))
                        (unless (or (zerop (mod (- x y) n))
                                    (zerop (mod (+ x y) n)))
                          ; Nein! Fertig!
                          (throw 'tryfactoring5 (gcd (- x y) n))
                 ))   ) )
                 (setq h (car l2)) ; Hilfsgleichung aus der Liste
                 ; Falls h und eqn an derselben Stelle ihre erste ungerade Zahl
                 ; haben, wird h zu eqn multipliziert:
                 (setq ha (car h)) (setq a (car eqn))
                 (if (do ((i 0 (1+ i)))
                         ((> i m) nil)
                       (cond ((and (oddp (aref ha i)) (evenp (aref a i)))
                              (return nil))
                             ((and (evenp (aref ha i)) (oddp (aref a i)))
                              ; muž eqn vor h einf�gen
                              (setq eqnlist (nreconc l1 (cons eqn l2)))
                              (return-from add-eqn nil))
                             ((and (oddp (aref ha i)) (oddp (aref a i)))
                              (return t))
                     ) )
                     (progn
                       (dotimes (i (1+ m))
                         (incf (aref a i) (aref ha i))
                       )
                       (setq eqn (cons a
                                 (cons (mod (* (second eqn) (second h)) n)
                                       (mod (* (cddr eqn) (cddr h)) n)
                     ) )         ))
                 )
                 (pop l2)
                 (push h l1)
           ) ) )
           (calc-factorbase (D) ; berechnet zu D=k*n eine passende Faktorbasis
             (do ((i 0)
                  (p 2 (1+ p)))
                 ((>= i m))
               (when (or (= p 2) (and (klprimtest p) (= (jacobi D p) 1)))
                     (setf (aref factorbase i) p)
                     (incf i)
           ) ) )
          )
      (let ((D N) ; D=k*n
            (Anfangszustand t)
            R R1 U U1 V V1 P P1 A S   ; Variablenbenennung wie in [Knuth]
            ; Mit den Bezeichnungen wie in [Wunderlich]:
            ; Es gilt f�r alle i=0,1,...:
            ; sqrt(D) = [q0,...,qi-1,(sqrt(D)+Pi)/Qi]
            ; qi = floor((sqrt(D)+Pi)/Qi) = floor((floor(sqrt(D))+Pi)/Qi),
            ; Q-1 = D - q0^2, P0=0, Q0=1, Pi+1 = qi * Qi - Pi,
            ; Qi+1 = (D - Pi+1)^2 / Qi = Qi-1 - qi * ( qi * Qi - 2 * Pi ),
            ; alle abs(Pi)<sqrt(D), alle Qi>0, alle qi>0, alle Variablen ganz.
            ; A-1=1, B-1=0, A0=q0, B0=1,
            ; Ai = qi * Ai-1 + Ai-2, Bi = qi * Bi-1 + Bi-2 f�r i>0,
            ; in Q(X) gilt (Ai+x*Ai-1)/(Bi+x*Bi-1) = [q0,...,qi-1,qi + x].
            ; F�r i>=0 ist Ai-1^2 - D * Bi-1^2 = (-1)^i * Qi
            ; F�r i>=0 ist Ai-1 * Ai - D * Bi-1 * Bi = (-1)^i * Pi+1
            ; Nach i Schleifendurchl„ufen ist
            ; R = floor(sqrt(D))
            ; R1 = 2*floor(sqrt(D))
            ; U = floor(sqrt(D)) + Pi+1
            ; U1 = floor(sqrt(D)) + Pi
            ; V = Qi+1
            ; V1 = Qi
            ; P = Ai mod N
            ; P1 = Ai-1 mod N
            ; A = qi
            ; S = (-1)^(i+1)
           )
        #'(lambda ()
            (catch 'tryfactoring5
              (do ((it 200) h1) ; Iterationenzahl und Hilfsvariable
                  ((<= it 0))
                (when Anfangszustand ; muž initialisieren
                  (calc-factorbase D) ; Faktorbasis berechnen.
                  (dotimes (i m)      ; Test auf Teilbarkeit in Faktorbasis.
                    (if (zerop (mod n (aref factorbase i)))
                        ; muž ein echter Teiler sein, weil n nicht prim
                        (throw 'tryfactoring5 (aref factorbase i))
                  ) )
                  (setq R (isqrt D))
                  (setq R1 (* 2 R))
                  (setq U R1)
                  (setq U1 R)
                  (setq V (- D (* R R)))
                  (if (zerop V)       ; k*N=D ist Quadratzahl = R^2
                      ; -> g:=ggT(R,N) teilt N.
                      ; Bei g=1 w„re ggT(R^2,N)=1 => R^2|k => N=1, Widerspruch.
                      ; Bei g=N w„re N|R => N^2|R^2=k*N => N|k, das passiert
                      ; zu unseren Lebzeiten nicht mehr.
                      (throw 'tryfactoring5 (gcd R n)) )
                  (setq V1 1)
                  (setq P R)
                  (setq P1 1)
                  (setq A R)
                  (setq S -1)
                  (setq Anfangszustand nil)
                  (decf it 2)
                )
                ; Hier gelten die ganzen Invarianten.
                ; (assert (zerop (mod (+ (* P1 P1) (* S V1)) n)))
                ; (assert (zerop (mod (- (* P P) (* S V)) n)))
                (if (setq h1 (base-factor (* S V)))
                    (add-eqn (cons (car h1) (cons P (cdr h1)))) )
                (when (= V 1) ; Periode des Kettenbruches erreicht ->
                  (incf D n)
                  (setq Anfangszustand t) ; neuer Anlauf mit gr”žerem k.
                )
                ; n„chste Kettenbruchiteration:
                (multiple-value-setq (A h1) (floor U V)) ; A:=qi+1
                (psetq  U1 U  U (- R1 h1)) ; U1:=R+Pi+1, U:=R+Pi+2
                (psetq  V1 V  V (- V1 (* A (- U U1)))) ; V1:=Qi+1, V:=Qi+2
                (setq S (- S)) ; S:=(-1)^(i+2)
                (psetq  P1 P  P (mod (+ (* A P) P1) n))
                  ; P1:=Ai mod N, P:=Ai+1 mod N
                ; jetzt ist ein Durchlauf fertig, i:=i+1
                (decf it)
) ) ) )   ) ) )

#|
; sechster Faktorisierungsalgorithmus: Quadratisches Sieben.
|#


; (cyclic-list x1 ... xn) mit n>0
; liefert die zyklische Liste #1=(x1 ... xn . #1#)
(defun cyclic-list (&rest args)
  (setf (cdr (last args)) args)
  args
)


; (verffz b) verfeinert die in b steckenden Faktorzerlegung von n = (fz-num b)
; mit bis zu k Iterationsschritten (k fehlt -> beliebig viele Schritte), solange
; bis eine bestimmte Bedingung (return-test b) zutrifft ( /= NIL liefert).
; Ergebnis ist die Anzahl der verbrauchten Iterationen.
(defun verffz (b &optional (k nearly-infinite)
                 &key (return-test #'(lambda (x)
                                       (declare (ignore x))
                                       nil) )
                 &aux (it 0))
   (if (eq (fz-cl b) 'prime) (return-from verffz it)) ; nichts zu tun.
   (if (eq (fz-cl b) 'probprime)
     (let ((n (fz-num b)))
        ; Jetzt muž der Primzahlverdacht erh„rtet werden.
        ; Daž n>2 ungerade ist, ist schon bekannt.
        (if (null (fz-other b)) ; noch keine wartenden Prozeduren
            (setf (fz-other b)
                  (cyclic-list
                    (tryprimeprover1 n)
                    (tryprimeprover2 n)
                    ; usw. (weitere Beweisversucher)
        )   )     )
        ; Jetzt m�ssen die wartenden Funktionen der Reihe nach aktiviert
        ; werden, bis sie ein Ergebnis liefern.
        (incf it (tryprimeprove b k))
        (return-from verffz it)
   ) )
   (if (eq (fz-cl b) 'notprime)
       ; erst einmal muž ein Faktor von n gefunden werden:
       (let ((n (fz-num b)))
         (if (null (fz-other b))
             (setf (fz-other b)
                   (cyclic-list
                     (tryfactor1 n)
                     (tryfactor2 n)
                     (tryfactor3 n)
                     ; usw.
         )   )     )
         (incf it (tryfactor b k))
   )   )
   (if (eq (fz-cl b) 'notprime)
       ; konnte noch nicht faktorisieren
       (return-from verffz it))
   (assert (eq (fz-cl b) 'factored))
   (if (fzl-recalc (fz-other b))
       (return-from verffz it)) ; Verfeinern unm”glich
   ; Jetzt wird die Anstrengung auf die einzelnen Faktoren verteilt.
   (let ((l nil) f) ; l = Liste noch zu verfeinernder Faktoren
     (loop
        (when (or (>= it k) (funcall return-test b)) (return))
        (if (null l) (setq l (fzl-factorlist (fz-other b)))) ; von neuem
        (setq f (first l))
        (if (not (integerp f))
          (progn
            (incf it (verffz f                 ; rekursiv verfeinern
                             (ceiling (max 0 (- k it)) (length l))))
            (when (eq (fz-cl f) 'factored)
                  ; in b wird f durch seine Faktoren ersetzt
                  (setf (fzl-factorlist (fz-other b))
                    (mergefz (remove f (fzl-factorlist (fz-other b))
                                     :count 1 :test #'eq)
                             (copy-list (fzl-factorlist (fz-other f)))
            )     ) )
                  ; (beachte: l bleibt hierbei die Liste der noch nicht
                  ; abgearbeiteten Faktoren, der Schwanz von
                  ; (fzl-factorlist (fz-other b)). )
            (when (member (fz-cl f) '(prime factored))
                  (if (fzl-recalc (fz-other b)) ; noch was zu tun?
                      (return-from verffz it))
            )
        ) )
        (pop l)
     )
     ; (fzl-factorlist (fz-other b)) ist jetzt zur Gen�ge verfeinert.
     it
) )


; (pfz n) liefert die Primfaktorzerlegung von n>0:
; alle Primfaktoren, in aufsteigender Reihenfolge, evtl. mehrfach vertreten.
(defun pfz (n)
  (if (< n *Kleinheit*)
      (klpfz n)
      (let ((b (pfz-table-lookup n)))
        (when (null b) ; n noch nicht in der Tabelle?
              (mr n) ; n in die Tabelle eintragen
              (setq b (pfz-table-lookup n))
        )
        (do ((i 5 (* 2 i))) ; i w„chst stark
            ((or (eq (fz-cl b) 'prime)
                 (and (eq (fz-cl b) 'factored) (fzl-allprimes (fz-other b)))
             ))
          (verffz b i) ; ewig verfeinern
        )
        (if (eq (fz-cl b) 'prime)
            (list n)
            (mapcar #'fz-key (fzl-factorlist (fz-other b)))
) )   ) )


; (setpfz n l) speichert in der Tabelle ein (von irgendwoher bekanntes)
; Primfaktorzerlegungsergebnis ab. Es sollte irgendwann l als Ergebnis von
; (pfz n) gekommen sein.
; Nach (setpfz n l) ist sichergestellt, daž in Zukunft bei (pfz n) ohne
; Rechnung sofort das Ergebnis l kommt.
(defun setpfz (n l)
  (if ; Plausibilit„tspr�fung, ob (pfz n) = l sein kann.
      (and (apply #'<= 2 l)
           (every #'(lambda (f) (or (>= f *Kleinheit*) (klprimtest f))) l)
           (= n (apply #'* l))
      )
      (progn
        ; erst die Primfaktoren als Primzahlen in die Tabelle abspeichern
        (setq l
          (mapcar
            #'(lambda (f)
                (if (< f *Kleinheit*)
                    f
                    (let ((b (pfz-table-lookup-insert f)))
                      (unless (and (eq (fz-cl b) 'prime)
                                   (integerp (fz-other b)))
                        (setf (fz-cl b) 'prime)
                        (setf (fz-other b) nil)
                      )
                      b
              ) )   )
            l
        ) )
        ; l ist die neue Faktorzerlegungsliste
        (when (and (>= n *Kleinheit*) (> (length l) 1))
          (let ((b (pfz-table-lookup-insert n)))
            (setf (fz-cl b) 'factored)
            (setf (fz-other b) (make-fzl :allprimes t :factorlist l))
        ) )
        n ; als Ergebnis
      )
      (error "Fehler: ~D hat nicht die Primfaktorzerlegung ~S.~%" n l)
) )


;-------------------------------------------------------------------------------
; Anwendungen (benutzen nur pfz):


; (pfzv n) ergibt die Primfaktorzerlegung von n>0 in der Gestalt
; ((p1 . e1) (p2 . e2) ... (pk . ek)) mit Primzahlen p1<...<pk und
; Exponenten e1,...,ek>0 mit n=p1^e1 ... pk^ek.
(defun pfzv (n)
  (do ((l1 (pfz n))
       (l2 nil))
      ((endp l1) (reverse l2))
    (let* ((p (first l1))
           (e (do ((i 0 (1+ i)))
                  ((or (null l1) (not (= p (first l1)))) i)
                (pop l1)
          ))  )
       (push (cons p e) l2)
) ) )


; (primteiler n) liefert die Primteiler von n>0, jeden nur einmal.
(defun primteiler (n) (mapcar #'car (pfzv n)))


; (pdivisors n) ergibt die positiven Teiler einer Zahl n>0, in
; irgendeiner Reihenfolge
(defun pdivisors (n)
  (reduce
    #'(lambda (l pe)          ; ergibt l, elementweise mit 1,p,...,p^e
        (let ((p (first pe))  ; multipliziert.
              (e (rest pe)))
          (mapcan
            #'(lambda (m)     ; ergibt die mit m multiplizierte Liste l
                (mapcar #'(lambda (x) (* m x)) l)
              )
            (mapcar #'(lambda (i) (expt p i)) (intlist 0 e))
      ) ) )
    (pfzv n)
    :initial-value '(1)
  )
)


; (divisors n) ergibt die ganzzahligen Teiler einer Zahl n /= 0
(defun divisors (n)
  (let ((pdiv (pdivisors (abs n))))
    (append pdiv (mapcar #'- pdiv))
) )


; (euler-phi n) ergibt f�r n>0 die Eulersche phi-Funktion von n,
; d.h. die Anzahl der zu n teilerfremden nat�rlichen Zahlen >0, <=n.
(defun euler-phi (n)
  (reduce #'*
          (mapcar #'(lambda (pe)
                      (let ((p (car pe)) ; Primzahl
                            (e (cdr pe)) ; Exponent, >0
                           )
                        (* (expt p (1- e)) (1- p))
                    ) )
                  (pfzv n)
) )       )


; (prim-von-bis a b) ergibt f�r eine Liste aller Primzahlen von a bis b
; (inklusive).
(defun prim-von-bis (a b)
  (remove-if-not #'isprime (intlist a b))
)