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.
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 Accuracy of Mathematical Functions in Single, Double, Double Extended, and Quadruple Precision - Inria - Institut national de recherche en sciences et technologies du numérique - 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.
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 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.
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 .
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.
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.
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.
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.
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:
current_prolog_flag(prefer_rationals, true) holds
current_prolog_flag(iso, true) holds
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.
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.
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
It’s been almost 4 years since this was last discussed, so I thought it was time to investigate the current state of play. To be clear, in all that time I’ve not come across any issues with the current strategy of explicit outward rounding of “TO_NEAREST” results from the various elementary function libraries, so I don’t think this is a priority issue. But It feels like a loose end that would be nice to tie up at some point.
Inria’s CORE-MATH appears to be the recognized “standard” approach and their source release now provide correct (in all rounding modes) C source for 64 bit float (and others) elementary functions (exp, log, sin, …) under the MIT license conditions. The intent is that these will eventually be incorporated into the various math libraries supported by the various organizations that build such things. But it’s unclear how quickly this is happening.
There is at least one Git repo (courtesy Bob Burger) that has bundled the 64 bit CORE-MATH code with a build system for generating a “crmath” static library and header file for general use on various platforms. I used this to build a MacOS Intel version of crmath.
With copious assistance from AI (Google Gemini), I spent a couple of days to build a version of SWI-Prolog that includes this crmath library and modified pl-arith.c to use it instead of the existing explicit rounding versions. And all this appears to work, producing slightly tighter intervals (in clpBNR) and comes with Inria’s guarantee of correctness.
The only problem with my exploratory work is that the performance is an order (or two) magnitude worse than the existing implementation. But Inria’s CORE-MATH performance numbers (see link above) indicate they should be comparable so I think I’m missing something (compiler options?) in my build environment that results in such a drastic loss of performance.
The other major issue, if it’s decided to pursue this further, is how to properly incorporate the CORE-MATH functions properly into the SWIP build process. I’m pretty sure my hack using Bob Burger’s repo isn’t how it should be done.
Comments/suggestions welcome.
P.S. I must commend Google Gemini on its patience; cmake is my personal programming nightmare.
That is not so clear. The best way to integrate is probably to add Bob Burger’s repo as a submodule and build the library as part of the build process. I guess the changes to pl-arith.c are limited to #define sin(x) cr_sin(x), etc. disabling the current safety margin.
All that is fairly easy. But, adding the source means 2.2Mb extra and the binary on AMD64 is 670Kbytes. That is quite a high price to pay (libswipl.so is only 2Mb). In particular because glibc already seems to be using these functions. If we can detect that we can disable our safety rounding for Linux 
In theory, probing the glibc version is enough.
I guess the hope that Microsoft and Apple update their C libraries any time soon is idle.
I don’t think adding 670K to the distributed binaries for MacOS and Windows can be justified, no?
Best effort is probably a CMake option -DCRMATH=ON that pulls the library source and builds as described above. That has fairly low impact on the system and allows building from source quite easily. In addition, if the claim on glibc is correct, test the version and disable the safety rounding.
If I’d understood submodules, I might have tried that approach. Is there an existing example that I can use as a template? I agree that the library should be built as part of the swipl build process, probably with source pulled directly from the Inria repo, i.e., I don’t see any advantage in using Bob Burger’s repo directly.
Yes the changes to pl-arith.c are fairly straight forward. The existing prototype conditionally builds either version depending on the existence of the crmath library. But given the performance disparity, I must be doing something wrong.
I’m not sure what I’m looking at but the libswipl.dylib I built (in build/src) is 2.4 MB. Just as Bob’s repo contains a subset of the CORE-MATH source, we can use a subset of the repo for building swipl. When I built a library with just the subset needed, the size of libcrmath.a is 479K (down from 786K). However the size of libswipl didn’t change at all, so I’m not sure what’s going on.
Having said that I don’t really understand the concern. The swipl frameworks are 54 MB. For example, libgmp.dynlib is 1.1MB. Isn’t a libcrmath more like one of these framework libraries? So a useful libcrmath is less than 1% of the total?
I don’t know how the Linix libraries are supporting the cr_ functions, but I’m happy to keep exploring these options given a bit more guidance, or existing example to work from, given my tenuous grasp on CMake.
Most likely something wrong. If libcrmath.a is statically linked to libswipl.dylib it should get bigger. Typically nm is your friend to see what is defined, imported and exported.
seven binary64 functions are integrated into the GNU libc up from release 2.43 (acosh, asinh, atanh, erf, erfc, lgamma, tgamma), and the following binary32 functions are integrated into the GNU libc up from release 2.42: acosf …
So, for these 7 we can test for glibc 2.43 and skip the safety rounding. It seems a bit strange selection though. Possible part of glibc already did proper rounding?
I’m still very reluctant to add 2.2 Mb external code to the source. We do not want to make our own selection but use some external repo as-is. It should support CMake based building. Similarly, growing the core shared library from 2 to 2.5Mb is IMO a high price. I have always disliked bloatware …
So, for now I think the only option I can agree with is to make CMake pull the repo and build it on e.g., a CRMATH option. I never did it, but I know CMake has libraries that can dynamically download additional source packages and include them into the build process.
The rocksDB add-on package is an example of using a submodule from GitHub. It gets automatically downloaded when a git clone --recurse is done.
For this situation, the standard build in the rocksDB repository isn’t suitable, so the build for rocksDB add-on package does that. I don’t see the point of downloading the package at build time – it’s easier to use git’s submodule feature and the source code will take up space either way. And the libedit source includes a submodule https://github.com/SWI-Prolog/winlibedit.git although it’s not a “foreign” submodule.
(PS: I didn’t make the rocksDb add-on package - @jan did that, including the rocksDB submodule, but I’ve modified the add-on package and and have updated the version of the rocksDB submodule.)
It is a different story. RocksDB needs this because the default version as it appears in package managers is often not suitable (incompatible allocator, no C++ RTTI). It is not part of the Prolog core anyway, so size matters less. Libedit is included because first we had a much modified version to support Windows, so we need it anyway. Using libedit removed a lot of legacy code, so the total impact on complexity is smaller than it looks. It is not small (640K source), but still a lot smaller than crmath. Next, BSD libedit has very poor Unicode support, so the current version is quite far of the BSD version 
crmath is quite big, really goes into the kernel (libedit is demand loaded in interactive sessions only) and improves clpBNR only marginally while providing pretty much no benefits for other users.
How can we tell? Only from glibc release notes? And is this the strategy all the platforms will eventually support? Bob Burger’s repo adapts a different strategy: a separate crmath library using the reserved alternative names for correctly rounded functions. To me, that makes more sense in the short-medium term, but who knows?
For the record, the subset of CORE-MATH currently required for SWIP is:
I agree the size of the source seems excessive, particularly since the the whole src directory of swipl-dev is 6.4 MB, but it is what it is. I just don’t know what would be acceptable numbers (greater than 0).
That makes sense; if only I knew how to that.
I don’t disagree. Particularly with current performance characteristics (which I don’t understand), I can’t justify any change to status quo. It works fine (until it doesn’t). The main purpose of the exercise was to evaluate the current state of play. Four years ago many of the 64bit elementary functions didn’t even exist.
Thanks Peter, working examples of how this might work are useful, even if not directly applicable. When it comes to CMake, I’m mostly floundering.
AI can help a lot. A quick good old search reveals Using Dependencies Guide — CMake 4.4.0 Documentation, where it seems you are looking for FetchContent. That should allow you to pull in the dependency. Then build it as static library and add it as dependency to libswipl. Most modern AI coding assistants can probably figure this out for you.
Claude code did 99% of the work to restructure the MacOS installer and allow for signing it. That is way more complicated
Status update:
I’ve formalized the performance test to compare the existing, explicitly rounded functions with the CORE-MATH versions. Each function is called (default rounding) with 1,000,000 random floats between -pi and pi. A base test is used to remove the overhead and the tests are run with the optimise flag on so I believe the net times below are just measuring just the time for the VM codes for the relevant function, e.g., a_func1(sin/1), etc. (The log test has an abs/1 as well to prevent negative input values). The results for the status quo (explicit rounding):
?- f_statistics.
Net time for base is 0.0025700000000004053 sec.
Net time for sin is 0.03979100000000013 sec.
Net time for cos is 0.03829900000000119 sec.
Net time for tan is 0.06262400000000046 sec.
Net time for exp is 0.025326999999999877 sec.
Net time for log is 0.05007800000000007 sec.
Average interval width for sin is 2.0 ULP's.
Average interval width for cos is 2.0 ULP's.
Average interval width for tan is 2.0 ULP's.
Average interval width for exp is 2.0 ULP's.
Average interval width for log is 2.0 ULP's.
true.
So, for example, the net execution time for a_func1(sin/1) would be ~0.04 microseconds. The interval widths are the difference between the to_positive and to_negative rounding in “units of least precision” (ULP’s); 2.0 is the expected value for explicit rounding.
Using the CORE-MATH implementations from Bob Burger’s repo:
?- f_statistics.
Net time for base is 0.002267999999997272 sec.
Net time for sin is 0.42767399999999967 sec.
Net time for cos is 0.4687599999999996 sec.
Net time for tan is 0.6317289999999964 sec.
Net time for exp is 0.07366000000000028 sec.
Net time for log is 0.19276999999999944 sec.
Average interval width for sin is 1.0 ULP's.
Average interval width for cos is 1.0 ULP's.
Average interval width for tan is 1.0 ULP's.
Average interval width for exp is 1.0 ULP's.
Average interval width for log is 1.0 ULP's.
true.
So the CORE-MATH trig functions take about 10 times as long, and exp and log 3 times as long. (Average interval width is 1.0 ULP’s as expected.) As mentioned before, I’m unpleasantly surprised by the performance gap given INRIA’s claims on its web site. It’s hard to lobby for CORE-MATH given these numbers.
On the build exercise, I’ve cobbled together (probably the right term) a build environment as follows:
USE_CRMATH to enable building with the CORE-MATH implementaions, default is OFF.USE_CRMATH=ON, download a current copy of Bob’s repo to the build directory. If successful, set the flag HAVE_CRMATH to 1, otherwise 0.HAVE_CRMATH and the local ‘crmath’ has not been built, build a local libcrmath.a. Add the library to the'LIBSWIPL_LIBRARIES so it’s included in the final libswipl.dylib. The target compile options for the crmath library are set to “-O3 -ffp-contract=on -fno-math-errno”.pl-crmath.h file to the src/ directory and initialize var HAVE_CRMATH_H in cmake/Config.cmake using check_include_file . (To this point, I haven’t found a use for the var).HAVE-CRMATH to the compile environment (#cmakedefine).pl-arith.c for either explicit rounding or CORE-MATH based on HAVE_CRMATH. (Note that all this is done only if O_ROUND_UP_DOWN is true as per current process.)The main cost in source is the new header file (55 lines) and the fairly modest increments to the 4 CMake files affected. Additional overhead for CORE-MATH generated in the build directory:
build/_deps directory containing the external repo (~3 MB)build/src/libcrmath.a (~ 1MB)build/src/libswipl.10.1.xx.dylib increases from 2.1 MB to 2.4MB.I don’t think there’s any need to retain the first two items after a successful build.
Hope this is a sensible way of doing things and addresses concerns raised in the previous posts. While I don’t see a short-medium term future for this work, I’d like to see it captured in a PR for review regardless of whether it gets applied or the USE_CRMATH option is ever enabled (other than by me).
Just create a PR. I think I’m happy to merge that based on what you write. But, I guess it remains future work, probably waiting for performance enhancements and merging into standard C libraries.
If you publish the performance evaluation code, I’m happy to check how this works out on Linux+GCC-16 on AMD.
I’ll do a PR shortly after a tidy up and some investigation into compiler flags for the crmath library. You should definitely review it because my expertise in this whole area is pretty limited.
Longer term we’ll see if the CORE-MATH performance issue can be sorted. If the standard C libraries catch up (I’m not holding my breath), this whole thing (including the status quo explicit outward rounding) can just be disabled using the O-ROUND_UP_DOWN flag.
That would be interesting; the source is below. I’ve commented out the interval width tests as that’s probably not of interest. But while doing this I saw something odd while running it on an Apple M3 laptop/SWIP 10.1.11. When initially loaded, the net times are very small, even negative. If a spy point is set and I skip/leap through the test, at some point it seems to “fix itself”, and that persists after turning debug off. An example below:
% first run - negative times for all functions !?
102 ?- f_statistics.
Net time for base is -0.026874999999999996 sec.
Net time for sin is -0.002577999999999983 sec.
Net time for cos is -0.0028840000000000116 sec.
Net time for tan is -0.0021999999999999936 sec.
Net time for exp is -0.004882000000000011 sec.
Net time for log is 0.002085000000000045 sec.
true.
% set spy point
103 ?- spy(f_timeit).
% New spy point on elemf_test:f_timeit/2
true.
% leap/skip through test. Notice times fixed on log test, but when this occurs varies
[debug] 104 ?- f_statistics.
* Call: (13) elemf_test:f_timeit(base, _56262) ? skip
* Exit: (13) elemf_test:f_timeit(base, 0.21563300000000007) ? leap
* Call: (13) elemf_test:f_timeit(base, _58302) ? skip
* Exit: (13) elemf_test:f_timeit(base, 0.21135199999999998) ? leap
Net time for base is -0.00428100000000009 sec.
* Call: (13) elemf_test:f_timeit(sin, _60304) ? skip
* Exit: (13) elemf_test:f_timeit(sin, 0.23986099999999988) ? leap
Net time for sin is 0.024227999999999805 sec.
* Call: (13) elemf_test:f_timeit(cos, _62306) ? skip
* Exit: (13) elemf_test:f_timeit(cos, 0.23684499999999997) ? leap
Net time for cos is 0.021211999999999898 sec.
* Call: (13) elemf_test:f_timeit(tan, _64308) ? skip
* Exit: (13) elemf_test:f_timeit(tan, 0.23745899999999986) ? leap
Net time for tan is 0.02182599999999979 sec.
* Call: (13) elemf_test:f_timeit(exp, _66310) ? skip
* Exit: (13) elemf_test:f_timeit(exp, 0.2342249999999999) ? leap
Net time for exp is 0.01859199999999983 sec.
* Call: (13) elemf_test:f_timeit(log, _68312) ? skip
* Exit: (13) elemf_test:f_timeit(log, 1.334591) ? leap
Net time for log is 1.1189580000000001 sec.
true.
[debug] 105 ?- nodebug.
true.
% after turning debug off all times look plausible
106 ?- f_statistics.
Net time for base is 0.030053999999999803 sec.
Net time for sin is 0.1210249999999995 sec.
Net time for cos is 0.11901799999999962 sec.
Net time for tan is 0.12725599999999915 sec.
Net time for exp is 0.09614999999999974 sec.
Net time for log is 0.15450799999999898 sec.
true.
I never observed this behaviour on my X86 Mac desktop. I profiled it to check port counts, but all looked normal. This is very repeatable on the M3 laptop, but hard to imagine what might be causing it.
:- module(elemf_test,
[
f_timeit/2,
f_precision/2,
f_statistics/0
]).
:- set_prolog_flag(optimise, true).
f_statistics :-
f_timeit(base,Tbase),
member(F,[base,sin,cos,tan,exp,log]),
f_timeit(F,T), DT is T-Tbase, format("Net time for ~w is ~w sec.\n",[F,DT]),
fail.
/*f_statistics :-
member(F,[base,sin,cos,tan,exp,log]),
f_precision(F,Average), format("Average interval width for ~w is ~w ULP's.\n",[F,Average]),
fail.
*/
f_statistics.
f_time(base) :- between(1,1000000,_), random_input(X), _ is X, fail.
f_time(sin) :- between(1,1000000,_), random_input(X), _ is sin(X), fail.
f_time(cos) :- between(1,1000000,_), random_input(X), _ is cos(X), fail.
f_time(tan) :- between(1,1000000,_), random_input(X), _ is tan(X), fail.
f_time(exp) :- between(1,1000000,_), random_input(X), _ is exp(X), fail.
f_time(log) :- between(1,1000000,_), random_input(X), _ is log(abs(X)), fail.
f_time(_).
f_timeit(F,T) :-
statistics(cputime,T0),
f_time(F),
statistics(cputime,T1),
T is T1-T0.
f_precision(F,Average) :-
f_precision_(0,1000000,F,0.0,Average).
f_precision_(Count,Count,_F,Acc,Average) :- !,
Average is Acc/Count.
f_precision_(InCount,Count,F,Acc,Average) :-
f_bounds_(F,L,H),
delta_(L,H,0,D),
NxtAcc is Acc+D,
NxtCount is InCount+1,
f_precision_(NxtCount,Count,F,NxtAcc,Average).
f_bounds_(sin,L,H) :-
random_input(X), L is roundtoward(sin(X),to_negative), H is roundtoward(sin(X),to_positive).
f_bounds_(cos,L,H) :-
random_input(X), L is roundtoward(cos(X),to_negative), H is roundtoward(cos(X),to_positive).
f_bounds_(tan,L,H) :-
random_input(X), L is roundtoward(tan(X),to_negative), H is roundtoward(tan(X),to_positive).
f_bounds_(exp,L,H) :-
random_input(X), L is roundtoward(exp(X),to_negative), H is roundtoward(exp(X),to_positive).
f_bounds_(log,L,H) :-
random_input(X), L is roundtoward(log(abs(X)),to_negative), H is roundtoward(log(abs(X)),to_positive).
delta_(B,B,D,D) :- !.
delta_(L,H,N,D) :-
NxtL is nexttoward(L,1.0Inf),
NxtN is N+1,
delta_(NxtL,H,NxtN,D).
random_input(X) :- X is 2*pi*random_float - pi.