The "failure" cases are a result of over-sensitivity to the boundary restrictions on asin(x) and depend on the current environment settings for floating point exception handling. From the glibc math library man page:
Domain error: x is outside the range [-1, 1]
errno is set to EDOM. An invalid floating-point exception (FE_INVALID) is raised.
This can be triggered by round-off error near the boundary. You wouldn't see the failure if your machine settings were less strict about floating-point exception handling. Gnuplot correctly catches that potential trap for the case x = 1+ε but does not check for x = -1-ε. I will fix that.
As to the fluctuating values in this range, it seems due to the limit of precision (1.e-08) provided by the asin routine in the glibc math library. I get the same fluctuation in values if I print directly from C rather than using the gnuplot routine. Tested with glibc version 2.36. I can't fix that.
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
If the calculation were carried out to infinite precision then it would be true that beta would strictly lie in the range [-1:1] and thus asin(beta) would satisfy the documented domain restriction. Given the 64-bit precision used here, beta may exceed this range by some amount ε on the order of 1.e-15. One can quibble that the math library itself should allow this much slop, but apparently it doesn't. So the gnuplot code can avoid the resulting floating point domain error by enforcing if (beta < -1) beta = -1 prior to the call to asin(beta).
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
BTW, this behaviour of glibc seems especially ridiculous since what I’m trying to visualize now¹⁾ needs calculations with 600 decimal places…
¹⁾ I’m trying to make an us-lowly-humans-digestible presentation of “hyper-asymptotics” and “hyperterminants” — and they become more “clean” when the “usual pesky terms of syper-asymptotics” become larger and larger.
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
If you need 600 decimal places of precision then you simply cannot use standard math operations and libraries in any language that I know of. Intrinsic variable types like double (64 bit) or even long double (128 bit) cannot store a value to that level of precision even before you get to the step of manipulating it or passing it to a function [1]. There are special "arbitrary precision" libraries that you might use instead of the normal arithmetic operators and math library functions , but gnuplot was not written to use any of them.
I much appreciate your effort in finding and reporting a fixable problem with gnuplot's asin() function (the cases where it returned NaN rather than -pi/2). The fix will be in 6.0.4 (next month).
[1] In the particular case you show as a test, modifying gnuplot's internal routine f_asin() to use long double intermediate variables and the corresponding glibc call asinl(beta) is sufficient to remove the visible artifacts in the resulting plot. I.e. 112 bits of precision rather than 52 bits. I tried that in the course of pinning down exactly where the greater precision is needed. There are still artifacts of course, but they are smaller than the visible resolution of the plot. If you were dumping numerical output rather than graphics, the finite precision would still be detectable. Generalizing this to replace double with long double everywhere would be a huge task. Furthermore the algorithms used for special functions may not be accurate to a matching level of precision.
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
If you need 600 decimal places of precision then you simply cannot use standard math operations and libraries in any language that I know of.
I suspect this is just a glitch in your memory. With any contemporary library “focussed for math”, anything ≪10⁵–10⁶ decimal places should not be a problem. ¹⁾ (And in many cases it is needed, since the standard formulas often include major cancellations — with the associated precision loss. Like the answer being ∼1, but to get it you need to subtract two 600-digits numbers.)
¹⁾ The only exception is solving ODEs. AFAIK, currently the “hard wall” is at about 1,000 decimal places.
I use gp/PARI, but when its plotting abilities are not flexible enough, I need to transfer results to gnuplot (via tables). (It is already about 15 years as it can efficiently integrate in ultra-high precision, so the last hurdle has been overcome.)
I think we may be talking past each other, probably because we're coming from quite different fields. You're obviously in math and have a different notion of "standard libraries and languages".
What I meant was that the native syntax of C/Fortran/R/whatever only provides variables holding IEEE 754 compatible 64- or 128- bit floating point representations. The corresponding standards like c17 only promise support for the respective basic math functions, asin() or asinl() in this case. Gnuplot is written in C and its internal representation of real or complex numbers is 64-bit. It does contain additional code and links to external libraries in order to support special functions, but the precision is still limited by the 64-bit variables that are being passed in and used for computation.
I am not familiar with gp/PARI, but I gather from a quick look I gather that PARI is indeed an arbitrary precision library, while gp provides a C-compatible interface to it. Great, but that's at least one step removed from what I normally deal with.
So OK, I'm with you as to your use of high-precision calculations to generate numbers that you then want to plot. But I don't know of a way to get floating point values into gnuplot that are not limited to the roughly 16 decimal places of precision that you get with IEEE 754 64-bit representation. I hope that isn't a fatal limitation for your intended use.
Anyhow, I believe I have fixed the code that produced the specific glitchy plots you report for asin(x) with x < -1. That was old, old code that probably dates back to the 1980s and predates the IEEE standards. I am guessing that noise-level deviation from the domain [-1:1] did not cause an FPE trap on older hardware. I found an equivalent glitch in the acos(x) code as well, also now fixed. I will have a look at the rest of the [old!] routines in that source file to see if there are other cases where treatment of real values as complex values that just happen to have zero imaginary component introduces extra computation that may produce noise-level variation in the low order bit[s].
Thanks again for finding and reporting this.
cheers,
Ethan
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
Sorry, my fault. I should have bracketed my remarks about “all libraries” by <pendantic>/</pendantic>!
As a developer of graphic interfaces to math facilities, I’m well aware of the inexorable constraints the maintainer is confined to. Obviously, I did not talk about libraries “easy to be used from gnuplot”. Still, I must have stated this explicitly!
Thanks for the fixes!
P.S. PARI is the C library — and an interpreted language, while gp is the CLI/graphing calculator with bridges to (95% of) capabilities of this library. It would not be hard to bridge things-withoug-side-effects from PARI to be useful from gnuplot (like I do in Math::Pari).
However, the bulk of work done by this module is to handle side effects (variables, defined functions etc.) transparently — and I do not see how to make it simple for gnuplot. And without this (and/or without representing the true polymorphism of the PARI data type) the functionality of this library is going to be severely limited.
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
Fix is now in 6.1, will be in 6.0.4
In addition to the avoidance of an FPE trap in the general case, an additional bit or two of precision can be obtained in the special case [obviously common!] of x purely real (no imaginary component). In that case for x outside the range [-1:1] we can skip the general calculation of beta(x) because it must be either -1 or +1. That avoids noise-level fluctuation in the low order bit(s) of beta(x).
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
The "failure" cases are a result of over-sensitivity to the boundary restrictions on asin(x) and depend on the current environment settings for floating point exception handling. From the glibc math library man page:
Domain error: x is outside the range [-1, 1]
errno is set to EDOM. An invalid floating-point exception (FE_INVALID) is raised.
This can be triggered by round-off error near the boundary. You wouldn't see the failure if your machine settings were less strict about floating-point exception handling. Gnuplot correctly catches that potential trap for the case
x = 1+εbut does not check forx = -1-ε. I will fix that.As to the fluctuating values in this range, it seems due to the limit of precision (1.e-08) provided by the asin routine in the glibc math library. I get the same fluctuation in values if I print directly from C rather than using the gnuplot routine. Tested with glibc version 2.36. I can't fix that.
3−1hardly looks like ε to me. 😌 So, frankly speaking, I do not know what you mean here…Why not? It is just one iteration of Newton’s method.
There is no "3-1' anywhere. There is in this case a calculation
If the calculation were carried out to infinite precision then it would be true that beta would strictly lie in the range [-1:1] and thus
asin(beta)would satisfy the documented domain restriction. Given the 64-bit precision used here, beta may exceed this range by some amount ε on the order of 1.e-15. One can quibble that the math library itself should allow this much slop, but apparently it doesn't. So the gnuplot code can avoid the resulting floating point domain error by enforcingif (beta < -1) beta = -1prior to the call to asin(beta).BTW, this behaviour of
glibcseems especially ridiculous since what I’m trying to visualize now¹⁾ needs calculations with 600 decimal places…If you need 600 decimal places of precision then you simply cannot use standard math operations and libraries in any language that I know of. Intrinsic variable types like
double(64 bit) or evenlong double(128 bit) cannot store a value to that level of precision even before you get to the step of manipulating it or passing it to a function [1]. There are special "arbitrary precision" libraries that you might use instead of the normal arithmetic operators and math library functions , but gnuplot was not written to use any of them.I much appreciate your effort in finding and reporting a fixable problem with gnuplot's asin() function (the cases where it returned NaN rather than -pi/2). The fix will be in 6.0.4 (next month).
[1] In the particular case you show as a test, modifying gnuplot's internal routine f_asin() to use
long doubleintermediate variables and the corresponding glibc callasinl(beta)is sufficient to remove the visible artifacts in the resulting plot. I.e. 112 bits of precision rather than 52 bits. I tried that in the course of pinning down exactly where the greater precision is needed. There are still artifacts of course, but they are smaller than the visible resolution of the plot. If you were dumping numerical output rather than graphics, the finite precision would still be detectable. Generalizing this to replacedoublewithlong doubleeverywhere would be a huge task. Furthermore the algorithms used for special functions may not be accurate to a matching level of precision.I suspect this is just a glitch in your memory. With any contemporary library “focussed for math”, anything ≪10⁵–10⁶ decimal places should not be a problem. ¹⁾ (And in many cases it is needed, since the standard formulas often include major cancellations — with the associated precision loss. Like the answer being ∼1, but to get it you need to subtract two 600-digits numbers.)
I use
gp/PARI, but when its plotting abilities are not flexible enough, I need to transfer results tognuplot(via tables). (It is already about 15 years as it can efficiently integrate in ultra-high precision, so the last hurdle has been overcome.)Last edit: Ilya Zakharevich 5 days ago
I think we may be talking past each other, probably because we're coming from quite different fields. You're obviously in math and have a different notion of "standard libraries and languages".
What I meant was that the native syntax of C/Fortran/R/whatever only provides variables holding IEEE 754 compatible 64- or 128- bit floating point representations. The corresponding standards like c17 only promise support for the respective basic math functions, asin() or asinl() in this case. Gnuplot is written in C and its internal representation of real or complex numbers is 64-bit. It does contain additional code and links to external libraries in order to support special functions, but the precision is still limited by the 64-bit variables that are being passed in and used for computation.
I am not familiar with gp/PARI, but I gather from a quick look I gather that PARI is indeed an arbitrary precision library, while gp provides a C-compatible interface to it. Great, but that's at least one step removed from what I normally deal with.
So OK, I'm with you as to your use of high-precision calculations to generate numbers that you then want to plot. But I don't know of a way to get floating point values into gnuplot that are not limited to the roughly 16 decimal places of precision that you get with IEEE 754 64-bit representation. I hope that isn't a fatal limitation for your intended use.
Anyhow, I believe I have fixed the code that produced the specific glitchy plots you report for asin(x) with x < -1. That was old, old code that probably dates back to the 1980s and predates the IEEE standards. I am guessing that noise-level deviation from the domain [-1:1] did not cause an FPE trap on older hardware. I found an equivalent glitch in the acos(x) code as well, also now fixed. I will have a look at the rest of the [old!] routines in that source file to see if there are other cases where treatment of real values as complex values that just happen to have zero imaginary component introduces extra computation that may produce noise-level variation in the low order bit[s].
Thanks again for finding and reporting this.
cheers,
Ethan
Sorry, my fault. I should have bracketed my remarks about “all libraries” by
<pendantic>/</pendantic>!As a developer of graphic interfaces to math facilities, I’m well aware of the inexorable constraints the maintainer is confined to. Obviously, I did not talk about libraries “easy to be used from
gnuplot”. Still, I must have stated this explicitly!Thanks for the fixes!
Fix is now in 6.1, will be in 6.0.4
In addition to the avoidance of an FPE trap in the general case, an additional bit or two of precision can be obtained in the special case [obviously common!] of x purely real (no imaginary component). In that case for x outside the range [-1:1] we can skip the general calculation of beta(x) because it must be either -1 or +1. That avoids noise-level fluctuation in the low order bit(s) of beta(x).