Rounding error in (/)/2

Help!
1

There seems to be also a rounding challenge in (/)/2. But I didn’t
provide a test case for that. My test case here was only for (**)/2.
What test cases do you have for (/)/2 ?

Rounding error in (**)/2
https://swi-prolog.discourse.group/t/rounding-error-in-2/5795/13

mpq_get_d_nearest adds some bits to the numerator to increase
precision. And uses modulo to determine whether HALF_EVEN can
round up and doesn’t need to consider the EVEN rule.

Roughly the algorithm does the following, so instead that it computes,
assuming a and b have already the same exponent

a // b

It computes:

a*2^54 // b

And from this it determines the float mantissa, and uses an adjusted
exponent. Ater this it does the aforenentioned rounding, Possibly the
result of a / 1 is indeed some float(a) with HALF_EVEN, but its utterly slow!

Edit 25.09.2022
Please be patient I will try to produce some test cases. In
the past ECLiPSe Prolog had the most precise rounding for (/)/2.
I already posted about this in comp.lang.prolog. But this was

a while ago. But I guess it is worth to investigate the same
again, and also do some counting testing, which gives a
statistical view of the precision. In my systems I use a

a / b which is bootstrapped from float(a) and float(b), those
float/1 implementations with HALF_EVEN rounding, and in
many cases it is on par with what ECLIPSe Prolog does.

Now lets look into SWI-Prolog as well.

2

Actually now I have the feeling ECLiPSe Prolog has degraded.
But there is now the issue, suppose there is a new SWI-Prolog
version with a new float/1 compatible with (**)/2,

do the results still hold?

?- aggregate_all(count, (between(1,1000000,N), 
(123456789*123456789*N)/123456789 =\= float(123456789*N)), C).

/* SWI-Prolog 8.5.17 */
C = 0

/* Jekejeke Prolog 1.5.4, both JDK 1.8 and JDK 19 */
C = 72226

/* ECLiPSe Prolog 7.0.61 */
C = 494175

Edit 25.09.2022
And of course there is the question of timing. SWI-Prolog,
the version 8.5.17 I have available, is rather behind, possibly rather
an issue of (/)/2, but could be an issue of both (/)/2

and float/1 in the future:

?- time((between(1,1000000,N), (123456789*123456789*N)/123456789
=\= float(123456789*N), fail; true)).

/* SWI-Prolog 8.5.17 */
% 2,000,002 inferences, 1.297 CPU in 1.285 seconds (101% CPU, 1542170 Lips)

/* Jekejeke Prolog 1.5.4, JDK 19 */
% Threads 328 ms, GC 3 ms, Up 333 ms (Current 09/25/22 16:01:26)

/* Jekejeke Prolog 1.5.4, JDK 1.8 */
% Threads 250 ms, GC 3 ms, Up 261 ms (Current 09/25/22 16:05:51)

/* ECLiPSe Prolog 7.0.61 */
% Yes (1.16s cpu)
3

I guess that is on Windows? GMP seems really slow on the Windows version. On Linux I get

% 2,000,003 inferences, 0.371 CPU in 0.371 seconds (100% CPU, 5394619 Lips)

I haven’t studied the problem. It could be that the Windows malloc implementation is much worse or GMP is compiled differently. As it both concerns GCC 11 on x86_64 one would assume similar results.

edit: We should suspect how GMP is compiled. Using a really simple malloc replacement that gets small allocations from a memory buffer we get practically the same results on both Linux and Windows. That is, the Windows versions is still about 4 times slower.

4

Its also lacking behind on Unix, via WSL2 on the same machine,
but this is not yet a version with some new float/1:

/* SWI-Prolog 8.4.2 */
% 2,000,001 inferences, 0.543 CPU in 0.543 seconds (100% CPU, 3685380 Lips)

JDK 1.8 does it in 261 ms. But I use a different (/)/2 algorithm.
Actually one could also introduce two float division algorithms, similar
to two power raise functions (**)/2 and (^)/2:

  • Fast but less precise (/)/2: Calling X is A/B will use a less precise algorithm,
    like for example float(A)/float(B). The benefit would be more speedier,
    disadvantage less precision. The error is not that bad, see above in the
    example I get only 72226 faults. Same idea could be applied to (**)/2.

  • Slow but more precise (rdiv)/2 and float/1 combo: Call X is float(A rdiv B) will
    use a more precise algorithm, for example divmod/4 and extra shifting
    based as in mpq_get_d_nearest . The benefit would be more precise,
    disadvantage less speed. Same idea could be appled to (^)/2.

Edit 25.09.2022
I also tested Trealla Prolog, Ciao Playground and Scryer Prolog, here precision:

/* Trealla Prolog 2.2.6 */
C = 72226

/* Ciao Playground 1.21.0 */
C = 72226

/* Scryer Prolog 0.9.0 */
C = 494423

Here the time measurements, Trealla Prolog nearly as fast as JDK 1.8:

/* Trealla Prolog 2.2.6 */
% Time elapsed 0.293s

/* Ciao Playground 1.21.0 */
% Wall time 931.0 ms

/* Scryer Prolog 0.9.0 */
%   % CPU time: 1.559s
5

One more question. Do you have a smallint fast path for A/B, or
does it always detour over some kind of mpq_get_d_nearest?
The advantage of float(A)/float(B) is, it is already a smallint

fastpath for A and B smallint. For testing the smallint fast
path would need another test case.

Edit 25.09.2022
Ok, this could serve as a test case for smallint fast path, but I did
not have a SWI-Prolog with some optimized float/1 available, so
figures might improve in the future:

?- between(1,1000000,N), (33*33*N)/33
=\= float(33*N).

Some results concering the claim float(A)/float(B) is a fast path
for smallint arguments, measuring the precision:

/* SWI-Prolog, Windows and Unix WSL2 */
C = 0

/* Jekejeke Prolog 1.5.4, both JDK 1.8 and JDK 19 */
C = 0

Some results concerning timing:

/* SWI-Prolog 8.5.17, Windows */
% 2,000,003 inferences, 0.297 CPU in 0.290 seconds (102% CPU, 6736852 Lips)

/* SWI-Prolog 8.4.2, Unix WSL2 */
% 2,000,003 inferences, 0.337 CPU in 0.337 seconds (100% CPU, 5929587 Lips)

/* Jekejeke Prolog 1.5.4, JDK 19 */
% Threads 297 ms, GC 3 ms, Up 290 ms (Current 09/25/22 18:38:00)

/* Jekejeke Prolog 1.5.4, JDK 1.8 */
% Threads 188 ms, GC 2 ms, Up 186 ms (Current 09/25/22 18:39:33)
6

Smalint fastpath results for (/)/2 and for other Prolog
systems, only the timing, otherwise they all agree C=0:

/* Trealla Prolog 2.2.6 */
% Time elapsed 0.151s

/* ECLiPSe Prolog 7.0.61 */
% Yes (0.75s cpu)

/* Ciao Playground 1.21.0 */
% Wall time 784.0 ms

/* Scryer Prolog 0.9.0 */
% CPU time: 1.314s

The only good result is by Trealla Prolog.

7

Just for completeness:

/* Ciao Prolog 1.21 */
% 0.197 seconds

8

So Ciao WASM, used by Ciao Playground, has 3-4 times slower arithmetic?
I get on my machine similar results, via WSL2:

/* Ciao Prolog 1.21.0, (2022-09-06 00:43:43 +0200) [LINUXx86_64] */
?- time((between(1,1000000,N), (123456789*123456789*N)/123456789
=\= float(123456789*N), fail; true)).
% 0.24308399999999997 seconds

?- time((between(1,1000000,N), (33*33*N)/33
=\= float(33*N), fail; true)).
% 0.20234699999999997 seconds

In SWI-Prolog its the other way around. For example I get the first test a little faster,
on the same machine, using this SWI-Prolog WASM shell:

/* SWI-Prolog (32 bits, version 8.5.17, WASM */
?- time((between(1,1000000,N), (123456789*123456789*N)/123456789
=\= float(123456789*N), fail; true)).
% 2,001,401 inferences, 1.079 CPU in 1.079 seconds (100% CPU, 1854867 Lips)
true.

?- time((between(1,1000000,N), (33*33*N)/33
=\= float(33*N), fail; true)).
% 2,001,408 inferences, 0.629 CPU in 0.629 seconds (100% CPU, 3181889 Lips)
true.

Maybe there is a difference in 32-bit GMP and 64-bit GMP,
or do I already test some new version of SWI-Prolog?

10

Recent patches in the bigint-to-float branch fix this. As you seem to have a Linux system, maybe you should checkout this branch :slight_smile:

12

No. Not planning for that either. As an open project it is hard to make such claims. Especially in this case, it is not the OS, but the C compiler and math libraries used to build SWI-Prolog that cause the differences (and sometimes possibly the CPU). As we do not control these it is hard to make any claims. The docs should point at such issues where this applies.

If someone wants it, the plunit framework provides hooks to create the reports that are also used to make it emit output that complies to various unit test frameworks. Use these and write the aggregator.

16

For some reason the continuous integration at SWI-Prolog Continuous integration status stopped at version 8.5.9

17

That isn’t updated very often. It has to be triggered by hand. I still do not know how to deal with this. All in all there are so many configurations and targets on which to build that it is practically impossible to cover them all. Often it is all not that useful as changes do not touch anything that may result in a different test result in different configurations or targets.

I’ve also looked with @ericzinda into using the GitHub CI services. If you do comprehensive tests this way you quickly exceed the amount of free CPU one is allowed to use, making it rather costly at not much gain. The one running on the dev server allows for both incremental and full builds, which makes it faster. But then, incremental builds do not 100% work; sometimes some cmake tweaks are needed to reconfigure certain parts.

For each release we go through compiling and testing on Ubuntu Linux, MacOS, Win32, Win64 and soon also WASM.

2 Likes
18

Was doing regression testing versus (/)/2. There is indeed an
improvement from version 8.5.17 to version 8.5.18. For
example I find the following difference between the two versions:

/* SWI-Prolog 8.5.17 */
?- X is 370370367037037036703703703671 / 123456789012345678901234567890.
X = 2.9999999999999996.

/* SWI-Prolog 8.5.18 */
?- X is 370370367037037036703703703671 / 123456789012345678901234567890.
X = 3.0000000000000004.

Unfortunately I found a ISO core standard passage, which defines
(/)/2 in some particular way, which gives my idea of having two (/)/2
a new meaning, one less precise and fast and one more precise

and slow. Namely the ISO core standard says this, page 116:

Unbenannt2

So the ISO core standard rather sees (/)/2 associated with the
less precise and fast and not with the more precise and slow.
Testing shows that on the other hand SWI-Prolog pursues

the more precise and slow:

/* SWI-Prolog 8.5.18 */
/* more precise and slow */
?- X is 370370367037037036703703703670 / 123456789012345678901234567890.
X = 3.

/* what the ISO core standard would want, if we take the standard literally */
/* ?- X is float(370370367037037036703703703670) /
           float(123456789012345678901234567890). */
X = 3.0000000000000004.
20

Unless the iso flag is set to true (default false):

?- set_prolog_flag(iso,true).
true.

?- X is 370370367037037036703703703670 / 123456789012345678901234567890.
X = 3.0000000000000004.
1 Like
23

Did I miss the test driver? It is not very clear what you are measuring. Input? Output? Division? MPQ → float conversion? And then, is the C runtime library involved, which may make the results compiler/platform dependent?

SWI-Prolog division using / performs integer division first. If the result is exact it returns an integer. Otherwise it does the usual float thing. Unless the flag prefer_rationals is true. In that case the result is a rational number. With an additional constraint that you can set a maximum size for rationals and what to do if this is exceeded (error or continue using floats).

I’ve always claimed ISO arithmetic is flawed. Using a dynamically typed language one should not be forced to return incorrect results when there are correct results. If you want floating point arithmetic, just introduce one argument to be float and all will be converted to float.

24

The tests results have nothing to do with ISO core
standard. I added the 10 and 20 failing test cases
manually to cases.p, since the repeat loop takes

around 100’000 iterations to find them, needle in the
haystack, and my test cases are only 100 long. Sorry,
didn’t post the test driver. I have started just now writing

test drivers, and they already land on my sources pages.
And I can post the SWI-Prolog test driver here, to be used
together with the previous cases.p, it is this file:

swi.p…log (2,1 KB)

Edit 14.10.2022
The tests results have nothing to do with ISO core
standard for the following reasons, its similar like when
I tested the parser. I am only cross testing with other systems.

Whereas you have ambiguity in the parser, you have different
rounding, maybe in the float numbers. Although somebody said
somewhere C has some default which is always the same.

You can read it as a disagreement matrice. For example the ISO
core standard doesn’t specify the rnd_F function. It only lists a
few laws it should satisfiy, like for example:

- rnd_F(x) = -rnd_F(x)

Which narrows down the possible rounding functions.
But I don’t find a particular float rounding function chosen.
It says “implementation defined”:

Unbenannt2

Thats why I though asking what “to_nearest” means
with respect to SWI-Prolog important. One could use this
information also to verify some “inner” inconsistencies.

25

For example this test case shows an “inner” inconsistency,
using the emulator float_half_even2/2 that provides a corrected
“to_nearest”, since the “to_nearest” in SWI-Prolog 8.5.18 has

still a few bugs, see case2/3 and case3/3 for float/1 bugs:

From an erroneous float/1 one can unlikely implement a correct (/)/2. So
I guess half of the case4/4 and case5/4 failures in SWI-Prolog 8.5.18
are maybe induced by case2/3 and case3/3 failures in SWI-Prolog 8.5.18.

But can’t say 100% whether there is this failure dependency, that
causes a cascading failure in the arithmetic, didn’t check the
source code of SWI-Prolog 8.5.18.

26

I recently submitted a PR to fix a bug in converting GMP bigints to doubles which should appear in 8.5.19. With this fix, using your driver with the last posted version of the test cases:

?- suite.
case, swi: 0
case2, swi: 0
case3, swi: 0
case4, swi: 0
case5, swi: 0
true.

so I think all the failing tests in 8.5.18 are due to the conversion bug.

2 Likes
27

Small test on my system (latest git version)

swipl swi.pl
Welcome to SWI-Prolog (threaded, 64 bits, version 8.5.18-36-g4c76db827)
...

3 ?- swi.
case, swi: 0
case2, swi: 0
case3, swi: 0
case4, swi: 0
case5, swi: 0
case6, swi: 0
true.

But, on Windows (same sources) … Again the sloppy pow() from the MinGW runtime I assume :frowning:

1 ?- swi.
case, swi: 0
case2, swi: 0
case3, swi: 0
case4, swi: 0
case5, swi: 0
case6, swi: 74

You can also see it running in SWISH at SWISH -- SWI-Prolog for SHaring
Note that setting the iso flag is not allowed by the sandbox (I tried in my local copy of SWISH, but SWISH doesn’t work in iso mode). This is solved by changing A is P/Q in to A is float(P)/Q, forcing the float promotion.