home *** CD-ROM | disk | FTP | other *** search
/ GEMini Atari / GEMini_Atari_CD-ROM_Walnut_Creek_December_1993.iso / files / language / examples / primzahl.lsp < prev    next >
Lisp/Scheme  |  1993-10-23  |  41KB  |  1,038 lines
  1. ; Bruno Haible 10.11.1987-12.12.1987, 25.2.1990
  2. ; Primfaktorzerleger für ganze Zahlen
  3.  
  4. (provide 'primzahl)
  5. ;(in-package 'primzahl)
  6. ;(export '(primteiler pfzv pdivisors divisors pfz prim-von-bis
  7. ;          *Kleinheit* prime probprime notprime factored *pfz-table* 
  8. ;          isprobprime isprime verffz setpfz))
  9. (require 'avl) ; AVL-Bäume
  10.  
  11.  
  12. (setq *print-circle* t) ; wegen der zyklischen Listen in fz-other
  13.  
  14.  
  15. ; (intlist a b) ergibt (a a+1 ... b-1 b), wenn a und b integers sind.
  16. (defun intlist (a b)
  17.   (do ((l nil (cons i l))
  18.        (i b (1- i)))
  19.       ((< i a) l)
  20. ) )
  21.  
  22.  
  23. ; exakter Quotient von Integers, schneller als / :
  24. #-CLISP
  25. (defun exquo (a b)
  26.   (multiple-value-bind (q r) (floor a b)
  27.     (unless (zerop r) (error "Quotient ~S/~S nicht exakt." a b))
  28.     q
  29. ) )
  30.  
  31.  
  32. (require 'jacobi)
  33. ; (jacobi a b) liefert für ganze Zahlen a,b mit ungeradem b>0 das Jacobi-Symbol
  34. ; (a / b). Es ist =0 genau dann, wenn a und b nicht teilerfremd sind.
  35.  
  36.  
  37. ; (pfz n) liefert die Primfaktorzerlegung von n>0, die Primfaktoren von n
  38. ; in aufsteigender Reihenfolge, evtl. mehrfach vertreten.
  39.  
  40. #|
  41. ; vorläufige Version des Primfaktorzerlegers
  42. (defun pfz (n)
  43.   (reverse
  44.     (let ((l nil))         ; l = Liste der bereits gefundenen Faktoren
  45.       (do () ((oddp n))
  46.         (setq n (exquo n 2))
  47.         (push 2 l))
  48.       (let ((p 3))         ; p = Test-teiler (ungerade)
  49.         (loop
  50.           (if (= n 1) (return l))
  51.           (setq s (isqrt n))
  52.           (do () ((or (> p s) (zerop (mod n p))))
  53.                  (incf p 2))
  54.           (if (> p s) (return (cons n l)))
  55.           (push p l)
  56.           (setq n (exquo n p))
  57. ) ) ) ) )
  58. |#
  59.  
  60.  
  61. ; (Miller-Rabin-test N) testet, ob eine gegebene Zahl n>1 wohl eine
  62. ; Primzahl ist.
  63. ; Falls das Ergebnis NIL ist, ist N garantiert zusammengesetzt,
  64. ; Eventuell kommt ein nichttrivialer Teiler von N als zweiter Wert heraus.
  65. ; Falls das Ergebnis T ist, ist N mit hoher Wahrscheinlichkeit prim.
  66. ; (Fehlerwahrscheinlichkeit < 1E-30 )
  67. ; Ist n gerade und >2, so ist das Ergebnis NIL.
  68. (defun miller-rabin-test (n)
  69.   ; Vorausschleife, die bereits 87% aller Zahlen faktorisiert
  70.   (dolist (p '(2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67))
  71.     (if (zerop (mod n p))
  72.         (return-from miller-rabin-test (if (= n p) t (values nil p)))))
  73.   ; erst jetzt geht's richtig los:
  74.   (let (blist n1 s t1)
  75.     (if (< n 2) (return-from miller-rabin-test nil))
  76.     (setq blist (cond ((< n 2000) '(2))
  77.                       ((< n 1300000) '(2 3))
  78.                       ((< n 25000000) '(2 3 5))
  79.                       ((< n 3200000000) '(2 3 5 7))
  80.                       (t '(2 3 5 7 11 13 17 19 23 29
  81.                            31 37 41 43 47 53 59 61 67 71
  82.                            73 79 83 89 97 101 103 107 109 113
  83.                            127 131 137 139 149 151 157 163 167 173
  84.                            179 181 191 193 197 199 211 223 227 229))))
  85.     ; mit 50 Primzahlen: P(falsch) < (1/4)^50 < 1E-30
  86.     (setq n1 (- n 1))
  87.     (setq s n1  t1 0)
  88.     (do () ((oddp s))
  89.            (setq s (exquo s 2)) (incf t1))
  90.     ; n-1 = s * 2^t1 mit ungeradem s zerlegt.
  91.     (dolist (b blist)
  92.       (if (zerop (mod n b))
  93.         (if (= n b) (return-from miller-rabin-test t)
  94.                     (return-from miller-rabin-test (values nil b))))
  95.       (do ((c (exptmod (mod b n) s n) (sqrmod c n))
  96.            (c1 n1 c)
  97.            (j 0 (1+ j)))
  98.           ((or (= j t1) (= c 1))
  99.            (cond ((not (= c 1)) ; sicher j=t1, aber c/=1
  100.                   (return-from miller-rabin-test nil))
  101.                  ((not (= c1 n1)) ; sicher c=1, c1/=n-1, also j>0
  102.                   ; ggt(c1-1,n) * ggt(c1+1,n) = ggt(c1^2-1,n) = ggt(c-1,n) = n
  103.                   (return-from miller-rabin-test (values nil (gcd n (- c1 1)))))
  104.       )   )) ; Eine Iteration über b beendet
  105.     )   ; blist abgearbeitet
  106.     t   ; n wahrscheinlich prim
  107. ) )
  108.  
  109. ; Hilfsfunktion für Miller-Rabin-Test:
  110. ; Quadrieren modulo n von x mit 0<=x<n
  111. (defun sqrmod (x n)
  112.   (mod (* x x) n))
  113.  
  114. ; Hilfsfunktion für Miller-Rabin-Test:
  115. ; Potenzieren modulo n von a mit 0<=a<n und b>0
  116. (defun exptmod (a b n)
  117.   (let ((c 1))
  118.     (loop    ; a^b*c mod n bleibt invariant
  119.       (if (oddp b) (setq c (mod (* a c) n)))
  120.       (if (= b 1) (return c))
  121.       (setq a (mod (* a a) n) b (floor b 2))
  122. ) ) )
  123.  
  124.  
  125. ; verifizierter Primzahltest für kleine Zahlen n>0
  126. (defun klprimtest (n)
  127.   (cond ((= n 1) nil)
  128.         ((evenp n) (= n 2))
  129.         (t (let ((s (isqrt n)))
  130.              (do ((i 3 (+ i 2)))
  131.                  ((> i s) t)
  132.                (if (zerop (mod n i)) (return nil))
  133. ) )     )  ) )
  134.  
  135.  
  136. ; Primfaktorzerlegung kleiner Zahlen n>0
  137. (defun klpfz (n)
  138.   (let ((l nil))
  139.     (do () ((oddp n)) (setq n (ash n -1)) (push 2 l)) ; Zweier heraus
  140.     (block Rest-prim
  141.       (do ((p 3))
  142.           ((= n 1))
  143.         (loop ; in dieser Schleife wird ein Teiler gesucht
  144.           (cond ((> (* p p) n) (return-from Rest-prim (push n l)))
  145.                 ((zerop (mod n p)) (push p l) (setq n (exquo n p)) (return))
  146.                 (t (incf p 2))
  147.     ) ) ) )
  148.     (nreverse l)
  149. ) )
  150.  
  151.  
  152. ; Datenstruktur: In *pfz-table* steht eine Tabelle von Primfaktorzerlegungs-
  153. ; ergebnissen. Zu jeder Zahl (sie ist der Schlüssel) steht drin:
  154. ; Ob sie Primzahl ist, eine eventuelle Faktorzerlegung bzw. evtl. eine
  155. ; Primitivwurzel.
  156.  
  157. (defparameter *Kleinheit* 10000
  158.   "Nur Zahlen >= *Kleinheit* werden in die PFZ-Tabelle aufgenommen.")
  159. ; *Kleinheit* sollte nur geändert werden, wenn gleichzeitig
  160. ; (setq *pfz-table* nil) ausgeführt wird!
  161. (assert (> *Kleinheit* 1))
  162.  
  163. (deftype fz-class () "Das Ergebnis einer Faktorzerlegung, ganz grob"
  164.   '(member              ; ein Symbol:
  165.             nil        ; (nur kurz nach der Initialisierung)
  166.             prime      ; n ist prim
  167.             probprime  ; n ist wahrscheinlich prim
  168.             notprime   ; n ist garantiert zusammengesetzt, Faktoren unbekannt
  169.             factored   ; n ist zusammengesetzt, schon faktorisiert
  170. )  )
  171.  
  172. (defstruct fz     "Eine Faktorisierung als Ganzes"
  173.   (num 2   :type integer :read-only t)  ; die zu faktorisierende Zahl, >1
  174.   (cl  nil :type fz-class) ; Klassifikation
  175.   (other nil) ; weitere Information:
  176.      ; bei cl = prime : NIL oder eine Primitivwurzel mod n
  177.      ; bei cl = probprime : Eine zyklische Liste von wartenden Funktionen.
  178.      ;                      Jede von ihnen liefert NIL, wenn sie wieder aufge-
  179.      ;                      rufen werden will, und (NIL . factor), falls sie
  180.      ;                      die Zahl als zusammengesetzt nachgewiesen hat, evtl.
  181.      ;                      einen nichttrivialen Faktor gefunden hat, und
  182.      ;                      (T . Prw), falls sie die Zahl als prim nachgewiesen
  183.      ;                      hat, mit evtl. einer Primitivwurzel Prw.
  184.      ;                      (Eventuell zu Beginn die leere Liste.)
  185.      ; bei cl = notprime : Eine zyklische Liste von wartenden Funktionen
  186.      ;                     (Closures), die mitten in der Suche nach einem
  187.      ;                     Faktor von n stecken. Jede von ihnen liefert NIL,
  188.      ;                     wenn sie wieder aufgerufen werden will, und einen
  189.      ;                     nichttrivialen Faktor, wenn sie einen gefunden hat.
  190.      ;                     (Eventuell zu Beginn die leere Liste.)
  191.      ; bei cl = factored : Ein Konstrukt vom Typ fzl, das die Faktorenzerlegung
  192.      ;                     von n widerspiegelt.
  193. )
  194.  
  195. (defstruct fzl "Faktorenzerlegung, Liste"
  196.   (allprimes nil  :type symbol) ; Flag, ob alle Faktoren von factorlist
  197.                                 ; bekanntermaßen Primzahlen sind. 
  198.   (factorlist nil :type list) ; Eine Liste (f1 ... fk) von Faktoren von n,
  199.         ; mit 1 < f1 <= ... <= fk und n = f1 * ... * fk. Die fj sind dabei
  200.         ; (Pointer auf) Konstrukte vom Typ fz, bzw. bei fj<*Kleinheit* ist fj
  201.         ; die Zahl selbst und Primzahl.
  202. )
  203.  
  204.  
  205. ; (fz-key f) ergibt zu einem Element einer Faktorzerlegungsliste die
  206. ; wirkliche Zahl.
  207. (defun fz-key (g)
  208.   (if (integerp g) g (fz-num g))
  209. )
  210.  
  211.  
  212. ; (fz-isprime f) stellt fest, ob ein Element einer Faktorzerlegungsliste
  213. ; prim ist.
  214. (defun fz-isprime (f)
  215.   (or (integerp f) (eq (fz-cl f) 'prime))
  216. )
  217.  
  218.  
  219. ; (mergefz l1 l2) mischt zwei Faktorzerlegungen l1 und l2 von n1 bzw. n2,
  220. ; so daß eine Faktorzerlegung von n1*n2 herauskommt. Vorsicht: Destruktiv!
  221. (defun mergefz (l1 l2)
  222.   (merge 'list l1 l2 #'< :key #'fz-key)
  223. )
  224.  
  225.  
  226. ; (fzl-recalc h) berechnet neu, ob alle Faktoren von h (h ist vom Typ fzl)
  227. ; prim sind, trägt diese Information in h ein und liefert sie als Ergebnis.
  228. (defun fzl-recalc (h)
  229.   (setf (fzl-allprimes h)
  230.     (reduce #'(lambda (allp f) (and allp (fz-isprime f)))
  231.             (fzl-factorlist h)
  232.             :initial-value t
  233. ) ) )
  234.   
  235.  
  236. (defun fz-comp (x y) (< (fz-num x) (fz-num y)))
  237.  
  238. (defun fz-eq-test (x y) (= (fz-num x) (fz-num y)))
  239.  
  240. (defvar *pfz-table* (avl:seq-to-avl nil #'fz-comp)
  241.   "Die globale Tabelle, die alle Primfaktorzerlegungsergebnisse enthält")
  242. ; Die in *pfz-table* abgespeicherten Werte sind alle vom Typ fz.
  243. ; Die Ordnungsrelation ist durch die Zahl selbst (fz-num *) gegeben.
  244.  
  245.  
  246. ; Ausgeben der globalen Tabelle:
  247. (defun pt (&optional (output-stream *standard-output*))
  248.   (pprint (avl:avl-to-seq *pfz-table*) output-stream)
  249. )
  250.  
  251.  
  252. (defun pfz-table-lookup (n)
  253. ; ergibt die in *pfz-table* stehende Faktorzerlegung von n
  254. ; 1. Version (langsamer):
  255. ; (avl:member (make-fz :num n) *pfz-table* #'fz-comp
  256. ;   #'(lambda (x y) (and (fz-eq-test x y) y)) ))
  257. ; 2. Version (schneller, aber weniger portabel):
  258.   (do ((tr *pfz-table*))
  259.       ((or (null tr) (= n (fz-num (avl::node-value tr))))
  260.        (cond (tr (avl::node-value tr))))
  261.     (cond ((< n (fz-num (avl::node-value tr))) (setq tr (avl::node-left tr)))
  262.           ((> n (fz-num (avl::node-value tr))) (setq tr (avl::node-right tr)))
  263.           (t (error "Unvergleichbarkeit zweier Zahlen!"))
  264. ) ) )
  265.  
  266.  
  267. ; (pfz-table-insert b) fügt zu *pfz-table* einen neuen Eintrag b (vom Typ fz)
  268. ; hinzu. Ergebnis ist b.
  269. (defun pfz-table-insert (b)
  270.   (setq *pfz-table* (avl:insert b *pfz-table* #'fz-comp #'fz-eq-test))
  271.   b
  272. )
  273.  
  274.  
  275. ; (pfz-table-lookup-insert n) sucht in *pfz-table* nach einem Eintrag zur Zahl
  276. ; n. Wird kein Eintrag gefunden, so wird ein neuer in die Tabelle eingefügt.
  277. ; Das Ergebnis ist stets der entsprechende Eintrag in der Tabelle *pfz-table*.
  278. (defun pfz-table-lookup-insert (n)
  279.   (assert (>= n *Kleinheit*))
  280.   (let ((b (pfz-table-lookup n)))
  281.     (unless b (pfz-table-insert (setq b (make-fz :num n))))
  282.     b
  283. ) )
  284.  
  285.  
  286. ; (recall-factors n) ergibt eine Liste aller Faktoren der natürlichen Zahl n>1,
  287. ; soweit sie bekannt sind, in einer sortierten Liste, im selben Format wie bei
  288. ; fzl. Der zweite Wert zeigt an, ob alles sicher Primzahlen sind.
  289. (defun recall-factors (n)
  290.   (if (< n *Kleinheit*)
  291.       (values (klpfz n) t)
  292.       (let ((b (pfz-table-lookup n)))
  293.         (if (null b)
  294.             ; nicht gefunden
  295.             (progn
  296.               (mr n)
  297.               (values (list (pfz-table-lookup n)) nil)
  298.             )
  299.             ; gefunden
  300.             (case (fz-cl b)
  301.               ((prime) (values (list b) t))
  302.               ((probprime notprime) (values (list b) nil))
  303.               ((factored)
  304.                (setq b (fz-other b))
  305.                (values (fzl-factorlist b) (fzl-allprimes b))
  306.               )
  307.               (t (error "Falsch aufgebauter fz-Struct!"))
  308. ) )   ) )   )
  309.  
  310.  
  311. ; (build-fzl l) ergibt die Faktorzerlegung vom Typ fzl von n, wenn bekannt ist,
  312. ; daß das Produkt aller Elemente von l (alles natürliche Zahlen >1) = n ist.
  313. (defun build-fzl (l)
  314.   (let ((l1 nil) (allp t))
  315.     (dolist (f l) (multiple-value-bind (f1 f2) (recall-factors f)
  316.                     (setq l1 (mergefz l1 (copy-list f1)))
  317.                     (setq allp (and allp f2))
  318.     )             )
  319.     (make-fzl :allprimes allp :factorlist l1)
  320. ) )
  321.  
  322.  
  323. ; Führt für eine Zahl >= *Kleinheit*, über die noch nichts bekannt ist,
  324. ; den Miller-Rabin-Test durch und speichert das Ergebnis in *pfz-table* ab.
  325. ; Bei Ergebnis t (also n wohl prim) ist n=2 oder n>2 ungerade.
  326. (defun mr (n)
  327.   (multiple-value-bind (pr fct) (miller-rabin-test n)
  328.     (cond (pr ; wohl prim
  329.               (pfz-table-insert (make-fz :num n :cl 'probprime))
  330.               t)
  331.           ; sonst sicher nicht prim
  332.           (fct ; Faktor gefunden
  333.               (pfz-table-insert
  334.                 (make-fz :num n :cl 'factored
  335.                          :other (build-fzl (list fct (exquo n fct)))))
  336.               nil)
  337.           (t ; kein Faktor gefunden
  338.               (pfz-table-insert (make-fz :num n :cl 'notprime))
  339.               nil)
  340. ) ) )
  341.  
  342.  
  343. ; (isprobprime n) stellt fest, ob die natürliche Zahl n wahrscheinlich eine
  344. ; Primzahl ist.
  345. ; Bei geraden Zahlen >2 liefert dies stets NIL.
  346. (defun isprobprime (n)
  347.   (if (< n *Kleinheit*)
  348.       (klprimtest n)
  349.       (let ((b (pfz-table-lookup n)))
  350.         (if b (member (fz-cl b) '(prime probprime)) ; bekannt
  351.               ; sonst muß gerechnet werden:
  352.               (mr n)
  353. ) )   ) )
  354.  
  355.  
  356. ; (isprime n) stellt fest, ob die natürliche Zahl n eine Primzahl ist.
  357. (defun isprime (n)
  358.   (if (< n *Kleinheit*)
  359.       (klprimtest n)
  360.       (let ((b (pfz-table-lookup n)))
  361.         (if (and b (member (fz-cl b) '(prime))) ; bekannt
  362.             (return-from isprime t))
  363.         (if (and b (member (fz-cl b) '(notprime factored))) ; bekannt
  364.             (return-from isprime nil))
  365.         (assert (or (null b) (eq (fz-cl b) 'probprime)))
  366.         (if (and (null b) (null (mr n))) (return-from isprime nil))
  367.         ; jetzt ist bekannt, daß n wohl prim ist, steht in *pfz-table*.
  368.         (if (null b) (setq b (pfz-table-lookup n)))
  369.         (assert (and b (member (fz-cl b) '(prime probprime))))
  370.         (if (eq (fz-cl b) 'prime) (return-from isprime t))
  371.         (verffz b) ; Primzahlverdacht erhärten
  372.         (assert (not (eq (fz-cl b) 'probprime)))
  373.         (eq (fz-cl b) 'prime)
  374. ) )   )
  375.  
  376.  
  377. (defparameter nearly-infinite 1000000000000000000
  378.   "Schier unendliche Zahl von Iterationen")
  379.  
  380.  
  381. ; Sei (fz-cl b) = 'probprime und (fz-other b) /= NIL.
  382. ; (tryprimeprove b) ruft die wartenden Funktionen in (fz-other b) der Reihe
  383. ; nach auf, bis ein Ergebnis kommt. Maximal k Iterationen (k fehlt->grenzenlos).
  384. ; Das Ergebnis ist die Anzahl der verbrauchten Iterationen.
  385. (defun tryprimeprove (b &optional (k nearly-infinite) &aux (it 0) erg)
  386.   (loop
  387.      (when (>= it k) (return it))
  388.      (incf it)
  389.      (when (setq erg (funcall (first (fz-other b))))
  390.        ; Ungewißheit beendet.
  391.        ; erg = (NIL) oder (NIL . factor) oder (T) oder (T . Prw)
  392.        (cond ((car erg) (setf (fz-cl b) 'prime)
  393.                         (setf (fz-other b) (cdr erg)) )
  394.              ((null (cdr erg)) (setf (fz-cl b) 'notprime)
  395.                                (setf (fz-other b) nil) )
  396.              (t (setq erg (cdr erg))
  397.                 (setf (fz-other b) nil)
  398.                 (setf (fz-cl b) 'factored)
  399.                 (setf (fz-other b)
  400.                       (build-fzl (list erg (exquo (fz-num b) erg))))
  401.        )     )
  402.        (return it)
  403.      )
  404.      (pop (fz-other b)) ; Funktionen-Schlange rotieren
  405. ) )
  406.  
  407.  
  408. ; 1. Primzahlbeweiser: Bei jedem Aufruf wird eine gewisse Anzahl von
  409. ; Teilern durchprobiert, maximal bis zur Wurzel.
  410. ; [Allgemein bekannt.]
  411. (defun tryprimeprover1 (n)
  412.   (let ((p 1)
  413.         (s (isqrt n)))
  414.     #'(lambda ()
  415.         (block tryprimeproving1
  416.           (dotimes (i (min 100 (ceiling (- s p) 2)))
  417.             (incf p 2)
  418.             (if (zerop (mod n p))
  419.                 (return-from tryprimeproving1 (cons nil p)))
  420.           )
  421.           (if (>= p s) (cons t nil) nil)
  422. ) )   ) )
  423.  
  424.  
  425. ; 2. Primzahlbeweiser: Aus einer hinreichend guten Faktorzerlegung von n-1
  426. ; wird eine Zahl der Ordnung n-1 mod n bestimmt, also ist n prim.
  427. ; [Donald E. Knuth, The art of computer programming, Band 2, Seminumerical
  428. ; algorithms, Abschnitt 4.5.4, Aufgabe 26, Seite 397.]
  429. (defun tryprimeprover2 (n)
  430.   (let ((enough nil)
  431.         f r plist
  432.         (n1 (1- n)))
  433.     (flet ((factored-enough (b) ; stellt fest, ob n-1 weit genug faktorisiert
  434.              (setq enough                                                 ; ist.
  435.                (and (= (fz-num b) n1)
  436.                     (eq (fz-cl b) 'factored)
  437.                     (progn
  438.                       ; f sammelt die primen Faktoren p von n,
  439.                       ; r sammelt die nichtprimen Faktoren von n,
  440.                       ; plist = Menge aller primen Faktoren von n.
  441.                       (setq f 1)
  442.                       (setq r 1)
  443.                       (setq plist nil)
  444.                       (dolist (p (fzl-factorlist (fz-other b)))
  445.                         (if (fz-isprime p)
  446.                             (progn
  447.                               (setq f (* f (fz-key p)))
  448.                               (setq plist (adjoin (fz-key p) plist))
  449.                             )
  450.                             (setq r (* r (fz-key p)))
  451.                       ) )
  452.                       ; alles ok, wenn r<=f+1
  453.                       (<= r (1+ f))
  454.            ) ) )    )
  455.            (tpp2-rest () ; führt den Primzahlbeweis aus, wenn f,r,plist da.
  456.              ; Bei r=1 könnte man sogar eine Primitivwurzel errechnen.
  457.              ; [D. E. Knuth, Band 2, Abschnitt 4.5.4, Aufgabe 10, Seite 395.]
  458.              (do ((x 2 (1+ x)))
  459.                  ((null plist) (cons t nil))
  460.                (unless (= (exptmod x n1 n) 1)
  461.                  (format t
  462.                    "~%Bei ~D hatte der Miller-Rabin-Test falsch geraten!~%"
  463.                    n)
  464.                  (return (cons nil nil))
  465.                )
  466.                (setq plist
  467.                  (remove-if #'(lambda (p)
  468.                                 (= (gcd (1- (exptmod x (exquo n1 p) n)) n) 1) )
  469.                             plist
  470.           )) ) ) )
  471.       (if (< n1 *Kleinheit*)
  472.           #'(lambda ()
  473.               (setq f n1)
  474.               (setq r 1)
  475.               (setq plist (primteiler n1))
  476.               (tpp2-rest))
  477.           (let ((b (pfz-table-lookup n1)))
  478.             (when (null b)
  479.                   (isprobprime n1) ; n1 in *pfz-table* eintragen
  480.                   (setq b (pfz-table-lookup n1)) )
  481.             #'(lambda ()
  482.                 (verffz b 1 :return-test #'factored-enough)
  483.                 ; Anmerkung: Immer nur eine Zeitscheibe ist wohl ungünstig,
  484.                 ; weil dann stets der kleinste mögliche Faktor von n1
  485.                 ; verfeinert wird.
  486.                 (when (or enough (factored-enough b))
  487.                       (tpp2-rest)
  488.               ) )
  489. ) ) ) )   )
  490.  
  491.  
  492. ; Sei (fz-cl b) = 'notprime und (fz-other b) /= NIL.
  493. ; (tryfactor b) ruft die wartenden Funktionen in (fz-other b) der Reihe nach
  494. ; auf, bis ein Ergebnis kommt. Maximal k Iterationen (k fehlt -> grenzenlos).
  495. ; Ergebnis ist die Anzahl der verbrauchten Iterationen.
  496. (defun tryfactor (b &optional (k nearly-infinite) &aux (it 0) erg)
  497.   (loop
  498.     (when (>= it k) (return it))
  499.     (incf it)
  500.     (when (setq erg (funcall (first (fz-other b))))
  501.       (setf (fz-other b) nil)
  502.       (setf (fz-cl b) 'factored)
  503.       (setf (fz-other b) (build-fzl (list erg (exquo (fz-num b) erg))))
  504.       (return it))
  505.     (pop (fz-other b)) ; Funktionen-Schlange rotieren
  506. ) )
  507.  
  508.  
  509. ; Alle Faktorisierungsalgorithmen bekommen eine ungerade Zahl n>1, die
  510. ; garantiert zusammengesetzt ist, und liefern eine parameterlose Funktion ab,
  511. ; die nach wiederholtem Aufruf einen Faktor von n liefert (und NIL, solange
  512. ; sie noch keinen Faktor gefunden hat).
  513.  
  514.  
  515. ; erster Faktorisierungsalgorithmus: Teilersuche bis zur Wurzel
  516. ; [Allgemein bekannt.]
  517. (defun tryfactor1 (n)
  518.   (let ((p 1)
  519.         (s (isqrt n)))
  520.     #'(lambda ()
  521.         (block tryfactoring1
  522.           (dotimes (i (min 100 (ceiling (- s p) 2)))
  523.             (incf p 2)
  524.             (if (zerop (mod n p)) (return-from tryfactoring1 p))
  525.           )
  526.           (if (>= p s) (error "Nicht-Primzahl ~D ohne Teiler." n))
  527.           nil
  528. ) )   ) )
  529.  
  530.  
  531. ; zweiter Faktorisierungsalgorithmus: Pollardsche Iteration, in der Version
  532. ; von R.P. Brent.
  533. ; [Richard P. Brent, An improved Monte Carlo factorization algorithm,
  534. ;  BIT 20 (1980), 176-184, Seite 182].
  535. (defun tryfactor2 (n &aux (m 100))
  536.   (let ((state 1) ; Zustand: 1 bei der Initialisierung,
  537.                   ;          2 bei den Iterationen ins Leere,
  538.                   ;          3 bei den Iterationen in Blöcken.
  539.         (s 1) ; Iterationsparameter
  540.         x     ; relative Anfangszahl
  541.         y     ; Iterationswert
  542.         r     ; wird immer wieder verdoppelt
  543.         k     ; Anzahl der bisher ausgeführten Iterationen (0<=k<=r)
  544.         q     ; bisher aufgelaufenes Produkt
  545.        ) ; m = Blocklänge
  546.     (flet ((brentiter (x)
  547.              (mod (+ (* x x) s) n)
  548.           ))
  549.       #'(lambda ()
  550.           (let ((it 100)) ; Noch verbleibende Iterationen
  551.             (loop
  552.               (if (<= it 0) (return nil))
  553.               (case state
  554.                 ((1) (setq y 0)
  555.                      (setq x y)
  556.                      (setq r 1)
  557.                      (setq k 0)
  558.                      (setq q 1)
  559.                      (setq state 2)
  560.                 )
  561.                 ((2) (let ((i (min it (- r k))))
  562.                        (dotimes (i1 i) (setq y (brentiter y)))
  563.                        (incf k i)
  564.                        (decf it i))
  565.                      (when (= k r)
  566.                        (setq k 0)
  567.                        (setq state 3)
  568.                 )    )
  569.                 ((3) (let ((ys y)
  570.                            (i (min it (- r k) m))
  571.                            g)
  572.                        (dotimes (i1 i)
  573.                          (setq y (brentiter y))
  574.                          (setq q (mod (* q (abs (- x y))) n))
  575.                        )
  576.                        (incf k i)
  577.                        (decf it i)
  578.                        (setq g (gcd q n))
  579.                        (cond ((> g 1) ; Habe wohl etwas gefunden
  580.                          (if (= g n) ; wirklich?
  581.                            (loop ; nachiterieren
  582.                              (setq ys (brentiter ys))
  583.                              (decf it)
  584.                              (setq g (gcd (- x ys) n))
  585.                              (if (> g 1) (return))
  586.                          ) )
  587.                          (if (< g n) ; ja!
  588.                            (return g))
  589.                          ; nein ...
  590.                          (incf s) ; ganz neuer Versuch
  591.                          (setq state 1)
  592.                         )
  593.                         ((= k r) ; state 3 zu Ende
  594.                          (setq x y)
  595.                          (setq r (* 2 r))
  596.                          (setq k 0)
  597.                          (setq state 2)
  598.                 )    ) ))
  599. ) ) )   ) ) ) )
  600.  
  601.  
  602. ; dritter Faktorisierungsalgorithmus: Fermats Methode.
  603. ; [Donald E. Knuth, The art of computer programming, Band 2, Seminumerical
  604. ;  algorithms, Abschnitt 4.5.4, Algorithmus C, Seite 371.]
  605. (defun tryfactor3 (n)
  606.   (let* ((s (isqrt n)) 
  607.          (x (1+ (* 2 s)))
  608.          (y 1)
  609.          (r (- (* s s) n)) )
  610.     #'(lambda ()
  611.         (dotimes (it 500)
  612.           (cond ((> r 0) (decf r y) (incf y 2))
  613.                 ((< r 0) (incf r x) (incf x 2))
  614.                 (t (assert (zerop r)) (return (exquo (- x y) 2)))
  615. ) )   ) ) )
  616.  
  617. #|
  618. ; vierter Faktorisierungsalgorithmus: Siebmethode.
  619. ; [Donald E. Knuth, The art of computer programming, Band 2, Seminumerical
  620. ;  algorithms, Abschnitt 4.5.4, Algorithmus D, Seite 373.]
  621. (defun tryfactor4 (n)
  622.   (let* ((x (isqrt n))
  623.          (mlist '#(3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61))
  624.          (r (length mlist))
  625. |#
  626.  
  627. ; fünfter Faktorisierungsalgorithmus: Kettenbruch.
  628. ; Auch Algorithmus von Kraitchik - Morrison - Brillhart genannt.
  629. ; [Donald E. Knuth, The art of computer programming, Band 2, Seminumerical
  630. ;  algorithms, Abschnitt 4.5.4, Algorithmus E, Seiten 380-384.]
  631. ; [Carl Pomerance: The quadratic sieve factoring algorithm, in:
  632. ;  Lecture Notes in Computer Science, Band 209, Advances in Cryptology,
  633. ;  Seiten 169-173.]
  634. ; [Marvin C. Wunderlich: A running time analysis of Brillhart's continued
  635. ;  fraction factoring method, in:
  636. ;  Lecture Notes in Mathematics, Band 751, Seiten 328-342.]
  637. (defun tryfactor5 (n)
  638.   (let* ((m (floor (exp (* 0.25 (sqrt (* (log n) (log (log n))))))))
  639.          ; man kann den Faktor 0.25 auch vergrößern, wenn der Platz reicht.
  640.          (factorbase (make-array (list m) :element-type 'integer))
  641.                        ; wird ein Array von m Primzahlen sein
  642.          (eqnlist nil) ; wird eine Liste von Gleichungen sein, wobei die Array-
  643.                        ; Komponenten in Gaußscher Treppennormalform über
  644.                        ; GF(2) sind.
  645.         )
  646.     ; (assert (>= m 1))
  647.     (flet ((base-factor (v) ; stellt fest, ob abs(v) ganz über der Faktorbasis
  648.              ; zerfällt. Wenn ja, also v = (-1)^e0 * ... * pm^em * v1^2, kommt
  649.              ; ( #(e0 ... em) . v1) als Ergebnis.
  650.              (block base-factoring
  651.                (let ((a (make-array (list (1+ m)) :initial-element 0)))
  652.                  (assert (not (zerop v)))
  653.                  (when (minusp v) ; Ziehe Faktor (-1) heraus
  654.                    (setq v (- v))
  655.                    (setf (aref a 0) 1)
  656.                  )
  657.                  (dotimes (i m) ; Ziehe kleine Faktoren heraus
  658.                    (let ((p (aref factorbase i)))
  659.                      (do ()
  660.                        ((cond ((zerop (mod v p))
  661.                                (setq v (exquo v p))
  662.                                (incf (aref a (1+ i)))
  663.                                nil)
  664.                               (t)
  665.                  ) ) ) ))
  666.                  (unless (= v 1) ; Ist der Rest ein Quadrat?
  667.                    (let ((v2 (isqrt v)))
  668.                      (if (= (* v2 v2) v)
  669.                          (setq v v2)
  670.                          (return-from base-factoring nil) ; nein -> wertlos
  671.                  ) ) )
  672.                  (cons a v)
  673.            ) ) )
  674.            (add-eqn (eqn) ; fügt eine neue Gleichung eqn = (#(e0..em) x . v)
  675.                           ; die x^2 = (-1)^e0 ... pm^em v^2 mod n bedeutet,
  676.                           ; zum Wissen in eqnlist hinzu.
  677.              ; Wird ein nichttrivialer Faktor von n gefunden, erfolgt ein THROW.
  678.              (block add-eqn
  679.                (do ((l1 nil)
  680.                     (l2 eqnlist)
  681.                     h ha a)
  682.                  ; stets (append (reverse l1) l2) = eqnlist
  683.                  ((endp l2) ; Liste eqnlist zu Ende untersucht.
  684.                   (if ; Sind in eqn alle Exponenten gerade?
  685.                       (some #'oddp (car eqn)) ; Nein: muß eqn behalten
  686.                       (setq eqnlist (nreconc l1 (list eqn)))
  687.                       ; Ja! Ist die Gleichung x^2 = y^2 mod n trivial ?
  688.                       (let ((x (second eqn))
  689.                             (y (let ((v (cddr eqn)))
  690.                                  (dotimes (i m)
  691.                                    (setq v (* v
  692.                                      (expt (aref factorbase i)
  693.                                            (exquo (aref (car eqn) (1+ i)) 2)
  694.                                  ) )))
  695.                                  v
  696.                            ))  )
  697.                         (assert (zerop (mod (- (* x x) (* y y)) n)))
  698.                         (unless (or (zerop (mod (- x y) n))
  699.                                     (zerop (mod (+ x y) n)))
  700.                           ; Nein! Fertig!
  701.                           (throw 'tryfactoring5 (gcd (- x y) n))
  702.                  ))   ) )
  703.                  (setq h (car l2)) ; Hilfsgleichung aus der Liste
  704.                  ; Falls h und eqn an derselben Stelle ihre erste ungerade Zahl
  705.                  ; haben, wird h zu eqn multipliziert:
  706.                  (setq ha (car h)) (setq a (car eqn))
  707.                  (if (do ((i 0 (1+ i)))
  708.                          ((> i m) nil)
  709.                        (cond ((and (oddp (aref ha i)) (evenp (aref a i)))
  710.                               (return nil))
  711.                              ((and (evenp (aref ha i)) (oddp (aref a i)))
  712.                               ; muß eqn vor h einfügen
  713.                               (setq eqnlist (nreconc l1 (cons eqn l2)))
  714.                               (return-from add-eqn nil))
  715.                              ((and (oddp (aref ha i)) (oddp (aref a i)))
  716.                               (return t))
  717.                      ) )
  718.                      (progn
  719.                        (dotimes (i (1+ m))
  720.                          (incf (aref a i) (aref ha i))
  721.                        )
  722.                        (setq eqn (cons a
  723.                                  (cons (mod (* (second eqn) (second h)) n)
  724.                                        (mod (* (cddr eqn) (cddr h)) n)
  725.                      ) )         ))
  726.                  )
  727.                  (pop l2)
  728.                  (push h l1)
  729.            ) ) )
  730.            (calc-factorbase (D) ; berechnet zu D=k*n eine passende Faktorbasis
  731.              (do ((i 0)
  732.                   (p 2 (1+ p)))
  733.                  ((>= i m))
  734.                (when (or (= p 2) (and (klprimtest p) (= (jacobi D p) 1)))
  735.                      (setf (aref factorbase i) p)
  736.                      (incf i)
  737.            ) ) )
  738.           )
  739.       (let ((D N) ; D=k*n
  740.             (Anfangszustand t)
  741.             R R1 U U1 V V1 P P1 A S   ; Variablenbenennung wie in [Knuth]
  742.             ; Mit den Bezeichnungen wie in [Wunderlich]:
  743.             ; Es gilt für alle i=0,1,...:
  744.             ; sqrt(D) = [q0,...,qi-1,(sqrt(D)+Pi)/Qi]
  745.             ; qi = floor((sqrt(D)+Pi)/Qi) = floor((floor(sqrt(D))+Pi)/Qi),
  746.             ; Q-1 = D - q0^2, P0=0, Q0=1, Pi+1 = qi * Qi - Pi,
  747.             ; Qi+1 = (D - Pi+1)^2 / Qi = Qi-1 - qi * ( qi * Qi - 2 * Pi ),
  748.             ; alle abs(Pi)<sqrt(D), alle Qi>0, alle qi>0, alle Variablen ganz.
  749.             ; A-1=1, B-1=0, A0=q0, B0=1,
  750.             ; Ai = qi * Ai-1 + Ai-2, Bi = qi * Bi-1 + Bi-2 für i>0,
  751.             ; in Q(X) gilt (Ai+x*Ai-1)/(Bi+x*Bi-1) = [q0,...,qi-1,qi + x].
  752.             ; Für i>=0 ist Ai-1^2 - D * Bi-1^2 = (-1)^i * Qi
  753.             ; Für i>=0 ist Ai-1 * Ai - D * Bi-1 * Bi = (-1)^i * Pi+1
  754.             ; Nach i Schleifendurchläufen ist
  755.             ; R = floor(sqrt(D))
  756.             ; R1 = 2*floor(sqrt(D))
  757.             ; U = floor(sqrt(D)) + Pi+1
  758.             ; U1 = floor(sqrt(D)) + Pi
  759.             ; V = Qi+1
  760.             ; V1 = Qi
  761.             ; P = Ai mod N
  762.             ; P1 = Ai-1 mod N
  763.             ; A = qi
  764.             ; S = (-1)^(i+1)
  765.            )
  766.         #'(lambda ()
  767.             (catch 'tryfactoring5
  768.               (do ((it 200) h1) ; Iterationenzahl und Hilfsvariable
  769.                   ((<= it 0))
  770.                 (when Anfangszustand ; muß initialisieren
  771.                   (calc-factorbase D) ; Faktorbasis berechnen.
  772.                   (dotimes (i m)      ; Test auf Teilbarkeit in Faktorbasis.
  773.                     (if (zerop (mod n (aref factorbase i)))
  774.                         ; muß ein echter Teiler sein, weil n nicht prim
  775.                         (throw 'tryfactoring5 (aref factorbase i))
  776.                   ) )
  777.                   (setq R (isqrt D))
  778.                   (setq R1 (* 2 R))
  779.                   (setq U R1)
  780.                   (setq U1 R)
  781.                   (setq V (- D (* R R)))
  782.                   (if (zerop V)       ; k*N=D ist Quadratzahl = R^2
  783.                       ; -> g:=ggT(R,N) teilt N.
  784.                       ; Bei g=1 wäre ggT(R^2,N)=1 => R^2|k => N=1, Widerspruch.
  785.                       ; Bei g=N wäre N|R => N^2|R^2=k*N => N|k, das passiert
  786.                       ; zu unseren Lebzeiten nicht mehr.
  787.                       (throw 'tryfactoring5 (gcd R n)) )
  788.                   (setq V1 1)
  789.                   (setq P R)
  790.                   (setq P1 1)
  791.                   (setq A R)
  792.                   (setq S -1)
  793.                   (setq Anfangszustand nil)
  794.                   (decf it 2)
  795.                 )
  796.                 ; Hier gelten die ganzen Invarianten.
  797.                 ; (assert (zerop (mod (+ (* P1 P1) (* S V1)) n)))
  798.                 ; (assert (zerop (mod (- (* P P) (* S V)) n)))
  799.                 (if (setq h1 (base-factor (* S V)))
  800.                     (add-eqn (cons (car h1) (cons P (cdr h1)))) )
  801.                 (when (= V 1) ; Periode des Kettenbruches erreicht ->
  802.                   (incf D n)
  803.                   (setq Anfangszustand t) ; neuer Anlauf mit größerem k.
  804.                 )
  805.                 ; nächste Kettenbruchiteration:
  806.                 (multiple-value-setq (A h1) (floor U V)) ; A:=qi+1
  807.                 (psetq  U1 U  U (- R1 h1)) ; U1:=R+Pi+1, U:=R+Pi+2
  808.                 (psetq  V1 V  V (- V1 (* A (- U U1)))) ; V1:=Qi+1, V:=Qi+2
  809.                 (setq S (- S)) ; S:=(-1)^(i+2)
  810.                 (psetq  P1 P  P (mod (+ (* A P) P1) n))
  811.                   ; P1:=Ai mod N, P:=Ai+1 mod N
  812.                 ; jetzt ist ein Durchlauf fertig, i:=i+1
  813.                 (decf it)
  814. ) ) ) )   ) ) )
  815.  
  816. #|
  817. ; sechster Faktorisierungsalgorithmus: Quadratisches Sieben.
  818. |#
  819.  
  820.  
  821. ; (cyclic-list x1 ... xn) mit n>0
  822. ; liefert die zyklische Liste #1=(x1 ... xn . #1#)
  823. (defun cyclic-list (&rest args)
  824.   (setf (cdr (last args)) args)
  825.   args
  826. )
  827.  
  828.  
  829. ; (verffz b) verfeinert die in b steckenden Faktorzerlegung von n = (fz-num b)
  830. ; mit bis zu k Iterationsschritten (k fehlt -> beliebig viele Schritte), solange
  831. ; bis eine bestimmte Bedingung (return-test b) zutrifft ( /= NIL liefert).
  832. ; Ergebnis ist die Anzahl der verbrauchten Iterationen.
  833. (defun verffz (b &optional (k nearly-infinite)
  834.                  &key (return-test #'(lambda (x)
  835.                                        (declare (ignore x))
  836.                                        nil) )
  837.                  &aux (it 0))
  838.    (if (eq (fz-cl b) 'prime) (return-from verffz it)) ; nichts zu tun.
  839.    (if (eq (fz-cl b) 'probprime)
  840.      (let ((n (fz-num b)))
  841.         ; Jetzt muß der Primzahlverdacht erhärtet werden.
  842.         ; Daß n>2 ungerade ist, ist schon bekannt.
  843.         (if (null (fz-other b)) ; noch keine wartenden Prozeduren
  844.             (setf (fz-other b)
  845.                   (cyclic-list
  846.                     (tryprimeprover1 n)
  847.                     (tryprimeprover2 n)
  848.                     ; usw. (weitere Beweisversucher)
  849.         )   )     )
  850.         ; Jetzt müssen die wartenden Funktionen der Reihe nach aktiviert
  851.         ; werden, bis sie ein Ergebnis liefern.
  852.         (incf it (tryprimeprove b k))
  853.         (return-from verffz it)
  854.    ) )
  855.    (if (eq (fz-cl b) 'notprime)
  856.        ; erst einmal muß ein Faktor von n gefunden werden:
  857.        (let ((n (fz-num b)))
  858.          (if (null (fz-other b))
  859.              (setf (fz-other b)
  860.                    (cyclic-list
  861.                      (tryfactor1 n)
  862.                      (tryfactor2 n)
  863.                      (tryfactor3 n)
  864.                      ; usw.
  865.          )   )     )
  866.          (incf it (tryfactor b k))
  867.    )   )
  868.    (if (eq (fz-cl b) 'notprime)
  869.        ; konnte noch nicht faktorisieren
  870.        (return-from verffz it))
  871.    (assert (eq (fz-cl b) 'factored))
  872.    (if (fzl-recalc (fz-other b))
  873.        (return-from verffz it)) ; Verfeinern unmöglich
  874.    ; Jetzt wird die Anstrengung auf die einzelnen Faktoren verteilt.
  875.    (let ((l nil) f) ; l = Liste noch zu verfeinernder Faktoren
  876.      (loop
  877.         (when (or (>= it k) (funcall return-test b)) (return))
  878.         (if (null l) (setq l (fzl-factorlist (fz-other b)))) ; von neuem
  879.         (setq f (first l))
  880.         (if (not (integerp f))
  881.           (progn
  882.             (incf it (verffz f                 ; rekursiv verfeinern
  883.                              (ceiling (max 0 (- k it)) (length l))))
  884.             (when (eq (fz-cl f) 'factored)
  885.                   ; in b wird f durch seine Faktoren ersetzt
  886.                   (setf (fzl-factorlist (fz-other b))
  887.                     (mergefz (remove f (fzl-factorlist (fz-other b))
  888.                                      :count 1 :test #'eq)
  889.                              (copy-list (fzl-factorlist (fz-other f)))
  890.             )     ) )
  891.                   ; (beachte: l bleibt hierbei die Liste der noch nicht
  892.                   ; abgearbeiteten Faktoren, der Schwanz von
  893.                   ; (fzl-factorlist (fz-other b)). )
  894.             (when (member (fz-cl f) '(prime factored))
  895.                   (if (fzl-recalc (fz-other b)) ; noch was zu tun?
  896.                       (return-from verffz it))
  897.             )
  898.         ) )
  899.         (pop l)
  900.      )
  901.      ; (fzl-factorlist (fz-other b)) ist jetzt zur Genüge verfeinert.
  902.      it
  903. ) )
  904.  
  905.  
  906. ; (pfz n) liefert die Primfaktorzerlegung von n>0:
  907. ; alle Primfaktoren, in aufsteigender Reihenfolge, evtl. mehrfach vertreten.
  908. (defun pfz (n)
  909.   (if (< n *Kleinheit*)
  910.       (klpfz n)
  911.       (let ((b (pfz-table-lookup n)))
  912.         (when (null b) ; n noch nicht in der Tabelle?
  913.               (mr n) ; n in die Tabelle eintragen
  914.               (setq b (pfz-table-lookup n))
  915.         )
  916.         (do ((i 5 (* 2 i))) ; i wächst stark
  917.             ((or (eq (fz-cl b) 'prime)
  918.                  (and (eq (fz-cl b) 'factored) (fzl-allprimes (fz-other b)))
  919.              ))
  920.           (verffz b i) ; ewig verfeinern
  921.         )
  922.         (if (eq (fz-cl b) 'prime)
  923.             (list n)
  924.             (mapcar #'fz-key (fzl-factorlist (fz-other b)))
  925. ) )   ) )
  926.  
  927.  
  928. ; (setpfz n l) speichert in der Tabelle ein (von irgendwoher bekanntes)
  929. ; Primfaktorzerlegungsergebnis ab. Es sollte irgendwann l als Ergebnis von
  930. ; (pfz n) gekommen sein.
  931. ; Nach (setpfz n l) ist sichergestellt, daß in Zukunft bei (pfz n) ohne
  932. ; Rechnung sofort das Ergebnis l kommt.
  933. (defun setpfz (n l)
  934.   (if ; Plausibilitätsprüfung, ob (pfz n) = l sein kann.
  935.       (and (apply #'<= 2 l)
  936.            (every #'(lambda (f) (or (>= f *Kleinheit*) (klprimtest f))) l)
  937.            (= n (apply #'* l))
  938.       )
  939.       (progn
  940.         ; erst die Primfaktoren als Primzahlen in die Tabelle abspeichern
  941.         (setq l
  942.           (mapcar
  943.             #'(lambda (f)
  944.                 (if (< f *Kleinheit*)
  945.                     f
  946.                     (let ((b (pfz-table-lookup-insert f)))
  947.                       (unless (and (eq (fz-cl b) 'prime)
  948.                                    (integerp (fz-other b)))
  949.                         (setf (fz-cl b) 'prime)
  950.                         (setf (fz-other b) nil)
  951.                       )
  952.                       b
  953.               ) )   )
  954.             l
  955.         ) )
  956.         ; l ist die neue Faktorzerlegungsliste
  957.         (when (and (>= n *Kleinheit*) (> (length l) 1))
  958.           (let ((b (pfz-table-lookup-insert n)))
  959.             (setf (fz-cl b) 'factored)
  960.             (setf (fz-other b) (make-fzl :allprimes t :factorlist l))
  961.         ) )
  962.         n ; als Ergebnis
  963.       )
  964.       (error "Fehler: ~D hat nicht die Primfaktorzerlegung ~S.~%" n l)
  965. ) )
  966.  
  967.  
  968. ;-------------------------------------------------------------------------------
  969. ; Anwendungen (benutzen nur pfz):
  970.  
  971.  
  972. ; (pfzv n) ergibt die Primfaktorzerlegung von n>0 in der Gestalt
  973. ; ((p1 . e1) (p2 . e2) ... (pk . ek)) mit Primzahlen p1<...<pk und
  974. ; Exponenten e1,...,ek>0 mit n=p1^e1 ... pk^ek.
  975. (defun pfzv (n)
  976.   (do ((l1 (pfz n))
  977.        (l2 nil))
  978.       ((endp l1) (reverse l2))
  979.     (let* ((p (first l1))
  980.            (e (do ((i 0 (1+ i)))
  981.                   ((or (null l1) (not (= p (first l1)))) i)
  982.                 (pop l1)
  983.           ))  )
  984.        (push (cons p e) l2)
  985. ) ) )
  986.  
  987.  
  988. ; (primteiler n) liefert die Primteiler von n>0, jeden nur einmal.
  989. (defun primteiler (n) (mapcar #'car (pfzv n)))
  990.  
  991.  
  992. ; (pdivisors n) ergibt die positiven Teiler einer Zahl n>0, in
  993. ; irgendeiner Reihenfolge
  994. (defun pdivisors (n)
  995.   (reduce
  996.     #'(lambda (l pe)          ; ergibt l, elementweise mit 1,p,...,p^e
  997.         (let ((p (first pe))  ; multipliziert.
  998.               (e (rest pe)))
  999.           (mapcan
  1000.             #'(lambda (m)     ; ergibt die mit m multiplizierte Liste l
  1001.                 (mapcar #'(lambda (x) (* m x)) l)
  1002.               )
  1003.             (mapcar #'(lambda (i) (expt p i)) (intlist 0 e))
  1004.       ) ) )
  1005.     (pfzv n)
  1006.     :initial-value '(1)
  1007.   )
  1008. )
  1009.  
  1010.  
  1011. ; (divisors n) ergibt die ganzzahligen Teiler einer Zahl n /= 0
  1012. (defun divisors (n)
  1013.   (let ((pdiv (pdivisors (abs n))))
  1014.     (append pdiv (mapcar #'- pdiv))
  1015. ) )
  1016.  
  1017.  
  1018. ; (euler-phi n) ergibt für n>0 die Eulersche phi-Funktion von n,
  1019. ; d.h. die Anzahl der zu n teilerfremden natürlichen Zahlen >0, <=n.
  1020. (defun euler-phi (n)
  1021.   (reduce #'*
  1022.           (mapcar #'(lambda (pe)
  1023.                       (let ((p (car pe)) ; Primzahl
  1024.                             (e (cdr pe)) ; Exponent, >0
  1025.                            )
  1026.                         (* (expt p (1- e)) (1- p))
  1027.                     ) )
  1028.                   (pfzv n)
  1029. ) )       )
  1030.  
  1031.  
  1032. ; (prim-von-bis a b) ergibt für eine Liste aller Primzahlen von a bis b
  1033. ; (inklusive).
  1034. (defun prim-von-bis (a b)
  1035.   (remove-if-not #'isprime (intlist a b))
  1036. )
  1037.  
  1038.