GEMini Atari/GEMini_Atari_CD-ROM_Walnut_Creek_December_1993.iso/files/language/examples/potr.lsp

Text File | 1993-10-23 | 39KB | 1,019 lines

;;;; Funktionen zum Arbeiten mit Potenzreihen, 11. 10. 1988
;;;; erweitert am 24.8.1989, 30.11.1989, 28.1.1990 von Bruno Haible
;;;; erweitert am 2.8.1990 von Michael Stoll und Bruno Haible
;;;; erweitert am 10.8.1990 von Bruno Haible

(provide 'potr)

;;; Potenzreihen sind funktionale Objekte, die mit einer nichtnegativen
;;; ganzen Zahl als Argument aufgerufen werden können und dann als
;;; Wert einen Vektor (adjustable, mit Fillpointer) zurückgeben,
;;; der die Koeffizienten der Potenzreihe enthält und dessen Länge
;;; mindestens so groß ist wie das angegebene Argument. Eine Potenzreihe
;;; benutzt immer denselben Vektor, so daß zurückgegebene Werte bei
;;; weiterer Berechnung mitgeführt werden.

;;; Variablenbezeichnungen:
;;; pvec ist stets der Vektor für die Koeffizienten,
;;; pcount ist die aktuelle Länge des Vektors (d.h. die Koeffizienten
;;;        bis pcount-1 sind bekannt, der nächste hat Index pcount)


;; |0| ist die Nullreihe
(defvar \0
  (let* ((pcount 0)
         (pvec (make-array 10 :fill-pointer pcount :adjustable t))
        )
    #'(lambda (n)
        (if (<= n pcount) pvec
          (dotimes (i (- n pcount) pvec)
            (vector-push-extend 0 pvec)
            (setq pcount (1+ pcount))
) )   ) ) )

;; |1| ist die Einsreihe
(defvar \1
  (let* ((pcount 1)
         (pvec (make-array 10 :fill-pointer pcount :adjustable t))
        )
    (setf (aref pvec 0) 1)
    #'(lambda (n)
        (if (<= n pcount) pvec
          (dotimes (i (- n pcount) pvec)
            (vector-push-extend 0 pvec)
            (setq pcount (1+ pcount))
) )   ) ) )

;; Mit defpot kann man Reihen wie exp definieren
;; kof-form ist eine Form für den einzusetzenden Koeffizienten,
;;  dabei ist der Index als pcount erreichbar.
;; Ist ein kof-init angegeben, dann wird eine Variable kof definiert
;;  und auf den angegebenen Wert gesetzt (sie kann dann in kof-form
;;  verwendet werden).
;; Ist kof-step (eine Form) angegeben, dann wird nach der Berechnung
;;  und Abspeicherung des Koeffizienten kof auf den Wert der Form
;;  kof-step gesetzt. pcount hat dabei den nächsten Index als Wert.
(defmacro defpot (kof-form &optional kof-init kof-step)
  `(let* ((pcount 0)
          (pvec (make-array 10 :fill-pointer pcount :adjustable t))
          ,@(if kof-init `((kof ,kof-init)) nil)
         )
     #'(lambda (n)
         (if (<= n pcount) pvec
           (dotimes (i (- n pcount) pvec)
             (vector-push-extend ,kof-form pvec)
             (setq pcount (1+ pcount)
                   ,@(if kof-step `(kof ,kof-step) nil)
)  )   ) ) ) )

(defvar exp
  (defpot kof 1 (/ kof pcount))
)

(defvar sin
  (defpot (if (oddp pcount)
            (if (evenp (floor pcount 2)) kof (- kof))
            0
          )
          1 (/ kof pcount)
) )

(defvar cos
  (defpot (if (evenp pcount)
            (if (evenp (floor pcount 2)) kof (- kof))
            0
          )
          1 (/ kof pcount)
) )

(defvar log
  (defpot (if (zerop pcount) 0
            (if (oddp pcount) (/ pcount) (- (/ pcount)))
) )       )

;; pot+ addiert Potenzreihen
(defun pot+ (&rest pots)
  (if (null pots) |0|
    (if (null (rest pots)) (first pots)
      (let* ((pcount 0)
             (pvec (make-array 10 :fill-pointer pcount :adjustable t))
             (summanden (mapcar #'(lambda (pot) (funcall pot 0)) pots))
            )
        #'(lambda (n)
            (if (<= n pcount) pvec
              (progn
                (dolist (pot pots) (funcall pot n))
                (dotimes (i (- n pcount) pvec)
                  (vector-push-extend
                    (do ((s summanden (rest s))
                         (sum 0 (+ sum (aref (first s) pcount)))
                        )
                        ((null s) sum)
                    )
                    pvec
                  )
                  (setq pcount (1+ pcount))
) ) ) )   ) ) ) )

;; pot- subtrahiert Potenzreihen
(defun pot- (pot1 &rest pots)
  (let* ((pcount 0)
         (pvec (make-array 10 :fill-pointer pcount :adjustable t))
         (erste (funcall pot1 0))
         (subtrahenden (mapcar #'(lambda (pot) (funcall pot 0)) pots))
        )
    #'(lambda (n)
        (if (<= n pcount) pvec
          (progn
            (funcall pot1 n)
            (dolist (pot pots) (funcall pot n))
            (dotimes (i (- n pcount) pvec)
              (vector-push-extend
                (do ((s subtrahenden (rest s))
                     (sum (aref erste pcount)
                          (- sum (aref (first s) pcount))
                    ))
                    ((null s) sum)
                )
                pvec
              )
              (setq pcount (1+ pcount))
) )   ) ) ) )

;; pot*skal multipliziert eine Potenzreihe mit einem skalaren Faktor
(defun pot*skal (pot skal)
  (let* ((pcount 0)
         (pvec (make-array 10 :fill-pointer pcount :adjustable t))
         (pvec1 (funcall pot 0))
        )
    #'(lambda (n)
        (if (<= n pcount) pvec
          (progn
            (funcall pot n)
            (dotimes (i (- n pcount) pvec)
              (vector-push-extend (* skal (aref pvec1 pcount)) pvec)
              (setq pcount (1+ pcount))
) )   ) ) ) )

;; potshift multipliziert eine Potenzreihe mit einer Potenz der Unbestimmten
(defun potshift (pot shift)
  (unless (integerp shift)
          (error "~: Der Shift muß eine nichtnegative ganze Zahl sein, nicht ~"
                 'potshift shift
  )       )
  (if (zerop shift)
    pot
    (let ((pvec1 (funcall pot (max 0 (- shift)))))
      (when (minusp shift)
        (dotimes (i (- shift))
          (unless (zerop (aref pvec1 i))
            (error "~: Potenzreihe ~ nicht durch x^~ teilbar"
                   'potshift pvec1 (- shift)
      ) ) ) )
      (let* ((pcount 0)
             (pvec (make-array 10 :fill-pointer pcount :adjustable t)))
        (when (plusp shift)
          (dotimes (i shift)
            (vector-push-extend 0 pvec) (setq pcount (1+ pcount))
        ) )
        #'(lambda (n)
            (if (<= n pcount) pvec
              (progn
                (funcall pot (- n shift))
                (dotimes (i (- n pcount) pvec)
                  (vector-push-extend (aref pvec1 (- pcount shift)) pvec)
                  (setq pcount (1+ pcount))
) ) ) )   ) ) ) )

;; pot*2 multipliziert zwei Potenzreihen
(defun pot*2 (pot1 pot2)
  (let* ((pcount 0)
         (pvec (make-array 10 :fill-pointer pcount :adjustable t))
         (factor1 (funcall pot1 0))
         (factor2 (funcall pot2 0))
        )
    #'(lambda (n)
        (if (<= n pcount) pvec
          (progn
            (funcall pot1 n) (funcall pot2 n)
            (dotimes (i (- n pcount) pvec)
              (vector-push-extend
                (do ((j 0 (1+ j))
                     (k pcount (1- k))
                     (sum 0 (+ sum (* (aref factor1 j) (aref factor2 k))))
                    )
                    ((minusp k) sum)
                )
                pvec
              )
              (setq pcount (1+ pcount))
) )   ) ) ) )

;; pot* multipliziert (beliebig viele) Potenzreihen
(defun pot* (&rest pots)
  (if (null pots) |1|
    (if (null (rest pots)) (first pots)
      (pot*2 (first pots) (apply #'pot* (rest pots)))
) ) )

;; potkehr bildet den Kehrwert einer invertierbaren Potenzreihe
(defun potkehr (pot)
  (let* ((pvec1 (funcall pot 1))
         (const (aref pvec1 0))
        )
    (when (zerop const)
          (error "~: Die Potenzreihe ~ ist nicht invertierbar"
                 'potkehr pvec1
    )     )
    (let* ((pcount 1)
           (pvec (make-array 10 :fill-pointer pcount :adjustable t))
           (const1 (- (/ const)))
          )
      (setf (aref pvec 0) (/ const))
      #'(lambda (n)
          (if (<= n pcount) pvec
            (progn
              (funcall pot n)
              (dotimes (i (- n pcount) pvec)
                (vector-push-extend
                  (do ((j 1 (1+ j))
                       (k (1- pcount) (1- k))
                       (sum 0 (+ sum (* (aref pvec1 j) (aref pvec k))))
                      )
                      ((minusp k) (* const1 sum))
                  )
                  pvec
                )
                (setq pcount (1+ pcount))
) ) )   ) ) ) )

;; pot/ dividiert zwei Potenzreihen, falls möglich
;; (Bem.: 0/0 führt auf Endlosschleife, das läßt sich nicht vermeiden)
; Methode:
; pot1 = x^offset*(a0+a1*x+...),
; pot2 = x^offset*(b0+b1*x+...) mit b0 /= 0,
; -> pot1/pot2 = c0+c1*x+... wobei
;    c[n] = (a[n]-(b[n]c[0]+...+b[1]c[n-1]))/b[0] für n=0,1,...
(defun pot/ (pot1 pot2)
  (let* ((pvec1 (funcall pot1 1))
         (pvec2 (funcall pot2 1))
         (offset 0)
        )
    (loop
      (unless (zerop (aref pvec2 offset)) (return))
      (unless (zerop (aref pvec1 offset))
              (error "~: Potenzreihe ~ läßt sich nicht durch ~ dividieren"
                     'pot/ pvec1 pvec2
      )       )
      (setq offset (1+ offset))
      (funcall pot1 (1+ offset))
      (funcall pot2 (1+ offset))
    )
    (let* ((const (/ (aref pvec2 offset)))
           (pcount 0)
           (pvec (make-array 10 :fill-pointer pcount :adjustable t))
          )
      #'(lambda (n)
          (if (<= n pcount) pvec
            (progn
              (funcall pot1 (+ offset n))
              (funcall pot2 (+ offset n))
              (dotimes (i (- n pcount) pvec)
                (vector-push-extend
                  (do ((j (1+ offset) (1+ j))
                       (k (1- pcount) (1- k))
                       (sum (aref pvec1 (+ offset pcount))
                            (- sum (* (aref pvec2 j) (aref pvec k)))
                      ))
                      ((minusp k) (* const sum))
                  )
                  pvec
                )
                (setq pcount (1+ pcount))
) ) )   ) ) ) )

(defvar bernoulli
  (potkehr (potshift (pot- exp \1) -1)) ; 1 / (e^x-1)/x
)
(defun bernoulli (n) ; n-te Bernoulli-Zahl
  (* (aref (funcall bernoulli (1+ n)) n) (! n))
)

;; potexpt_Z erhebt eine Potenzreihe in die n-te Potenz (n ganze Zahl)
(defun potexpt_Z (pot n)
  (if (minusp n) (potexpt_Z (potkehr pot) (- n))
    (if (zerop n) |1|
      (if (oddp n) (pot*2 pot (potexpt_Z pot (1- n)))
        (let ((pot1 (potexpt_Z pot (floor n 2))))
          (pot*2 pot1 pot1)
) ) ) ) )

; exp(Potenzreihe)
; Methode:
; exp(a1*x+a2*x^2+...) = b0+b1*x+b2*x^2+... wobei
; b[0]=1 und
; 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.
(defun potexp (pot)
  (let ((pvec1 (funcall pot 1)))
    (unless (zerop (aref pvec1 0))
      (error "Nur Potenzreihen ohne Absolutglied können exponenziert werden.")
    )
    (let* ((pcount 1)
           (pvec (make-array 10 :fill-pointer pcount :adjustable t)))
      (setf (aref pvec 0) 1)
      #'(lambda (n)
          (if (<= n pcount) pvec
            (progn
              (funcall pot n)
              (loop
                ; Berechne b[pcount]:
                (vector-push-extend
                  (do ((i 1 (1+ i)) (sum 0))
                      ((> i pcount) (/ sum pcount))
                    (incf sum (* i (aref pvec1 i) (aref pvec (- pcount i))))
                  )
                  pvec
                )
                (setq pcount (1+ pcount))
                (when (<= n pcount) (return pvec))
) ) )   ) ) ) )

; Potenzreihe ^ e (e Zahl)
; Methode:
; f=a0+a1*x+a2*x^2+... mit a0=1 -> g:=f^e=b0+b1*x+b2*x^2+... mit
; g'/g = e * f'/f, also f g' = e f' g. Darin Koeffizient von x^(n-1):
; a[0] n b[n] + ... + a[n-1] 1 b[1] = e * ( 1 a[1] b[n-1] + ... + n a[n] b[0] ).
; b[0]=1 und
; 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.
(defun potexpt_Q (pot e &aux (e1 (+ e 1)))
  (let ((pvec1 (funcall pot 1)))
    (unless (= (aref pvec1 0) 1)
      (error "Nur Potenzreihen mit Absolutglied 1 können potenziert werden.")
    )
    (let* ((pcount 1)
           (pvec (make-array 10 :fill-pointer pcount :adjustable t)))
      (setf (aref pvec 0) (aref pvec1 0))
      #'(lambda (n)
          (if (<= n pcount) pvec
            (progn
              (funcall pot n)
              (loop
                ; Berechne b[pcount]:
                (vector-push-extend
                  (do ((i 1 (1+ i))
                       (sum 0))
                      ((> i pcount) (/ sum pcount))
                    (incf sum (* (- (* e1 i) pcount) (aref pvec1 i) (aref pvec (- pcount i))))
                  )
                  pvec
                )
                (setq pcount (1+ pcount))
                (when (<= n pcount) (return pvec))
) ) )   ) ) ) )

; log(Potenzreihe)
; Methode:
; log(a0+a1*x+a2*x^2+...) = b1*x+b2*x^2+... wobei
; a[0]=1 und
; b[n] = (n a[n] - 1 a[n-1] b[1] - ... - (n-1) a[1] b[n-1]) / n für n>0.
(defun potlog (pot)
  (let ((pvec1 (funcall pot 1)))
    (unless (= (aref pvec1 0) 1)
      (error "Nur Potenzreihen mit Absolutglied 1 können logarithmiert werden.")
    )
    (let* ((pcount 1)
           (pvec (make-array 10 :fill-pointer pcount :adjustable t)))
      (setf (aref pvec 0) 0)
      #'(lambda (n)
          (if (<= n pcount) pvec
            (progn
              (funcall pot n)
              (loop
                ; Berechne b[pcount]:
                (vector-push-extend
                  (do ((i 1 (1+ i))
                       (sum (* pcount (aref pvec1 pcount))))
                      ((>= i pcount) (/ sum pcount))
                    (decf sum (* i (aref pvec1 (- pcount i)) (aref pvec i)))
                  )
                  pvec
                )
                (setq pcount (1+ pcount))
                (when (<= n pcount) (return pvec))
) ) )   ) ) ) )

;; potderiv leitet eine Potenzreihe ab.
(defun potderiv (pot)
  (let* ((pcount 0)
         (pvec (make-array 10 :fill-pointer pcount :adjustable t))
         (pvec1 (funcall pot 0))
        )
    #'(lambda (n)
        (if (<= n pcount) pvec
          (progn
            (funcall pot (1+ n))
            (dotimes (i (- n pcount) pvec)
              (vector-push-extend
                (* (1+ pcount) (aref pvec1 (1+ pcount)))
                pvec
              )
              (setq pcount (1+ pcount))
) )   ) ) ) )

;; potint integriert eine Potenzreihe.
;; Optional: Vorgabe des konstanten Gliedes.
(defun potint (pot &optional (const 0))
  (let* ((pcount 1)
         (pvec (make-array 10 :fill-pointer pcount :adjustable t))
         (pvec1 (funcall pot 0))
        )
    (setf (aref pvec 0) const)
    #'(lambda (n)
        (if (<= n pcount) pvec
          (progn
            (funcall pot (1- n))
            (dotimes (i (- n pcount) pvec)
              (vector-push-extend (/ (aref pvec1 (1- pcount)) pcount) pvec)
              (setq pcount (1+ pcount))
) )   ) ) ) )

;; pot*hadamard bildet das Hadamard-Produkt (koeffizientenweise Produkt)
;; zweier Potenzreihen
(defun pot*hadamard (pot1 pot2)
  (let* ((pcount 0)
         (pvec (make-array 10 :fill-pointer pcount :adjustable t))
         (factor1 (funcall pot1 0))
         (factor2 (funcall pot2 0))
        )
    #'(lambda (n)
        (if (<= n pcount) pvec
          (progn
            (funcall pot1 n) (funcall pot2 n)
            (dotimes (i (- n pcount) pvec)
              (vector-push-extend
                (* (aref factor1 pcount) (aref factor2 pcount))
                pvec
              )
              (setq pcount (1+ pcount))
) )   ) ) ) )

;; potcopy liefert eine Kopie der in einem Symbol enthaltenen Potenzreihe
(defmacro potcopy (sym)
  `(potcopy1 ',sym)
)
(defun potcopy1 (sym)
  (let* ((pcount 0)
         (pvec (make-array 10 :fill-pointer pcount :adjustable t))
         (pvec1 nil)
        )
    #'(lambda (n)
        (if (<= n pcount) pvec
          (progn
            (setq pvec1 (funcall (symbol-value sym) n))
            (dotimes (i (- n pcount) pvec)
              (vector-push-extend (aref pvec1 pcount) pvec)
              (setq pcount (1+ pcount))
) )   ) ) ) )
; Anwendungsbeispiel: Lösung von Fixpunktgleichungen der Form y=F(y),
; wobei F streng kontrahierend ist (d.h. für y1-y2 = O(x^n) gilt
; F(y1)-F(y2) = O(x^(n+1)) ).
; Z.B. Löse   y = 1 + x y
;      durch  (setq y (pot+ |1| (potshift (potcopy y) 1)))

;; potlindgl löst eine homogene lineare Differentialgleichung
;; mit Potenzreihenansatz.
;; Input: Gleichung ((a00 a01 a02 ...) (a10 a11 a12 ...) ...), was
;;          0 = (a00+a01*x+a02*x^2+...) * D^0 F(x)
;;              + (a10+a11*x+a12*x^2+...) * D^1 F(x)
;;              + ...
;;              + (ar0+ar1*x+ar2*x^2+...) * D^r F(x)
;;          bedeutet
;;          (wobei die Koeffizienten Polynome oder Potenzreihen sind),
;;        ausreichend viele Anfangskoeffizienten (b0 b1 ...) für F(x),
;;        wobei Differentialoperator D = d/dx.
;; Output: Lösungspotenzreihe F(x).
; Methode:
; Inspiziere für n=0,1,... den Koeffizienten von x^n in dieser Gleichung:
; sum(i=0..r, sum(j>=0, sum(k>=0, k-i+j=n, a[i,j]*k...(k-i+1)*b[k] ))) = 0.
; Sammle dabei die Koeffizienten aller b[k] (k=n+r, n+r-1, ...) und
; überprüfe, ob die entstehende Gleichung nach genau einem b[k] aufgelöst
; werden kann, der bis dato noch unbekannt ist.
(defun potlindgl (gleichung anfang)
  (let* ((r (1- (length gleichung)))
         (pcount 0)
         (pvec (make-array 10 :fill-pointer pcount :adjustable t)))
    (dolist (b anfang)
      (vector-push-extend b pvec)
      (setq pcount (1+ pcount))
    )
    (labels
          ((teilsumme (k L) ; liefert zu einer Liste L=(c[k] c[k-1] ...)
             (let ((sum 0)) ; die Teilsumme c[k]*b[k]+c[k-1]*b[k-1]+...
               (dolist (c L) (incf sum (* c (aref pvec k))) (decf k))
               sum
           ) )
           (inspiziere (n) ; Koeffizienten von x^n inspizieren
             (let ((bk-koeffs nil) ; Liste aller b[k]-Koeffs
                   (max-k (+ n r)) ; höchstes auftretendes k
                   (faktoren (reverse gleichung))) ; (a[r,*],...,a[0,*])
               (do ((k max-k (1- k)))
                   ((minusp k))
                 ; Koeffizient von b[k] bestimmen:
                 (let ((nochwas nil) (koeff 0))
                   (do ((i r (1- i))
                        (faktr faktoren (cdr faktr))) ; (a[i,*],...,a[0,*])
                       ((minusp i))
                     (let ((j (+ (- n k) i)) ; j mit j-i=n-k
                           (prod ; k...(k-i+1)
                             (do* ((x 1 (* x (- k y))) (y 0 (1+ y)))
                                  ((= y i) x)
                          )) )
                       (when (minusp j) (return)) ; j<0 -> für dieses k fertig
                       ; Koeffizient von x^j in a[i,*] = (car faktr) :
                       (if (listp (car faktr))
                         ; Polynom
                         (when (car faktr) ; Polynom zu Ende -> nichts zu tun
                           (setq nochwas t) ; sonst Koeffizienten verbrauchen
                           (incf koeff (* prod (pop (car faktr))))
                         )
                         ; Potenzreihe
                         (progn
                           (setq nochwas t)
                           (incf koeff (* prod (aref (funcall (car faktr) (1+ j)) j)))
                   ) ) ) )
                   (unless nochwas ; dieser und alle weiteren Koeffs von b[k]
                     (return)      ; = 0 -> Schleife beenden
                   )
                   (push koeff bk-koeffs)
               ) )
               (setq bk-koeffs (nreverse bk-koeffs))
               ; bk-koeffs = (Koeff von b[max-k],
               ;              Koeff von b[max-k - 1],
               ;              ...)
               (loop
                 (when (or (null bk-koeffs) (not (zerop (car bk-koeffs))))
                   (return)
                 )
                 (decf max-k) (pop bk-koeffs)
               )
               (unless (null bk-koeffs)
                 ; Die Gleichung bestimmt b[max-k] eindeutig
                 (if (< max-k pcount)
                   ; Gleichung nur überprüfen
                   (unless (zerop (teilsumme max-k bk-koeffs))
                     (error "Anfangskoeffizienten ~S erfüllen Differentialgleichung nicht."
                            anfang
                   ) )
                   ; Gleichung ausnutzen
                   (progn
                     (when (> max-k pcount)
                       (warn "Koeffizienten von x^~S bis x^~S werden willkürlich auf 0 gesetzt."
                             pcount (1- max-k)
                       )
                       (loop
                         (vector-push-extend 0 pvec)
                         (setq pcount (1+ pcount))
                         (when (= pcount max-k) (return))
                     ) )
                     ; pcount = max-k
                     (vector-push-extend
                       (- (/ (teilsumme (1- max-k) (cdr bk-koeffs))
                             (car bk-koeffs)
                       )  )
                       pvec
                     )
                     (setq pcount (1+ pcount))
          )) ) ) ) )
      (let ((m 0)) ; alle x^n mit n<m sind bereits inspiziert
        #'(lambda (n)
            (if (<= n pcount) pvec
              (loop
                (inspiziere m)
                (incf m)
                (when (<= n pcount) (return pvec))
          ) ) )
) ) ) )

;; liefert zu einer Potenzreihe G mit G(0)=0 die Substitution ( F -> F o G ).
;; Dies ist - genau genommen - eine Funktion, die zu einem n>=0 einen Array
;; der Länge n+1 liefert, dessen Element k (für 0<=k<=n) der Koeffizient
;; von x^n/n! in G(x)^k/k! ist.
(defun make-subst (pot)
  (let ((pvec (funcall pot 1)))
    (unless (zerop (aref pvec 0))
      (error "Nur Potenzreihen mit Absolutglied 0 können eingesetzt werden.")
    )
    (let ((tables (make-array 10 :fill-pointer 0 :adjustable t)))
      (vector-push-extend '#(1) tables)
      #'(lambda (n)
          (when (>= n (fill-pointer tables))
            (funcall pot (1+ n))
            (do ((m (fill-pointer tables) (1+ m)))
                ((> m n))
              ; Vektor der Koeffizienten von x^m/m! in G(x)^k/k! berechnen:
              (let ((Vm (make-array (1+ m)))
                    (k 0))
                (vector-push-extend Vm tables)
                (setf (aref Vm k) 0) ; Koeff. von x^m/m! in 1 ist 0
                (loop
                  (when (= k m) (return))
                  ; Um G(x)^(k+1)/(k+1)! zu bekommen, muß man G(x)^k/k! mit
                  ; G(x) multiplizieren und durch (k+1) dividieren. Der
                  ; Koeffizient von x^m/m! ist also
                  ; = sum(i=k..m, m!/i! * (Koeff. von x^i/i! in G(x)^k/k!)
                  ;                     * (Koeff. von x^(m-i) in G(x)) )
                  ;   / (k+1)
                  (let ((index m) (prod 1) (sum 0))
                    (loop
                      (incf sum (* prod (aref (aref tables index) k)
                                        (aref pvec (- m index))
                      )         )
                      (when (= index k) (return))
                      (setq prod (* prod index))
                      (decf index)
                    )
                    (incf k)
                    (setf (aref Vm k) (/ sum k))
          ) ) ) ) )
          (aref tables n)
) ) )   )

;; potsubst setzt eine Potenzreihe in eine andere Potenzreihe ein: F o G,
;; wobei G(0)=0 ist. G muß in Form einer Substitution vorliegen.
(defun potsubst (pot subst)
  (let* ((pvec1 (funcall pot 0))
         (pcount 0)
         (pvec (make-array 10 :fill-pointer pcount :adjustable t)))
    #'(lambda (n)
        (if (<= n pcount) pvec
          (progn
            (funcall pot n)
            (dotimes (i (- n pcount) pvec)
              (vector-push-extend
                ; Koeffizient von x^pcount ist =
                ; sum(i=0..pcount, i!/pcount!*pot[i]*(Koeff.von x^pcount/pcount! in G(x)^i/i!))
                (let ((V (funcall subst pcount))
                      (index pcount) (prod 1) (sum 0))
                  (loop
                    (incf sum (* prod (aref pvec1 index) (aref V index)))
                    (when (zerop index) (return))
                    (setq prod (/ prod index))
                    (setq index (1- index))
                  )
                  sum
                )
                pvec
              )
              (setq pcount (1+ pcount))
) )   ) ) ) )

;; eine Funktion zum Ausdrucken:
(defun potprint (pot index
                 &key (var "T") (shift 0) (printer nil) (width 79)
                      (stream (if printer *printer-output* *standard-output*))
                )
  ;; pot = Potenzreihe
  ;; index = Grenze für den Ausdruck der Koeffizienten (ausschließlich)
  ;; stream = Output-Stream
  ;; var = String oder Symbol, Bezeichnung der Unbestimmten
  ;; shift = Exponent des ersten Gliedes
  ;; printer = Flag, ob Ausgabe auf Drucker (Exponenten hochstellen)
  ;; width = Zeilenlänge
  (let* ((pvec (funcall pot index)) ; Vektor mit den Koeffizienten
         (start-flag t)  ; Flag für ersten Koeffizienten /= 0
         (hochstring1 (if printer (coerce '(#\Escape #\S #\Null) 'string) "^"))
             ; String vor dem Exponenten (für EPSON LQ-800)
         (hochstring2 (if printer (coerce '(#\Escape #\T) 'string) ""))
             ; String nach dem Exponenten (für EPSON LQ-800)
         (hochsummand (if printer 0 1))
             ; druckende Länge der beiden Strings
         (vstring (string var))  ; Name der Variablen als String
         (vlen (length vstring)) ; und seine Länge
         (plusstring " + ")      ; Vorzeichenstrings mit Längen
         (pluslen 3)
         (minusstring " - ")
         (minuslen 3)
         (charcount 0)   ; Zähler für verbrauchte Zeichen in der Zeile
        )
    (macrolet ((out (string stringwidth) ; einen String gegebener Druckbreite
                                         ; ausgeben
                 `(progn
                    (setq charcount (+ charcount ,stringwidth)) ; Zeichenzähler weiterzählen
                    ; wenn's nicht mehr hinpaßt, neue Zeile anfangen:
                    (when (> charcount width) (terpri stream) (setq charcount ,stringwidth))
                    (write-string ,string stream)
                  )
              ))
      (flet ((monom (kof expo) ; monom gibt ein Glied (mit Vorzeichen) aus
               (unless (zerop kof) ; Monom mit Koeffizient=0 wird übergangen
                 (let (str strlen) ; str = bisherige Zeichen fürs Monom,
                                   ; strlen = druckende Länge von str
                   (if (plusp kof)
                     (progn ; Koeffizient >0
                       (if start-flag ; evtl. Pluszeichen
                         (setq start-flag nil str "" strlen 0)
                         (setq str plusstring strlen pluslen)
                     ) )
                     (progn ; Koeffizient <0
                       (setq start-flag nil)
                       (setq str minusstring strlen minuslen) ; Minuszeichen
                       (setq kof (abs kof))
                   ) )
                   (unless (and (= kof 1) (not (zerop expo)))
                     (let ((kofstring (princ-to-string kof)))
                       (setq str (string-concat str kofstring)
                             strlen (+ strlen (length kofstring))
                   ) ) )
                   (unless (zerop expo)
                     (if (= expo 1)
                       (setq str (string-concat str vstring)
                             strlen (+ strlen vlen)
                       )
                       (let ((expostring (princ-to-string expo)))
                         (setq str (string-concat str vstring hochstring1 expostring hochstring2)
                               strlen (+ strlen vlen hochsummand (length expostring))
                   ) ) ) )
                   (out str strlen)
            )) ) )
        (terpri stream)
        (dotimes (i index)
          (monom (aref pvec i) (+ i shift))
        )
        (out " + ..." 6)
        (values)
) ) ) )

;; potsolve erzeugt aus einer Gleichung X = F(X,T) mit F Polynom,
;; F = O(X^2) + O(T) eine Potenzreihe X = P(T).
;; F ist in der Form ( { (i j Aij) } ) anzugeben, wenn
;; F = SUMME(i,j; Aij*X^i*T^j).
(defun potsolve (polyxt)
  (let* ((pcount 1)
         (pvec (make-array 10 :fill-pointer pcount :adjustable t))
         glieder
        )
    (setf (aref pvec 0) 0)
    (labels ((result (n)
               (if (<= n pcount) pvec
                 (let ((summanden
                         (mapcar #'(lambda (pot) (funcall pot n)) glieder)
                      ))
                   (dotimes (i (- n pcount) pvec)
                     (vector-push-extend
                       (do ((s summanden (rest s))
                            (sum 0 (+ sum (aref (first s) pcount)))
                           )
                           ((null s) sum)
                       )
                       pvec
                     )
                     (setq pcount (1+ pcount))
             ) ) ) )
             (make-glied (koeff &aux (i (first koeff)) (j (second koeff))
                                     (a (third koeff))
                         )
               (unless (or (>= i 2) (>= j 1))
                       (error "~: Polynom enthält falsches Glied ~"
                              'potsolve koeff
               )       )
               (pot*skal (potshift (potexpt_Z #'result i) j) a)
            ))
      (setq glieder (mapcar #'make-glied polyxt))
      #'result
) ) )

;; Liefert die Dirichlet-Faltung zweier Folgen (a1,a2,a3,...), die
;; als Potenzreihen a1+a2*X+a3*X^2+... dargestellt sind.
(defun pot*dirichlet (folge1 folge2)
  (let* ((pcount 0)
         (pvec (make-array 10 :fill-pointer pcount :adjustable t))
         (pvec1 (funcall folge1 0))
         (pvec2 (funcall folge2 0))
        )
    #'(lambda (n)
        (if (<= n pcount) pvec
          (progn
            (funcall folge1 n) (funcall folge2 n)
            (dotimes (i (- n pcount) pvec)
              (vector-push-extend
                (let ((m (1+ pcount)))
                  ; Koeffizient von X^(m-1) ist c[m] = sum(p*q=m, a[p]*b[q])
                  (do ((p m (1- p))
                       (sum 0))
                      ((zerop p) sum)
                    (multiple-value-bind (q r) (floor m p)
                      (when (zerop r)
                        (incf sum
                              (* (aref pvec1 (1- p)) (aref pvec2 (1- q)))
                ) ) ) ) )
                pvec
              )
              (setq pcount (1+ pcount))
) )   ) ) ) )

;; Liefert den Quotienten bezüglich Dirichlet-Faltung
;; zweier Folgen (a1,a2,a3,...), die
;; als Potenzreihen a1+a2*X+a3*X^2+... dargestellt sind.
(defun pot/dirichlet (folge1 folge2)
  ; (c1,c2,c3,...) := (a1,a2,a3,...) / (b1,b2,b3,...)
  ; bedeutet  (a1,a2,a3,...) = (b1,b2,b3,...) * (c1,c2,c3,...) , also
  ; c[m] = (a[m] - sum(p*q=m,p>1,q<m, b[p]*c[q]))/b[1] für m=1,2,3,...
  (let* ((pcount 0)
         (pvec (make-array 10 :fill-pointer pcount :adjustable t))
         (pvec1 (funcall folge1 0))
         (pvec2 (funcall folge2 1))
         (const (/ (aref pvec2 0)))
        )
    #'(lambda (n)
        (if (<= n pcount) pvec
          (progn
            (funcall folge1 n) (funcall folge2 n)
            (dotimes (i (- n pcount) pvec)
              (vector-push-extend
                (let ((m (1+ pcount)))
                  ; Koeffizient von X^(m-1) ist c[m] = (a[m]-...)/b[1]
                  (do ((q (1- m) (1- q))
                       (sum (aref pvec1 (1- m))))
                      ((zerop q) (* sum const))
                    (multiple-value-bind (p r) (floor m q)
                      (when (zerop r)
                        (decf sum
                              (* (aref pvec2 (1- p)) (aref pvec (1- q)))
                ) ) ) ) )
                pvec
              )
              (setq pcount (1+ pcount))
) )   ) ) ) )

;; Möbius-Funktion als Folge ist das Dirichlet-Faltungs-Inverse der
;; konstanten Folge (defpot 1).
(defvar möbius (pot/dirichlet \1 (defpot 1)))

;; potlogdiff bildet die logarithmische Ableitung einer Potenzreihe /=0.
(defun potlogderiv (pot)
  (pot/ (potderiv pot) pot)
)

;; potkanprodukt entwickelt eine Potenzreihe aus 1+Z[[X]] in ein
;; "kanonisches Produkt"  (1-X)^c1*(1-X^2)^c2*...  und liefert (c1,c2,...).
(defun potkanprodukt (pot)
  (let ((pvec (funcall pot 1)))
    (unless (= (aref pvec 0) 1)
      (error "Nur Potenzreihen mit Absolutglied 1 sind in ein kanonisches Produkt entwickelbar.")
  ) )
  ; Ziel: pot als prod(i>=1, (1-X^i)^c[i]) schreiben.
  ; 1. Schritt: logarithmische Ableitung, liefert
  ;             sum(i>=1, c[i] * d/dX log(1-X^i) )
  ;             = sum(i>=1, -c[i]*i*X^(i-1)/(1-X^i) )
  ;             = 1/X * sum(i>=1, -c[i]*i*sum(j>=1,X^(i*j)) )
  ;             = 1/X * sum(n>=1, sum(i*j=n,-c[i]*i) * X^n)
  ; 2. Schritt: Dies entspricht der Folge ( sum(i*j=n,-c[i]*i) )[n>=1]
  ;             was die Dirichlet-Faltung der Folgen  ( -1 )[j>=1]
  ;             und  ( c[i]*i )[i>=1]  ist. Division bzgl. Dirichlet-Faltung
  ;             liefert also die Folge  ( c[i]*i )[i>=1]  als Reihe
  ;             sum(i>=1, c[i]*i*X^(i-1) ).
  ; 3. Schritt: Integrieren liefert  sum(i>=1, c[i]*X^i) .
  ;             Dies entspricht der Folge (0,c1,c2,...).
  ; 4. Schritt: Division durch X, und man bekommt die Folge (c1,c2,...).
  (potshift (potint (pot/dirichlet (potlogderiv pot) (defpot -1))) -1)
)

;; binom liefert die Potenzreihe (1+T)^e.
(defun potbinom (e)
  (let* ((pcount 0)
         (pvec (make-array 10 :fill-pointer pcount :adjustable t))
         (koeff 1)
         (h e)
        )
    #'(lambda (n)
        (if (<= n pcount) pvec
          (dotimes (i (- n pcount) pvec)
            (vector-push-extend koeff pvec)
            (setq pcount (1+ pcount)
                  koeff (* koeff (/ h pcount))
                  h (1- h)
) )   ) ) ) )

(defun rn_n (r) ; asymptotische Entwicklung von (rn)!/n!^r, x=1/n
  (potexp
    (defpot (if (zerop pcount)
              0
              (* (/ (bernoulli (1+ pcount)) pcount (1+ pcount))
                 (- (/ (expt r pcount)) r)
            ) )
) ) )

(defun n_nr (r) ; asymptotische Entwicklung von n!/(n/r)!^r, x=1/n
  (potexp
    (defpot (if (zerop pcount)
              0
              (* (/ (bernoulli (1+ pcount)) pcount (1+ pcount))
                 (- 1 (expt r (1+ pcount)))
            ) )
) ) )

; wandelt eine Reihe sum(n>=0, c[n+r]/n!*x^n) in eine Reihe c[0]+c[1]*x+...,
; wobei die ersten r Koeffizienten angegeben werden müssen.
(defun asympt-Reihe (pot &optional (anfang '(1)))
  (let* ((r (length anfang))
         (pvec1 (funcall pot 0))
         (pcount 0)
         (pvec (make-array 10 :fill-pointer pcount :adjustable t)))
    (dolist (c anfang)
      (vector-push-extend c pvec)
      (setq pcount (1+ pcount))
    )
    #'(lambda (n)
        (if (<= n pcount) pvec
          (progn
            (funcall pot (- n r))
            (dotimes (i (- n pcount) pvec)
              (vector-push-extend
                (let ((m (- pcount r))) ; Koeff. von x^m ist c[m+r]/m!
                  (* (aref pvec1 m) (! m)) ; mal m! -> liefert c[m+r]=c[pcount]
                )
                pvec
              )
              (setq pcount (1+ pcount))
) )   ) ) ) )


(setq hn2 ; asymptotische Entwicklung von Hn(2)
  (asympt-Reihe ; F(x)=sum(n>=0,c[n+1]/n!) löst die Differentialgleichung
    (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
      (let ((e^x exp) (e^2x (pot* exp exp)))
        (list (pot+ (pot*skal e^2x 2) (pot*skal e^x 2) (pot*skal |1| 3))
              (pot+ (pot*skal e^2x 2) (pot*skal e^x -2) (pot*skal |1| -8))
              (pot+ (pot*skal e^x -4) (pot*skal |1| 4))
      ) )
      (list 1/4)
) ) )

(require 'stirling) ; Funktionen stirling1-table, stirling2-table
; Verwende  #'stirling1-table  für die Substitution  x=log(1+t)
; und       #'stirling2-table  für die Substitution  t=e^x-1 .

; Stellt fest, ob eine asymptotische Reihe von einem Zp-Maß herkommt:
; Asymptotische Reihe sum(k>=0,c[k]/n^k) = exp( sum(k>=1,b[k]/k/n^k) )
; kommt genau dann von einem Maß, wenn sum(k>=1,b[k]/k!*x^k) in Zp[[e^x-1]]
; liegt. Wir liefern die Potenzreihe, bei der hierin t=e^x-1 substituiert
; ist. Es kann leicht gesehen werden, für welche p sie in Zp[[t]] liegt.
(defun maßtest (asym)
  (potsubst
    (pot*hadamard
      (potlog asym) ; sum(k>=1,b[k]/k*x^k)
      (potshift exp 1) ; sum(k>=1,1/(k-1)!*x^k)
    ) ; Hadamard-Produkt = sum(k>=1,b[k]/k!*x^k)
    #'stirling1-table ; darin x=log(1+t) einsetzen
) )

#|

;; Modulfunktionen, insbesondere Eisensteinreihen:

; (sigma s n) = sum(d|n,d^s)
(defun sigma (s n)
  (let ((n^1/2 (isqrt n)))
    (do ((d1 1 (1+ d1))
         (sum 0))
        ((> d1 n^1/2) sum)
      (when (zerop (mod n d1))
        (incf sum (expt d1 s)) ; Teiler <=n^1/2 berücksichtigen
        (let ((d2 (floor n d1)))
          (when (> d2 d1)
            (incf sum (expt d2 s)) ; Teiler >n^1/2 berücksichtigen
) ) ) ) ) )

; (eisenstein k) liefert die Potenzreihe Ek(q) = 1/(2*zeta(2k)) * Gk(q),
; Gk(q) = 2*zeta(2k) + 2*(-1)^k*(2pi)^2k/(2k-1)!*sum(n=1..inf,sigma(2k-1,n)*q^n)
; mit zeta(2k) = 2^(2k-1)*pi^2k*abs(bernoulli(2k))/(2k)!

; 8k*(-1)^k abs(bernoulli(2k))
(defun eisenstein (k &aux (2k-1 (- (* 2 k) 1)))
  (let* ((pcount 1)
         (pvec (make-array 1 :fill-pointer pcount :adjustable t))
         (koeff (/ (* 4 k (expt -1 k)) (abs (bernoulli (* 2 k)))))
        )
    (setf (aref pvec 0) 1)
    #'(lambda (n)
        (if (<= n pcount) pvec
          (dotimes (i (- n pcount) pvec)
            (vector-push-extend (* koeff (sigma 2k-1 pcount)) pvec)
            (setq pcount (1+ pcount))
) )   ) ) )

; Delta-Reihe:  = g2^3-27*g3^2 = (2pi)^12/(2^6*3^3)*(E2^3-E3^2)
;                = (2pi)^12*Delta
(setq Delta
  (pot*skal (pot- (potexpt_Z (eisenstein 2) 3) (potexpt_Z (eisenstein 3) 2))
            1/1728
) )

; DeltaX(q) = Prod(n=1..inf,1-q^n)^X
(setq Delta24 (potshift Delta -1))
(setq Delta3 (potexpt_Q Delta24 1/8))
(setq Delta1 (potexpt_Q Delta24 1/24))

|#