An alternative for the present inaccurate Math.ln(x)
Hi all, 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? 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. To see the inaccuracy just calculate ln(5.0) with a calculator or on www.wolframalpha.com 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. 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. 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. PROCEDURE ln* (x: REAL): REAL; (* hk 18 May 2024 *) CONST c1 = 0.70710678; (* 1/ sqrt(2) *) c2 = 0.69314718; (* ln(2) *) c3 = 0.89554059; c4 = 1.82984424; c5 = 1.65677798; VAR e: INTEGER; y: REAL; BEGIN ASSERT(x > 0.0); UNPK(x, e); y := x * 0.5; INC(e); IF y < c1 THEN y := 2.0*y; DEC(e) END; y := (y-1.0)/(y+1.0); y := c2*FLT(e) + y*(c3 + c4/(c5 - y*y)) RETURN y END ln; PROCEDURE ln*(x: REAL): REAL; (* INCORRECT version *) CONST c1 = 0.70710680; (* 1/sqrt(2) *) c2 = 0.69314720; (* ln(2) *) p0 = -9.01746917E1; p1 = 9.34639006E1; p2 = -1.83278704E1; q0 = -4.50873458E1; q1 = 6.76106560E1; q2 = -2.07334879E1; VAR e: INTEGER; xx, y: REAL; BEGIN ASSERT(x > 0.0); UNPK(x, e); IF x < c1 THEN x := x*2.0; e := e-1 END ; x := (x-1.0)/(x+1.0); xx := x; y := c2*FLT(e) + x*((p2*xx + p1)*xx + p0) / (((xx + q2)*xx + q1)*xx + q0); RETURN y END ln; 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. 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. You can find them on my GitHub pages, along with a TestMath module. (https://github.com/hansklav/Oberon-07-Math.ln/tree/main) There you can also find the (imho) missing arctan() function. (https://github.com/hansklav/Oberon-07) As far as I could test them the other functions in module Math (sqrt, exp, sin, cos, arctan) are accurate enough. 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. With the new Math.ln() function the H fractal looks much prettier ;-) Regards, Hans Klaver
Hans Thx for pointing out. Find a very minor improvement below: PROCEDURE ln*(x: REAL): REAL; (* ln(x) = 2*arctanh( (x-1)/(x+1) ) around 0, arctan() is almost linear with slope 1 *) CONST c1 = 1.4142135; (* sqrt(2) *) c2 = 0.6931472; (* ln(2) *) c3 = 0.89554059; c4 = 1.82984424; c5 = 1.65677798; VAR e: INTEGER; BEGIN ASSERT(x > 0.0); UNPK(x, e); (* x in 1 .. 2 *) IF x > c1 THEN x := x*0.5; INC(e) END; (* x in 0.7 .. 1.4) x := (x - 1.0)/(x + 1.0); (* x in -0.17 .. 0.17 *) x := FLT(e)*c2 + x*(c3 + c4/(c5 - x*x)) RETURN x END ln; Jörg
Am 19.05.2024 um 21:45 schrieb Hans Klaver <hklaver@dds.nl>:
Hi all,
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?
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.
<Schermafbeelding 2024-05-17 om 14.58.06.png>
To see the inaccuracy just calculate ln(5.0) with a calculator or on www.wolframalpha.com <http://www.wolframalpha.com/> 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.
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.
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.
PROCEDURE ln* (x: REAL): REAL; (* hk 18 May 2024 *) CONST c1 = 0.70710678; (* 1/ sqrt(2) *) c2 = 0.69314718; (* ln(2) *) c3 = 0.89554059; c4 = 1.82984424; c5 = 1.65677798; VAR e: INTEGER; y: REAL; BEGIN ASSERT(x > 0.0); UNPK(x, e); y := x * 0.5; INC(e); IF y < c1 THEN y := 2.0*y; DEC(e) END; y := (y-1.0)/(y+1.0); y := c2*FLT(e) + y*(c3 + c4/(c5 - y*y)) RETURN y END ln;
PROCEDURE ln*(x: REAL): REAL; (* INCORRECT version *) CONST c1 = 0.70710680; (* 1/sqrt(2) *) c2 = 0.69314720; (* ln(2) *) p0 = -9.01746917E1; p1 = 9.34639006E1; p2 = -1.83278704E1; q0 = -4.50873458E1; q1 = 6.76106560E1; q2 = -2.07334879E1; VAR e: INTEGER; xx, y: REAL; BEGIN ASSERT(x > 0.0); UNPK(x, e); IF x < c1 THEN x := x*2.0; e := e-1 END ; x := (x-1.0)/(x+1.0); xx := x; y := c2*FLT(e) + x*((p2*xx + p1)*xx + p0) / (((xx + q2)*xx + q1)*xx + q0); RETURN y END ln;
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.
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.
You can find them on my GitHub pages, along with a TestMath module. (https://github.com/hansklav/Oberon-07-Math.ln/tree/main) There you can also find the (imho) missing arctan() function. (https://github.com/hansklav/Oberon-07)
As far as I could test them the other functions in module Math (sqrt, exp, sin, cos, arctan) are accurate enough.
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.
With the new Math.ln() function the H fractal looks much prettier ;-)
Regards,
Hans Klaver
<Schermafbeelding 2024-05-18 om 14.46.07.png>
-- Oberon@lists.inf.ethz.ch mailing list for ETH Oberon and related systems https://lists.inf.ethz.ch/mailman/listinfo/oberon
Hi Jörg, You wrote:
Thx for pointing out. Find a very minor improvement below:
PROCEDURE ln*(x: REAL): REAL; (* ln(x) = 2*arctanh( (x-1)/(x+1) ) around 0, arctan() is almost linear with slope 1 *) CONST c1 = 1.4142135; (* sqrt(2) *) c2 = 0.6931472; (* ln(2) *) c3 = 0.89554059; c4 = 1.82984424; c5 = 1.65677798; VAR e: INTEGER; BEGIN ASSERT(x > 0.0); UNPK(x, e); (* x in 1 .. 2 *) IF x > c1 THEN x := x*0.5; INC(e) END; (* x in 0.7 .. 1.4) x := (x - 1.0)/(x + 1.0); (* x in -0.17 .. 0.17 *) x := FLT(e)*c2 + x*(c3 + c4/(c5 - x*x)) RETURN x END ln;
Beautiful! Thanks. You managed to find an even more elegant form for ln(). But please note that there is an unclosed comment in the above code. You mention the arctan() function, which reminds me of one episode in which I corresponded with Niklaus Wirth, about typos in his online book 'Programming in Oberon (a Tutorial)', later called 'Programming, A Tutorial'. I asked him if he had left out the function procedure arctan deliberately in his definition of Math.Mod (p. 46). He replied: "I did it intentionally, because I find that arctan is rarely used anyway". I was shocked. To me it was as if someone left out the icosahedron from the list of Platonic solids "because it is rarely used anyway"! The list of standard mathematical functions in programming languages has a long tradition (going back to FORTRAN, and Algol-60), and all these early languages were thought out by mathematicians or physicists, so I thought there was a good *mathemetical* reason to have arctan in the list. Wirth was a pragmatic programmer, and for him obviously it was no problem at all to leave arctan out. But he was not extremely thorough, because arctan is still mentioned in 'Programming, A Tutorial', on p. 15 ;-) I'm not a mathematician, but in my humble opinion Math.Mod without arctan is lacking. And maybe there are good computer science reasons to include more mathematical functions in it (as proposed in the Oakwood Guidelines). Regards, Hans
Hans The calculation of ln() is based on arctanh() (not arctan(), sorry for the typo) Unfortunately, arctan() is not as “easy” as ln(). The range reduction gets more complex. br Jörg
Am 23.05.2024 um 13:43 schrieb Hans Klaver <hklaver@dds.nl>:
Hi Jörg,
You wrote:
Thx for pointing out. Find a very minor improvement below:
PROCEDURE ln*(x: REAL): REAL; (* ln(x) = 2*arctanh( (x-1)/(x+1) ) around 0, arctan() is almost linear with slope 1 *) CONST c1 = 1.4142135; (* sqrt(2) *) c2 = 0.6931472; (* ln(2) *) c3 = 0.89554059; c4 = 1.82984424; c5 = 1.65677798; VAR e: INTEGER; BEGIN ASSERT(x > 0.0); UNPK(x, e); (* x in 1 .. 2 *) IF x > c1 THEN x := x*0.5; INC(e) END; (* x in 0.7 .. 1.4) x := (x - 1.0)/(x + 1.0); (* x in -0.17 .. 0.17 *) x := FLT(e)*c2 + x*(c3 + c4/(c5 - x*x)) RETURN x END ln;
Beautiful!
Thanks. You managed to find an even more elegant form for ln().
But please note that there is an unclosed comment in the above code.
You mention the arctan() function, which reminds me of one episode in which I corresponded with Niklaus Wirth, about typos in his online book 'Programming in Oberon (a Tutorial)', later called 'Programming, A Tutorial'. I asked him if he had left out the function procedure arctan deliberately in his definition of Math.Mod (p. 46). He replied: "I did it intentionally, because I find that arctan is rarely used anyway".
I was shocked. To me it was as if someone left out the icosahedron from the list of Platonic solids "because it is rarely used anyway"!
The list of standard mathematical functions in programming languages has a long tradition (going back to FORTRAN, and Algol-60), and all these early languages were thought out by mathematicians or physicists, so I thought there was a good *mathemetical* reason to have arctan in the list.
Wirth was a pragmatic programmer, and for him obviously it was no problem at all to leave arctan out. But he was not extremely thorough, because arctan is still mentioned in 'Programming, A Tutorial', on p. 15 ;-)
I'm not a mathematician, but in my humble opinion Math.Mod without arctan is lacking. And maybe there are good computer science reasons to include more mathematical functions in it (as proposed in the Oakwood Guidelines).
Regards,
Hans
-- Oberon@lists.inf.ethz.ch mailing list for ETH Oberon and related systems https://lists.inf.ethz.ch/mailman/listinfo/oberon
Hello: I am re-reading the Digital Circuit Design by NW. I got interested in Lola-2 (not the original Lola). Looking at the RISC5Top.Lola I am seeing a reference to PROM, declared as PROM = MODULE (IN clk: BIT; IN adr: [9] BIT; OUT data: WORD) ^; This module is missing from the Lola page https://people.inf.ethz.ch/wirth/Lola/index.html. It this omission intentional? Does anyone have the source of PROM.Lola? A related question: how does the Lola-2 compiler import the Lola-2 modules? I am not seeing the export / import marks like in Oberon. I am not seeing a project file either, as it would be the case for ISE or Vivado projects. Does the compiler find the separate module files on its own? Is anyone using Lola-2 for sizeable projects of the complexity similar to RISC5 System on Chip? If so, what was the practical experience? Does Lola-2 help, or the opposite? Thank you, Wojtek
On Mon, Jun 3, 2024 at 9:40 AM Skulski, Wojciech <skulski@pas.rochester.edu> wrote:
Hello:
I am re-reading the Digital Circuit Design by NW. I got interested in Lola-2 (not the original Lola). Looking at the RISC5Top.Lola I am seeing a reference to PROM, declared as
PROM = MODULE (IN clk: BIT; IN adr: [9] BIT; OUT data: WORD) ^;
This module is missing from the Lola page https://people.inf.ethz.ch/wirth/Lola/index.html. It this omission intentional? Does anyone have the source of PROM.Lola?
A related question: how does the Lola-2 compiler import the Lola-2 modules? I am not seeing the export / import marks like in Oberon. I am not seeing a project file either, as it would be the case for ISE or Vivado projects. Does the compiler find the separate module files on its own?
Is anyone using Lola-2 for sizeable projects of the complexity similar to RISC5 System on Chip? If so, what was the practical experience? Does Lola-2 help, or the opposite?
Thank you, Wojtek
-- Oberon@lists.inf.ethz.ch mailing list for ETH Oberon and related systems https://lists.inf.ethz.ch/mailman/listinfo/oberon
The Lola-2 compiler translates a single Lola-2 module to a corresponding Verilog file. Call it (in any sequence) for each of the Lola-2 modules that you want to translate. You can then create an ISE or Vivado Project to reference the resulting individual Verilog modules. The Verilog source file PROM.v file that you also need to include, is located in the Project Oberon -> SourcesVerilog folder. -- Regards, Chris Burrows CFB Software https://www.astrobe.com/RISC5
Chris:
The Lola-2 compiler translates a single Lola-2 module to a corresponding Verilog file. Call it (in any sequence) for each of the Lola-2 modules that you want to translate. You can then create an ISE or Vivado Project to reference the resulting individual Verilog modules. The Verilog source file PROM.v file that you also need to include, is located in the Project Oberon -> SourcesVerilog folder.
Thank you! I thought that the missing PROM was a mistake. What you are saying sounds different from the other NW work where he always stressed the consistency checks among files and types. I like the way NW handled the registers and tri-state buffers. The VHDL and Verilog way is too verbose. Lola-2 is providing the essence. I miss IF THEN ELSE. The ternary operator notation looks kind of awful. Lola-3, if it ever materializes, could provide IF THEN ELSE, the generate statement, and the COMPONENT keyword in addition to the MODULE, meaning the same. I wonder if substantial projects can be handled with Lola-2. Wojtek
This procedure is using REAL type which is now 64-bit (double precision). But constants c3,c4,c5 have only 8 fractional digits (old single precision REAL type?) What are these constants and how to make them having more digits to get double precision? Zdenek ---------- Původní e-mail ---------- Od: Jörg Straube <joerg.straube@iaeth.ch> Komu: ETH Oberon and related systems <oberon@lists.inf.ethz.ch> Datum: 22. 5. 2024 7:49:14 Předmět: Re: [Oberon] An alternative for the present inaccurate Math.ln(x) " Hans Thx for pointing out. Find a very minor improvement below: PROCEDURE ln*(x: REAL): REAL; (* ln(x) = 2*arctanh( (x-1)/(x+1) ) around 0, arctan() is almost linear with slope 1 *) CONST c1 = 1.4142135; (* sqrt(2) *) c2 = 0.6931472; (* ln(2) *) c3 = 0.89554059; c4 = 1.82984424; c5 = 1.65677798; VAR e: INTEGER; BEGIN ASSERT(x > 0.0); UNPK(x, e); (* x in 1 .. 2 *) IF x > c1 THEN x := x*0.5; INC(e) END; (* x in 0.7 .. 1.4) x := (x - 1.0)/(x + 1.0); (* x in -0.17 .. 0.17 *) x := FLT(e)*c2 + x*(c3 + c4/(c5 - x*x)) RETURN x END ln; Jörg " Am 19.05.2024 um 21:45 schrieb Hans Klaver <hklaver@dds.nl>: Hi all, 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? 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. <Schermafbeelding 2024-05-17 om 14.58.06.png> To see the inaccuracy just calculate ln(5.0) with a calculator or on www. wolframalpha.com(http://www.wolframalpha.com/) 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. 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. 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. PROCEDURE ln* (x: REAL): REAL; (* hk 18 May 2024 *) CONST c1 = 0.70710678; (* 1/ sqrt(2) *) c2 = 0.69314718; (* ln(2) *) c3 = 0.89554059; c4 = 1.82984424; c5 = 1.65677798; VAR e: INTEGER; y: REAL; BEGIN ASSERT(x > 0.0); UNPK(x, e); y := x * 0.5; INC(e); IF y < c1 THEN y := 2.0*y; DEC(e) END; y := (y-1.0)/(y+1.0); y := c2*FLT(e) + y*(c3 + c4/(c5 - y*y)) RETURN y END ln; PROCEDURE ln*(x: REAL): REAL; (* INCORRECT version *) CONST c1 = 0.70710680; (* 1/sqrt(2) *) c2 = 0.69314720; (* ln(2) *) p0 = -9.01746917E1; p1 = 9.34639006E1; p2 = -1.83278704E1; q0 = -4.50873458E1; q1 = 6.76106560E1; q2 = -2.07334879E1; VAR e: INTEGER; xx, y: REAL; BEGIN ASSERT(x > 0.0); UNPK(x, e); IF x < c1 THEN x := x*2.0; e := e-1 END ; x := (x-1.0)/(x+1.0); xx := x; y := c2*FLT(e) + x*((p2*xx + p1)*xx + p0) / (((xx + q2)*xx + q1)*xx + q 0); RETURN y END ln; 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. 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. You can find them on my GitHub pages, along with a TestMath module. ( https://github.com/hansklav/Oberon-07-Math.ln/tree/main (https://github.com/hansklav/Oberon-07-Math.ln/tree/main)) There you can also find the (imho) missing arctan() function. (https:// github.com/hansklav/Oberon-07(https://github.com/hansklav/Oberon-07)) As far as I could test them the other functions in module Math (sqrt, exp, sin, cos, arctan) are accurate enough. 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. With the new Math.ln() function the H fractal looks much prettier ;-) Regards, Hans Klaver <Schermafbeelding 2024-05-18 om 14.46.07.png> -- Oberon@lists.inf.ethz.ch mailing list for ETH Oberon and related systems https://lists.inf.ethz.ch/mailman/listinfo/oberon " -- Oberon@lists.inf.ethz.ch mailing list for ETH Oberon and related systems https://lists.inf.ethz.ch/mailman/listinfo/oberon "
On Sun, May 26, 2024 at 4:58 PM <oberon-eth@email.cz> wrote:
This procedure is using REAL type which is now 64-bit (double precision). But constants c3,c4,c5 have only 8 fractional digits (old single precision REAL type?) What are these constants and how to make them having more digits to get double precision?
Zdenek
The procedure quoted is intended for use in Project Oberon where the REAL / LONGREAL type is 32-bit.
A similar procedure was implemented using 64-bit LONGREALs in ETH Native Oberon System3. The source code can be downloaded from here: https://github.com/rsdoiel/NativeOberon/blob/main/MathL.Mod.txt The use of hexadecimal numbers for the constants circumvents rounding issues that might occur when different compilers convert them from text to binary. -- Regards, Chris Burrows CFB Software http://www.astrobe.com/RISC5
participants (5)
-
Chris Burrows -
Hans Klaver -
Jörg Straube -
oberon-eth＠email.cz -
Skulski, Wojciech