OpenEuphoria

1. Phix, numbers, printf, log10

• Posted by RobertS in February
• 1178 views

Phix 64-bit, but it similarily happens with 32-bit:

printf(1, "%1.2e", 2.25)
2.26e+0

printf(1, "%1.2e", 2.48)
2.49e+0

printf(1, "%1.2e", 2.428)
2.42e+0 truncating instead of rounding

It happens with all kinds of numbers, and any number of decimal digits — the last digit tends to be hit or miss.

Also, though a bit less annoying, %f has a "5" rounding issue, this too happens with all kinds of numbers and any number of decimal digits (also happens with Euphoria, as far as I can remember):

printf(1, "%0.1f", 2.75)
2.8 rounding up

printf(1, "%0.1f", 2.85)
2.8 rounding down

Something else: what is the best way to find the magnitude (exponent in scientific notation) of x? I've been using floor(log10(x)), but, probably due to how log10 is calculated, it's not entirely reliable:

? floor(log10(1000))
3 correct

? floor(log10(1000 - 1e-16))
2 correct

? floor(log10(1000 + 1e-16))
2 badly wrong

printf("%e", x) seems to get the exponent more reliably, but using sprintf("%e", x) and then reading the number after the "e" from the string is rather clumsy, and rounding can become an issue, too.

I think I've managed to work around these issues (thought this before, though), but things might be easier if sprint were able to return more than 10 significant digits (I'd be content with 15).

Thanks for any help or clarification where I may be confused,
Robert

3. Re: Phix, numbers, printf, log10

• Posted by petelomax in February
• 1109 views

Ah, the problem here is that 2.85 cannot be held exactly, and the nearest IEEE-754 representation is ~2.849999999999999999 which is correctly rounded down.
You might be able to do something with rationals, otherwise that's a fundamental and unsolveable hardware issue. [At least now I can start on the first three.]

PS I stumbled on what might be a fairly nifty magnitude() function, but it needs the range of values it must cope with, from you, and probably a bit of polishing, before I can post it.

Update: I've now got a fix for the first part ("%e" was not updating f the same way "%f" does as it printed digits):

            digit = floor(f) 
--21/2/24: 
            f = (f-digit)*10 

and

--21/2/24 
--                  f = (f-digit)*10 
                    digit = floor(f) 
                    f = (f-digit)*10 

Many thanks. I'd recommend running "p -test" after editing that in and running "p -cp", to make sure nothing got broken.

5. Re: Phix, numbers, printf, log10

• Posted by RobertS in February
• 1063 views

Sorry, maybe it's my fault and I did it wrong, but now I get:

printf(1, "%e", 12)
1.220000e+1

It works for numbers < 10, gets the last digit right, but when there is more than 1 digit before the decimal point, the digit before the decimal point gets repeated:

printf(1, "%e", 123.456)
1.2334563e+2

But I may have misunderstood how to edit the file. I now have

        else 
            digit = floor(f) 
--21/2/24:  
            f = (f-digit)*10  
        end if 
        result &= digit+'0' 
        expadj = 1 
        if precision>0 then 

and

--30/10/21 
                    fadj = fwk 
                else 
--21/2/24 
                    digit = floor(f) 
                    f = (f-digit)*10 
                    fadj = 0 
                end if 
--12/7/16: 
                if digit=10 then exit end if 

Anyway, when this is fixed, it will help a lot, thank you!

Robert

7. Re: Phix, numbers, printf, log10

• Posted by petelomax in February
• 1028 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

9. Re: Phix, numbers, printf, log10

• Posted by petelomax in February
• 975 views

OK, I've aligned the "fwk" side with the "(f-digit)*10" side and got rid of fadj, and pushed all this to https://raw.githubusercontent.com/petelomax/Phix/master/builtins/VM/pprntfN.e which you may find much easier:

@ -218,20 +218,26 @@ end for 
                fwk += epwr 
                digit += 1 
            end while 
--24/2/24: 
            f -= fwk 
        else 
            digit = floor(f) 
--21/2/24: 
            f = (f-digit)*10 
        end if 
        result &= digit+'0' 
        expadj = 1 
        if precision>0 then 
            result &= '.' 
            dotdone = 1 
--24/2/24: 
--          atom fadj = 0 
--1/11/15 
--          for i=1 to precision do 
            for i=1 to precision-(charflag='g') do 
                if ewk>0 then 
--24/2/24: 
--                  f -= fwk 
                    ewk -= 1 
                    epwr = power(10,ewk) 
                    fwk = 0--epwr 
@ -240,19 +246,25 @@ end for 
                        fwk += epwr 
                        digit += 1 
                    end while 
--24/2/24: 
--30/10/21 
--                  fadj = fwk 
                    f -= fwk 
                else 
--21/2/24 
--                  f = (f-digit)*10 
                    digit = floor(f) 
                    f = (f-digit)*10 
--24/2/24: 
--                  fadj = 0 
                end if 
--12/7/16: 
                if digit=10 then exit end if 
                result &= digit+'0' 
                expadj += 1 
            end for 
--24/2/24: 
--          f -= fadj 
        else 
            if ewk>0 then 
                f -= fwk 

Clearly "p -test" is missing a whole bunch of "%e"-related stuff.

11. Re: Phix, numbers, printf, log10

• Posted by RobertS in February
• 932 views
• Last edited in March

Thank you for the %e fix, I'll try it out tomorrow.

I am doing this, rounding the value to the specified number of decimal digits before calling sprintf():

-- atom r is the number to be printed 
-- integer fmt is the number of decimal digits 
-- string rf is a temporary variable 
-- string rr is the formatted result (yes, bad habit of mine, to use minimalistic variable names) 
 
r = round(r, power(10, fmt)) 
rf = "%1." & sprint(fmt) & "f" 
rr = trim(sprintf(rf, r)) 

I'm not sure if it's 100%, but I think it is a significant improvement.

I've also written a little function to work around the fact that power() uses logarithms, and therefore doesn't always give integer results for power(a, b) even if a and b are integer — seems to work fine (updated March 9):

function power_int(atom a, atom b, integer maxb = 59) 
 
	atom p 
	p = power(a, b) 
	if integer(a) and integer(b) and b > 0 then           -- if result should be integer ... 
		if not integer(p) then                        -- ... but isn't 
			if abs(p) < 1e18 and b <= maxb then   -- maxb limit to avoid using too much time for this, not sure if it is necessary 
				p = a                            -- power(2, 60) is > 1e18, therefore 59 is the highest possible exponent 
				for j = 2 to b do 
					p = p * a 
				end for 
			end if 
		end if 
	end if 
	return p 
end function 

? power(3, 18)
? integer(power(3, 18))
? power_int(3, 18)
? integer(power_int(3, 18))
387420489.0
0
387420489
1

What I still find a bit troubling, though, is that 8.98 * 100 is not an integer. In the upcoming 64-bit version of Hypatia I'll use a function (which can be disabled) which rounds non-integer values to 17 significant digits, and if that is integer, then it replaces the original value — meaning, everything very close to integer is considered to be intended to be integer — this also resolves the 8.98* 100 issue well enough for me.

I think that, with %e working correctly, I have all I need for all practical purposes. Still have to be wary, though.

Robert

13. Re: Phix, numbers, printf, log10

• Posted by petelomax in February
• 896 views

I missed a *10, new version at same link

15. Re: Phix, numbers, printf, log10

• Posted by petelomax in February
• 805 views

Thanks for letting me know. Entirely independently of all this I had also previously (8/2/24)

Bugfix: couple of sprint() in pDiagN.e still passing maxlen to asChar. 

So that does not happen here on the dev system, and I assume it is/was something mid-exception handling.
Should you need it, grab yourself an updated copy of pDiagN.e, and that should fix it. Or just leave it till 1.0.5 ships.

17. Re: Phix, numbers, printf, log10

• Posted by RobertS in April
• 653 views

Happy end, Hypatia 4.0 is here, with 18 digits integers and 15 digits floating-point accuracy. I wrote my own replacements for power(), sprintf(), to_number() and round(), and I think rounding is working the way it should now. In case anyone should be interested, we can talk about it. The "%e" fix was essential (thanks again, Pete!), and also never to forget to use / power(10, n) instead of * power(10, -n).

hypatia-rpn.net

Btw, I have no idea how this happens, but in March Hypatia had 1387 downloads. No feedback whatever, though.

Robert

19. Re: Phix, numbers, printf, log10

• Posted by axtens_bruce in April
• 636 views

bugmagnet@LAPTOP-H6HBEGA9:~/hypatia/source$ ~/phix/p hy.exw 
Hypatia RPN Calculator 3.3 
? ? 
? 
For help use HELP command 
 
I/O file location: program folder /home/bugmagnet/phix/ 
External editor: 
Current screen width not available under Linux 
Result format: decimal 
Auto include $: ON 
Angle mode: radians 
Zero threshold: 1e-13 
Max. number of passes in loop: 100000 
Echo mode: OFF 
Copy all results to clipboard: OFF 
Logging: OFF 
Buffer mode: OFF 

21. Re: Phix, numbers, printf, log10

• Posted by axtens_bruce in April
• 635 views

Hmmm ... no command line editing. And apparently then no history of commands. Or is it just that Hypatia hasn't been clued in to the relevant keystrokes?

-Bruce

23. Re: Phix, numbers, printf, log10

• Posted by axtens_bruce in April
• 621 views

And it looks like it's a problem for everyone: https://unix.stackexchange.com/questions/714209/get-cursor-location-in-tty-without-reading-from-stdin-need-help-with-redirectio

25. Re: Phix, numbers, printf, log10

• Posted by petelomax in April
• 509 views

I think the latest docs say everything that needs to be said: http://phix.x10.mx/docs/html/get_position.htm

The correct/easier way to properly solve this would be to write a test harness with (nb 1st upper) Print(), Printf(), Qu(), and Clear_Screen(), Position(), Put_Screen_Char(), and Display_Text_Image() shims, (that obviously pre-process their inputs to update some central position vars, with no real need for any low-level grubbiness, before invoking the real thing) with a Get_Position() that checks it would yield the same values as get_position(). It makes much more sense for someone working on Linux to write that on Linux (just as it makes much more sense for the person actually dealing with all that to actually care), and of course in Euphoria rather than Phix, since the whole point is to craft/test the code that would need to be added to the Phix builtins. Candidates with no formal qualifications in unicode and ncurses need not bother applying.