Dear all,
I am hitting a wall with the attached code in CL and C versions. For those interested, this is code lifted from Knuth's pages on random number generation (https://www-cs-faculty.stanford.edu/~knuth/programs.html#rng).
The problem is that there is a discrepancy in the "problem" loop between the CL and the C and I cannot figure out why, although the culprit may be that CL does not have "proper" 64 bits ints.
I am running on LW 8.x on an Intel Mac. I have not tried on other CL implementations/platforms. I compile the C code with clang (Apple clang version 15.0.0 (clang-1500.0.40.1)
If you run the code, you will see that the discrepancy appears at 'j = 14' in the loop.
What gives?
Thanks
Marco
On 9 Nov 2023, at 21:21, Marco Antoniotti marco.antoniotti@unimib.it wrote:
<problem-loop.lisp>
From the start, it looks like the ulp is more precise in C:
sbcl: 2.220446049250313d-16 clang: 2.2204460492503130808e-16
(using %.20g instead of %.20f)
Or perhaps it’s only the display procedure that truncates in lisp?
Hi
Thanks Pascal.
For LW on Intel (Mac) the ULP seems the same. With SBCL you should actually be able to peek at the actual bits making up the double float. Can you do something similar with LM?
Just curious: has anybody tried this on a M*/Arm Mac? Or, with LW, on your smartphone? :)
Cheers
MA
On Fri, Nov 10, 2023 at 8:29 AM Pascal Bourguignon (as pjb at informatimago dot com) lisp-hug@lispworks.com wrote:
On 9 Nov 2023, at 21:21, Marco Antoniotti marco.antoniotti@unimib.it wrote:
<problem-loop.lisp>
From the start, it looks like the ulp is more precise in C:
sbcl: 2.220446049250313d-16 clang: 2.2204460492503130808e-16
(using %.20g instead of %.20f)
Or perhaps it’s only the display procedure that truncates in lisp?
-- __Pascal J. Bourguignon__
Actually, ULP is not the same on LispWorks, which is why it diverges at 14. You can't see it when printing 20 decimal places.
The reason is that LispWorks calculates (expt 2.0d0 -52) using natural logarithms, so loses a few bits of precision. If you change it to use (/ 1d0 (ash 1 30) (ash 1 22)) or (float (expt 2 -52) 1d0) then LispWorks also appears to diverge at 21 like on other Lisps.
Looking at the binary representation of the floats shows that this divergence at 21 is just in the printing. The C library prints some extra digits that are not needed.
Thank you Martin.
Meanwhile, I fell in a rabbit hole (your fault guys!) and I produced this thing. I am sure it is code that has been written a godzillion times before. I will use it for the next phase of debugging.
Please poke holes into it.
All the best
Marco
On Fri, Nov 10, 2023 at 2:45 PM Martin Simmons martin@lispworks.com wrote:
Actually, ULP is not the same on LispWorks, which is why it diverges at 14. You can't see it when printing 20 decimal places.
The reason is that LispWorks calculates (expt 2.0d0 -52) using natural logarithms, so loses a few bits of precision. If you change it to use (/ 1d0 (ash 1 30) (ash 1 22)) or (float (expt 2 -52) 1d0) then LispWorks also appears to diverge at 21 like on other Lisps.
Looking at the binary representation of the floats shows that this divergence at 21 is just in the printing. The C library prints some extra digits that are not needed.
-- Martin Simmons LispWorks Ltd http://www.lispworks.com/
On Fri, 10 Nov 2023 09:55:38 +0100, Marco Antoniotti said:
Hi
Thanks Pascal.
For LW on Intel (Mac) the ULP seems the same. With SBCL you should actually be able to peek at the actual bits making up the double float. Can you do something similar with LM?
Just curious: has anybody tried this on a M*/Arm Mac? Or, with LW, on
your
smartphone? :)
Cheers
MA
On Fri, Nov 10, 2023 at 8:29 AM Pascal Bourguignon (as pjb at
informatimago
dot com) lisp-hug@lispworks.com wrote:
On 9 Nov 2023, at 21:21, Marco Antoniotti marco.antoniotti@unimib.it wrote:
<problem-loop.lisp>
From the start, it looks like the ulp is more precise in C:
sbcl: 2.220446049250313d-16 clang: 2.2204460492503130808e-16
(using %.20g instead of %.20f)
Or perhaps it’s only the display procedure that truncates in lisp?
-- __Pascal J. Bourguignon__
Lisp Hug - the mailing list for LispWorks users lisp-hug@lispworks.com http://www.lispworks.com/support/lisp-hug.html