Correct floating point arithmetic

[jan]

@ridgeworks

I strongly suspect very few people care.Otherwise, there would be correct, precise C math libraries in common usage.

What I care about is that the float rounding modes correctly bound any precise real values. The underlying assumption is that compiled arithmetic is IEEE compliant for all rounding modes and the functions in math libraries are correct within 1 ulp in default “to_nearest” rounding mode. Largely that seems to be true, but rigorous testing is hard and there are no guarantees.

Yes. I found crlibm again, and recall the main issue: it is GPL, not LGPL. That makes it usable only as an optional install for people building from sources as (effectively) changing changing SWI-Prolog to GPL is a no go. So, if we want really correct rounding we need to find another library or stick with your trick assuming the result is correctly rounded within 1 ulp and then step one up/down depending on the rounding mode being up or down. We could use crlibm to do the type of testing @j4n_bur53 is performing to validate the C runtime library. So far, only the MinGW library seems poor. This implies we need to live with it, use another compiler or help the MinGW fix their math library.

Hmmm. A search for crlibm alternatives finds https://openlibm.org/, quoting:

OpenLibm is an effort to have a high quality, portable, standalone C mathematical library (libm). It can be used standalone in applications and programming language implementations.

And

This project is maintained by the Julia Project

That could be worth evaluating. It seems to be a standard package in Ubuntu as well as Macports and uses a BSD license.

[ridgeworks]

@jan

A search for crlibm alternatives finds https://openlibm.org/, quoting:

OpenLibm is an effort to have a high quality, portable, standalone C mathematical library (libm). It can be used standalone in applications and programming language implementations.

And

This project is maintained by the Julia Project

That could be worth evaluating. It seems to be a standard package in Ubuntu as well as Macports and uses a BSD license.

Provided there are no build or licensing issues, this would appear to provide SWIP with a consistent, platform independent arithmetic library; obviously a good thing. The fact that Julia uses and supports it is a definite plus.

However it doesn’t claim to be provably correct under all rounding modes as crlibm is. So the current assumptions made about the underlying math library would continue to be in effect. I guess this is also a good thing, since it should not entail any SWIP source changes to move to OpenLibm - just a library replacement. (Is that right?)

As an aside, I see Julia also has a wrapper for crlibm (GitHub - JuliaIntervals/CRlibm.jl: Correctly-rounded mathematical functions for Julia)

The CRlibm library is state-of-the-art as regards correctly-rounded functions of Float64 arguments. It is required for our interval arithmetic library, ValidatedNumerics. Having gone to the trouble of wrapping it, it made sense to release it separately; for example, it could be used to test the quality of the OpenLibm functions.

It also claims:

CRlibm is distributed under the GNU Lesser General Public License (LGPL).

and when I look at

GitHub - taschini/crlibm: A mirror of the CRLibm project from INRIA Forge

a sampling of the copyright notices does seem to indicate that LGPL licensing is in effect. So I’m not sure what to believe.

[jan]

@ridgeworks

just a library replacement. (Is that right?)

Probably

@ridgeworks

CRlibm is distributed under the GNU Lesser General Public License (LGPL).

That seems right. The strange thing is that the repo contains both the GPL and LGPL license documents (missed the latter). The relevant source files refer to the LGPL, so I guess that is it. As we already have GMP that is acceptable (although I still want to have a GMP replacement to provide full functionality under the BSD license)

But then, crlibm does not seem to come with a configuration to build a dynamic library. If we do not fix that it effectively has more license consequences. I took a short look at the configuration (configure.ac) and it merely lists a series of CPUs, OSes and compilers. That means it gets outdated and we’ll have to do the necessary research and possible modifications for new CPUs, such as the ARM (and Apple’s M1). Together with a library that is basically unmaintained, that seems a no-go.

Another option might be mpfr At least this library is widely supported. The MPFR page lists a large set of other candidates. The aims to keep in mind:

• License (high preference for permissive open source (BSD,MIT,Apache,…), optionally LGPL).
• Properly maintained and cross-platform (OS, compiler and CPU)
• Performance
• Size

On the MPFR page I also find LibBF, which claims correct IEEE 754 rounding and both unbounded integers and floats.

Any volunteer willing to make a table?

(assuming correctness is fulfilled; without that we just as well stick with whatever C runtime libm that comes with the system).

@j4n_bur53

I guess this is something for a new thread. Or is this more specifically (**)/2
related?

You latest cases6 test obviously depends on float **, which seems pretty much ok on Linux and MacOS, but broken on MinGW on Windows. @ridgeworks interval arithmetic (thus) won’t work on this platform.

[ridgeworks]

@jan

Another option might be mpfr At least this library is widely supported.

I had briefly looked at mpfr before but dismissed it because it appeared to be a multi-precision float library. That was a bit hasty because it looks like you can set a mpfr_t’s precision to 53 bits for compatibility with IEEE doubles.

The relevant functions appear to calculate results using extra precision and then round according to rounding mode argument. So that seems to meet my requirements.

We already have cr_FN wrappers to standard elementary functions in pl-arith.c (exp,log, pow, sin, …), so it should be fairly straight froward to re-implement these to use the mpfr API instead. For example, currently we have:

// #define general form of unary function with no limits on result
#define CR_FUNC(func) \
static double cr_##func(double in) \
{ double result; \
  int roundMode = fegetround(); \
  if ( roundMode != FE_TONEAREST ) fesetround(FE_TONEAREST); \
  result = func(in); \
  if (isfinite(result)) switch( roundMode ) \
  { case FE_UPWARD     : result = nexttoward(result, INFINITY); break; \
    case FE_DOWNWARD   : result = nexttoward(result,-INFINITY); break; \
    case FE_TOWARDZERO : result = nexttoward(result, 0); break; \
  } \
  if ( roundMode != FE_TONEAREST ) fesetround(roundMode); \
  return result; \
}

using mfpr instead:

# include <mpfr.h>
// #define general form of unary function with no limits on result
#define CR_FUNC(func) \
static double cr_##func(double in) \
{ double result; \
  mpfr_t arg = mpfr_init_set_d(arg, in, MPFR_RNDN); \
  mpfr_rnd_t rnd;  \
  switch (fegetround()) \
  { case FE_TONEAREST : rnd = MPFR_RNDN; \
    case FE_UPWARD : rnd = MPFR_RNDU; \
    case FE_DOWNWARD : rnd = MPFR_RNDD; \
    case FE_TOWARDZERO : rnd = MPFR_RNDZ; \
  }; \
  mpfr_##func(arg, arg, rnd); \
  result = mpfr_get_d(arg, rnd); \
  mpfr_clear(arg); \
  return result; \
}

A better version would use a per-thread static variable for arg so it wouldn’t have to be initialized and cleared on each call. As currently structured it’s straight forward to swap in different API’s by redefining a few macros.

@jan

On the MPFR page I also find LibBF,

Also a possibility, but mfpr would seem to have the upper edge on the support side.

@jan

(assuming correctness is fulfilled; without that we just as well stick with whatever C runtime libm that comes with the system).

Agreed. Assuming that mpfr uses extra precision bits and then rounds as requested, it should meet the correctness requirement unless there are bugs. That possibility of serious bugs seems to be fairly remote as it’s been under active development for over 20 years, currently at release 4.1. Also performance should be decent given its lengthy history, although this needs to be verified. One would hope converting between doubles and 53-bit mfpr_t’s would be fairly efficient.

So at this point, mpfr seems to tick the boxes - licensing, support, functionality, platform independence; certainly worth an in-depth look. I’d be happy to work on the required changes to pl-arith.c if someone else can handle the build (library and app) side.

@jan

depends on float **, which seems pretty much ok on Linux and MacOS, but broken on MinGW on Windows.

Can someone elaborate? Is it a matter of a ulp difference in the results, or is it something more serious?

[ridgeworks]

mpfr background: https://members.loria.fr/PZimmermann/talks/mpfr18.pdf

[jan]

@ridgeworks

@jan

On the MPFR page I also find LibBF,

Also a possibility, but mfpr would seem to have the upper edge on the support side.

I had a closer look at LibBF. It is a very attractive library. It is pure C with a demo that you get a significantly faster version by enabling some CPU specific optimizations. It is a small pure C library using the MIT license, so we no longer need conditional compilation (with/without GMP). What makes it attractive is that it provides a single representation for computing with binary floats, decimal floats and big integers (bigints are represented as infinite precision floats). The next great point is that it is stateless. Functions take arguments to deal with precision and rounding (all modes) and a little context argument passes in the memory management functions. Functions properly return failure due to allocation failure. This means we can have a per-thread allocation context with and we do not need the hacky stuff to protect the process from aborting that we need with GMP (guess size, when in reasonable boundaries try the GMP, but wrap it using setjmp/longjmp as you may not return from an allocation attempt and do some tricky allocation wrapping).

I’ve implemented sin(X) and mul(X,Y) for bigints. sin(X) is about 10 times slower than Linux libm, so I did the same using mpfr, which is about 5 times slower. Mul wasn’t integrated as it eventually needs to be, but the performance seems not too far away from GMP. The package claims it looses on small numbers, but actually wins against GMP on (some) really big numbers. In our case, we can avoid a lot of the overhead due to the cleaner interface, providing hope for “small” big integers.

The main issue is the lack of support for rational numbers. That would need to be resolved. I do not know enough about these issues to judge how hard that would be. Of course, a rational number is just a pair of big integers, but how easy is the required normalization? Are there more clever tricks than high school math to deal with the basic multiplication, addition, etc? Probably some things are missing (e.g., random numbers). Some should be easy to replace.

Some float functions like sinh() are missing. We could wrap the libm versions for that. Possibly we could do that for all libm functions that we use to provide a runtime switch between quick and dirty or precise floats.

Another issue is how one should make rounding, precision and binary vs. decimal floats accessible to the Prolog user?

All in all, I think LibBF should be a serious candidate. Copying big numbers to the stack is comparable to GMP. The uniform representation, the fact that we do not need conditional compilation (#ifdef) all over the place and proper failure codes on allocation failures should make the implementation a lot cleaner.

Portability seems no issue. The library has been extensively tested against mpfr. I tried sin(X), comparing native, LibBF and mpfr on random floats. About 1 in every few thousand attempts the native one disagrees, while LibBF and mpfr agree with each other. Using wcalc (based on mpfr), computing at higher precision I could validate that LibBF and mpfr were correct (by hand, so not many verifications).

Only, this is a serious amount of work. First access whether it is really a good idea, then implement. To get this done, we need programmers and/or money.

@ridgeworks

Can someone elaborate? Is it a matter of a ulp difference in the results, or is it something more serious?

I don’t recall the details, but I think it is more serious than that on some functions.

[ridgeworks]

@jan

Only, this is a serious amount of work.

My agenda is much smaller than yours, specifically to replace the unreliable, platform dependent calls to libm with calls to a correctly rounded, platform independent library. From this perspective, there is no need to replace GMP or change anything else that the Prolog user would notice (i.e., decimal floats, precision, etc.). I wouldn’t think this would be a serious amount of work.

GMP may have its issues but replacing it would be a significantly large project and, as a Prolog user, I don’t see much (any?) payback.

It remains an open question whether mpfr or libBf would be preferable. libBf would definitely be the preferred choice if GMP was being removed (mpfr requires GMP); otherwise my reading is that mpfr would win based on support and common usage.

Maybe this can be split into two projects. The short-medium term would use GMP/mpfr to provide correctly rounded functions and require a fairly small effort (my agenda). A second project could evaluate/implement a replacement for GMP/mpfr with some alternative, e.g., libBf, and would require significantly more resources (?? agenda).

[jan]

@j4n_bur53

So LibBF could deliver better results for (**)/2?

I think it pretty much guarantees accurate results for all the functions it supports.

@ridgeworks

I wouldn’t think this would be a serious amount of work.

Probably not. If you replace the core functions in the ar_() functions with cr_(), e.g. cr_sin(), cr_pow(), etc (as is demanded by the standard anyway) and create a new file pl-mpfr.c that defines these functions defined in mpfr, the impact is small. Here, I assume that mpfr does its allocations on top of gmp and thus the complicated stuff that must prevent crashes and memory leaks using gmp still works.

We might want a flag or some other mechanism to switch to the libm implementations. In my measurement sin(X) of mpfr is about 4 times slower than libm’s sin(). This still brings me back to an older plan: have a struct/array with function pointers such that we can set the float functions through the foreign API. In that case we can make a loadable module to provide correct rounding vs. libm arithmetic and you can do :- use_module(library(mpfr)). if you need it.

@ridgeworks

It remains an open question whether mpfr or libBf would be preferable. libBf would definitely be the preferred choice if GMP was being removed (mpfr requires GMP); otherwise my reading is that mpfr would win based on support and common usage.

Right. Going libBf only makes sense if we use it to also get rid of GMP. The main reasons to do so are that we can make libBF unconditionally part of SWI-Prolog while we cannot do so for GMP for license reasons. Using LibBF for the whole code base is indeed a considerable amount of work, but should overall imply a serious simplification of the code. Simpler code means less bugs Probably some things will get slower and some will get faster. We’d need a more serious comparison of the performance.

@ridgeworks

Maybe this can be split into two projects. The short-medium term would use GMP/mpfr to provide correctly rounded functions and require a fairly small effort (my agenda). A second project could evaluate/implement a replacement for GMP/mpfr with some alternative, e.g., libBf, and would require significantly more resources (?? agenda).

I’m not a big fan of that as it means adding and later removing code and dependencies. It is probably the most practical though … Playing around with the precision is probably not a good idea as it would add another type to the already complicated mix of types in arithmetic. That is the beauty of LibBF: it removes a type.

[ridgeworks]

@jan

If you replace the core functions in the ar_() functions with cr_(), e.g. cr_sin(), cr_pow(), etc (as is demanded by the standard anyway) and create a new file pl-mpfr.c that defines these functions defined in mpfr, the impact is small. Here, I assume that mpfr does its allocations on top of gmp and thus the complicated stuff that must prevent crashes and memory leaks using gmp still works.

Note that we’re already “replacing” the core functions in ar_* functions using macros that are consolidated in ~200 lines of C in pl-arith.c to explicitly outward round the libm elementary functions. These would just have to be rewritten to use the correctly rounded library of choice. Taking this route, there’s no need to create a pl-mpfr.c. The only requirement is to include the necessary header (<mpfr.h> or <libbf.h>).

Also note that I’m not talking about exposing any of the Prolog machinery to instances of mpfr_t (or bf_t) or providing any general purpose API to the library. The lifetime of an mpfr_t (or bf_t) instance is very short (lifetime of a cr_* function call). If that’s the case, that should mean that none of the existing GMP “complexities” should be affected independent of which library is used.

@jan

We might want a flag or some other mechanism to switch to the libm implementations. In my measurement sin(X) of mpfr is about 4 times slower than libm’s sin().

Is that the raw C function, or as measured using Prolog arithmetic?

@jan

Using LibBF for the whole code base is indeed a considerable amount of work, but should overall imply a serious simplification of the code. Simpler code means less bugs

Except for rationals. Wouldn’t these have to be completely re-invented (or dropped?!).

@jan

I’m not a big fan of that as it means adding and later removing code and dependencies.

The code itself could be conditional just like the current version of cr_* functions (flag O_ROUND_UP_DOWN). Obviously there’s a dependance on an extra library (mpfr or libbf) but is that a big issue? Again, we’re talking about a very localized usage (200 lines in 1 file?).

Longer term we seem to be talking about removing dependancies, at least. But I don’t have a lot of confidence the long term will ever happen, which just leaves the status quo. That’s not bad except for Windows and the “stated assumption”.

@jan

That is the beauty of LibBF: it removes a type.

Which type? I’m only talking about using the correctly rounded functions on doubles, not about introducing new, or removing existing arithmetic, types. The libraries being discussed are overkill for just this requirement, but they happen to do the job for doubles.

For “my agenda” either mpfr or libBf seem up to the task. libBf is missing some functions, but these could just be left using the explicit rounding of libm already in place.

[jan]

@ridgeworks

These would just have to be rewritten to use the correctly rounded library of choice.

@ridgeworks

Also note that I’m not talking about exposing any of the Prolog machinery to instances of mpfr_t (or bf_t) or providing any general purpose API to the library. The lifetime of an mpfr_t (or bf_t) instance is very short (lifetime of a cr_* function call).

That is also my understanding. I think I prefer a new source file though. pl-arith.c is bad enough as it is Anyway, that is minor. The #define pow cr_pow is apparently not allowed, so it may be a good idea to get rid of it as the same time and call the thing called by ar_pow() directly cr_pow.

@ridgeworks

Is that the raw C function, or as measured using Prolog arithmetic?

Including Prolog. So, the draw performance difference even larger. Of course, sin(X) is just one (tested with random angles between 0 and 2pi).

@ridgeworks

Longer term we seem to be talking about removing dependancies, at least. But I don’t have a lot of confidence the long term will ever happen, which just leaves the status quo. That’s not bad except for Windows and the “stated assumption”.

I don’t know. Depends on the interest in high precision floats, decimal floats (has been asked for in the past) and big integers in the version with only permissive licenses. In addition to my desire to simplify the code rather than making it uglier when possible

@ridgeworks

Which type?

mpz_t. LibBF only has floats and considers bigints a special case of infinite precision floats. Rather funny perspective. Seems to work well though. In addition, as LibBF is a dependency that can always be present, we could consider removing the triple integer representation (tagged, indirect 64 bit and MPZ) and move to a dual model (tagged, bf_t).

Anyway, using mpfr is a serious option for the short term. Also note that libBF and mpfr (which agree) quickly proved that assumption that libm is correctly rounding to nearest to be wrong (at least for sin(X).

But, I wonder what makes more people happy: fast approximate floats or 100% correct slow floats? And, not considering the numbers, does making sin(X) 4 times slower cause serious trouble for some applications?

[ridgeworks]

@jan

Anyway, using mpfr is a serious option for the short term. Also note that libBF and mpfr (which agree) quickly proved that assumption that libm is correctly rounding to nearest to be wrong (at least for sin(X).

But, I wonder what makes more people happy: fast approximate floats or 100% correct slow floats? And, not considering the numbers, does making sin(X) 4 times slower cause serious trouble for some applications?

Implementation details aside, I think you understand why a short term solution is helpful, at least to me. As for speed vs. correctness, I suspect there would be little performance impact for the vast majority of applications given your measurements. But at the risk of introducing yet another arithmetic environment flag, it could be made a runtime option which would direct the cr_* functions to the appropriate strategy.

[jan]

@ridgeworks

I suspect there would be little performance impact for the vast majority of applications given your measurements. But at the risk of introducing yet another arithmetic environment flag, it could be made a runtime option which would direct the cr_* functions to the appropriate strategy.

I agree that the affected float functions are unlikely to have a lot of impact on 99.9% of the users. I typically would like to guarantee that existing programs do not slow down more than a few percent by new versions. The limit for a “reasonable” slowdown of course depends on what is gained. Moving to soft floats is rather dubious in that sense. They are a lot slower and we do not give the benefits back (e.g., arbitrary precision).

Having to load a library to get soft and precise floats is probably the cleanest solution. An implementation hack could be for that library to merely set some flag. I think a global setting is fine here. Enabling soft floats for your solver makes the semantics more precise, which should not negatively affect anyone. Other code could be negatively affected by the performance loss, but as you claim, that should not affect too many programs.

In that light, would it not be a better to provide predicates for your correctly rounding float functions?, e.g., a library that provides sin(RoundMode, Angle, Result) (not sure about the argument order)?

Please do be aware that adding dependencies is not a choice to make easily. It means updating the build scripts for the various platforms, many of which are not maintained by us, but e.g, by Linux distro package managers.

[ridgeworks]

@jan

Moving to soft floats is rather dubious in that sense. They are a lot slower and we do not give the benefits back (e.g., arbitrary precision).

Absolutely agree. I do not want soft and (arbitrarily) precise floats, IEEE 64 bit doubles are more than adequate for my purposes. What I do want is correctly rounded arithmetic for all rounding modes and the only current deficiency is lack of correctly rounded elementary functions (current hack notwithstanding), further exacerbated by platform dependant implementations of libm.

It appears the only correctly rounded library for doubles is crlibm and it only supports a partial set of the required functions. Furthermore, there are build and (possibly) licensing issues that have so far excluded it from consideration.

The other option on the table is to use soft floats with arbitrary precision and correct rounding support. libBF is small and simple (no state) but only supports a subset (~50%) of the required functions. mpfr is GMP based, more complicated to use (has some state), but is complete (100%).

Whichever library is used, it entails an additional dependancy (but possibly removes the dependancy on platform specific versions of libm) and will impact performance (elementary functions only). On that front, given:

?- time((between(1,1000000,_),_ is sin(1.0),fail)).
% 2,000,001 inferences, 0.215 CPU in 0.223 seconds (96% CPU, 9297833 Lips)
false.

?- time((between(1,1000000,_),fail)).
% 1,000,001 inferences, 0.078 CPU in 0.079 seconds (98% CPU, 12895086 Lips)
false.

I conclude the cost of the sin evaluation in the 8.5.20 is around 150 ns. If I apply a factor of 4 to this based on your earlier testing, that’s still well under a microsecond which I would find more than acceptable. So unless there’s some other evidence to the contrary, I don’t think performance is an issue.

Bottom line, my opinion: no to soft floats/arbitrary precision; yes to correctly rounded elementary functions using whatever library is up to the task.

[jan]

@j4n_bur53

The ISO core standard allows that already for pow() since it has (**)/2 and (^)/2.

? Doesn’t the ISO standard simply define **/2 in the same way as the C runtime library defines pow() and ^ as its integer counterpart, i.e. (not SWI-Prolog)

?- Y is 2**5.
Y = 32.0.
?- Y is 2^5.
Y = 32.

I’ve always considered this weird. We don’t have A is 2+2 to return 4.0 either. Prolog is a dynamically typed language and logical language. As equivalence is big thing in logic we’d better be precise when possible. So, both are considered equivalent to SWI-Prolog and the system will always try to produce a precise result (notably if prefer_rationals is set to true). If you do not want that, cast something in the chain of computations to float and all subsequent operations will be executed using floating point arithmetic. Good thing is that that practice is transparent and the code works in ISO compliant systems as well.

Good thing is we can decide when to do the transformation, i.e. we can (now) get a perfectly rounded float from integer power using X is float(A**B) which will do the computation using integers or rationals when possible and use our reliable rounding (thanks to @ridgeworks) to get the correct result. Or we do X is float(A)**B to perform the computation using C pow() with (as we now know) not always perfect results.

[j4n_bur53]

@jan

? Doesn’t the ISO standard simply define **/2 in the same way as the C runtime library defines pow() and ^ as its integer counterpart, i.e. (not SWI-Prolog)

?- Y is 2**5.
Y = 32.0.
?- Y is 2^5.
Y = 32.

The difference between the two notions “less precise and more
speed” verus “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.

[jan]

@ridgeworks

Bottom line, my opinion: no to soft floats/arbitrary precision; yes to correctly rounded elementary functions using whatever library is up to the task.

Sorry to say, but I see a bit too much “for my purposes”. Dropping float performance by a factor 4 and including a new dependency is a rather high price for a few ulps Precise rounding is one improvement, arbitrary precision or decimals might be something that is on the wish list of others. I have had a discussion with a commercial user about supporting decimal floats. If the redesign would also simplify other areas (like LibBF simplifying license handling and memory management) that is another argument.

The more I think about it, the more I’m attracted to put this stuff outside Prolog’s core. I see two ways for that, one is to provide a foreign interface that allows something like this, where (of course), the signature of my_very_nice_sin_function() must be the same as the C runtime function (double → double). The implementation would simply change a function pointer.

PL_set_function("sin", my_very_nice_sin_function);

A sensible alternative would be (as claimed before) to provide a foreign library with predicates for these functions and use that.

Note that, as far as I can tell, it is only the MinGW runtime that is (sometimes) rather far away from the correct result. We could also seek for a solution to that problem only.

[ridgeworks]

@jan

Sorry to say, but I see a bit too much “for my purposes”. Dropping float performance by a factor 4 and including a new dependency is a rather high price for a few ulps

Not really about “a few ulps” and more about correctly bounded results “for my purposes”. The better the libraries, the tighter the bounds. And I’m not advocating “Dropping float performance by a factor 4” in general, only on those elementary functions under discussion.

@jan

Precise rounding is one improvement, arbitrary precision or decimals might be something that is on the wish list of others.

No problem, but they’ve yet to appear on this thread.

@jan

A sensible alternative would be (as claimed before) to provide a foreign library with predicates for these functions and use that.

I’ve no interest in new predicates or external modules.

It’s really just about proper support for IEEE floats (correctly rounded results for all arithmetic operations on doubles). The current system provides all this modulo the assumption regarding 1 ulp precision on nearest rounding for elementary function. Your call if and how SWIP installations satisfy that assumption including:

@jan

Note that, as far as I can tell, it is only the MinGW runtime that is (sometimes) rather far away from the correct result. We could also seek for a solution to that problem only.

and just to repeat, I’m concerned about bounding, not “a few ulps”, but it can be difficult to separate these issues.

[jan]

Bottom line is that I very much dislike the idea of these complications and performance degradation to the core of the system. It would of course be great if all was fast and correct, but it appears that all alternatives to achieve this have their problems in terms of license, portability, support, etc.

I think the simplest is to stay at the square we are. Add some code to check for Windows (current_prolog_flag(windows, true)) and tell such users to use WSL2 with the Ubuntu version.

P.s. Why exactly is correct rounding important for this library? Does the computation fail or doesn’t it terminate if this demand is not met or does it just return wrong results?

[ridgeworks]

@jan

I think the simplest is to stay at the square we are. Add some code to check for Windows (current_prolog_flag(windows, true)) and tell such users to use WSL2 with the Ubuntu version.

And I have no problem with that. Hopefully the platform implementations of libm will improve over time, but I’m not holding my breath.

@jan

P.s. Why exactly is correct rounding important for this library? Does the computation fail or doesn’t it terminate if this demand is not met or does it just return wrong results?

My requirement is for constructing mathematically correct reals (not floats) which can be represented as a containing interval with lower and upper bound values, usually imprecise floats. Intervals are initially wide and narrow as additional constraints are defined. But to be correct the real value must always be contained in the interval, which can be done by outward rounding of any floating point calculation on the interval bounds. Conservative bounds are more important than precision to preserve the correctness properties of real arithmetic. 64 bits is more than adequate for most (all?) purposes as long as the the answer is properly bounded.

Failure is expected when the constraints are inconsistent. Non-termination (in fact, just appearance of non-termination since intervals cannot get wider) is not an issue because internal throttling mechanisms are in place. Results should never be “wrong” if the set of constraints is adequately defined (subject to interpretation).

[jan]

Thanks for the clarification. To me, that sort of closes the issue. As is, all seems to work well except for the Windows version. Anyone who wants this correctness on Windows can use WSL and hopefully libm of notably MinGW will improve over time. They too are probably faced with the fact that there is a lack of good and portable libraries for floating point arithmetic with license conditions they can accept.