OpenEuphoria

Re: Phix, numbers, printf, log10

• Posted by petelomax in February
• 1018 views

Agh, yep, my bad: it goes wrong on a properly built 64-bit - I was using a built-by-32-bit thing,
which happened to go wrong in completely different way, that just happened to undermine that test...
(I've also just realised that means I have always shipped a duff 64-bit-for-Linux executable...)
Anyway, it simply turns out to be exceeding the accuracy limits of log10:

?log10(1000 - 1e-16)-3 -- -2.168404345e-19 (right sign-wise) 
?log10(1000 + 1e-16)-3 -- -2.168404345e-19 (wrong sign-wise) 
?log10(1000 - 1e-15)-3 -- -4.33680869e-19 (right sign-wise) 
?log10(1000 + 1e-15)-3 -- 2.168404345e-19 (right sign-wise) 
?atom_to_float80(log10(1000 - 1e-16))&-1 -- {255,255,255,255,255,255,255,191,0,64'@',-1} 
?atom_to_float80(log10(1000 + 1e-16))&-1 -- {255,255,255,255,255,255,255,191,0,64'@',-1} 
?atom_to_float80(log10(1000 - 1e-15))&-1 -- {254,255,255,255,255,255,255,191,0,64'@',-1} 
?atom_to_float80(log10(1000 + 1e-15))&-1 -- {1,0,0,0,0,0,0,192,0,64'@',-1} 

So log10 of 1000 +/- 1e-16 returns the exact same value, not much printf() etc can do about that.

You quite sure about that? ?apply({99800000000000000,9.98 * power(10, 16),9.98e16},integer) gets me {1,0,0}.
Take a look at ptok.e/completeFloat(), line 1603. There's a glaring error, try ?{1e1000,1e2000}, which I just fixed by:

--22/2/24: 
--      if exponent>308 then 
        if exponent>308 and machine_bits()=32 then 
            -- rare case: avoid power() overflow 

There may be a case for "exponent>4932" under 64 bit, and/or similar on the -ve bounds, and rather than capping it
at 1000, use an is_inf() [which almost certainly did not exist when I wrote that code] check inside the loop.

As you can see, otherwise I modified that 2, 9, and 14 years ago, ie not very often, and the last time looks suspiciously similar to the case you just uncovered. The matching entry in the readme was:

28/09/2022: Bugfix: 1e18 was not generating an integer on 64bit. Rather than 
            (further) mess with %opPow in VM\pPower.e, I simply added a *10 
            loop in ptok.e/completeFloat() for (decimal) exponents 0..20. 
            Aka: integer i=1e18 <== "type error (storing float in integer)". 

In fact I suspect that test should [now] be elsif exponent=18 [and machine_bits()=64] then, or probably
better, line 321 in VM\pPower.e should be cmp rcx,18 and that branch/loop in completeFloat() removed.
(I've already made the latter two of those changes, and tested them, on the development version.)

It should be possible to defer the TokN += fraction/dec statement until after we also have exponent,
resurrect the "end if / exponent -= ndp / if exponent then" and do something a little better/smarter.
At that mid-point, everything that can be an integer should still be one. Using the 32-bit 1e308 limit for clarity,
though the same applies to the 64-bit 1e4932 limit, one thing we must not do is temporarily create 4.99e309 on our
way to creating 4.99e307, iyswim, but otherwise we can re-assemble our 4, 99, and 307 any way we like, such as
(4*100+99)*power(10,305) aka (TokN*power(10,ndp)+fraction)*[/]power([-](exponent-ndp)), to keep things integer.
Please keep a list of anything that goes wrong, however briefly, so I can later add them all to "p -test".
[No longer entirely sure, but I think the /= power(10,-exponent) was put there for accuracy reasons, which may
only matter for numbers very close to said limits, or very nearly zero.]

Anyhow, feel free to have a play with completeFloat(), performance is not even slightly an issue, my only concern
would be that 32-bit does not get crippled, and your only concern should be having a backup of ptok.e and p[w].exe
in case things go south, you also have "p -test" to keep you out of trouble. FYI, I run my editor etc off a manual
copy of pw.exe, so I don't have to keep on shutting things down to rebuild it. It is of course entirely your choice
to play with that in the compiler itself, or take a copy to play with outside, or even incorporate into Hypatia, but
rest assured the first of those sounds far scarier than it actually is.

I am now thinking that %e is in fact the best way to get the magnitude:

function magnitude(atom x)  --, dflt=x|0) 
--  if is_inf(x) or is_nan(x) then return x|0|dlft end if   -- ?? 
    string s = sprintf("%.20e",x) 
    integer e = find('e',s), 
            n = to_number(s[e+1..$]) 
--  atom chk = power(10,n)              -- ?? 
--  if x>=chk*10 then n += 1            -- ?? 
--  elsif x<chk then n -= 1 end if      -- ?? 
    return n 
end function 

You are of course quite right about several edge cases being cause for concern, though, and certainly enough to prevent me from even thinking of adding something like that as a new builtin, at least just yet.

Don't worry, this stuff is all very welcome!
Should you feel there is something missing from the docs, try writing something I can paste right in,
and tell me precisely where you think that would best go.

Oops, I'm on it...

PS for anyone concerned about that "duff 64-bit-for-Linux executable":

printf(1,"%.20f\n",PI) 
-- 3.14159265358979311600 -- as downloaded    (1.2246e-16 out) 
-- 3.14159265358979323851 -- after "/.p -c p" (0.0005e-16 out) 

I have improved the build script so that will no longer occur in 1.0.5 and later