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