real operations imprecise?

[H. S. Teoh]

On Tue, May 05, 2020 at 01:44:18PM +0000, WebFreak001 via Digitalmars-d-learn wrote:
> I was dumping the full PI value on my machine with the highest
> precision it could get and got:
> 
> $ rdmd --eval='printf("%.70llf\n", PI)'
> 3.1415926535897932385128089594061862044327426701784133911132812500000000

Whoa, hold your horses right there!  What does `pragma(msg, real.dig);`
output on your machine?

On my machine, it's 18, i.e., `real` is capable of holding only ~18
digits of meaningful information. 70 digits is WAY beyond anything
that's actually represented in a `real`. And if you check online for the
actual digits of pi, you'll see that after about the 19th digit of your
output above, the rest of the digits are just pure garbage.  What you
see is just meaningless output from the typesetting algorithm from
information that isn't actually there.


> now this all looks good, but when I tried to print PI_2 I got
> 
> $ rdmd --eval='printf("%.70llf\n", PI_2)'
> 1.5707963267948965579989817342720925807952880859375000000000000000000000
> 
> note how many more 0's there are.
> 
> When manually writing down the identifier as it is defined for PI but
> with the exponent reduced by one it looks good again:
> $ rdmd --eval='printf("%.70llf\n",
> cast(real)0x1.921fb54442d18469898cc51701b84p+0L)'
> 1.5707963267948966192564044797030931022163713350892066955566406250000000

Again, what's the value of `real.dig` on your machine?  Don't be
deceived by the number of trailing zeroes; most of the digits that come
before it are complete garbage long before it dwindled to zero.  After
about the 17th digit, the two printouts above have already diverged.
Looks to me like it's a difference of just 1 ulp or so in the actual
representation, assuming you're on x86 where `real` has about 18 digits
of precision.


> I would expect that the manifest constant PI_2 being defined as PI/2
> would simply modify the exponent and not the actual value in such a
> drastic manner.
> 
> While I don't need it myself right now, is there some compiler switch to
> make real operations like this division more precise to not lose all these
> bits of information? (LDC 1.20.1)

At the most, the above difference is only *one* bit of information. Most
of your trailing digits are meaningless garbage because they don't
actually exist in the representation of `real`.  Unless you have a
256-bit representation of `real`, you can't expect to get that many
digits out of it!!


> Sure, 1^-50 precision might not be everyones typical use-case but I
> think there is potential improvement either in the compiler or
> std.math to be done here :p

Um, no, if you need 1^-50 precision maybe you should be looking at
arbitrary-precision float libraries, like MPFR (but be prepared for a
big performance hit once you move away from hardware floats). The
hardware `real` type simply does not have that many bits to store that
many digits. You're asking for more digits (*way* more) than are
actually stored in the type, so your test results are invalid.  D's
floating-point types have a .dig property for a reason.  Use it! ;-)


T

-- 
The problem with the world is that everybody else is stupid.