<html><body><div>This procedure is using REAL type which is now 64-bit (double precision).</div><div>But constants c3,c4,c5 have only 8 fractional digits (old single precision REAL type?)</div><div>What are these constants and how to make them having more digits to get double precision?</div><div><br></div><div>Zdenek <br></div><div><br></div><aside>---------- Původní e-mail ----------<br>Od: Jörg Straube <joerg.straube@iaeth.ch><br>Komu: ETH Oberon and related systems <oberon@lists.inf.ethz.ch><br>Datum: 22. 5. 2024 7:49:14<br>Předmět: Re: [Oberon] An alternative for the present inaccurate Math.ln(x)</aside><br><blockquote data-email="joerg.straube@iaeth.ch"><p class="-wm-MsoNormal" style="font-size:16px;margin:0cm;font-family:Aptos,sans-serif;caret-color:rgb(33,33,33);color:rgb(33,33,33)"><span lang="EN-US" style="font-size:11pt;line-height:1.2">Hans</span></p><p class="-wm-MsoNormal" style="font-size:16px;margin:0cm;font-family:Aptos,sans-serif;caret-color:rgb(33,33,33);color:rgb(33,33,33)"><span lang="EN-US" style="font-size:11pt;line-height:1.2"> </span></p><p class="-wm-MsoNormal" style="font-size:16px;margin:0cm;font-family:Aptos,sans-serif;caret-color:rgb(33,33,33);color:rgb(33,33,33)"><span lang="EN-US" style="font-size:11pt;line-height:1.2">Thx for pointing out.</span></p><p class="-wm-MsoNormal" style="font-size:16px;margin:0cm;font-family:Aptos,sans-serif;caret-color:rgb(33,33,33);color:rgb(33,33,33)"><span lang="EN-US" style="font-size:11pt;line-height:1.2">Find a very minor improvement below:</span></p><p class="-wm-MsoNormal" style="font-size:16px;margin:0cm;font-family:Aptos,sans-serif;caret-color:rgb(33,33,33);color:rgb(33,33,33)"><span lang="EN-US" style="font-size:11pt;line-height:1.2"> </span></p><p class="-wm-MsoNormal" style="font-size:16px;margin:0cm;font-family:Aptos,sans-serif;caret-color:rgb(33,33,33);color:rgb(33,33,33)"><span lang="EN-US">PROCEDURE ln*(x: REAL): REAL;</span></p><p class="-wm-MsoNormal" style="font-size:16px;margin:0cm 0cm 0cm 35.4pt;font-family:Aptos,sans-serif;caret-color:rgb(33,33,33);color:rgb(33,33,33)"><span lang="EN-US">(* ln(x) = 2*arctanh( (x-1)/(x+1) )<br> around 0, arctan() is almost linear with slope 1</span></p><p class="-wm-MsoNormal" style="font-size:16px;margin:0cm;font-family:Aptos,sans-serif;caret-color:rgb(33,33,33);color:rgb(33,33,33);text-indent:35.4pt"><span lang="EN-US">*)</span></p><p class="-wm-MsoNormal" style="font-size:16px;margin:0cm;font-family:Aptos,sans-serif;caret-color:rgb(33,33,33);color:rgb(33,33,33);text-indent:35.4pt"><span lang="EN-US">CONST</span></p><p class="-wm-MsoNormal" style="font-size:16px;margin:0cm 0cm 0cm 35.4pt;font-family:Aptos,sans-serif;caret-color:rgb(33,33,33);color:rgb(33,33,33);text-indent:35.4pt"><span lang="EN-US">c1 = 1.4142135; (* sqrt(2) *)</span></p><p class="-wm-MsoNormal" style="font-size:16px;margin:0cm 0cm 0cm 35.4pt;font-family:Aptos,sans-serif;caret-color:rgb(33,33,33);color:rgb(33,33,33);text-indent:35.4pt"><span lang="EN-US">c2 = 0.6931472; (* ln(2) *)</span></p><p class="-wm-MsoNormal" style="font-size:16px;margin:0cm 0cm 0cm 35.4pt;font-family:Aptos,sans-serif;caret-color:rgb(33,33,33);color:rgb(33,33,33);text-indent:35.4pt"><span lang="EN-US">c3 = 0.89554059;</span></p><p class="-wm-MsoNormal" style="font-size:16px;margin:0cm 0cm 0cm 35.4pt;font-family:Aptos,sans-serif;caret-color:rgb(33,33,33);color:rgb(33,33,33);text-indent:35.4pt"><span lang="EN-US">c4 = 1.82984424;</span></p><p class="-wm-MsoNormal" style="font-size:16px;margin:0cm 0cm 0cm 35.4pt;font-family:Aptos,sans-serif;caret-color:rgb(33,33,33);color:rgb(33,33,33);text-indent:35.4pt"><span lang="EN-US">c5 = 1.65677798;</span></p><p class="-wm-MsoNormal" style="font-size:16px;margin:0cm 0cm 0cm 35.4pt;font-family:Aptos,sans-serif;caret-color:rgb(33,33,33);color:rgb(33,33,33)"><span lang="EN-US">VAR e: INTEGER;</span></p><p class="-wm-MsoNormal" style="font-size:16px;margin:0cm 0cm 0cm 35.4pt;font-family:Aptos,sans-serif;caret-color:rgb(33,33,33);color:rgb(33,33,33)"><span lang="EN-US">BEGIN</span></p><p class="-wm-MsoNormal" style="font-size:16px;margin:0cm 0cm 0cm 35.4pt;font-family:Aptos,sans-serif;caret-color:rgb(33,33,33);color:rgb(33,33,33);text-indent:35.4pt">ASSERT(x > 0.0); UNPK(x, e); (* x in 1 .. 2 *)</p><p class="-wm-MsoNormal" style="font-size:16px;margin:0cm 0cm 0cm 35.4pt;font-family:Aptos,sans-serif;caret-color:rgb(33,33,33);color:rgb(33,33,33);text-indent:35.4pt"><span lang="EN-US">IF x > c1 THEN x := x*0.5; INC(e) END; (* x in 0.7 .. 1.4)</span></p><p class="-wm-MsoNormal" style="font-size:16px;margin:0cm 0cm 0cm 35.4pt;font-family:Aptos,sans-serif;caret-color:rgb(33,33,33);color:rgb(33,33,33);text-indent:35.4pt">x := (x - 1.0)/(x + 1.0); (* x in -0.17 .. 0.17 *)</p><p class="-wm-MsoNormal" style="font-size:16px;margin:0cm 0cm 0cm 35.4pt;font-family:Aptos,sans-serif;caret-color:rgb(33,33,33);color:rgb(33,33,33);text-indent:35.4pt">x := FLT(e)*c2 + x*(c3 + c4/(c5 - x*x))</p><p class="-wm-MsoNormal" style="font-size:16px;margin:0cm 0cm 0cm 35.4pt;font-family:Aptos,sans-serif;caret-color:rgb(33,33,33);color:rgb(33,33,33)"><span lang="EN-US">RETURN x END ln;</span></p><p class="-wm-MsoNormal" style="font-size:16px;margin:0cm;font-family:Aptos,sans-serif;caret-color:rgb(33,33,33);color:rgb(33,33,33)"><span lang="EN-US" style="font-size:11pt;line-height:1.2"> </span></p><p class="-wm-MsoNormal" style="font-size:16px;margin:0cm;font-family:Aptos,sans-serif;caret-color:rgb(33,33,33);color:rgb(33,33,33)"><span lang="EN-US" style="font-size:11pt;line-height:1.2">Jörg</span></p><div><div dir="auto" style="caret-color:rgb(0,0,0);color:rgb(0,0,0);letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;text-decoration:none;word-wrap:break-word;line-break:after-white-space"><div><br></div></div><br class="-wm-Apple-interchange-newline">
</div>
<div><br><blockquote><div>Am 19.05.2024 um 21:45 schrieb Hans Klaver <hklaver@dds.nl>:</div><br class="-wm-Apple-interchange-newline"><div><div style="word-wrap:break-word;line-break:after-white-space"><div dir="auto" style="word-wrap:break-word;line-break:after-white-space"><div dir="auto" style="word-wrap:break-word;line-break:after-white-space">Hi all,<div><br></div><div>Did anyone notice that the Project Oberon 2013 Math.ln(x) function gives very sloppy results, only the first decimal or first few decimals being correct? </div><div><br></div><div>While trying to produce a fractal on letter H (a nice programming task when your name is Hans ;-) and getting the result in figure 1, I first blamed my home grown graphics routines, but it soon became clear that Math.ln was the culprit. I use that function in my implementation of Math.power() which I need to calculate a contraction factor for the fractal.</div><div><br></div><div><span id="-wm-cid:DB9BCBE8-D6FF-41CD-8E60-B4F72D299A8D"><Schermafbeelding 2024-05-17 om 14.58.06.png></span></div><div><br></div><div>To see the inaccuracy just calculate ln(5.0) with a calculator or on <a href="http://www.wolframalpha.com/">www.wolframalpha.com</a> and compare it with the result of Math.ln(5.0): the correct result to 6 decimal places should be 1.609438E+00 but Math.ln gives 1.621313E+00.</div><div><br></div><div>First I hoped that a typo was the cause: in the PO book (Applications, p. 83) constant q1 = 6.176106560E1, whereas in the source q1 = 6.76106560E1. But changing that value gave little improvement.</div><div><br></div><div>With the help of two older Oberon implementations I was able to figure out a correct version of Math.ln(x) which also uses the UNPK(x, e) standard procedure, and doesn't need to import SYSTEM. The new correct version turns out to be even simpler than the present one.</div><div><br></div><div> PROCEDURE ln* (x: REAL): REAL; (* hk 18 May 2024 *)<br> CONST<br> c1 = 0.70710678; (* 1/ sqrt(2) *)<br> c2 = 0.69314718; (* ln(2) *)<br> c3 = 0.89554059;<br> c4 = 1.82984424;<br> c5 = 1.65677798;<br> VAR e: INTEGER; y: REAL;<br> BEGIN ASSERT(x > 0.0); UNPK(x, e);<br> y := x * 0.5; INC(e);<br> IF y < c1 THEN y := 2.0*y; DEC(e) END;<br> y := (y-1.0)/(y+1.0);<br> y := c2*FLT(e) + y*(c3 + c4/(c5 - y*y))<br> RETURN y<br> END ln;</div><div><br></div><div><div> PROCEDURE ln*(x: REAL): REAL; (* INCORRECT version *)</div><div> CONST c1 = 0.70710680; (* 1/sqrt(2) *)</div><div> c2 = 0.69314720; (* ln(2) *)</div><div> p0 = -9.01746917E1;</div><div> p1 = 9.34639006E1;</div><div> p2 = -1.83278704E1;</div><div> q0 = -4.50873458E1;</div><div> q1 = 6.76106560E1;</div><div> q2 = -2.07334879E1;</div><div> VAR e: INTEGER; xx, y: REAL;</div><div> BEGIN ASSERT(x > 0.0); UNPK(x, e);</div><div> IF x < c1 THEN x := x*2.0; e := e-1 END ;</div><div> x := (x-1.0)/(x+1.0);</div><div> xx := x;</div><div> y := c2*FLT(e) + x*((p2*xx + p1)*xx + p0) / (((xx + q2)*xx + q1)*xx + q0);</div><div> RETURN y</div><div> END ln;</div><br></div><div>The two older versions (one from 1995 by Michael Griebling and another by an unknown author from the 1999 ETH Oberon System 3 release 2.3) give the exact same correct results, but could not use UNPK and therefore are much more and somewhat more convoluted. </div><div><br></div><div>The procedures UNPK and PACK prove very useful to implement these mathematical functions. In Griebling's version the functionality of UNPK is given by two functions fraction() and exponent() and that of PACK by a function scale(); they need SYSTEM.VAL() and SET. Also in the 1999 ETH version SYSTEM.VAL() and SET is used.</div><div><br></div><div>You can find them on my GitHub pages, along with a TestMath module. (<a href="https://github.com/hansklav/Oberon-07-Math.ln/tree/main">https://github.com/hansklav/Oberon-07-Math.ln/tree/main</a>)</div><div>There you can also find the (imho) missing arctan() function. (<a href="https://github.com/hansklav/Oberon-07">https://github.com/hansklav/Oberon-07</a>)</div><div><br></div><div>As far as I could test them the other functions in module Math (sqrt, exp, sin, cos, arctan) are accurate enough.</div><div><br></div><div>For those who are puzzled by procedures UNPK and PACK it may be informative to know that in most other languages their functionality is given by functions named frexp() and ldexp() or scalbn(). These routines are used to extract the binary fraction and exponent of 2 from a floating-point value, and to synthesize a float from these components. They provide a way to get at the internal representation of a float without doing direct bit manipulation. Wirth only used them in module Math.</div><div><br></div><div>With the new Math.ln() function the H fractal looks much prettier ;-)</div><div><br></div><div>Regards,</div><div><br></div><div>Hans Klaver</div><div><br></div><div><span id="-wm-cid:CDFA0705-72E7-4A6B-822E-69C25C0574BA"><Schermafbeelding 2024-05-18 om 14.46.07.png></span></div><div><br></div><div><br></div></div></div></div>--<br>Oberon@lists.inf.ethz.ch mailing list for ETH Oberon and related systems<br>https://lists.inf.ethz.ch/mailman/listinfo/oberon<br></div></blockquote></div><br>--
<br>Oberon@lists.inf.ethz.ch mailing list for ETH Oberon and related systems
<br>https://lists.inf.ethz.ch/mailman/listinfo/oberon
<br></blockquote></body></html>