Logarithmen berechnen

Das Verfahren von R. P. Kelisky und T. J. Rivlin basiert auf rekursiven Definitionen für den Zähler und den Nenner eines Bruchs, dessen Wert eine Näherung für den natürlichen Logarithmus ist.

Die ersten beiden Zähler sind 2(z-1) und 4(z-1)(z+1).

(defmethod logarithm-numerator (z steps) (= steps 1)
  (* 2 (- z 1)))

(defmethod logarithm-numerator (z steps) (= steps 2)
  (* 4 (- z 1) (+ z 1)))

Jeder weitere Zähler wird rekursiv aus seinen beiden Vorgängern, hier als p und pp bezeichnet, berechnet.

(defmethod logarithm-numerator (z p pp steps) (and (> steps 2) (integer? steps))
  (- (* 2 (+ z 1) p)
     (* (- 1 z) (- 1 z) pp)
     (if
       (even? steps)
       0
       (* 2 (expt (- 1 z) steps) (/ -2 steps (- steps 2))))))

(defmethod logarithm-numerator (z steps) (and (> steps 2) (integer? steps))
  (let
    ((a 0)
     (p (logarithm-numerator z 2))
     (pp (logarithm-numerator z 1)))
    (dotimes (i steps p)
      (setq a (logarithm-numerator z p pp (+ i 3)))
      (setq pp p)
      (setq p a))))

Eine etwas einfachere Rekursionsvorschrift dient zur Berechnung der Nenner. Die ersten beiden sind z+1 und 1+6z+z².

(defmethod logarithm-denominator (z steps) (= steps 1)
  (+ 1 z))

(defmethod logarithm-denominator (z steps) (= steps 2)
  (+ 1 (* 6 z) (* z z)))

(defmethod logarithm-denominator (z p pp steps) (and (> steps 2) (integer? steps))
  (- (* 2 (+ 1 z) p)
     (* (- 1 z) (- 1 z) pp)))

(defmethod logarithm-denominator (z steps) (and (> steps 2) (integer? steps))
  (let
    ((a 0)
     (p (logarithm-denominator z 2))
     (pp (logarithm-denominator z 1)))
    (dotimes (i steps p)
      (setq a (logarithm-denominator z p pp (+ i 3)))
      (setq pp p)
      (setq p a))))

Das oben beschriebene Näherungsverfahren konvergiert nur gut für Argumente zwischen 1 und 2.

(defproc small-logarithm (z steps)
  (/ (logarithm-numerator z steps)
     (logarithm-denominator z steps)))

Größere Argumente werden anhand einer additiven Zerlegung auf kleinere zurückgeführt.

(defproc big-logarithm (z h steps)
  (if
    (> z 2)
    (big-logarithm (/ z 2) (+ h 1) (+ steps 1))
    (+ (small-logarithm z steps)
       (* h (small-logarithm 2 steps)))))

Argumente kleiner als 1 werden durch eine Transformation auf Argumente größer als 1 abgebildet. Um eine Genauigkeit von etwa d Dezimalstellen zu erreichen, muss das Näherungsverfahren 5/4 d Schritte durchlaufen.

(defproc logarithm (z digits)
  (if
    (< z 1)
    (- (logarithm (/ z) digits))
    (let
      ((steps (floor (* 5/4 digits))))
      (approximate
        (if
          (> z 2)
          (big-logarithm z 0 steps)
          (small-logarithm z steps))
        (expt 1/10 digits)))))

Die Hilfsfunktion even? liefert den Wert t, wenn es sich bei ihrem Argument um eine gerade Zahl handelt.

(defproc even? (n) (integer? (/ n 2)))

Die Hilfsfunktion expt berechnet die Potenz einer Zahl (siehe Potenzen berechnen).

Für den natürlichen Logarithmus von 2 ergibt sich eine Näherung auf 100 Stellen von

0.693147180559945309417232121458176568075500134360
2552541206800094933936219696947156058633269964186875.


Quelle

R. P. Kelisky, T. J. Rivlin
"A Rational Approximation to the Logarithm"
IBM Watson Research Center
Yorktown Heights, New York
1967


logarithm.sheet.xml