Some suspicion about the random number skewing.
Mostlikely #define MAX_URANDOMM_ITER is too low. See
here what GMP does. And it then falls back to low quality result:

if (count == 0)
/* Too many iterations; return result mod n == result - n */
mpn_sub_n (rp, rp, np, size);
https://salsa.debian.org/science-team/gmp/-/blob/master/mpz/urandomm.c

Note that the example interval [0…2^n+1) is a very nasty non power
of 2 interval. That I was using a nasty interval was rather not planned,
when writing the fuzzer, just some accident.

In my implementation I have no MAX_URANDOMM_ITER limit.
My suspiciion, therefore it generates a better probability shaping
on my side, which can be even measured,

already for 100 bits? Is this the reason? Or its some other bug
in mpz_urandomm() caused by the last limb cmp juggling?

I didn’t test the random/1 evaluable function in Trealla and Cliao yet.
Cannot run the same test case in Trealla Prolog or Ciao Prolog,
too much stuff missing:

- Trealla Prolog: Didn’t find bigint random.

- Ciao Prolog: Found some random, with an incompatible
protocol to SWI-Prolog, different meaning of the parameters,
but there is so much other stuff missing, like
msb/1, aggregate_all/3, etc…

Mostowski Collapse schrieb am Mittwoch, 12. Oktober 2022 um 19:46:58 UTC+2:
> Some suspicion about the random number skewing.
> Mostlikely #define MAX_URANDOMM_ITER is too low. See
> here what GMP does. And it then falls back to low quality result:
>
> if (count == 0)
> /* Too many iterations; return result mod n == result - n */
> mpn_sub_n (rp, rp, np, size);
> https://salsa.debian.org/science-team/gmp/-/blob/master/mpz/urandomm.c
>
> Note that the example interval [0…2^n+1) is a very nasty non power
> of 2 interval. That I was using a nasty interval was rather not planned,
> when writing the fuzzer, just some accident.
>
> In my implementation I have no MAX_URANDOMM_ITER limit.
> My suspiciion, therefore it generates a better probability shaping
> on my side, which can be even measured,
>
> already for 100 bits? Is this the reason? Or its some other bug
> in mpz_urandomm() caused by the last limb cmp juggling?

On Ciao Prolog its:
random(0,M,...)

On SWI-Prolog its:
random(0,M+1,...)

Should there be a discussion about random/3, its protocoll?
I don’t know, what are the plans of Ciao Prolog and SWI-Prolog?
I don’t know? Is there a road map? I don’t know.

Recently SWI-Prolog deprecated random/3 and introduced
random_between/3 with a different protocoll. Maybe this can
point the way forward? I am more interested about the

skewness here, and not the protcoll of some random API.
Skewed shouldn’t mean screwed. Its quite possible that a
skewed random generator can be fixed.

Mostowski Collapse schrieb am Mittwoch, 12. Oktober 2022 um 19:48:05 UTC+2:
> I didn’t test the random/1 evaluable function in Trealla and Cliao yet.
> Cannot run the same test case in Trealla Prolog or Ciao Prolog,
> too much stuff missing:
>
> - Trealla Prolog: Didn’t find bigint random.
>
> - Ciao Prolog: Found some random, with an incompatible
> protocol to SWI-Prolog, different meaning of the parameters,
> but there is so much other stuff missing, like
> msb/1, aggregate_all/3, etc…
> Mostowski Collapse schrieb am Mittwoch, 12. Oktober 2022 um 19:46:58 UTC+2:
> > Some suspicion about the random number skewing.
> > Mostlikely #define MAX_URANDOMM_ITER is too low. See
> > here what GMP does. And it then falls back to low quality result:
> >
> > if (count == 0)
> > /* Too many iterations; return result mod n == result - n */
> > mpn_sub_n (rp, rp, np, size);
> > https://salsa.debian.org/science-team/gmp/-/blob/master/mpz/urandomm.c
> >
> > Note that the example interval [0…2^n+1) is a very nasty non power
> > of 2 interval. That I was using a nasty interval was rather not planned,
> > when writing the fuzzer, just some accident.
> >
> > In my implementation I have no MAX_URANDOMM_ITER limit.
> > My suspiciion, therefore it generates a better probability shaping
> > on my side, which can be even measured,
> >
> > already for 100 bits? Is this the reason? Or its some other bug
> > in mpz_urandomm() caused by the last limb cmp juggling?

So maybe the requirement is eclectic? A non skewed
random ganerator for bigints. Often power of 2 is only
required or smallint is enough.

Ha Ha, a Casino website could be interested that their
random generators are not skewed. Or maybe they even
go for a particular skewed random generator?

What if the users of the Casino website figure out the
skewing and play accordingly? I am thinking of games
like roulette or card games. So its in their own interest?

Legitimacy & Regulation of the RNG
There are a number of independent companies which specialise in iGaming testing.
https://www.casinoguardian.co.uk/casino-guide/random-number-generator/

Mostowski Collapse schrieb am Mittwoch, 12. Oktober 2022 um 22:22:29 UTC+2:
> On Ciao Prolog its:
> random(0,M,...)
>
> On SWI-Prolog its:
> random(0,M+1,...)
>
> Should there be a discussion about random/3, its protocoll?
> I don’t know, what are the plans of Ciao Prolog and SWI-Prolog?
> I don’t know? Is there a road map? I don’t know.
>
> Recently SWI-Prolog deprecated random/3 and introduced
> random_between/3 with a different protocoll. Maybe this can
> point the way forward? I am more interested about the
>
> skewness here, and not the protcoll of some random API.
> Skewed shouldn’t mean screwed. Its quite possible that a
> skewed random generator can be fixed.
> Mostowski Collapse schrieb am Mittwoch, 12. Oktober 2022 um 19:48:05 UTC+2:
> > I didn’t test the random/1 evaluable function in Trealla and Cliao yet.
> > Cannot run the same test case in Trealla Prolog or Ciao Prolog,
> > too much stuff missing:
> >
> > - Trealla Prolog: Didn’t find bigint random.
> >
> > - Ciao Prolog: Found some random, with an incompatible
> > protocol to SWI-Prolog, different meaning of the parameters,
> > but there is so much other stuff missing, like
> > msb/1, aggregate_all/3, etc…
> > Mostowski Collapse schrieb am Mittwoch, 12. Oktober 2022 um 19:46:58 UTC+2:
> > > Some suspicion about the random number skewing.
> > > Mostlikely #define MAX_URANDOMM_ITER is too low. See
> > > here what GMP does. And it then falls back to low quality result:
> > >
> > > if (count == 0)
> > > /* Too many iterations; return result mod n == result - n */
> > > mpn_sub_n (rp, rp, np, size);
> > > https://salsa.debian.org/science-team/gmp/-/blob/master/mpz/urandomm.c
> > >
> > > Note that the example interval [0…2^n+1) is a very nasty non power
> > > of 2 interval. That I was using a nasty interval was rather not planned,
> > > when writing the fuzzer, just some accident.
> > >
> > > In my implementation I have no MAX_URANDOMM_ITER limit.
> > > My suspiciion, therefore it generates a better probability shaping
> > > on my side, which can be even measured,
> > >
> > > already for 100 bits? Is this the reason? Or its some other bug
> > > in mpz_urandomm() caused by the last limb cmp juggling?

2 out of 3 Prolog system cannot generate uniform random
bigints? Don't know yet, have to dig deeper. Strange that
the below result has always 100.0:

The error in Ciao Prolog is much larger than in SWI-Prolog,
i.e. the distance to the expected value. Was using this query,
taking into account the different random/3 protocoll, and

manually compensating for the missing aggregate_all/3
and msb/1 evaluable function:

$ ciaosh
Ciao 1.22.0 (2022-09-28 21:20:35 +0200) [LINUXx86_64]
?- use_module(library(between)).

?- ['aggregate.p'].

?- E is (1<<100), M is (1<<101), K is (1<<20),
between(1,10,_), aggregate_all(sum(X),
(between(1,K,_), random(0,M,X)), S),
H is abs(S-E*K), msb(H,J), N is J-20, write(N), nl, fail; true.
100.0
100.0
100.0
100.0
100.0
100.0
100.0
100.0
100.0
100.0

And this helper code (file aggregate.p):

:- use_module(library(numlists)).
:- use_module(library(random)).

aggregate_all(sum(X), G, S) :-
findall(X, G, L),
sum_list(L, S).

msb(X, Y) :- Y is log(X)/log(2).

Mostowski Collapse schrieb am Donnerstag, 13. Oktober 2022 um 14:53:46 UTC+2:
> So maybe the requirement is eclectic? A non skewed
> random ganerator for bigints. Often power of 2 is only
> required or smallint is enough.
>
> Ha Ha, a Casino website could be interested that their
> random generators are not skewed. Or maybe they even
> go for a particular skewed random generator?
>
> What if the users of the Casino website figure out the
> skewing and play accordingly? I am thinking of games
> like roulette or card games. So its in their own interest?
>
> Legitimacy & Regulation of the RNG
> There are a number of independent companies which specialise in iGaming testing.
> https://www.casinoguardian.co.uk/casino-guide/random-number-generator/

Does the C floating point API have HALF_EVEN default?
Maybe, maybe not. Is there a test case to figure out
what rounding a system, Prolog or not, does, for its

floating point numbers, independent of some bigint
conversion issues? Was trying the folllowing test and got
some surprise. There was first quite some agreement:

/* SWI-Prolog 8.5.18 */
?- X is 665723957508124892951320985600.0+26388279066624.0.
X = 6.65723957508125e+29.

/* Ciao Prolog 1.22.0 */
?- X is 665723957508124892951320985600.0+26388279066624.0.
X = 6.65723957508125e29 ?

/* Trealla Prolog 2.4.3 */
?- X is 665723957508124892951320985600.0+26388279066624.0.
X = 6.65723957508125e+29.

/* Jekejeke Prolog 1.5.5 */
?- X is 665723957508124892951320985600.0+26388279066624.0.
X = 6.65723957508125E29.

But then suddently:

/* Scryer Prolog 0.9.0 */
?- X is 665723957508124892951320985600.0+26388279066624.0.
X = 6.657239575081248e29.

/* ECLiPSe Prolog 7.0.61 */
?- X is 665723957508124892951320985600.0+26388279066624.0.
X = 6.6572395750812482e+29

?- X is 6.6572395750812482e+29-6.65723957508125E29.
X = -140737488355328.0.

Maybe this also explains the divergence in bigint to float conversion.
Maybe they have by accident a misconfigured FPU? Some init call
to the runtime system missing? Or a cheap floating point library?

But one is WSL2 and the other Windows!

Mostowski Collapse schrieb am Donnerstag, 13. Oktober 2022 um 20:16:02 UTC+2:
> Does the C floating point API have HALF_EVEN default?
> Maybe, maybe not. Is there a test case to figure out
> what rounding a system, Prolog or not, does, for its
>
> floating point numbers, independent of some bigint
> conversion issues? Was trying the folllowing test and got
> some surprise. There was first quite some agreement:
>
> /* SWI-Prolog 8.5.18 */
> ?- X is 665723957508124892951320985600.0+26388279066624.0.
> X = 6.65723957508125e+29.
>
> /* Ciao Prolog 1.22.0 */
> ?- X is 665723957508124892951320985600.0+26388279066624.0.
> X = 6.65723957508125e29 ?
>
> /* Trealla Prolog 2.4.3 */
> ?- X is 665723957508124892951320985600.0+26388279066624.0.
> X = 6.65723957508125e+29.
>
> /* Jekejeke Prolog 1.5.5 */
> ?- X is 665723957508124892951320985600.0+26388279066624.0.
> X = 6.65723957508125E29.

Have to dig deeper, could be a parsing problem.
For the first summand, I get:

/* SWI-Prolog */
?- X is 665723957508124892951320985600.0.
X = 6.65723957508125e+29.

/* ECLiPSe Prolog */
X is 665723957508124892951320985600.0.
X = 6.6572395750812482e+29

I am not yet at a test case, which would give me
some information about rounding.

Mostowski Collapse schrieb am Donnerstag, 13. Oktober 2022 um 20:18:04 UTC+2:
> But then suddently:
>
> /* Scryer Prolog 0.9.0 */
> ?- X is 665723957508124892951320985600.0+26388279066624.0.
> X = 6.657239575081248e29.
>
> /* ECLiPSe Prolog 7.0.61 */
> ?- X is 665723957508124892951320985600.0+26388279066624.0.
> X = 6.6572395750812482e+29
>
> ?- X is 6.6572395750812482e+29-6.65723957508125E29.
> X = -140737488355328.0.
>
> Maybe this also explains the divergence in bigint to float conversion.
> Maybe they have by accident a misconfigured FPU? Some init call
> to the runtime system missing? Or a cheap floating point library?
>
> But one is WSL2 and the other Windows!
> Mostowski Collapse schrieb am Donnerstag, 13. Oktober 2022 um 20:16:02 UTC+2:
> > Does the C floating point API have HALF_EVEN default?
> > Maybe, maybe not. Is there a test case to figure out
> > what rounding a system, Prolog or not, does, for its
> >
> > floating point numbers, independent of some bigint
> > conversion issues? Was trying the folllowing test and got
> > some surprise. There was first quite some agreement:
> >
> > /* SWI-Prolog 8.5.18 */
> > ?- X is 665723957508124892951320985600.0+26388279066624.0.
> > X = 6.65723957508125e+29.
> >
> > /* Ciao Prolog 1.22.0 */
> > ?- X is 665723957508124892951320985600.0+26388279066624.0.
> > X = 6.65723957508125e29 ?
> >
> > /* Trealla Prolog 2.4.3 */
> > ?- X is 665723957508124892951320985600.0+26388279066624.0.
> > X = 6.65723957508125e+29.
> >
> > /* Jekejeke Prolog 1.5.5 */
> > ?- X is 665723957508124892951320985600.0+26388279066624.0.
> > X = 6.65723957508125E29.

Looks like my system belongs to the 2 system that
cannot generate good random numbers? I guess the
system with a smaller value is possible the less good

random generator, because some periodicity was hit.
Java uses a Lehmer algorithm which has a very short
period and uses very few bits in its state. On the other

hand the Twister uses quite some bits I guess, and has
a much larger periodicity. Twister would be therefore the
better randum number generator.

Mostowski Collapse schrieb am Donnerstag, 13. Oktober 2022 um 20:36:51 UTC+2:
> Have to dig deeper, could be a parsing problem.
> For the first summand, I get:
>
> /* SWI-Prolog */
> ?- X is 665723957508124892951320985600.0.
> X = 6.65723957508125e+29.
>
> /* ECLiPSe Prolog */
> X is 665723957508124892951320985600.0.
> X = 6.6572395750812482e+29
>
> I am not yet at a test case, which would give me
> some information about rounding.
> Mostowski Collapse schrieb am Donnerstag, 13. Oktober 2022 um 20:18:04 UTC+2:
> > But then suddently:
> >
> > /* Scryer Prolog 0.9.0 */
> > ?- X is 665723957508124892951320985600.0+26388279066624.0.
> > X = 6.657239575081248e29.
> >
> > /* ECLiPSe Prolog 7.0.61 */
> > ?- X is 665723957508124892951320985600.0+26388279066624.0.
> > X = 6.6572395750812482e+29
> >
> > ?- X is 6.6572395750812482e+29-6.65723957508125E29.
> > X = -140737488355328.0.
> >
> > Maybe this also explains the divergence in bigint to float conversion.
> > Maybe they have by accident a misconfigured FPU? Some init call
> > to the runtime system missing? Or a cheap floating point library?
> >
> > But one is WSL2 and the other Windows!
> > Mostowski Collapse schrieb am Donnerstag, 13. Oktober 2022 um 20:16:02 UTC+2:
> > > Does the C floating point API have HALF_EVEN default?
> > > Maybe, maybe not. Is there a test case to figure out
> > > what rounding a system, Prolog or not, does, for its
> > >
> > > floating point numbers, independent of some bigint
> > > conversion issues? Was trying the folllowing test and got
> > > some surprise. There was first quite some agreement:
> > >
> > > /* SWI-Prolog 8.5.18 */
> > > ?- X is 665723957508124892951320985600.0+26388279066624.0.
> > > X = 6.65723957508125e+29.
> > >
> > > /* Ciao Prolog 1.22.0 */
> > > ?- X is 665723957508124892951320985600.0+26388279066624.0.
> > > X = 6.65723957508125e29 ?
> > >
> > > /* Trealla Prolog 2.4.3 */
> > > ?- X is 665723957508124892951320985600.0+26388279066624.0.
> > > X = 6.65723957508125e+29.
> > >
> > > /* Jekejeke Prolog 1.5.5 */
> > > ?- X is 665723957508124892951320985600.0+26388279066624.0.
> > > X = 6.65723957508125E29.

Now turning my attention to (/)/2. SWI-Prolog wants me
to use set_prolog_flag(iso, true), to get the semantics
of the ISO core standard, page 116.

But the fuzzer still finds the needle in the haystack.
Using more bits to the numerator, and I find:

/* Jekejeke Prolog 1.5.5 */
?- X is -15915286457581575168 / -238455449856.
X = 6.674322800000001E7.

/* SWI-Prolog 8.5.18 */
?- current_prolog_flag(iso, X).
X = true.

?- X is -15915286457581575168 / -238455449856.
X = 66743227.99999999.

Now I have test cases so nasty, even Ciao Prolog and
SICStus Prolog are stumbling. Interestingly Ciao Prolog
and SICStus Prolog stumble at the exact same test cases.

/* Ciao Prolog 1.22.0 */
?- case4(N, P, Q, Y), Z is P/Q, Y =\= Z.
N = 3,
P = -138732791312937296168,
Q = 70775453,
Y = -1960182315087.9346,
Z = -1960182315087.9343 ? ;
no

/* SICStus Prolog 4.7.1 */
?- case4(N, P, Q, Y), Z is P/Q, Y =\= Z.
N = 3,
P = -138732791312937296168,
Q = 70775453,
Y = -1960182315087.9346,
Z = -1960182315087.9343 ? ;
no

Mostowski Collapse schrieb am Freitag, 14. Oktober 2022 um 10:41:45 UTC+2:
> Now turning my attention to (/)/2. SWI-Prolog wants me
> to use set_prolog_flag(iso, true), to get the semantics
> of the ISO core standard, page 116.
>
> But the fuzzer still finds the needle in the haystack.
> Using more bits to the numerator, and I find:
> /* Jekejeke Prolog 1.5.5 */
> ?- X is -15915286457581575168 / -238455449856.
> X = 6.674322800000001E7.
>
> /* SWI-Prolog 8.5.18 */
> ?- current_prolog_flag(iso, X).
> X = true.
>
> ?- X is -15915286457581575168 / -238455449856.
> X = 66743227.99999999.

Last chapter of the testing campaign. Testing (**)/2 now.

BTW, this is a little cringe what the fuzzer found for Dogelog
Player on Python. One can try Python 3.11 and will get:

Python 3.11.0rc1 (main, Aug 8 2022, 11:30:54)
>>> 396** -1
0.002525252525252525
>>> 1/396
0.0025252525252525255

Same test case that makes SICStus Prolog tumbling. Such
surprises lower the pain of testing, make it rather addictive.

Mostowski Collapse schrieb am Freitag, 14. Oktober 2022 um 10:43:32 UTC+2:
> Now I have test cases so nasty, even Ciao Prolog and
> SICStus Prolog are stumbling. Interestingly Ciao Prolog
> and SICStus Prolog stumble at the exact same test cases.
> /* Ciao Prolog 1.22.0 */
> ?- case4(N, P, Q, Y), Z is P/Q, Y =\= Z.
> N = 3,
> P = -138732791312937296168,
> Q = 70775453,
> Y = -1960182315087.9346,
> Z = -1960182315087.9343 ? ;
> no
>
> /* SICStus Prolog 4.7.1 */
> ?- case4(N, P, Q, Y), Z is P/Q, Y =\= Z.
> N = 3,
> P = -138732791312937296168,
> Q = 70775453,
> Y = -1960182315087.9346,
> Z = -1960182315087.9343 ? ;
> no
> Mostowski Collapse schrieb am Freitag, 14. Oktober 2022 um 10:41:45 UTC+2:
> > Now turning my attention to (/)/2. SWI-Prolog wants me
> > to use set_prolog_flag(iso, true), to get the semantics
> > of the ISO core standard, page 116.
> >
> > But the fuzzer still finds the needle in the haystack.
> > Using more bits to the numerator, and I find:
> > /* Jekejeke Prolog 1.5.5 */
> > ?- X is -15915286457581575168 / -238455449856.
> > X = 6.674322800000001E7.
> >
> > /* SWI-Prolog 8.5.18 */
> > ?- current_prolog_flag(iso, X).
> > X = true.
> >
> > ?- X is -15915286457581575168 / -238455449856.
> > X = 66743227.99999999.

Next challenge, how to keep my comparison matrice
small. There are so many JavaScript platforms, nodeJS,
browsers, and WTF they differ:

/* Chrome Windows Nok */
?- X is 1889**29.
X = 1.0250068783051206E95.

/* FireFox Windows Ok */
?- X is 1889**29.
X = 1.0250068783051207E95.

Whats the Novacore bigint versus float mantra? The difference
between the two notions “less precise and more speed” versus
“more precise and less speed” can be maybe capture as follows:

- Less Precise and More Speed: Even if the arguments are bigint, and
have potentially more than 52 bit mantissa, they are always first
converted to float. So that the operation is only performed with 52
bit arguments, allow for faster processing.

- More Precise and Less Speed: If the arguments are bigint or if
the arguments are rational numbers or arbitrary floats, all bits are
taken into account. So that the operation is performed with more bits
than only 52 bits, resulting in slower processing.

Theoretically in the intersection, when the arguments are already
52 bit mantissa, there should be no difference in speed and also
no difference in accuracy, except if the second mode is

also required to return more bits. So we also have for the
first mode the result is **always** required to be only a float.
Then there is mostlikely always a speed difference, and the

first mode is always faster and less accurate.

Mostowski Collapse schrieb am Mittwoch, 19. Oktober 2022 um 16:42:37 UTC+2:
> Whats the Novacore bigint versus float mantra? The difference
> between the two notions “less precise and more speed” versus
> “more precise and less speed” can be maybe capture as follows:
>
> - Less Precise and More Speed: Even if the arguments are bigint, and
> have potentially more than 52 bit mantissa, they are always first
> converted to float. So that the operation is only performed with 52
> bit arguments, allow for faster processing.
>
> - More Precise and Less Speed: If the arguments are bigint or if
> the arguments are rational numbers or arbitrary floats, all bits are
> taken into account. So that the operation is performed with more bits
> than only 52 bits, resulting in slower processing.

Ok did some performance measurement. Basically used this harness.
Was reusing the test data, that was used for accuracy testing, so its
based on the same cases.p file:

% swi_perf
swi_perf :-
write('case, swi: '), time((between(1,20000,_),
case(_, X, Y), Y =\= float(X), fail; true)),
write('case2, swi: '), time((between(1,20000,_),
case2(_, X, Y), Y =\= float(X), fail; true)),
write('case3, swi: '), time((between(1,20000,_),
case3(_, X, Y), Y =\= float(X), fail; true)),
write('case4, swi: '), time((between(1,20000,_),
case4(_, P, Q, Y), Y =\= P/Q, fail; true)),
write('case5, swi: '), time((between(1,20000,_),
case5(_, P, Q, Y), Y =\= P/Q, fail; true)),
write('case6, swi: '), time((between(1,20000,_),
case6(_, P, Q, Y), Y =\= P**Q, fail; true)).

Mostowski Collapse schrieb am Mittwoch, 19. Oktober 2022 um 16:42:37 UTC+2:
> Whats the Novacore bigint versus float mantra? The difference
> between the two notions “less precise and more speed” versus
> “more precise and less speed” can be maybe capture as follows:
>
> - Less Precise and More Speed: Even if the arguments are bigint, and
> have potentially more than 52 bit mantissa, they are always first
> converted to float. So that the operation is only performed with 52
> bit arguments, allow for faster processing.
>
> - More Precise and Less Speed: If the arguments are bigint or if
> the arguments are rational numbers or arbitrary floats, all bits are
> taken into account. So that the operation is performed with more bits
> than only 52 bits, resulting in slower processing.

Woa! My Jekejeke Prolog (was using JDK 1.8) is quite a number cruncher,
on the other hand Dogelog Player (was using nodeJS) is lacking behind.
Inbetween the result of Prolog systems that I tested, results in milliseconds.

System case case2 case3 case4 case5 case6
jekejeke 270 266 255 305 361 539
ciao 344 343 343 386 386 544
eclipse 414 388 396 518 541 939
trealla 553 547 549 569 565 605
swi 989 941 979 528 532 861
scryer 889 857 872 1177 1188 1147
sicstus 1087 1083 1104 1221 1218 1226
dogelog 1171 1163 1171 1241 1300 1330

Mostowski Collapse schrieb am Mittwoch, 19. Oktober 2022 um 16:51:45 UTC+2:
> Ok did some performance measurement. Basically used this harness.
> Was reusing the test data, that was used for accuracy testing, so its
> based on the same cases.p file:
>
> % swi_perf
> swi_perf :-
> write('case, swi: '), time((between(1,20000,_),
> case(_, X, Y), Y =\= float(X), fail; true)),
> write('case2, swi: '), time((between(1,20000,_),
> case2(_, X, Y), Y =\= float(X), fail; true)),
> write('case3, swi: '), time((between(1,20000,_),
> case3(_, X, Y), Y =\= float(X), fail; true)),
> write('case4, swi: '), time((between(1,20000,_),
> case4(_, P, Q, Y), Y =\= P/Q, fail; true)),
> write('case5, swi: '), time((between(1,20000,_),
> case5(_, P, Q, Y), Y =\= P/Q, fail; true)),
> write('case6, swi: '), time((between(1,20000,_),
> case6(_, P, Q, Y), Y =\= P**Q, fail; true)).
> Mostowski Collapse schrieb am Mittwoch, 19. Oktober 2022 um 16:42:37 UTC+2:
> > Whats the Novacore bigint versus float mantra? The difference
> > between the two notions “less precise and more speed” versus
> > “more precise and less speed” can be maybe capture as follows:
> >
> > - Less Precise and More Speed: Even if the arguments are bigint, and
> > have potentially more than 52 bit mantissa, they are always first
> > converted to float. So that the operation is only performed with 52
> > bit arguments, allow for faster processing.
> >
> > - More Precise and Less Speed: If the arguments are bigint or if
> > the arguments are rational numbers or arbitrary floats, all bits are
> > taken into account. So that the operation is performed with more bits
> > than only 52 bits, resulting in slower processing.

Code inspection sometimes helps find defects:

/* Jekejeke Prolog Ok */
?- X is (2^53+1)*2^64+1, Y is float(X).
X = 166153499473114502559719956244594689, Y = 1.6615349947311452E35.

/* SWI-Prolog 8.5.20 Ok */
?- X is float(166153499473114502559719956244594689).
X = 1.6615349947311452e+35.

/* Ciao Prolog Nok */
?- X is float(166153499473114502559719956244594689).
X = 1.661534994731145e35 ?

Pitty my fuzzer didn't find it. But I guess I will manually add it.

Mostowski Collapse schrieb am Mittwoch, 19. Oktober 2022 um 16:52:41 UTC+2:
> Woa! My Jekejeke Prolog (was using JDK 1.8) is quite a number cruncher,
> on the other hand Dogelog Player (was using nodeJS) is lacking behind.
> Inbetween the result of Prolog systems that I tested, results in milliseconds.
>
> System case case2 case3 case4 case5 case6
> jekejeke 270 266 255 305 361 539
> ciao 344 343 343 386 386 544
> eclipse 414 388 396 518 541 939
> trealla 553 547 549 569 565 605
> swi 989 941 979 528 532 861
> scryer 889 857 872 1177 1188 1147
> sicstus 1087 1083 1104 1221 1218 1226
> dogelog 1171 1163 1171 1241 1300 1330
> Mostowski Collapse schrieb am Mittwoch, 19. Oktober 2022 um 16:51:45 UTC+2:
> > Ok did some performance measurement. Basically used this harness.
> > Was reusing the test data, that was used for accuracy testing, so its
> > based on the same cases.p file:
> >
> > % swi_perf
> > swi_perf :-
> > write('case, swi: '), time((between(1,20000,_),
> > case(_, X, Y), Y =\= float(X), fail; true)),
> > write('case2, swi: '), time((between(1,20000,_),
> > case2(_, X, Y), Y =\= float(X), fail; true)),
> > write('case3, swi: '), time((between(1,20000,_),
> > case3(_, X, Y), Y =\= float(X), fail; true)),
> > write('case4, swi: '), time((between(1,20000,_),
> > case4(_, P, Q, Y), Y =\= P/Q, fail; true)),
> > write('case5, swi: '), time((between(1,20000,_),
> > case5(_, P, Q, Y), Y =\= P/Q, fail; true)),
> > write('case6, swi: '), time((between(1,20000,_),
> > case6(_, P, Q, Y), Y =\= P**Q, fail; true)).
> > Mostowski Collapse schrieb am Mittwoch, 19. Oktober 2022 um 16:42:37 UTC+2:
> > > Whats the Novacore bigint versus float mantra? The difference
> > > between the two notions “less precise and more speed” versus
> > > “more precise and less speed” can be maybe capture as follows:
> > >
> > > - Less Precise and More Speed: Even if the arguments are bigint, and
> > > have potentially more than 52 bit mantissa, they are always first
> > > converted to float. So that the operation is only performed with 52
> > > bit arguments, allow for faster processing.
> > >
> > > - More Precise and Less Speed: If the arguments are bigint or if
> > > the arguments are rational numbers or arbitrary floats, all bits are
> > > taken into account. So that the operation is performed with more bits
> > > than only 52 bits, resulting in slower processing.

Explanation of the defect:

Mostly based on code from Torbjorn Granlund, Johan Andersson,
and Mats Carlsson (Ref: Knuth vol. 2 sec. 4.3.1)

flt64_t bn_to_float(bignum_t *bn) {
while (i > 1) {
f = f*norm2 + u;
https://github.com/ciao -lang/ciao/blob/master/core/engine/eng_bignum.c#L1043

I think a defect could be if HALF_EVEN prematurely
rounds down the value of f, and then later a value of u,
would indicate that this rounding down was not legit,

but float addition cannot help it.

Mostowski Collapse schrieb am Donnerstag, 20. Oktober 2022 um 03:55:47 UTC+2:
> Code inspection sometimes helps find defects:
>
> /* Jekejeke Prolog Ok */
> ?- X is (2^53+1)*2^64+1, Y is float(X).
> X = 166153499473114502559719956244594689, Y = 1.6615349947311452E35.
>
> /* SWI-Prolog 8.5.20 Ok */
> ?- X is float(166153499473114502559719956244594689).
> X = 1.6615349947311452e+35.
>
> /* Ciao Prolog Nok */
> ?- X is float(166153499473114502559719956244594689).
> X = 1.661534994731145e35 ?
>
> Pitty my fuzzer didn't find it. But I guess I will manually add it.
> Mostowski Collapse schrieb am Mittwoch, 19. Oktober 2022 um 16:52:41 UTC+2:
> > Woa! My Jekejeke Prolog (was using JDK 1.8) is quite a number cruncher,
> > on the other hand Dogelog Player (was using nodeJS) is lacking behind.
> > Inbetween the result of Prolog systems that I tested, results in milliseconds.
> >
> > System case case2 case3 case4 case5 case6
> > jekejeke 270 266 255 305 361 539
> > ciao 344 343 343 386 386 544
> > eclipse 414 388 396 518 541 939
> > trealla 553 547 549 569 565 605
> > swi 989 941 979 528 532 861
> > scryer 889 857 872 1177 1188 1147
> > sicstus 1087 1083 1104 1221 1218 1226
> > dogelog 1171 1163 1171 1241 1300 1330
> > Mostowski Collapse schrieb am Mittwoch, 19. Oktober 2022 um 16:51:45 UTC+2:
> > > Ok did some performance measurement. Basically used this harness.
> > > Was reusing the test data, that was used for accuracy testing, so its
> > > based on the same cases.p file:
> > >
> > > % swi_perf
> > > swi_perf :-
> > > write('case, swi: '), time((between(1,20000,_),
> > > case(_, X, Y), Y =\= float(X), fail; true)),
> > > write('case2, swi: '), time((between(1,20000,_),
> > > case2(_, X, Y), Y =\= float(X), fail; true)),
> > > write('case3, swi: '), time((between(1,20000,_),
> > > case3(_, X, Y), Y =\= float(X), fail; true)),
> > > write('case4, swi: '), time((between(1,20000,_),
> > > case4(_, P, Q, Y), Y =\= P/Q, fail; true)),
> > > write('case5, swi: '), time((between(1,20000,_),
> > > case5(_, P, Q, Y), Y =\= P/Q, fail; true)),
> > > write('case6, swi: '), time((between(1,20000,_),
> > > case6(_, P, Q, Y), Y =\= P**Q, fail; true)).
> > > Mostowski Collapse schrieb am Mittwoch, 19. Oktober 2022 um 16:42:37 UTC+2:
> > > > Whats the Novacore bigint versus float mantra? The difference
> > > > between the two notions “less precise and more speed” versus
> > > > “more precise and less speed” can be maybe capture as follows:
> > > >
> > > > - Less Precise and More Speed: Even if the arguments are bigint, and
> > > > have potentially more than 52 bit mantissa, they are always first
> > > > converted to float. So that the operation is only performed with 52
> > > > bit arguments, allow for faster processing.
> > > >
> > > > - More Precise and Less Speed: If the arguments are bigint or if
> > > > the arguments are rational numbers or arbitrary floats, all bits are
> > > > taken into account. So that the operation is performed with more bits
> > > > than only 52 bits, resulting in slower processing.

I added some hand picked test cases, from the
Ciao Prolog blooper, including some negative bigints.
Interestingly other Prolog systems are also pray to the

same blooper. In total I found that 3 Prolog systems have
this blooper, whereas Trealla only for 50%. Here are now
the number of failures from the new test cases:

% case3, trealla: 5
% case3, ciao: 10
% case3, sicstus: 10

Mostowski Collapse schrieb am Donnerstag, 20. Oktober 2022 um 03:58:38 UTC+2:
> Explanation of the defect:
>
> Mostly based on code from Torbjorn Granlund, Johan Andersson,
> and Mats Carlsson (Ref: Knuth vol. 2 sec. 4.3.1)
>
> flt64_t bn_to_float(bignum_t *bn) {
> while (i > 1) {
> f = f*norm2 + u;
> https://github.com/ciao -lang/ciao/blob/master/core/engine/eng_bignum.c#L1043
>
> I think a defect could be if HALF_EVEN prematurely
> rounds down the value of f, and then later a value of u,
> would indicate that this rounding down was not legit,
>
> but float addition cannot help it.
> Mostowski Collapse schrieb am Donnerstag, 20. Oktober 2022 um 03:55:47 UTC+2:
> > Code inspection sometimes helps find defects:
> >
> > /* Jekejeke Prolog Ok */
> > ?- X is (2^53+1)*2^64+1, Y is float(X).
> > X = 166153499473114502559719956244594689, Y = 1.6615349947311452E35..
> >
> > /* SWI-Prolog 8.5.20 Ok */
> > ?- X is float(166153499473114502559719956244594689).
> > X = 1.6615349947311452e+35.
> >
> > /* Ciao Prolog Nok */
> > ?- X is float(166153499473114502559719956244594689).
> > X = 1.661534994731145e35 ?
> >
> > Pitty my fuzzer didn't find it. But I guess I will manually add it.
> > Mostowski Collapse schrieb am Mittwoch, 19. Oktober 2022 um 16:52:41 UTC+2:
> > > Woa! My Jekejeke Prolog (was using JDK 1.8) is quite a number cruncher,
> > > on the other hand Dogelog Player (was using nodeJS) is lacking behind..
> > > Inbetween the result of Prolog systems that I tested, results in milliseconds.
> > >
> > > System case case2 case3 case4 case5 case6
> > > jekejeke 270 266 255 305 361 539
> > > ciao 344 343 343 386 386 544
> > > eclipse 414 388 396 518 541 939
> > > trealla 553 547 549 569 565 605
> > > swi 989 941 979 528 532 861
> > > scryer 889 857 872 1177 1188 1147
> > > sicstus 1087 1083 1104 1221 1218 1226
> > > dogelog 1171 1163 1171 1241 1300 1330
> > > Mostowski Collapse schrieb am Mittwoch, 19. Oktober 2022 um 16:51:45 UTC+2:
> > > > Ok did some performance measurement. Basically used this harness.
> > > > Was reusing the test data, that was used for accuracy testing, so its
> > > > based on the same cases.p file:
> > > >
> > > > % swi_perf
> > > > swi_perf :-
> > > > write('case, swi: '), time((between(1,20000,_),
> > > > case(_, X, Y), Y =\= float(X), fail; true)),
> > > > write('case2, swi: '), time((between(1,20000,_),
> > > > case2(_, X, Y), Y =\= float(X), fail; true)),
> > > > write('case3, swi: '), time((between(1,20000,_),
> > > > case3(_, X, Y), Y =\= float(X), fail; true)),
> > > > write('case4, swi: '), time((between(1,20000,_),
> > > > case4(_, P, Q, Y), Y =\= P/Q, fail; true)),
> > > > write('case5, swi: '), time((between(1,20000,_),
> > > > case5(_, P, Q, Y), Y =\= P/Q, fail; true)),
> > > > write('case6, swi: '), time((between(1,20000,_),
> > > > case6(_, P, Q, Y), Y =\= P**Q, fail; true)).
> > > > Mostowski Collapse schrieb am Mittwoch, 19. Oktober 2022 um 16:42:37 UTC+2:
> > > > > Whats the Novacore bigint versus float mantra? The difference
> > > > > between the two notions “less precise and more speed” versus
> > > > > “more precise and less speed” can be maybe capture as follows:
> > > > >
> > > > > - Less Precise and More Speed: Even if the arguments are bigint, and
> > > > > have potentially more than 52 bit mantissa, they are always first
> > > > > converted to float. So that the operation is only performed with 52
> > > > > bit arguments, allow for faster processing.
> > > > >
> > > > > - More Precise and Less Speed: If the arguments are bigint or if
> > > > > the arguments are rational numbers or arbitrary floats, all bits are
> > > > > taken into account. So that the operation is performed with more bits
> > > > > than only 52 bits, resulting in slower processing.

Woa! My testing campaign has not yet ended. Now I
found a test case, where the Unix platform also fails.
Was increasing the bits in the case6/3 fuzzer:

pow2(B, E, X) :- E < 0, integer(E), !, X is 1/(B^(-E)).
pow2(B, E, X) :- X is B^E.

?- current_prolog_flag(iso, X).
X = false.

/* SWI-Prolog 8.5.20 WLS2 Nok */
?- repeat, L is -(1<<15), H is (1<<15)+1,
M is -(1<<3), J is (1<<3)+1, random(L, H, P),
random(M, J, Q), pow2(P, Q, A), B is P**Q, A=\=B.
P = -18573,
Q = -7,
A = -1.3116730299820406e-30,
B = -1.3116730299820408e-30 .

I have SWIPL on Unix platform only since weekend,
could not do this testing beforehand, so I guess I
need to update case6/3 test cases,

because they showed Unix platform completely
passing, which seems to be wrong. I used in the old
fuzzer 12 bits + 6 bits, the above uses 15 bits + 3 bits,

and new nasty test cases pop up. Its a test case where
JDK 19 non-strict also fails. Cool! I couldn’t believe that
it is that good. It was too good to be true.

So its also not that good!

/* Jekejeke Prolog 1.5.5, JDK 19 non-strict Nok */
?- X is -18573** -7.
X = -1.3116730299820408E-30.

Mostowski Collapse schrieb am Dienstag, 25. Oktober 2022 um 14:56:11 UTC+2:
> Woa! My testing campaign has not yet ended. Now I
> found a test case, where the Unix platform also fails.
> Was increasing the bits in the case6/3 fuzzer:
>
> pow2(B, E, X) :- E < 0, integer(E), !, X is 1/(B^(-E)).
> pow2(B, E, X) :- X is B^E.
>
> ?- current_prolog_flag(iso, X).
> X = false.
>
> /* SWI-Prolog 8.5.20 WLS2 Nok */
> ?- repeat, L is -(1<<15), H is (1<<15)+1,
> M is -(1<<3), J is (1<<3)+1, random(L, H, P),
> random(M, J, Q), pow2(P, Q, A), B is P**Q, A=\=B.
> P = -18573,
> Q = -7,
> A = -1.3116730299820406e-30,
> B = -1.3116730299820408e-30 .

Nope the JDK 19 is correct here. Bug in long division?

/* SWI-Prolog 8.5.20 WSL2 */
?- current_prolog_flag(iso, X).
X = false.

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

?- X is float(1 rdiv (18573^7)).
X = 1.3116730299820408e-30.

the Mostowski Collapse schrieb am Dienstag, 25. Oktober 2022 um 14:59:11 UTC+2:
> I have SWIPL on Unix platform only since weekend,
> could not do this testing beforehand, so I guess I
> need to update case6/3 test cases,
>
> because they showed Unix platform completely
> passing, which seems to be wrong. I used in the old
> fuzzer 12 bits + 6 bits, the above uses 15 bits + 3 bits,
>
> and new nasty test cases pop up. Its a test case where
> JDK 19 non-strict also fails. Cool! I couldn’t believe that
> it is that good. It was too good to be true.
>
> So its also not that good!
>
> /* Jekejeke Prolog 1.5.5, JDK 19 non-strict Nok */
> ?- X is -18573** -7.
> X = -1.3116730299820408E-30.
> Mostowski Collapse schrieb am Dienstag, 25. Oktober 2022 um 14:56:11 UTC+2:
> > Woa! My testing campaign has not yet ended. Now I
> > found a test case, where the Unix platform also fails.
> > Was increasing the bits in the case6/3 fuzzer:
> >
> > pow2(B, E, X) :- E < 0, integer(E), !, X is 1/(B^(-E)).
> > pow2(B, E, X) :- X is B^E.
> >
> > ?- current_prolog_flag(iso, X).
> > X = false.
> >
> > /* SWI-Prolog 8.5.20 WLS2 Nok */
> > ?- repeat, L is -(1<<15), H is (1<<15)+1,
> > M is -(1<<3), J is (1<<3)+1, random(L, H, P),
> > random(M, J, Q), pow2(P, Q, A), B is P**Q, A=\=B.
> > P = -18573,
> > Q = -7,
> > A = -1.3116730299820406e-30,
> > B = -1.3116730299820408e-30 .

Enhanced test cases, new viewpoint. Now WASM is
the more accurate platform, both for SWI and Ciao.
JDK 19 is still ultra accurate. But I cant say that it

is always correctly rounded. For non-integer arguments
I found a faulty example for (**)/2, but the below test
cases only test integer arguments of (**)/2:

System (=:=)/2 (/)/2 (**)/2 Total Variant
jekejeke 0 0 0 0 JDK 19
dogelog 0 0 2 2 PyPy
swi 0 0 4 4 WASM
sicstus 10 0 2 12 Windows
trealla 5 0 11 16 WSL2
ciao 10 0 10 20 WASM
scryer 33 46 11 90 WSL2
eclipse 33 45 68 146 Windows

Mostowski Collapse schrieb am Dienstag, 25. Oktober 2022 um 17:36:46 UTC+2:
> Nope the JDK 19 is correct here. Bug in long division?
>
> /* SWI-Prolog 8.5.20 WSL2 */
> ?- current_prolog_flag(iso, X).
> X = false.
> ?- X is 1/(18573^7).
> X = 1.3116730299820406e-30.
>
> ?- X is float(1 rdiv (18573^7)).
> X = 1.3116730299820408e-30.
> the Mostowski Collapse schrieb am Dienstag, 25. Oktober 2022 um 14:59:11 UTC+2:
> > I have SWIPL on Unix platform only since weekend,
> > could not do this testing beforehand, so I guess I
> > need to update case6/3 test cases,
> >
> > because they showed Unix platform completely
> > passing, which seems to be wrong. I used in the old
> > fuzzer 12 bits + 6 bits, the above uses 15 bits + 3 bits,
> >
> > and new nasty test cases pop up. Its a test case where
> > JDK 19 non-strict also fails. Cool! I couldn’t believe that
> > it is that good. It was too good to be true.
> >
> > So its also not that good!
> >
> > /* Jekejeke Prolog 1.5.5, JDK 19 non-strict Nok */
> > ?- X is -18573** -7.
> > X = -1.3116730299820408E-30.
> > Mostowski Collapse schrieb am Dienstag, 25. Oktober 2022 um 14:56:11 UTC+2:
> > > Woa! My testing campaign has not yet ended. Now I
> > > found a test case, where the Unix platform also fails.
> > > Was increasing the bits in the case6/3 fuzzer:
> > >
> > > pow2(B, E, X) :- E < 0, integer(E), !, X is 1/(B^(-E)).
> > > pow2(B, E, X) :- X is B^E.
> > >
> > > ?- current_prolog_flag(iso, X).
> > > X = false.
> > >
> > > /* SWI-Prolog 8.5.20 WLS2 Nok */
> > > ?- repeat, L is -(1<<15), H is (1<<15)+1,
> > > M is -(1<<3), J is (1<<3)+1, random(L, H, P),
> > > random(M, J, Q), pow2(P, Q, A), B is P**Q, A=\=B.
> > > P = -18573,
> > > Q = -7,
> > > A = -1.3116730299820406e-30,
> > > B = -1.3116730299820408e-30 .

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.