home *** CD-ROM | disk | FTP | other *** search
/ GEMini Atari / GEMini_Atari_CD-ROM_Walnut_Creek_December_1993.iso / files / language / examples / potr.lsp < prev    next >
Encoding:
Text File  |  1993-10-23  |  38.1 KB  |  1,019 lines
  1. ;;;; Funktionen zum Arbeiten mit Potenzreihen, 11. 10. 1988
  2. ;;;; erweitert am 24.8.1989, 30.11.1989, 28.1.1990 von Bruno Haible
  3. ;;;; erweitert am 2.8.1990 von Michael Stoll und Bruno Haible
  4. ;;;; erweitert am 10.8.1990 von Bruno Haible
  5.  
  6. (provide 'potr)
  7.  
  8. ;;; Potenzreihen sind funktionale Objekte, die mit einer nichtnegativen
  9. ;;; ganzen Zahl als Argument aufgerufen werden können und dann als
  10. ;;; Wert einen Vektor (adjustable, mit Fillpointer) zurückgeben,
  11. ;;; der die Koeffizienten der Potenzreihe enthält und dessen Länge
  12. ;;; mindestens so groß ist wie das angegebene Argument. Eine Potenzreihe
  13. ;;; benutzt immer denselben Vektor, so daß zurückgegebene Werte bei
  14. ;;; weiterer Berechnung mitgeführt werden.
  15.  
  16. ;;; Variablenbezeichnungen:
  17. ;;; pvec ist stets der Vektor für die Koeffizienten,
  18. ;;; pcount ist die aktuelle Länge des Vektors (d.h. die Koeffizienten
  19. ;;;        bis pcount-1 sind bekannt, der nächste hat Index pcount)
  20.  
  21.  
  22. ;; |0| ist die Nullreihe
  23. (defvar \0
  24.   (let* ((pcount 0)
  25.          (pvec (make-array 10 :fill-pointer pcount :adjustable t))
  26.         )
  27.     #'(lambda (n)
  28.         (if (<= n pcount) pvec
  29.           (dotimes (i (- n pcount) pvec)
  30.             (vector-push-extend 0 pvec)
  31.             (setq pcount (1+ pcount))
  32. ) )   ) ) )
  33.  
  34. ;; |1| ist die Einsreihe
  35. (defvar \1
  36.   (let* ((pcount 1)
  37.          (pvec (make-array 10 :fill-pointer pcount :adjustable t))
  38.         )
  39.     (setf (aref pvec 0) 1)
  40.     #'(lambda (n)
  41.         (if (<= n pcount) pvec
  42.           (dotimes (i (- n pcount) pvec)
  43.             (vector-push-extend 0 pvec)
  44.             (setq pcount (1+ pcount))
  45. ) )   ) ) )
  46.  
  47. ;; Mit defpot kann man Reihen wie exp definieren
  48. ;; kof-form ist eine Form für den einzusetzenden Koeffizienten,
  49. ;;  dabei ist der Index als pcount erreichbar.
  50. ;; Ist ein kof-init angegeben, dann wird eine Variable kof definiert
  51. ;;  und auf den angegebenen Wert gesetzt (sie kann dann in kof-form
  52. ;;  verwendet werden).
  53. ;; Ist kof-step (eine Form) angegeben, dann wird nach der Berechnung
  54. ;;  und Abspeicherung des Koeffizienten kof auf den Wert der Form
  55. ;;  kof-step gesetzt. pcount hat dabei den nächsten Index als Wert.
  56. (defmacro defpot (kof-form &optional kof-init kof-step)
  57.   `(let* ((pcount 0)
  58.           (pvec (make-array 10 :fill-pointer pcount :adjustable t))
  59.           ,@(if kof-init `((kof ,kof-init)) nil)
  60.          )
  61.      #'(lambda (n)
  62.          (if (<= n pcount) pvec
  63.            (dotimes (i (- n pcount) pvec)
  64.              (vector-push-extend ,kof-form pvec)
  65.              (setq pcount (1+ pcount)
  66.                    ,@(if kof-step `(kof ,kof-step) nil)
  67. )  )   ) ) ) )
  68.  
  69. (defvar exp
  70.   (defpot kof 1 (/ kof pcount))
  71. )
  72.  
  73. (defvar sin
  74.   (defpot (if (oddp pcount)
  75.             (if (evenp (floor pcount 2)) kof (- kof))
  76.             0
  77.           )
  78.           1 (/ kof pcount)
  79. ) )
  80.  
  81. (defvar cos
  82.   (defpot (if (evenp pcount)
  83.             (if (evenp (floor pcount 2)) kof (- kof))
  84.             0
  85.           )
  86.           1 (/ kof pcount)
  87. ) )
  88.  
  89. (defvar log
  90.   (defpot (if (zerop pcount) 0
  91.             (if (oddp pcount) (/ pcount) (- (/ pcount)))
  92. ) )       )
  93.  
  94. ;; pot+ addiert Potenzreihen
  95. (defun pot+ (&rest pots)
  96.   (if (null pots) |0|
  97.     (if (null (rest pots)) (first pots)
  98.       (let* ((pcount 0)
  99.              (pvec (make-array 10 :fill-pointer pcount :adjustable t))
  100.              (summanden (mapcar #'(lambda (pot) (funcall pot 0)) pots))
  101.             )
  102.         #'(lambda (n)
  103.             (if (<= n pcount) pvec
  104.               (progn
  105.                 (dolist (pot pots) (funcall pot n))
  106.                 (dotimes (i (- n pcount) pvec)
  107.                   (vector-push-extend
  108.                     (do ((s summanden (rest s))
  109.                          (sum 0 (+ sum (aref (first s) pcount)))
  110.                         )
  111.                         ((null s) sum)
  112.                     )
  113.                     pvec
  114.                   )
  115.                   (setq pcount (1+ pcount))
  116. ) ) ) )   ) ) ) )
  117.  
  118. ;; pot- subtrahiert Potenzreihen
  119. (defun pot- (pot1 &rest pots)
  120.   (let* ((pcount 0)
  121.          (pvec (make-array 10 :fill-pointer pcount :adjustable t))
  122.          (erste (funcall pot1 0))
  123.          (subtrahenden (mapcar #'(lambda (pot) (funcall pot 0)) pots))
  124.         )
  125.     #'(lambda (n)
  126.         (if (<= n pcount) pvec
  127.           (progn
  128.             (funcall pot1 n)
  129.             (dolist (pot pots) (funcall pot n))
  130.             (dotimes (i (- n pcount) pvec)
  131.               (vector-push-extend
  132.                 (do ((s subtrahenden (rest s))
  133.                      (sum (aref erste pcount)
  134.                           (- sum (aref (first s) pcount))
  135.                     ))
  136.                     ((null s) sum)
  137.                 )
  138.                 pvec
  139.               )
  140.               (setq pcount (1+ pcount))
  141. ) )   ) ) ) )
  142.  
  143. ;; pot*skal multipliziert eine Potenzreihe mit einem skalaren Faktor
  144. (defun pot*skal (pot skal)
  145.   (let* ((pcount 0)
  146.          (pvec (make-array 10 :fill-pointer pcount :adjustable t))
  147.          (pvec1 (funcall pot 0))
  148.         )
  149.     #'(lambda (n)
  150.         (if (<= n pcount) pvec
  151.           (progn
  152.             (funcall pot n)
  153.             (dotimes (i (- n pcount) pvec)
  154.               (vector-push-extend (* skal (aref pvec1 pcount)) pvec)
  155.               (setq pcount (1+ pcount))
  156. ) )   ) ) ) )
  157.  
  158. ;; potshift multipliziert eine Potenzreihe mit einer Potenz der Unbestimmten
  159. (defun potshift (pot shift)
  160.   (unless (integerp shift)
  161.           (error "~: Der Shift muß eine nichtnegative ganze Zahl sein, nicht ~"
  162.                  'potshift shift
  163.   )       )
  164.   (if (zerop shift)
  165.     pot
  166.     (let ((pvec1 (funcall pot (max 0 (- shift)))))
  167.       (when (minusp shift)
  168.         (dotimes (i (- shift))
  169.           (unless (zerop (aref pvec1 i))
  170.             (error "~: Potenzreihe ~ nicht durch x^~ teilbar"
  171.                    'potshift pvec1 (- shift)
  172.       ) ) ) )
  173.       (let* ((pcount 0)
  174.              (pvec (make-array 10 :fill-pointer pcount :adjustable t)))
  175.         (when (plusp shift)
  176.           (dotimes (i shift)
  177.             (vector-push-extend 0 pvec) (setq pcount (1+ pcount))
  178.         ) )
  179.         #'(lambda (n)
  180.             (if (<= n pcount) pvec
  181.               (progn
  182.                 (funcall pot (- n shift))
  183.                 (dotimes (i (- n pcount) pvec)
  184.                   (vector-push-extend (aref pvec1 (- pcount shift)) pvec)
  185.                   (setq pcount (1+ pcount))
  186. ) ) ) )   ) ) ) )
  187.  
  188. ;; pot*2 multipliziert zwei Potenzreihen
  189. (defun pot*2 (pot1 pot2)
  190.   (let* ((pcount 0)
  191.          (pvec (make-array 10 :fill-pointer pcount :adjustable t))
  192.          (factor1 (funcall pot1 0))
  193.          (factor2 (funcall pot2 0))
  194.         )
  195.     #'(lambda (n)
  196.         (if (<= n pcount) pvec
  197.           (progn
  198.             (funcall pot1 n) (funcall pot2 n)
  199.             (dotimes (i (- n pcount) pvec)
  200.               (vector-push-extend
  201.                 (do ((j 0 (1+ j))
  202.                      (k pcount (1- k))
  203.                      (sum 0 (+ sum (* (aref factor1 j) (aref factor2 k))))
  204.                     )
  205.                     ((minusp k) sum)
  206.                 )
  207.                 pvec
  208.               )
  209.               (setq pcount (1+ pcount))
  210. ) )   ) ) ) )
  211.  
  212. ;; pot* multipliziert (beliebig viele) Potenzreihen
  213. (defun pot* (&rest pots)
  214.   (if (null pots) |1|
  215.     (if (null (rest pots)) (first pots)
  216.       (pot*2 (first pots) (apply #'pot* (rest pots)))
  217. ) ) )
  218.  
  219. ;; potkehr bildet den Kehrwert einer invertierbaren Potenzreihe
  220. (defun potkehr (pot)
  221.   (let* ((pvec1 (funcall pot 1))
  222.          (const (aref pvec1 0))
  223.         )
  224.     (when (zerop const)
  225.           (error "~: Die Potenzreihe ~ ist nicht invertierbar"
  226.                  'potkehr pvec1
  227.     )     )
  228.     (let* ((pcount 1)
  229.            (pvec (make-array 10 :fill-pointer pcount :adjustable t))
  230.            (const1 (- (/ const)))
  231.           )
  232.       (setf (aref pvec 0) (/ const))
  233.       #'(lambda (n)
  234.           (if (<= n pcount) pvec
  235.             (progn
  236.               (funcall pot n)
  237.               (dotimes (i (- n pcount) pvec)
  238.                 (vector-push-extend
  239.                   (do ((j 1 (1+ j))
  240.                        (k (1- pcount) (1- k))
  241.                        (sum 0 (+ sum (* (aref pvec1 j) (aref pvec k))))
  242.                       )
  243.                       ((minusp k) (* const1 sum))
  244.                   )
  245.                   pvec
  246.                 )
  247.                 (setq pcount (1+ pcount))
  248. ) ) )   ) ) ) )
  249.  
  250. ;; pot/ dividiert zwei Potenzreihen, falls möglich
  251. ;; (Bem.: 0/0 führt auf Endlosschleife, das läßt sich nicht vermeiden)
  252. ; Methode:
  253. ; pot1 = x^offset*(a0+a1*x+...),
  254. ; pot2 = x^offset*(b0+b1*x+...) mit b0 /= 0,
  255. ; -> pot1/pot2 = c0+c1*x+... wobei
  256. ;    c[n] = (a[n]-(b[n]c[0]+...+b[1]c[n-1]))/b[0] für n=0,1,...
  257. (defun pot/ (pot1 pot2)
  258.   (let* ((pvec1 (funcall pot1 1))
  259.          (pvec2 (funcall pot2 1))
  260.          (offset 0)
  261.         )
  262.     (loop
  263.       (unless (zerop (aref pvec2 offset)) (return))
  264.       (unless (zerop (aref pvec1 offset))
  265.               (error "~: Potenzreihe ~ läßt sich nicht durch ~ dividieren"
  266.                      'pot/ pvec1 pvec2
  267.       )       )
  268.       (setq offset (1+ offset))
  269.       (funcall pot1 (1+ offset))
  270.       (funcall pot2 (1+ offset))
  271.     )
  272.     (let* ((const (/ (aref pvec2 offset)))
  273.            (pcount 0)
  274.            (pvec (make-array 10 :fill-pointer pcount :adjustable t))
  275.           )
  276.       #'(lambda (n)
  277.           (if (<= n pcount) pvec
  278.             (progn
  279.               (funcall pot1 (+ offset n))
  280.               (funcall pot2 (+ offset n))
  281.               (dotimes (i (- n pcount) pvec)
  282.                 (vector-push-extend
  283.                   (do ((j (1+ offset) (1+ j))
  284.                        (k (1- pcount) (1- k))
  285.                        (sum (aref pvec1 (+ offset pcount))
  286.                             (- sum (* (aref pvec2 j) (aref pvec k)))
  287.                       ))
  288.                       ((minusp k) (* const sum))
  289.                   )
  290.                   pvec
  291.                 )
  292.                 (setq pcount (1+ pcount))
  293. ) ) )   ) ) ) )
  294.  
  295. (defvar bernoulli
  296.   (potkehr (potshift (pot- exp \1) -1)) ; 1 / (e^x-1)/x
  297. )
  298. (defun bernoulli (n) ; n-te Bernoulli-Zahl
  299.   (* (aref (funcall bernoulli (1+ n)) n) (! n))
  300. )
  301.  
  302. ;; potexpt_Z erhebt eine Potenzreihe in die n-te Potenz (n ganze Zahl)
  303. (defun potexpt_Z (pot n)
  304.   (if (minusp n) (potexpt_Z (potkehr pot) (- n))
  305.     (if (zerop n) |1|
  306.       (if (oddp n) (pot*2 pot (potexpt_Z pot (1- n)))
  307.         (let ((pot1 (potexpt_Z pot (floor n 2))))
  308.           (pot*2 pot1 pot1)
  309. ) ) ) ) )
  310.  
  311. ; exp(Potenzreihe)
  312. ; Methode:
  313. ; exp(a1*x+a2*x^2+...) = b0+b1*x+b2*x^2+... wobei
  314. ; b[0]=1 und
  315. ; b[n] = (1 a[1] b[n-1] + ... + (n-1) a[n-1] b[1] + n a[n] b[0]) / n für n>0.
  316. (defun potexp (pot)
  317.   (let ((pvec1 (funcall pot 1)))
  318.     (unless (zerop (aref pvec1 0))
  319.       (error "Nur Potenzreihen ohne Absolutglied können exponenziert werden.")
  320.     )
  321.     (let* ((pcount 1)
  322.            (pvec (make-array 10 :fill-pointer pcount :adjustable t)))
  323.       (setf (aref pvec 0) 1)
  324.       #'(lambda (n)
  325.           (if (<= n pcount) pvec
  326.             (progn
  327.               (funcall pot n)
  328.               (loop
  329.                 ; Berechne b[pcount]:
  330.                 (vector-push-extend
  331.                   (do ((i 1 (1+ i)) (sum 0))
  332.                       ((> i pcount) (/ sum pcount))
  333.                     (incf sum (* i (aref pvec1 i) (aref pvec (- pcount i))))
  334.                   )
  335.                   pvec
  336.                 )
  337.                 (setq pcount (1+ pcount))
  338.                 (when (<= n pcount) (return pvec))
  339. ) ) )   ) ) ) )
  340.  
  341. ; Potenzreihe ^ e (e Zahl)
  342. ; Methode:
  343. ; f=a0+a1*x+a2*x^2+... mit a0=1 -> g:=f^e=b0+b1*x+b2*x^2+... mit
  344. ; g'/g = e * f'/f, also f g' = e f' g. Darin Koeffizient von x^(n-1):
  345. ; a[0] n b[n] + ... + a[n-1] 1 b[1] = e * ( 1 a[1] b[n-1] + ... + n a[n] b[0] ).
  346. ; b[0]=1 und
  347. ; b[n] = (((e+1)*1-n) a[1] b[n-1] + ... + ((e+1)*(n-1)-n) a[n-1] b[1] + (n-(e+1)*n) a[n] b[0]) / n für n>0.
  348. (defun potexpt_Q (pot e &aux (e1 (+ e 1)))
  349.   (let ((pvec1 (funcall pot 1)))
  350.     (unless (= (aref pvec1 0) 1)
  351.       (error "Nur Potenzreihen mit Absolutglied 1 können potenziert werden.")
  352.     )
  353.     (let* ((pcount 1)
  354.            (pvec (make-array 10 :fill-pointer pcount :adjustable t)))
  355.       (setf (aref pvec 0) (aref pvec1 0))
  356.       #'(lambda (n)
  357.           (if (<= n pcount) pvec
  358.             (progn
  359.               (funcall pot n)
  360.               (loop
  361.                 ; Berechne b[pcount]:
  362.                 (vector-push-extend
  363.                   (do ((i 1 (1+ i))
  364.                        (sum 0))
  365.                       ((> i pcount) (/ sum pcount))
  366.                     (incf sum (* (- (* e1 i) pcount) (aref pvec1 i) (aref pvec (- pcount i))))
  367.                   )
  368.                   pvec
  369.                 )
  370.                 (setq pcount (1+ pcount))
  371.                 (when (<= n pcount) (return pvec))
  372. ) ) )   ) ) ) )
  373.  
  374. ; log(Potenzreihe)
  375. ; Methode:
  376. ; log(a0+a1*x+a2*x^2+...) = b1*x+b2*x^2+... wobei
  377. ; a[0]=1 und
  378. ; b[n] = (n a[n] - 1 a[n-1] b[1] - ... - (n-1) a[1] b[n-1]) / n für n>0.
  379. (defun potlog (pot)
  380.   (let ((pvec1 (funcall pot 1)))
  381.     (unless (= (aref pvec1 0) 1)
  382.       (error "Nur Potenzreihen mit Absolutglied 1 können logarithmiert werden.")
  383.     )
  384.     (let* ((pcount 1)
  385.            (pvec (make-array 10 :fill-pointer pcount :adjustable t)))
  386.       (setf (aref pvec 0) 0)
  387.       #'(lambda (n)
  388.           (if (<= n pcount) pvec
  389.             (progn
  390.               (funcall pot n)
  391.               (loop
  392.                 ; Berechne b[pcount]:
  393.                 (vector-push-extend
  394.                   (do ((i 1 (1+ i))
  395.                        (sum (* pcount (aref pvec1 pcount))))
  396.                       ((>= i pcount) (/ sum pcount))
  397.                     (decf sum (* i (aref pvec1 (- pcount i)) (aref pvec i)))
  398.                   )
  399.                   pvec
  400.                 )
  401.                 (setq pcount (1+ pcount))
  402.                 (when (<= n pcount) (return pvec))
  403. ) ) )   ) ) ) )
  404.  
  405. ;; potderiv leitet eine Potenzreihe ab.
  406. (defun potderiv (pot)
  407.   (let* ((pcount 0)
  408.          (pvec (make-array 10 :fill-pointer pcount :adjustable t))
  409.          (pvec1 (funcall pot 0))
  410.         )
  411.     #'(lambda (n)
  412.         (if (<= n pcount) pvec
  413.           (progn
  414.             (funcall pot (1+ n))
  415.             (dotimes (i (- n pcount) pvec)
  416.               (vector-push-extend
  417.                 (* (1+ pcount) (aref pvec1 (1+ pcount)))
  418.                 pvec
  419.               )
  420.               (setq pcount (1+ pcount))
  421. ) )   ) ) ) )
  422.  
  423. ;; potint integriert eine Potenzreihe.
  424. ;; Optional: Vorgabe des konstanten Gliedes.
  425. (defun potint (pot &optional (const 0))
  426.   (let* ((pcount 1)
  427.          (pvec (make-array 10 :fill-pointer pcount :adjustable t))
  428.          (pvec1 (funcall pot 0))
  429.         )
  430.     (setf (aref pvec 0) const)
  431.     #'(lambda (n)
  432.         (if (<= n pcount) pvec
  433.           (progn
  434.             (funcall pot (1- n))
  435.             (dotimes (i (- n pcount) pvec)
  436.               (vector-push-extend (/ (aref pvec1 (1- pcount)) pcount) pvec)
  437.               (setq pcount (1+ pcount))
  438. ) )   ) ) ) )
  439.  
  440. ;; pot*hadamard bildet das Hadamard-Produkt (koeffizientenweise Produkt)
  441. ;; zweier Potenzreihen
  442. (defun pot*hadamard (pot1 pot2)
  443.   (let* ((pcount 0)
  444.          (pvec (make-array 10 :fill-pointer pcount :adjustable t))
  445.          (factor1 (funcall pot1 0))
  446.          (factor2 (funcall pot2 0))
  447.         )
  448.     #'(lambda (n)
  449.         (if (<= n pcount) pvec
  450.           (progn
  451.             (funcall pot1 n) (funcall pot2 n)
  452.             (dotimes (i (- n pcount) pvec)
  453.               (vector-push-extend
  454.                 (* (aref factor1 pcount) (aref factor2 pcount))
  455.                 pvec
  456.               )
  457.               (setq pcount (1+ pcount))
  458. ) )   ) ) ) )
  459.  
  460. ;; potcopy liefert eine Kopie der in einem Symbol enthaltenen Potenzreihe
  461. (defmacro potcopy (sym)
  462.   `(potcopy1 ',sym)
  463. )
  464. (defun potcopy1 (sym)
  465.   (let* ((pcount 0)
  466.          (pvec (make-array 10 :fill-pointer pcount :adjustable t))
  467.          (pvec1 nil)
  468.         )
  469.     #'(lambda (n)
  470.         (if (<= n pcount) pvec
  471.           (progn
  472.             (setq pvec1 (funcall (symbol-value sym) n))
  473.             (dotimes (i (- n pcount) pvec)
  474.               (vector-push-extend (aref pvec1 pcount) pvec)
  475.               (setq pcount (1+ pcount))
  476. ) )   ) ) ) )
  477. ; Anwendungsbeispiel: Lösung von Fixpunktgleichungen der Form y=F(y),
  478. ; wobei F streng kontrahierend ist (d.h. für y1-y2 = O(x^n) gilt
  479. ; F(y1)-F(y2) = O(x^(n+1)) ).
  480. ; Z.B. Löse   y = 1 + x y
  481. ;      durch  (setq y (pot+ |1| (potshift (potcopy y) 1)))
  482.  
  483. ;; potlindgl löst eine homogene lineare Differentialgleichung
  484. ;; mit Potenzreihenansatz.
  485. ;; Input: Gleichung ((a00 a01 a02 ...) (a10 a11 a12 ...) ...), was
  486. ;;          0 = (a00+a01*x+a02*x^2+...) * D^0 F(x)
  487. ;;              + (a10+a11*x+a12*x^2+...) * D^1 F(x)
  488. ;;              + ...
  489. ;;              + (ar0+ar1*x+ar2*x^2+...) * D^r F(x)
  490. ;;          bedeutet
  491. ;;          (wobei die Koeffizienten Polynome oder Potenzreihen sind),
  492. ;;        ausreichend viele Anfangskoeffizienten (b0 b1 ...) für F(x),
  493. ;;        wobei Differentialoperator D = d/dx.
  494. ;; Output: Lösungspotenzreihe F(x).
  495. ; Methode:
  496. ; Inspiziere für n=0,1,... den Koeffizienten von x^n in dieser Gleichung:
  497. ; sum(i=0..r, sum(j>=0, sum(k>=0, k-i+j=n, a[i,j]*k...(k-i+1)*b[k] ))) = 0.
  498. ; Sammle dabei die Koeffizienten aller b[k] (k=n+r, n+r-1, ...) und
  499. ; überprüfe, ob die entstehende Gleichung nach genau einem b[k] aufgelöst
  500. ; werden kann, der bis dato noch unbekannt ist.
  501. (defun potlindgl (gleichung anfang)
  502.   (let* ((r (1- (length gleichung)))
  503.          (pcount 0)
  504.          (pvec (make-array 10 :fill-pointer pcount :adjustable t)))
  505.     (dolist (b anfang)
  506.       (vector-push-extend b pvec)
  507.       (setq pcount (1+ pcount))
  508.     )
  509.     (labels
  510.           ((teilsumme (k L) ; liefert zu einer Liste L=(c[k] c[k-1] ...)
  511.              (let ((sum 0)) ; die Teilsumme c[k]*b[k]+c[k-1]*b[k-1]+...
  512.                (dolist (c L) (incf sum (* c (aref pvec k))) (decf k))
  513.                sum
  514.            ) )
  515.            (inspiziere (n) ; Koeffizienten von x^n inspizieren
  516.              (let ((bk-koeffs nil) ; Liste aller b[k]-Koeffs
  517.                    (max-k (+ n r)) ; höchstes auftretendes k
  518.                    (faktoren (reverse gleichung))) ; (a[r,*],...,a[0,*])
  519.                (do ((k max-k (1- k)))
  520.                    ((minusp k))
  521.                  ; Koeffizient von b[k] bestimmen:
  522.                  (let ((nochwas nil) (koeff 0))
  523.                    (do ((i r (1- i))
  524.                         (faktr faktoren (cdr faktr))) ; (a[i,*],...,a[0,*])
  525.                        ((minusp i))
  526.                      (let ((j (+ (- n k) i)) ; j mit j-i=n-k
  527.                            (prod ; k...(k-i+1)
  528.                              (do* ((x 1 (* x (- k y))) (y 0 (1+ y)))
  529.                                   ((= y i) x)
  530.                           )) )
  531.                        (when (minusp j) (return)) ; j<0 -> für dieses k fertig
  532.                        ; Koeffizient von x^j in a[i,*] = (car faktr) :
  533.                        (if (listp (car faktr))
  534.                          ; Polynom
  535.                          (when (car faktr) ; Polynom zu Ende -> nichts zu tun
  536.                            (setq nochwas t) ; sonst Koeffizienten verbrauchen
  537.                            (incf koeff (* prod (pop (car faktr))))
  538.                          )
  539.                          ; Potenzreihe
  540.                          (progn
  541.                            (setq nochwas t)
  542.                            (incf koeff (* prod (aref (funcall (car faktr) (1+ j)) j)))
  543.                    ) ) ) )
  544.                    (unless nochwas ; dieser und alle weiteren Koeffs von b[k]
  545.                      (return)      ; = 0 -> Schleife beenden
  546.                    )
  547.                    (push koeff bk-koeffs)
  548.                ) )
  549.                (setq bk-koeffs (nreverse bk-koeffs))
  550.                ; bk-koeffs = (Koeff von b[max-k],
  551.                ;              Koeff von b[max-k - 1],
  552.                ;              ...)
  553.                (loop
  554.                  (when (or (null bk-koeffs) (not (zerop (car bk-koeffs))))
  555.                    (return)
  556.                  )
  557.                  (decf max-k) (pop bk-koeffs)
  558.                )
  559.                (unless (null bk-koeffs)
  560.                  ; Die Gleichung bestimmt b[max-k] eindeutig
  561.                  (if (< max-k pcount)
  562.                    ; Gleichung nur überprüfen
  563.                    (unless (zerop (teilsumme max-k bk-koeffs))
  564.                      (error "Anfangskoeffizienten ~S erfüllen Differentialgleichung nicht."
  565.                             anfang
  566.                    ) )
  567.                    ; Gleichung ausnutzen
  568.                    (progn
  569.                      (when (> max-k pcount)
  570.                        (warn "Koeffizienten von x^~S bis x^~S werden willkürlich auf 0 gesetzt."
  571.                              pcount (1- max-k)
  572.                        )
  573.                        (loop
  574.                          (vector-push-extend 0 pvec)
  575.                          (setq pcount (1+ pcount))
  576.                          (when (= pcount max-k) (return))
  577.                      ) )
  578.                      ; pcount = max-k
  579.                      (vector-push-extend
  580.                        (- (/ (teilsumme (1- max-k) (cdr bk-koeffs))
  581.                              (car bk-koeffs)
  582.                        )  )
  583.                        pvec
  584.                      )
  585.                      (setq pcount (1+ pcount))
  586.           )) ) ) ) )
  587.       (let ((m 0)) ; alle x^n mit n<m sind bereits inspiziert
  588.         #'(lambda (n)
  589.             (if (<= n pcount) pvec
  590.               (loop
  591.                 (inspiziere m)
  592.                 (incf m)
  593.                 (when (<= n pcount) (return pvec))
  594.           ) ) )
  595. ) ) ) )
  596.  
  597. ;; liefert zu einer Potenzreihe G mit G(0)=0 die Substitution ( F -> F o G ).
  598. ;; Dies ist - genau genommen - eine Funktion, die zu einem n>=0 einen Array
  599. ;; der Länge n+1 liefert, dessen Element k (für 0<=k<=n) der Koeffizient
  600. ;; von x^n/n! in G(x)^k/k! ist.
  601. (defun make-subst (pot)
  602.   (let ((pvec (funcall pot 1)))
  603.     (unless (zerop (aref pvec 0))
  604.       (error "Nur Potenzreihen mit Absolutglied 0 können eingesetzt werden.")
  605.     )
  606.     (let ((tables (make-array 10 :fill-pointer 0 :adjustable t)))
  607.       (vector-push-extend '#(1) tables)
  608.       #'(lambda (n)
  609.           (when (>= n (fill-pointer tables))
  610.             (funcall pot (1+ n))
  611.             (do ((m (fill-pointer tables) (1+ m)))
  612.                 ((> m n))
  613.               ; Vektor der Koeffizienten von x^m/m! in G(x)^k/k! berechnen:
  614.               (let ((Vm (make-array (1+ m)))
  615.                     (k 0))
  616.                 (vector-push-extend Vm tables)
  617.                 (setf (aref Vm k) 0) ; Koeff. von x^m/m! in 1 ist 0
  618.                 (loop
  619.                   (when (= k m) (return))
  620.                   ; Um G(x)^(k+1)/(k+1)! zu bekommen, muß man G(x)^k/k! mit
  621.                   ; G(x) multiplizieren und durch (k+1) dividieren. Der
  622.                   ; Koeffizient von x^m/m! ist also
  623.                   ; = sum(i=k..m, m!/i! * (Koeff. von x^i/i! in G(x)^k/k!)
  624.                   ;                     * (Koeff. von x^(m-i) in G(x)) )
  625.                   ;   / (k+1)
  626.                   (let ((index m) (prod 1) (sum 0))
  627.                     (loop
  628.                       (incf sum (* prod (aref (aref tables index) k)
  629.                                         (aref pvec (- m index))
  630.                       )         )
  631.                       (when (= index k) (return))
  632.                       (setq prod (* prod index))
  633.                       (decf index)
  634.                     )
  635.                     (incf k)
  636.                     (setf (aref Vm k) (/ sum k))
  637.           ) ) ) ) )
  638.           (aref tables n)
  639. ) ) )   )
  640.  
  641. ;; potsubst setzt eine Potenzreihe in eine andere Potenzreihe ein: F o G,
  642. ;; wobei G(0)=0 ist. G muß in Form einer Substitution vorliegen.
  643. (defun potsubst (pot subst)
  644.   (let* ((pvec1 (funcall pot 0))
  645.          (pcount 0)
  646.          (pvec (make-array 10 :fill-pointer pcount :adjustable t)))
  647.     #'(lambda (n)
  648.         (if (<= n pcount) pvec
  649.           (progn
  650.             (funcall pot n)
  651.             (dotimes (i (- n pcount) pvec)
  652.               (vector-push-extend
  653.                 ; Koeffizient von x^pcount ist =
  654.                 ; sum(i=0..pcount, i!/pcount!*pot[i]*(Koeff.von x^pcount/pcount! in G(x)^i/i!))
  655.                 (let ((V (funcall subst pcount))
  656.                       (index pcount) (prod 1) (sum 0))
  657.                   (loop
  658.                     (incf sum (* prod (aref pvec1 index) (aref V index)))
  659.                     (when (zerop index) (return))
  660.                     (setq prod (/ prod index))
  661.                     (setq index (1- index))
  662.                   )
  663.                   sum
  664.                 )
  665.                 pvec
  666.               )
  667.               (setq pcount (1+ pcount))
  668. ) )   ) ) ) )
  669.  
  670. ;; eine Funktion zum Ausdrucken:
  671. (defun potprint (pot index
  672.                  &key (var "T") (shift 0) (printer nil) (width 79)
  673.                       (stream (if printer *printer-output* *standard-output*))
  674.                 )
  675.   ;; pot = Potenzreihe
  676.   ;; index = Grenze für den Ausdruck der Koeffizienten (ausschließlich)
  677.   ;; stream = Output-Stream
  678.   ;; var = String oder Symbol, Bezeichnung der Unbestimmten
  679.   ;; shift = Exponent des ersten Gliedes
  680.   ;; printer = Flag, ob Ausgabe auf Drucker (Exponenten hochstellen)
  681.   ;; width = Zeilenlänge
  682.   (let* ((pvec (funcall pot index)) ; Vektor mit den Koeffizienten
  683.          (start-flag t)  ; Flag für ersten Koeffizienten /= 0
  684.          (hochstring1 (if printer (coerce '(#\Escape #\S #\Null) 'string) "^"))
  685.              ; String vor dem Exponenten (für EPSON LQ-800)
  686.          (hochstring2 (if printer (coerce '(#\Escape #\T) 'string) ""))
  687.              ; String nach dem Exponenten (für EPSON LQ-800)
  688.          (hochsummand (if printer 0 1))
  689.              ; druckende Länge der beiden Strings
  690.          (vstring (string var))  ; Name der Variablen als String
  691.          (vlen (length vstring)) ; und seine Länge
  692.          (plusstring " + ")      ; Vorzeichenstrings mit Längen
  693.          (pluslen 3)
  694.          (minusstring " - ")
  695.          (minuslen 3)
  696.          (charcount 0)   ; Zähler für verbrauchte Zeichen in der Zeile
  697.         )
  698.     (macrolet ((out (string stringwidth) ; einen String gegebener Druckbreite
  699.                                          ; ausgeben
  700.                  `(progn
  701.                     (setq charcount (+ charcount ,stringwidth)) ; Zeichenzähler weiterzählen
  702.                     ; wenn's nicht mehr hinpaßt, neue Zeile anfangen:
  703.                     (when (> charcount width) (terpri stream) (setq charcount ,stringwidth))
  704.                     (write-string ,string stream)
  705.                   )
  706.               ))
  707.       (flet ((monom (kof expo) ; monom gibt ein Glied (mit Vorzeichen) aus
  708.                (unless (zerop kof) ; Monom mit Koeffizient=0 wird übergangen
  709.                  (let (str strlen) ; str = bisherige Zeichen fürs Monom,
  710.                                    ; strlen = druckende Länge von str
  711.                    (if (plusp kof)
  712.                      (progn ; Koeffizient >0
  713.                        (if start-flag ; evtl. Pluszeichen
  714.                          (setq start-flag nil str "" strlen 0)
  715.                          (setq str plusstring strlen pluslen)
  716.                      ) )
  717.                      (progn ; Koeffizient <0
  718.                        (setq start-flag nil)
  719.                        (setq str minusstring strlen minuslen) ; Minuszeichen
  720.                        (setq kof (abs kof))
  721.                    ) )
  722.                    (unless (and (= kof 1) (not (zerop expo)))
  723.                      (let ((kofstring (princ-to-string kof)))
  724.                        (setq str (string-concat str kofstring)
  725.                              strlen (+ strlen (length kofstring))
  726.                    ) ) )
  727.                    (unless (zerop expo)
  728.                      (if (= expo 1)
  729.                        (setq str (string-concat str vstring)
  730.                              strlen (+ strlen vlen)
  731.                        )
  732.                        (let ((expostring (princ-to-string expo)))
  733.                          (setq str (string-concat str vstring hochstring1 expostring hochstring2)
  734.                                strlen (+ strlen vlen hochsummand (length expostring))
  735.                    ) ) ) )
  736.                    (out str strlen)
  737.             )) ) )
  738.         (terpri stream)
  739.         (dotimes (i index)
  740.           (monom (aref pvec i) (+ i shift))
  741.         )
  742.         (out " + ..." 6)
  743.         (values)
  744. ) ) ) )
  745.  
  746. ;; potsolve erzeugt aus einer Gleichung X = F(X,T) mit F Polynom,
  747. ;; F = O(X^2) + O(T) eine Potenzreihe X = P(T).
  748. ;; F ist in der Form ( { (i j Aij) } ) anzugeben, wenn
  749. ;; F = SUMME(i,j; Aij*X^i*T^j).
  750. (defun potsolve (polyxt)
  751.   (let* ((pcount 1)
  752.          (pvec (make-array 10 :fill-pointer pcount :adjustable t))
  753.          glieder
  754.         )
  755.     (setf (aref pvec 0) 0)
  756.     (labels ((result (n)
  757.                (if (<= n pcount) pvec
  758.                  (let ((summanden
  759.                          (mapcar #'(lambda (pot) (funcall pot n)) glieder)
  760.                       ))
  761.                    (dotimes (i (- n pcount) pvec)
  762.                      (vector-push-extend
  763.                        (do ((s summanden (rest s))
  764.                             (sum 0 (+ sum (aref (first s) pcount)))
  765.                            )
  766.                            ((null s) sum)
  767.                        )
  768.                        pvec
  769.                      )
  770.                      (setq pcount (1+ pcount))
  771.              ) ) ) )
  772.              (make-glied (koeff &aux (i (first koeff)) (j (second koeff))
  773.                                      (a (third koeff))
  774.                          )
  775.                (unless (or (>= i 2) (>= j 1))
  776.                        (error "~: Polynom enthält falsches Glied ~"
  777.                               'potsolve koeff
  778.                )       )
  779.                (pot*skal (potshift (potexpt_Z #'result i) j) a)
  780.             ))
  781.       (setq glieder (mapcar #'make-glied polyxt))
  782.       #'result
  783. ) ) )
  784.  
  785. ;; Liefert die Dirichlet-Faltung zweier Folgen (a1,a2,a3,...), die
  786. ;; als Potenzreihen a1+a2*X+a3*X^2+... dargestellt sind.
  787. (defun pot*dirichlet (folge1 folge2)
  788.   (let* ((pcount 0)
  789.          (pvec (make-array 10 :fill-pointer pcount :adjustable t))
  790.          (pvec1 (funcall folge1 0))
  791.          (pvec2 (funcall folge2 0))
  792.         )
  793.     #'(lambda (n)
  794.         (if (<= n pcount) pvec
  795.           (progn
  796.             (funcall folge1 n) (funcall folge2 n)
  797.             (dotimes (i (- n pcount) pvec)
  798.               (vector-push-extend
  799.                 (let ((m (1+ pcount)))
  800.                   ; Koeffizient von X^(m-1) ist c[m] = sum(p*q=m, a[p]*b[q])
  801.                   (do ((p m (1- p))
  802.                        (sum 0))
  803.                       ((zerop p) sum)
  804.                     (multiple-value-bind (q r) (floor m p)
  805.                       (when (zerop r)
  806.                         (incf sum
  807.                               (* (aref pvec1 (1- p)) (aref pvec2 (1- q)))
  808.                 ) ) ) ) )
  809.                 pvec
  810.               )
  811.               (setq pcount (1+ pcount))
  812. ) )   ) ) ) )
  813.  
  814. ;; Liefert den Quotienten bezüglich Dirichlet-Faltung
  815. ;; zweier Folgen (a1,a2,a3,...), die
  816. ;; als Potenzreihen a1+a2*X+a3*X^2+... dargestellt sind.
  817. (defun pot/dirichlet (folge1 folge2)
  818.   ; (c1,c2,c3,...) := (a1,a2,a3,...) / (b1,b2,b3,...)
  819.   ; bedeutet  (a1,a2,a3,...) = (b1,b2,b3,...) * (c1,c2,c3,...) , also
  820.   ; c[m] = (a[m] - sum(p*q=m,p>1,q<m, b[p]*c[q]))/b[1] für m=1,2,3,...
  821.   (let* ((pcount 0)
  822.          (pvec (make-array 10 :fill-pointer pcount :adjustable t))
  823.          (pvec1 (funcall folge1 0))
  824.          (pvec2 (funcall folge2 1))
  825.          (const (/ (aref pvec2 0)))
  826.         )
  827.     #'(lambda (n)
  828.         (if (<= n pcount) pvec
  829.           (progn
  830.             (funcall folge1 n) (funcall folge2 n)
  831.             (dotimes (i (- n pcount) pvec)
  832.               (vector-push-extend
  833.                 (let ((m (1+ pcount)))
  834.                   ; Koeffizient von X^(m-1) ist c[m] = (a[m]-...)/b[1]
  835.                   (do ((q (1- m) (1- q))
  836.                        (sum (aref pvec1 (1- m))))
  837.                       ((zerop q) (* sum const))
  838.                     (multiple-value-bind (p r) (floor m q)
  839.                       (when (zerop r)
  840.                         (decf sum
  841.                               (* (aref pvec2 (1- p)) (aref pvec (1- q)))
  842.                 ) ) ) ) )
  843.                 pvec
  844.               )
  845.               (setq pcount (1+ pcount))
  846. ) )   ) ) ) )
  847.  
  848. ;; Möbius-Funktion als Folge ist das Dirichlet-Faltungs-Inverse der
  849. ;; konstanten Folge (defpot 1).
  850. (defvar möbius (pot/dirichlet \1 (defpot 1)))
  851.  
  852. ;; potlogdiff bildet die logarithmische Ableitung einer Potenzreihe /=0.
  853. (defun potlogderiv (pot)
  854.   (pot/ (potderiv pot) pot)
  855. )
  856.  
  857. ;; potkanprodukt entwickelt eine Potenzreihe aus 1+Z[[X]] in ein
  858. ;; "kanonisches Produkt"  (1-X)^c1*(1-X^2)^c2*...  und liefert (c1,c2,...).
  859. (defun potkanprodukt (pot)
  860.   (let ((pvec (funcall pot 1)))
  861.     (unless (= (aref pvec 0) 1)
  862.       (error "Nur Potenzreihen mit Absolutglied 1 sind in ein kanonisches Produkt entwickelbar.")
  863.   ) )
  864.   ; Ziel: pot als prod(i>=1, (1-X^i)^c[i]) schreiben.
  865.   ; 1. Schritt: logarithmische Ableitung, liefert
  866.   ;             sum(i>=1, c[i] * d/dX log(1-X^i) )
  867.   ;             = sum(i>=1, -c[i]*i*X^(i-1)/(1-X^i) )
  868.   ;             = 1/X * sum(i>=1, -c[i]*i*sum(j>=1,X^(i*j)) )
  869.   ;             = 1/X * sum(n>=1, sum(i*j=n,-c[i]*i) * X^n)
  870.   ; 2. Schritt: Dies entspricht der Folge ( sum(i*j=n,-c[i]*i) )[n>=1]
  871.   ;             was die Dirichlet-Faltung der Folgen  ( -1 )[j>=1]
  872.   ;             und  ( c[i]*i )[i>=1]  ist. Division bzgl. Dirichlet-Faltung
  873.   ;             liefert also die Folge  ( c[i]*i )[i>=1]  als Reihe
  874.   ;             sum(i>=1, c[i]*i*X^(i-1) ).
  875.   ; 3. Schritt: Integrieren liefert  sum(i>=1, c[i]*X^i) .
  876.   ;             Dies entspricht der Folge (0,c1,c2,...).
  877.   ; 4. Schritt: Division durch X, und man bekommt die Folge (c1,c2,...).
  878.   (potshift (potint (pot/dirichlet (potlogderiv pot) (defpot -1))) -1)
  879. )
  880.  
  881. ;; binom liefert die Potenzreihe (1+T)^e.
  882. (defun potbinom (e)
  883.   (let* ((pcount 0)
  884.          (pvec (make-array 10 :fill-pointer pcount :adjustable t))
  885.          (koeff 1)
  886.          (h e)
  887.         )
  888.     #'(lambda (n)
  889.         (if (<= n pcount) pvec
  890.           (dotimes (i (- n pcount) pvec)
  891.             (vector-push-extend koeff pvec)
  892.             (setq pcount (1+ pcount)
  893.                   koeff (* koeff (/ h pcount))
  894.                   h (1- h)
  895. ) )   ) ) ) )
  896.  
  897. (defun rn_n (r) ; asymptotische Entwicklung von (rn)!/n!^r, x=1/n
  898.   (potexp
  899.     (defpot (if (zerop pcount)
  900.               0
  901.               (* (/ (bernoulli (1+ pcount)) pcount (1+ pcount))
  902.                  (- (/ (expt r pcount)) r)
  903.             ) )
  904. ) ) )
  905.  
  906. (defun n_nr (r) ; asymptotische Entwicklung von n!/(n/r)!^r, x=1/n
  907.   (potexp
  908.     (defpot (if (zerop pcount)
  909.               0
  910.               (* (/ (bernoulli (1+ pcount)) pcount (1+ pcount))
  911.                  (- 1 (expt r (1+ pcount)))
  912.             ) )
  913. ) ) )
  914.  
  915. ; wandelt eine Reihe sum(n>=0, c[n+r]/n!*x^n) in eine Reihe c[0]+c[1]*x+...,
  916. ; wobei die ersten r Koeffizienten angegeben werden müssen.
  917. (defun asympt-Reihe (pot &optional (anfang '(1)))
  918.   (let* ((r (length anfang))
  919.          (pvec1 (funcall pot 0))
  920.          (pcount 0)
  921.          (pvec (make-array 10 :fill-pointer pcount :adjustable t)))
  922.     (dolist (c anfang)
  923.       (vector-push-extend c pvec)
  924.       (setq pcount (1+ pcount))
  925.     )
  926.     #'(lambda (n)
  927.         (if (<= n pcount) pvec
  928.           (progn
  929.             (funcall pot (- n r))
  930.             (dotimes (i (- n pcount) pvec)
  931.               (vector-push-extend
  932.                 (let ((m (- pcount r))) ; Koeff. von x^m ist c[m+r]/m!
  933.                   (* (aref pvec1 m) (! m)) ; mal m! -> liefert c[m+r]=c[pcount]
  934.                 )
  935.                 pvec
  936.               )
  937.               (setq pcount (1+ pcount))
  938. ) )   ) ) ) )
  939.  
  940.  
  941. (setq hn2 ; asymptotische Entwicklung von Hn(2)
  942.   (asympt-Reihe ; F(x)=sum(n>=0,c[n+1]/n!) löst die Differentialgleichung
  943.     (potlindgl ; (2*e^2x+2*e^x+3)*F(x) + (2*e^2x-2*e^x-8)*F'(x) + (-4*e^x+4)*F''(x) = 0
  944.       (let ((e^x exp) (e^2x (pot* exp exp)))
  945.         (list (pot+ (pot*skal e^2x 2) (pot*skal e^x 2) (pot*skal |1| 3))
  946.               (pot+ (pot*skal e^2x 2) (pot*skal e^x -2) (pot*skal |1| -8))
  947.               (pot+ (pot*skal e^x -4) (pot*skal |1| 4))
  948.       ) )
  949.       (list 1/4)
  950. ) ) )
  951.  
  952. (require 'stirling) ; Funktionen stirling1-table, stirling2-table
  953. ; Verwende  #'stirling1-table  für die Substitution  x=log(1+t)
  954. ; und       #'stirling2-table  für die Substitution  t=e^x-1 .
  955.  
  956. ; Stellt fest, ob eine asymptotische Reihe von einem Zp-Maß herkommt:
  957. ; Asymptotische Reihe sum(k>=0,c[k]/n^k) = exp( sum(k>=1,b[k]/k/n^k) )
  958. ; kommt genau dann von einem Maß, wenn sum(k>=1,b[k]/k!*x^k) in Zp[[e^x-1]]
  959. ; liegt. Wir liefern die Potenzreihe, bei der hierin t=e^x-1 substituiert
  960. ; ist. Es kann leicht gesehen werden, für welche p sie in Zp[[t]] liegt.
  961. (defun maßtest (asym)
  962.   (potsubst
  963.     (pot*hadamard
  964.       (potlog asym) ; sum(k>=1,b[k]/k*x^k)
  965.       (potshift exp 1) ; sum(k>=1,1/(k-1)!*x^k)
  966.     ) ; Hadamard-Produkt = sum(k>=1,b[k]/k!*x^k)
  967.     #'stirling1-table ; darin x=log(1+t) einsetzen
  968. ) )
  969.  
  970. #|
  971.  
  972. ;; Modulfunktionen, insbesondere Eisensteinreihen:
  973.  
  974. ; (sigma s n) = sum(d|n,d^s)
  975. (defun sigma (s n)
  976.   (let ((n^1/2 (isqrt n)))
  977.     (do ((d1 1 (1+ d1))
  978.          (sum 0))
  979.         ((> d1 n^1/2) sum)
  980.       (when (zerop (mod n d1))
  981.         (incf sum (expt d1 s)) ; Teiler <=n^1/2 berücksichtigen
  982.         (let ((d2 (floor n d1)))
  983.           (when (> d2 d1)
  984.             (incf sum (expt d2 s)) ; Teiler >n^1/2 berücksichtigen
  985. ) ) ) ) ) )
  986.  
  987. ; (eisenstein k) liefert die Potenzreihe Ek(q) = 1/(2*zeta(2k)) * Gk(q),
  988. ; Gk(q) = 2*zeta(2k) + 2*(-1)^k*(2pi)^2k/(2k-1)!*sum(n=1..inf,sigma(2k-1,n)*q^n)
  989. ; mit zeta(2k) = 2^(2k-1)*pi^2k*abs(bernoulli(2k))/(2k)!
  990.  
  991. ; 8k*(-1)^k abs(bernoulli(2k))
  992. (defun eisenstein (k &aux (2k-1 (- (* 2 k) 1)))
  993.   (let* ((pcount 1)
  994.          (pvec (make-array 1 :fill-pointer pcount :adjustable t))
  995.          (koeff (/ (* 4 k (expt -1 k)) (abs (bernoulli (* 2 k)))))
  996.         )
  997.     (setf (aref pvec 0) 1)
  998.     #'(lambda (n)
  999.         (if (<= n pcount) pvec
  1000.           (dotimes (i (- n pcount) pvec)
  1001.             (vector-push-extend (* koeff (sigma 2k-1 pcount)) pvec)
  1002.             (setq pcount (1+ pcount))
  1003. ) )   ) ) )
  1004.  
  1005. ; Delta-Reihe:  = g2^3-27*g3^2 = (2pi)^12/(2^6*3^3)*(E2^3-E3^2)
  1006. ;                = (2pi)^12*Delta
  1007. (setq Delta
  1008.   (pot*skal (pot- (potexpt_Z (eisenstein 2) 3) (potexpt_Z (eisenstein 3) 2))
  1009.             1/1728
  1010. ) )
  1011.  
  1012. ; DeltaX(q) = Prod(n=1..inf,1-q^n)^X
  1013. (setq Delta24 (potshift Delta -1))
  1014. (setq Delta3 (potexpt_Q Delta24 1/8))
  1015. (setq Delta1 (potexpt_Q Delta24 1/24))
  1016.  
  1017. |#
  1018.  
  1019.