Rounding error in (**)/2

Its really a combination of two problems. float_half_even/2
doesn’t help that much. I could make it run via the mpz_get_d/2
simulation here. And then I can measure:

?- aggregate_all(sum(S), (between(1,1000000,N),
    S is abs(N^10-float(N)**10)), T).
T = 9.901018060037038e+48.

?- aggregate_all(sum(S), (between(1,1000000,N),
    A is N^10, float_half_even(A,B), S is abs(B-float(N)**10)), T).
T = 7.871954492015358e+48.

Whereas JDK 19 gives, aggregate error is 1000 times smaller:

?- aggregate_all(sum(S), (between(1,1000000,N),
  S is abs(N^10-float(N)**10)), T).
T = 7.00914689466423E45.

So on the windows platform it must be some low quality math.h
library. But I am not sure whether ar_pow() even uses some math.h
library. One value that has a problem is for example this here:

/* SWI-Prolog 8.5.17 */
?- X is 978400.0**10.0.
X = 8.038304251093646e+59.

But if it were to use math.h, it would be:

#include <math.h>
printf("X=%.15e", pow(978400.0, 10.0));
X=8.038304251093644e+59

https://www.onlinegdb.com/online_c_compiler

I can use this to validate my float_half_even/2, seems to be ok:

?- X is 978400^10, Y is float(X).
X = 803830425109364419968771652984040429977600000000000000000000,
Y = 8.038304251093643e+59.

?- X is 978400^10, float_half_even(X,Y).
X = 803830425109364419968771652984040429977600000000000000000000,
Y = 8.038304251093644e+59.

So the value is already a next after from 43 to 44. How the heck
does it even next after one more time from 44 to 46?