Correct floating point arithmetic

[jan]

@ridgeworks

That’s the kind of thing than should probably be done using the foreign language interface and separated from the core SWIP (along the lines of @jan 's suggestion earlier in the thread).

That was not my proposal If we hide mpfr or LibBF behind an interface that exchanges doubles we could consider having the core of the arithmetic using a array of function pointers. So, ar_sin() (the function implementing Prolog sin/1) does not call libm’s sin(), but instead calls something like funcs[FUNC_SIN](). Now we can add some API that allows a foreign library to replace funcs[FUNC_SIN]. That way we can avoid this whole discussion from affecting the core. The default are the libm functions and we can provide foreign libraries that do things differently. Now we can create mpfr wrappers or plugin crlibm. As anyone can contribute such libraries and distribute them as appropriate, people can choose between license, performance and correctness. I can live with that.

In contrast, adding mpfr, crlibm, LibBF to the core locks us to one of these choices and adds license and/or dependency issues to the core. I’m not against considering that, but if we go that far I’d also like to see the other benefits come to SWI-Prolog. This requires more work though as we need change the simple representation of floats as doubles. That might not be too bad as besides the arithmetic, read/write and term ordering, very little of the system is even aware floats exist.

Considering most of the system is agnostic about numbers anyway, it might be possible to define a clean abstract interface that allows anyone to add numeric types, their conversions and their functions.

[ridgeworks]

I was thinking more about:

@jan

In that light, would it not be a better to provide predicates for your correctly rounding float functions?, e.g., a library that provides sin(RoundMode, Angle, Result) (not sure about the argument order)?

which sounded more like an FLI solution to me.

I did see a C2x draft (somewhere?) that reserved all the cr_* function names so maybe that will be the final resolution if and when it happens.

Anything more, like "change the simple representation of floats as doubles. " isn’t of much interest to me, but maybe some other stakeholder will step up.

[ridgeworks]

@ridgeworks

I did see a C2x draft (somewhere?) that reserved all the cr_* function names so maybe that will be the final resolution if and when it happens.

For the record:
https://www.open-std.org/jtc1/sc22/wg14/www/docs/n3054.pdf (cf 7.33.8, item 4, page 457)

[jan]

Good. Eventually your problem may be solved with little work from our side These things tend to take long, but if there is a standard I suspect there is sooner or later a proper free implementation.

[ridgeworks]

@jan

Good. Eventually your problem may be solved with little work from our side These things tend to take long, but if there is a standard I suspect there is sooner or later a proper free implementation.

Hopefully the CORE-MATH project will help accelerate the process:

core-math-arith2022.pdf

From the Conclusion slide:

• We present for the first time implementations of correctly rounded routines as fast as in the best current math libraries or even faster.
• Full set of C99 binary32 functions ready for integration (either as expf or cr_expf).
• Full set of C99 binary64 functions planned for end of 2023.
• Time for IEEE-754 to require correct rounding!
• Not yet another libm, but aimed at integration in existing libms.

[peter.ludemann]

BTW, you might want to also add “cospi” and “sinpi” (and other radian-oriented trig functions?) - apparently they are a “hidden” part of the standard “scientific” package that Python provides, and provide better accuracy for those who do analytic geometry using radians.

- X is cos(pi/2).
X = 6.123233995736766e-17.

https://mail.python.org/archives/list/scipy-dev@python.org/thread/ES6FDLI5I4JF7GBQKJZIUZYJZRN2UT7I/

Also, such packages might also have other suggestions for maintaining precision.

[ridgeworks]

@j4n_bur53

It already fails somehow in SWI-Prolog 8.5.20, Windows differs from
Unix, …

What are your expectations? In general, programming languages depend on the math libraries (libm) supported by platform vendors (and others), and today these implementations do no purport to support correct rounding. For details of some C libraries see https://hal.inria.fr/hal-03141101 - recommended to read at least the first page.

There are some signs of progress in producing correctly rounded elementary functions, but I think it’s going to take a couple of years minimum.

[j4n_bur53]

@ridgeworks

There are some signs of progress in producing
correctly rounded elementary functions, but I think
it’s going to take a couple of years minimum.

In the case of (**)/2 and for smallint arguments, you
cannot blame low quality ar_pow() function with float
arguments. Its obviously a bug in the (**)/2 implementation

of SWI-Prolog 8.5.20. You can check yourself. Positive
exponents work fine:

/* SWI-Prolog 8.5.20 Ok */
?- X is float(1889^29).
X = 1.0250068783051207e+95.

/* Jekejeke Prolog 1.5.5 Ok */
?- X is 1889**29.
X = 1.0250068783051207E95

But there is a bug for negative exponents, it takes a float
fast path, which makes it less precise:

/* SWI-Prolog 8.5.20 Nok */
?- X is 1889** -29.
X = 9.756032092716607e-96.

/* Jekejeke Prolog 1.5.5 Ok */
?- X is 1889** (-29).
X = 9.756032092716584E-96.

This could be easily fixed, by using another computation
path for negative exponents, that is as exact as what
it does for positive exponents. This computation path is

also available in SWI-Prolog, I find:

?- set_prolog_flag(iso, false).
true

/* SWI-Prolog 8.5.20 Ok */
?- X is 1/1889^29.
X = 9.756032092716584e-96.

Edit 25.10.2022
So I guess by changing the computation path in SWI-Prolog,
the precision of negative exponents could be increased, without
waiting for cr_XXX libraries. I showed this bug already in the

(**)/2 thread (lookup multiplicative inverse), but it wasn’t
fixed in release 8.5.19 or release 8.5.20. One might also consider
to separate float fast path into a function (**)/2 and bigint based

power, slower but more precise, even for negative exponents into
a function (^)/2. So to give the application programmer a choice
between cr_XXX functions and bigint based approaches.

This way end-user would not need to wait for more precise cr_XXX
lib, and could already enjoy more precise from (^)/2 function. If
I would open a GitHub issue you could tag it won’t fix, I

wouldn’t mind. But I get always told I should use SWI-Prolog
discourse and not GitHub issues, so for a year or so, already
I don’t use GitHub issues anymore. But on SWI-Prolog discourse,

I didn’t get a feedback yet about the roadmap of SWI-Prolog
concerning some very particular issues of (**)/2 and (^)/2.

[ridgeworks]

@j4n_bur53

But there is a bug for negative exponents, it takes a float
fast path, which makes it less precise:

/* SWI-Prolog 8.5.20 Nok */
?- X is 1889** -29.
X = 9.756032092716607e-96.

/* Jekejeke Prolog 1.5.5 Ok */
?- X is 1889** (-29).
X = 9.756032092716584E-96.

It still looks like a libm or GMP issue to me; On Mac-Intel, 8.5.20:

?- X is 1889** -29.
X = 9.756032092716584e-96.

[j4n_bur53]

@ridgeworks

It still looks like a libm or GMP issue to me; On Mac-Intel, 8.5.20:

?- X is 1889** -29.
X = 9.756032092716584e-96.

You didn’t test thoroughly, I now also found Unix failures:

/* SWI-Prolog 8.5.20 WSL2 Nok */
?- X is -18573** -7.
X = -1.3116730299820408e-30.

?- X is 1/(-18573**7).
X = -1.3116730299820406e-30.

This was just announced 1hour ago, your post is from 28min ago?
It was announced here, but affects all threads dealing with
correct power functions, also this thread here and the libBf thread:
https://swi-prolog.discourse.group/t/rounding-error-in-2/5795/67

I will reopen my testing compaign and produce a new case6/3,
but hopefully after that, it should be clear that the Unix platform
is no exception, and also subject to the “nearly rounded” phaenomenom.

Thats just a false rumor that only Windows is affected, which was
sustained because nobody did some more extended testing. I couldn’t
test before now, because I didn’t have a Unix SWIPL installation,

WSL2 building didn’t work on my side, had to fall back to PPA.
The new example is again an example where the arguments are
still smallint and are exactly representable as float as well.

[jan]

You can keep testing float operations at infinitum, but all you are testing is libm and hardware/compiler support for the float computations. We already know that none of these are perfect. Even your JDK might give guarantees under strict math if this is defined by the language, but otherwise you probably also depend on the choice of C/C++ compiler and runtime library used to build the JDK runtime system.

With the work of @ridgeworks (on comments by you), translation of bigints and rationals to floats is hopefully correct. So, do the math using rationals and you get a correct float result Otherwise results are going to be wrong, MinGW a lot worse than Linux or MacOS.

[ridgeworks]

@j4n_bur53

I thought with iso=false it never
uses float in the division, in the sense that arguments would
be converted to float, and would give the same as float( rdiv )
combo?

Not quite. Integer division will only produce a (non-integer) rational when the prefer_rationals global flag is also true. See SWI-Prolog -- Rational number examples .

@j4n_bur53

The early float conversion introduces a rounding error

If the iso flag is true or the prefer_rationals is false, and the integer division is not exact, the two arguments to the ar_divide function are converted to C doubles and then are divided via normal compiled C arithmetic (at least as I understand the source). If I do this explicitly, I get:

?- X is 1.0/float(18573^7).
X = 1.3116730299820406e-30.

But if everything is done using rational arithmetic:

?- set_prolog_flag(iso,false), set_prolog_flag(prefer_rationals,true).
true.

?- X is 1/(18573^7),F is float(X).
X = 1r762385119722780192080867194597,
F = 1.3116730299820408e-30.

I appreciate that this is a little confusing and it definitely complicates comparing results across different platforms. But I don’t think there are any mysteries here.

[jan]

@ridgeworks

Not quite. Integer division will only produce a (non-integer) rational when the prefer_rationals global flag is also true. See SWI-Prolog – Manual .

Actually, this is quite confusing as I noticed. See below. The problem is that rational exponentiation performs rational arithmetic, while integer performs float arithmetic. As integers are rationals and values are always kept in canonical form, the second is considered integer arithmetic.

?- A is 7r5 ^ -5.
A = 3125r16807.

?- A is 7r1 ^ -5.
A = 5.9499018266198606e-5.

@ridgeworks

I appreciate that this is a little confusing

This, I think, is not confusing. It is the same as A = float(X*Y) vs A is X*float(Y) and works at is is supposed to in ISO Prolog.

@j4n_bur53

+Expr1 / +Expr2
Otherwise (both arguments are integer) the result is integer
if the division is exact and float otherwise.

That is pretty outdated I’ll update that. edit You’re only picking a little fragment of the long documentation. I think the documentation is accurate, be it a bit hard to read due to all the conditions. Possibly should be a table? If anyone wants to send a PR with a clearer description, please do.

[jan]

@j4n_bur53

Could SWI-Prolog also directly produce the correctly rounded one,
when invoking the evaluable function (/)/2 and when iso=false
and prefer_rationals=false ?

If you run it on a system that has a good pow() function, yes. If not, you need to set prefer_rationals to true and you get

?- A is float(2**100 / 18573**7).
A = 1.6627431037599143.

The rules for integer division (and int ** -int) are not that hard:

• If current_prolog_flag(prefer_rationals, true) holds
  • Integer division returns a rational number.
• else if current_prolog_flag(iso, true) holds
  • Integer division is performed after converting both operands to a float (and thus returns a float).
• else
  • Integer division returns an integer if the division is exact. Otherwise Integer division is performed after converting both operands to a float (and thus returns a float).

If you want proper arithmetic, leave iso false, set prefer_rationals and I also like rational_syntax set to natural, so it reads and writes n/m as a rational number. Initially I could not convince Joachim this made sense. Later he agreed it makes as much sense as that -n reads as a single token. This should have been the ISO standard I’m still tempted to make this the default. But yes, it will break some applications and will make some slow. I use it on daily basis though and experience very few problems with it. The main problem is while using Prolog as a calculator and you get this rather uninformative answer

?- A is 431/42.
A = 431/42.

[jan]

@j4n_bur53

Could SWI-Prolog also directly produce the correctly rounded one,
when invoking the evaluable function (/)/2 and when iso=false
and prefer_rationals=false ?

Finally I see what you mean. Well, division is defined to convert to float first, in this setting also for SWI-Prolog. So, the result is the same as

?- X is float(1267650600228229401496703205376)/float(762385119722780192080867194597).
X = 1.662743103759914.

Now, the rounding is correct. The one from 762385119722780192080 looked a bit suspicious, but is correct; the first is indeed closer.

?- X is float(762385119722780192080867194597).
X = 7.623851197227803e+29.

?- X is nexttoward(float(762385119722780192080867194597), 0).
X = 7.623851197227801e+29.

So, that is the expected result. And yes, it can be precise if we define integer division as rational number division and conversion to float. A quick check says this is about 2 times slower on these specific numbers. I don’t know whether this is a good idea. After all, we have a good way to deal with this using prefer_rationals. Probably it is a good idea if people can trade performance and precision as for 99% of the applications using floats a few ULPs is irrelevant.

[j4n_bur53]

Here is a trick to force Python compatibility mode, I guess
works for SWI-Prolog 8.5.20 which has improved rational
number rounding? Thats the incorrectly rounded result:

/* Incorrectly rounded of the real number quotient */
?- X is 1267650600228229401496703205376 /
           762385119722780192080867194597.
X = 1.662743103759914.

Now do switch on Python compatibility mode, and
watch the result getting correctly rounded:

/* Use this combo for Python compatibility:
             prefer_rationals=true,
             max_rational_size=0,
             max_rational_size_action=float */
?- set_prolog_flag(prefer_rationals, true).
true.

?- set_prolog_flag(max_rational_size, 0).
true.

?- set_prolog_flag(max_rational_size_action, float).
true.

/* Correctly rounded of the real number quotient, via
   detour over rational numbers, which are correctly
   rounded since release 8.5.20 of SWI-Prolog? */
?- X is 1267650600228229401496703205376 /
            762385119722780192080867194597.
X = 1.6627431037599143.

Edit 26.10.2022
BTW the real number quotient is this real number:

1267650600228229401496703205376 /
762385119722780192080867194597 =
1.66274310375991431383116365048...

Was using microsoft calculator desk accessory. But
maybe should do more testing of this new Python
mode trick. Nasty test cases might lure everywhere.

[jan]

Cute B.t.w., if you want arbitrary precision decimal representations, you can use format/2 as below. If the argument to ~f is an expression it is evaluated and if the result is a rational it uses GMP to output a decimal approximation rather than first converting the expression result to a float. Of course you must make sure it is evaluated the way you want, e.g., prefer_rationals should be true.

?- format('~50f', [1267650600228229401496703205376/762385119722780192080867194597]).
1.66274310375991431383116365048697983300087128704386

[j4n_bur53]

@jan

?- format(‘~50f’, …)

Pretty sure helpful for debugging. BTW, the Python mode completely
solves the pow() problem for integer arguments on all platforms and
irrespective of libm. I get this here:

/* SWI-Prolog 8.5.20 WSL2 Nok */
?- X is 28639** -2.
X = 1.2192262404758602e-9.

And then this here:

?- set_prolog_flag(prefer_rationals, true).
true.

?- set_prolog_flag(max_rational_size, 0).
true.

?- set_prolog_flag(max_rational_size_action, float).
true.

/* SWI-Prolog 8.5.20 WSL2 Ok */
?- X is 28639** -2.
X = 1.2192262404758604e-9.

Edit 26.10.2022
Somehow the Python mode can be also combined with
iso=true, and then (/)/2 is not anymore Python.

Here is a test harness:

swi.p.log (4,1 KB)

Now I get, on Windows:

?- swi.
case3, swi: 0
case4, swi: 0
case6, swi: 43

?- swi_python.
case3, swi: 0
case4, swi: 0
case6, swi: 0

A little take away, test cases for (**)/2 are correct.