Wurzelapproximation

17.01.2025, 07:16

Dopap

Wurzelapproximation

Neulich fand ich in meinem mathe/physik Karton den edlen ledergebundenen simplen solaren
Klapptaschenrehner im Scheckkartenformat aus den 80' Jahren wieder.
Lediglich die Wurzeltaste funktionierte nicht mehr.
Eine memorierbare Funktion sollte jetzt Abhilfe schaffen, wofür aber nur Polynome samt Hornerschema infrage kamen.
Ein Approximatios-Polynom z.B eines 4.grades ist dafür aber ziemlich ungeeignet.
Klar, im entscheidenden Intervall von 1 bis 99 kann man auch mit einem ad hoc Schätzwert $\begin{eqnarray*} x_0\approx \sqrt {y^*} \end{eqnarray*}$ samt anschließender
Verbesserung nach Newton arbeiten... Die Funktion sollte aber ohne Schätzwert auskommen und Rechenstabenauigkeit haben!

Idee : Bruch von Polynomen mit ganzzahligen Koeffizienten
$\begin{eqnarray*} \sqrt{(y)}\approx W(y)=P(y)/Q(y) \end{eqnarray*}$ zum Beispiel

$\begin{eqnarray*} x:= W(y)= \frac {(A + By+Cy^2)}{(1+by+cy^2)} \end{eqnarray*}$ mit den 5 Koeffizienten { A B C b c}.
Mit der Wahl von 5 Interpolationspunkten $\begin{eqnarray*} (x_i,y_i) = (1,1) ,(2,4), (5,25), (7,49) ,(9,81) \end{eqnarray*}$ Danach ist das GS
$\begin{eqnarray*} x_i\cdot (1+by_i+cy_i^2)=A+By_i+Cy_i^2 \quad i=1..5 \end{eqnarray*}$ zu lösen.

Die so numerisch berechneten Koeffizienten sind unhandlich und auch nicht ganzzahlig.
Mit etwas Runden und "Feilen" konnte ich dann W(y) zu

$\begin{eqnarray*} x_0 := W(y^*)=\frac{(630 + 768y^*+24y^{*2})}{(1231+208y^*+y^{*2})}\approx \sqrt {y^*} \end{eqnarray*}$

bestimmen, ohne den Gesamtfehler $\begin{eqnarray*} [\epsilon \epsilon]= \sum\limits_{k=1}^{99} {\left(\sqrt k -W(k)\right)} ^2 = 0.1240 \end{eqnarray*}$ wesentlich erhöht zu haben.
Die Graphen von $\begin{eqnarray*} W(x) \text{ und } \sqrt x \end{eqnarray*}$ in einem Koordinatensystem:

man beachte das wohlwollende Verhalten von W(x) ab x>99 sowie die Asymptoden y=24 und x=-6.1.
Dahingegen kann ein Polynom mit seinen "Wellen" einpacken.
Ein evtl. gewünschter anschließender Newton ist schnell getippt:

$\begin{eqnarray*} \frac {x_0^2+y^*}{\underbrace {2x_0}_{\text{Ableitung von y=x^2 (!)}}}=:x_1 \end{eqnarray*}$ mit ca. 5 bis 6 gültigen Ziffern

bleibt noch die Frage :

wie könnten die optimalen 5 Stützpunkte mit einem Minimum von $\begin{eqnarray*} [\epsilon \epsilon] \end{eqnarray*}$ aussehen?

17.01.2025, 09:00

HAL 9000

Auf diesen Beitrag antworten »

Statt 1,2,5,7,9 mit dieser zu groß scheinenden Lücke zwischen 2 und 5 scheint mir eine Wahl der $\begin{eqnarray*} x \end{eqnarray*}$ -Punkte sinnvoller, wo die Differenzenfolge monoton wächst, also z.B. 1,2,4,6,9, zumal ja Kurvenverlauf der Quadrat- bzw. Wurzelfunktion im hinteren Verlauf weniger gekrümmt ist als vorn, und damit einer geringeren Stützstellendichte bedarf - heuristisch betrachtet.

---------------------------------------------------------------------------

Ich hab mal folgende Alternative ausprobiert:

1) Skalieren der Floating-Point-Zahl $\begin{eqnarray*} x \end{eqnarray*}$ auf das Intervall $\begin{eqnarray*} [0.5,2) \end{eqnarray*}$ :

Das ist relativ einfach, da im Floating-Point-Format $\begin{eqnarray*} x = m\cdot 2^e \end{eqnarray*}$ der Exponent $\begin{eqnarray*} e \end{eqnarray*}$ direkt vorliegt. Wir setzen nun $\begin{eqnarray*} f = \left\lfloor \frac{e+1}{2}\right\rfloor \end{eqnarray*}$ und betrachten statt $\begin{eqnarray*} x \end{eqnarray*}$ den Radikand $\begin{eqnarray*} x' = x\cdot 2^{-2f} = m\cdot 2^{e-2f} \end{eqnarray*}$ , dann haben wir das gewünschte $\begin{eqnarray*} x'\in [0.5,2) \end{eqnarray*}$ . Anzumerken ist, dass bei dieser Operation die Mantisse unverändert bleibt, rechentechnisch beschränkt sich das auf ein paar superschnelle einfache Bitschiebereien.

2) Durchführen von nur zwei Newton-Schritten ausgehend von Startwert 1, d.h. $\begin{eqnarray*} z_1 = \frac{1+x'}{2} \end{eqnarray*}$ und $\begin{eqnarray*} z_2 = \frac{1}{2}\left(z_1+\frac{x'}{z_1}\right) \end{eqnarray*}$

3) Wieder hochskalieren $\begin{eqnarray*} y = z_2\cdot 2^f \end{eqnarray*}$ (wiederum nur eine Exponentenanpassung mit unveränderter Mantisse von $\begin{eqnarray*} z_2 \end{eqnarray*}$ ).

code:
1: 2: 3: 4: 5: 6: 7:	double halsqrt(double x) { int f = (ilogb(x) + 1) >> 1; x = ldexp(x, -(f << 1)); double y = ldexp(1+x, -1); return ldexp(y + x/y, f-1); }

Hab mal nachgerechnet: Ergibt für deine Testwerte $\begin{eqnarray*} \sum_{k=1}^{99} \left(\sqrt{k}-y(k)\right)^2 \approx 0.000575815 \end{eqnarray*}$ , und was Rechenoperationen betrifft sehr sparsam: In Floating-Point zwei Additionen und eine Division - der Rest sind nur Additionen,Subtraktionen und Shiftoperationen von Integer-Zahlen (und zwar der Floatingpoint-Exponenten). Und das Verfahren ist nicht auf $\begin{eqnarray*} x\in [1,99] \end{eqnarray*}$ beschränkt, sondern für den gesamten positiven Wertebereich der Floating-Point-Zahlen geeignet.

P.S.: Kann natürlich sein, dass du in deinem Museumsstück keinen Zugriff auf die Binärdarstellung der Zahlen hast, womöglich verwenden die auch gar nicht das IEEE754-Floatingpoint-Format.

20.01.2025, 02:17

Dopap

Auf diesen Beitrag antworten »

nein, das Museumsstück hat nur Grundrechenarten + Wurzel + 1/x und wurde rein
geschäftlich genutzt.
Bitschiebereien kann nur der HP 50g und dort funktioniert die Wurzel noch.

habe jetzt deine Punktliste ausprobiert, was exakt gerechnet zu
$\begin{eqnarray*} W(x)=\frac{432 + 596x+22x^2}{876+ 173x +x^2 } \end{eqnarray*}$ und zu $\begin{eqnarray*} [\epsilon \epsilon]= 0.0313 \end{eqnarray*}$ führt.
Im Vergleich zu meiner unüberlegten Liste mit $\begin{eqnarray*} [\epsilon \epsilon]= 0.0154 \end{eqnarray*}$ verwirrt

Aber egal, meistens wird eh' noch der schnelle Newton angehängt und optisch nicht zu unterscheiden

20.01.2025, 08:02

HAL 9000

Auf diesen Beitrag antworten »

Sich widersprechende Angaben
Was gilt denn nun für deine Orignialfunktion: Das

Zitat:

Original von Dopap
ohne den Gesamtfehler $\begin{eqnarray*} [\epsilon \epsilon]= \sum\limits_{k=1}^{99} {\left(\sqrt k -W(k)\right)} ^2 = 0.1240 \end{eqnarray*}$ wesentlich erhöht zu haben.

oder das

Zitat:

Original von Dopap
Im Vergleich zu meiner unüberlegten Liste mit $\begin{eqnarray*} [\epsilon \epsilon]= 0.0154 \end{eqnarray*}$

----------------------------------------------------------------------------------

Da Vertrauen in so einem Fall dann nicht mehr angebracht ist, habe ich selbst nachgerechnet:

Für $\begin{eqnarray*} W(y^*)=\frac{(630 + 768y^*+24y^{*2})}{(1231+208y^*+y^{*2})} \end{eqnarray*}$ bekomme ich $\begin{eqnarray*} \sum\limits_{k=1}^{99} {\left(\sqrt k -W(k)\right)} ^2 = 0.3531 \end{eqnarray*}$ heraus.

Das hingegen

Zitat:

Original von Dopap
$\begin{eqnarray*} W(x)=\frac{432 + 596x+22x^2}{876+ 173x +x^2 } \end{eqnarray*}$ und zu $\begin{eqnarray*} [\epsilon \epsilon]= 0.0313 \end{eqnarray*}$ führt.

kann ich bestätigen.

EDIT: Für $\begin{eqnarray*} W(y^*)=\frac{(630 + 7{\color{red}86}y^*+24y^{*2})}{(1231+208y^*+y^{*2})}{y^*} \end{eqnarray*}$ stimmt das $\begin{eqnarray*} \sum\limits_{k=1}^{99} {\left(\sqrt k -W(k)\right)} ^2 = 0.0154 \end{eqnarray*}$ . Ein doppeltes Forum Kloppe

für den Zahlendreher sowie die initial falsche Quadratsummenberechnung.

[attach]58071[/attach]

Um die Güte der Approximation im Bild bewerten zu können, plottet man besser nicht die Wurzel- und die Approximationsfunktion, sondern die Differenz der beiden: deine Approximation (blau) und meine "Bitschieber"-basierte (orange).

20.01.2025, 17:48

Dopap

Auf diesen Beitrag antworten »

tut mir leid. Der Zahlendreher liegt wohl an der deutschen Eigentümlichkeit der Endziffern.
Weiterhin konnte ich feststellen, dass mein Programm in der ersten post
am Ende noch eine Wurzel gezogen hatte also $\begin{eqnarray*} \sqrt{[\epsilon \epsilon]} \end{eqnarray*}$ unglücklich

und zudem klemmen nach 8 Jahren manche Tasten am Laptop ...

-----------------------------

Informative Grafik !

20.01.2025, 22:49

HAL 9000

Auf diesen Beitrag antworten »

Zitat:

Original von Dopap
Der Zahlendreher liegt wohl an der deutschen Eigentümlichkeit der Endziffern.

Originelle Ausrede - sowas muss einem erstmal einfallen. Chapeau! Teufel

Neue Frage »

Antworten »

Wurzelapproximation

Verwandte Themen