One more question. Do SWI-Prolog 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. That
is a smallint fastpath, seen that this has count C=0:

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

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)

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.

Mostowski Collapse schrieb:
> One more question. Do SWI-Prolog 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. That
> is a smallint fastpath, seen that this has count C=0:
>
> ?- between(1,1000000,N), (33*33*N)/33
> =\= float(33*N).
>
> 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)

Here is a 100% pure Prolog implementation of HALF_EVEN, without
any materialization, it assumes that the current float/1 does only mpz_get_d:

float_half_even(X,Y) :-
(getbit(X, msb(X)-53) =:= 0;
getbit(X, msb(X)-52) =:= 0, msb(X)-53 =:= lsb(X)), !, Y is float(X).
float_half_even(X,Y) :-
Y is nexttoward(float(X), 1E300).

Was using 1E300 for infinity, is there something else that would work?
float_half_even/2 is just a translation of the following into Prolog:

public double doubleValue() {
https://hg.openjdk.java.net/jdk8/jdk8/jdk/file/tip/src/share/classes/java/math/BigInteger.java#l3871

Here some test whether it agrees with BigInteger.
doubleValue(), first testing the old version of
SWI-Prolog, where float is mpz_get_d:

/* SWI-Prolog 8.5.17 */
?- between(1,10,N), X is 10^20+N*2^13,
float_half_even(X,Y), write(Y), nl, fail; true.
1.0e+20
1.0000000000000002e+20
1.0000000000000003e+20
1.0000000000000003e+20
1.0000000000000003e+20
1.0000000000000005e+20
1.0000000000000007e+20
1.0000000000000007e+20
1.0000000000000007e+20
1.0000000000000008e+20

Comparing with a Prolog system that calls
BigInteger.doubleValue():

/* Jekejeke Prolog 1.5.4 */
?- between(1,10,N), X is 10^20+N*2^13,
Y is float(X), write(Y), nl, fail; true.
1.0E20
1.0000000000000002E20
1.0000000000000003E20
1.0000000000000003E20
1.0000000000000003E20
1.0000000000000005E20
1.0000000000000007E20
1.0000000000000007E20
1.0000000000000007E20
1.0000000000000008E20

Yeah, they agree! Riddle solved, the misterious
veil of the secret souce behind an ultrafast
float/1 for bigint has been lifted.

Mostowski Collapse wrote:
> Here is a 100% pure Prolog implementation of HALF_EVEN, without
> any materialization, it assumes that the current float/1 does only mpz_get_d:
>
> float_half_even(X,Y) :-
> (getbit(X, msb(X)-53) =:= 0;
> getbit(X, msb(X)-52) =:= 0, msb(X)-53 =:= lsb(X)), !, Y is float(X).
> float_half_even(X,Y) :-
> Y is nexttoward(float(X), 1E300).
>
> Was using 1E300 for infinity, is there something else that would work?
> float_half_even/2 is just a translation of the following into Prolog:
>
> public double doubleValue() {
> https://hg.openjdk.java.net/jdk8/jdk8/jdk/file/tip/src/share/classes/java/math/BigInteger.java#l3871
>

Quite some disagreement concerning (**)/2:

/* Ciao Playground 1.21.0 */
?- Z is 51**10.
Z = 119042423827613000.0 ?

/* ECLiPSe Prolog 7.0.61 */
?- X is 51**10.
X = 1.1904242382761301e+17

Guess which one corresponds to HALF_EVEN ?

Whats also funny when a Prolog system disagrees internally:

/* SWI-Prolog 8.5.17, Windows 10 */
?- X is 51^10, Z is float(X), Y is nexttoward(Z, 1E300).
X = 119042423827613001,
Z = 1.1904242382761301e+17,
Y = 1.1904242382761302e+17.

?- Z is float(51^10), Y is nexttoward(Z, 1E300).
Z = 1.19042423827613e+17,
Y = 1.1904242382761301e+17.

The above "bug" makes it hard for me to test float_half_even/2.

Mostowski Collapse schrieb am Mittwoch, 28. September 2022 um 12:33:22 UTC+2:
> Quite some disagreement concerning (**)/2:
> /* Ciao Playground 1.21.0 */
> ?- Z is 51**10.
> Z = 119042423827613000.0 ?
> /* ECLiPSe Prolog 7.0.61 */
> ?- X is 51**10.
> X = 1.1904242382761301e+17
>
> Guess which one corresponds to HALF_EVEN ?

If I do nothing, I cannot make mistakes. Could be
the motto that Ciao Prolog doesn't support (^)/2
as an evaluable predicate.

But among those Prolog systems that take the
risk (clap clap!!) and support (^)/2 the situation
is a little immature, I find:

Scryer Prolog and Trealla Prolog, they both do wrong
float/1 rounding (1.19042423827613e+17
instead 1.1904242382761301E17):

/* Scryer Prolog 0.9.0 */
?- X is 51^10, Z is float(X).
X = 119042423827613001, Z = 1.19042423827613e17.
?- Z is float(51^10).
Z = 1.19042423827613e17.

But Trealla Prolog has a problem with (^)/2, wrong
result of bigint pow() (i.e. 119042423827613008 instead
119042423827613001):

/* Trealla Prolog 2.2.10 */
?- X is 51^10, Z is float(X).
X = 119042423827613008, Z = 1.19042423827613e+17.
?- Z is float(51^10).
Z = 1.19042423827613e+17.

Mostowski Collapse schrieb am Mittwoch, 28. September 2022 um 12:34:44 UTC+2:
> Whats also funny when a Prolog system disagrees internally:
>
> /* SWI-Prolog 8.5.17, Windows 10 */
> ?- X is 51^10, Z is float(X), Y is nexttoward(Z, 1E300).
> X = 119042423827613001,
> Z = 1.1904242382761301e+17,
> Y = 1.1904242382761302e+17.
>
> ?- Z is float(51^10), Y is nexttoward(Z, 1E300).
> Z = 1.19042423827613e+17,
> Y = 1.1904242382761301e+17.
>
> The above "bug" makes it hard for me to test float_half_even/2.
> Mostowski Collapse schrieb am Mittwoch, 28. September 2022 um 12:33:22 UTC+2:
> > Quite some disagreement concerning (**)/2:
> > /* Ciao Playground 1.21.0 */
> > ?- Z is 51**10.
> > Z = 119042423827613000.0 ?
> > /* ECLiPSe Prolog 7.0.61 */
> > ?- X is 51**10.
> > X = 1.1904242382761301e+17
> >
> > Guess which one corresponds to HALF_EVEN ?

To judge that ECLiPSe Prolog is correct, one
needs to think in binary representation. We find,
using a precise float to bigint conversion:

?- X is integer(1.1904242382761301E17).
X = 119042423827613008.

?- X is integer(1.19042423827613E17).
X = 119042423827612992.

Now which one is closer?

?- X is 119042423827613001-119042423827613008.
X = -7.

?- X is 119042423827613001-119042423827612992.
X = 9.

So 1.1904242382761301E17 is the correct float.

Mostowski Collapse schrieb am Mittwoch, 28. September 2022 um 12:33:22 UTC+2:
> Quite some disagreement concerning (**)/2:
> /* Ciao Playground 1.21.0 */
> ?- Z is 51**10.
> Z = 119042423827613000.0 ?
> /* ECLiPSe Prolog 7.0.61 */
> ?- X is 51**10.
> X = 1.1904242382761301e+17
>
> Guess which one corresponds to HALF_EVEN ?

Be careful, integer/1 isn't in the ISO core standard, thats
my take in one system to make it a precise conversion from
float to bigint. Maybe not reproducable in other Prolog systems.

It works also in Dogelog player:

/* Dogelog Player, JavaScript */
?- X is integer(1.1904242382761301E17).
X = 119042423827613008.
?- X is integer(1.19042423827613E17).
X = 119042423827612992.

/* Dogelog Player, Python */
?- X is integer(1.1904242382761301E17).
X = 119042423827613008.
?- X is integer(1.19042423827613E17).
X = 119042423827612992.

Mostowski Collapse schrieb am Mittwoch, 28. September 2022 um 12:56:30 UTC+2:
> To judge that ECLiPSe Prolog is correct, one
> needs to think in binary representation. We find,
> using a precise float to bigint conversion:
>
> ?- X is integer(1.1904242382761301E17).
> X = 119042423827613008.
>
> ?- X is integer(1.19042423827613E17).
> X = 119042423827612992.
>
> Now which one is closer?
>
> ?- X is 119042423827613001-119042423827613008.
> X = -7.
>
> ?- X is 119042423827613001-119042423827612992.
> X = 9.
>
> So 1.1904242382761301E17 is the correct float.
> Mostowski Collapse schrieb am Mittwoch, 28. September 2022 um 12:33:22 UTC+2:
> > Quite some disagreement concerning (**)/2:
> > /* Ciao Playground 1.21.0 */
> > ?- Z is 51**10.
> > Z = 119042423827613000.0 ?
> > /* ECLiPSe Prolog 7.0.61 */
> > ?- X is 51**10.
> > X = 1.1904242382761301e+17
> >
> > Guess which one corresponds to HALF_EVEN ?

Warning, possibly this ticket can be closed. Its not necessarily
a display problem. For example many Math.pow(double,double)
implementations cannot do the problem correctly.

They cannot or do not round to the theoretical optimum, they
have inherently a higher error than only one Unit in the last place (ULP).
So a certain error is to expect. You see this here:

/* Dogelog Player, JavaScript */
?- X is 51**10.
X = 1.19042423827613E17.

?- X is float(51^10).
X = 119042423827613010.0.

In the above the (**)/2 is Math.pow(double,double) and the
(^)/2 is exact bigint exponentiation. I took the pair (51,10)
from a list of pairs where (**)/2 and (^)/2 often disagree

based on Math.pow(double,double) typically being weak.

Mostowski Collapse schrieb am Mittwoch, 28. September 2022 um 12:33:22 UTC+2:
> Quite some disagreement concerning (**)/2:
> /* Ciao Playground 1.21.0 */
> ?- Z is 51**10.
> Z = 119042423827613000.0 ?
> /* ECLiPSe Prolog 7.0.61 */
> ?- X is 51**10.
> X = 1.1904242382761301e+17
>
> Guess which one corresponds to HALF_EVEN ?