Incorrect results for a couple of real cpp_bin_float functions:ℝ→ℝ #264
Comments
|
I should have sorted EDIT: I fixed this in OP |
|
I would like to look at a few instances of spurious results for LogGamma[x]. But I have a few questions.
I get a complex value for that one. |
MPFR can dynamically increase its precision (and the digits in the type at hand) when calculating near singularities. I think most Boost.Math special functions use fixed-width from start to end. So I believe calculations often could get lossy near singularities, unless something like a Taylor series could be found for some of the ranges you suspect. |
yes, in this test I only work with real valued arguments and results. And thanks, of course lgamma can't work when tgamma is negative! I will edit the code. (EDIT: OK, I fixed the code in the opening post)
yes, I'm aware of that. We can make these tests pass just by increasing the acceptable ULP error, and this solution is enough for me, because (as I see it) And also: all other |
|
Many thanks for the test cases, working on some error plots now.... |
you are welcome. In a few days I will add a couple more |
cosurgi
changed the title
|
OK, some ULP plots coming up.... acos, asin, atan all look good to me, though if we're being picky atan near zero might be improvable: acos:
asin:
atan:
|
|
sin/cos/tan exhibit the usual high errors near the roots, I do wonder if higher than normal precision for the calculation of N * Pi would be enough to (mostly) squash these: cos:
sin:
tan:
|
|
Log looks good, though I might have hoped for smaller errors on something that's so crazy-well-conditioned:
exp as you would expect gets tougher for larger arguments:
|
|
tgamma behaves much as you would expect, especially for larger arguments:
erf looks good, though I wonder if there is low hanging fruit near zero:
erfc has high error for large arguments, but it's ~ e^-x^2 in that domain, so we have the same ill-conditioned issues as exp, but more so:
|
|
Wow, I'm glad to see these plots work well! (Well, other than the couple bugs you had to fix in #417) @cosurgi : Note how these plots give a totally different way of understanding accuracy; @jzmaddock notes that erf at zero is inaccurate, but it's only inaccurate relative to the condition number; in fact in absolute terms it is pretty good. Now, there is another philosophy, adopted by MPFR, that for a representable
This is sort of half-way between the philosophy of MPFR and the backward stability philosophy. I mean, no one would complain if things got a bit more accurate, but at the same time, what is the claimed accuracy in the docs? I believe that the claimed accuracy of |
|
Thank you very much for your reply!
yes, thank you. As I see it, the horizontal axis is the x argument and vertical axis is the ULP error drawn with blue dots. What is the green line though?
I fully agree and understand the problems involved.
You mean for example
Which fixed-point type do you have in mind? I search multiprecision docs but I see no mention of a fixed point type.
yes. That's the important point. Write what kind of accuracy to expect, and all problems are solved. Maybe add in docs a more mention about how to obtain the error estimate in a small domain range in which the user is interested in (not the (-100,100) which I tested here ;). I'm sorry, if I missed this part if it's in docs already. |
The green line is +-|xf'(x)/f(x)|; the condition number of function evaluation.
Hopefully when you work out the analysis you don't necessarily use twice the number, just
There is no fixed point class currently; I'm just stating the mechanism that is used to solve this problem. I've been agitating for fixed point for a while; but also have expended zero effort towards making it happen, so currently I have what I deserve, which is nothing. |
heheh. I also would like to see a fixed point type. We can muster only a finite amount of time though, so we do what we can! |
In fact extending the EDIT: not useful - #264 (comment) , also page 46 in Kahan's paper and page 4 in another Kahan's paper |
Maybe we could exploit symmetry in acos (a shifted oddness in fact, acos(x)=-acos(-x)+pi), because the errors are large close to 1, but are small close to -1.
I guess so. |
|
I'm wondering if Automatic Differentiation in: https://www.boost.org/doc/libs/master/libs/math/doc/html/math_toolkit/autodiff.html |
There was a legitimate effort to move forward on fixed-point in GSoC15. I would like to re-activate this work some day. John McFarlane has done a lot of work on fixed-point and other UDTs in CNL. As far as I know, this work does not specifically target Boost as backend technology, but the author is Boost-aware. |
|
From: Nick <[email protected]>
Sent: 10 August 2020 17:31
To: boostorg/multiprecision <[email protected]>
Cc: Subscribed <[email protected]>
Subject: Re: [boostorg/multiprecision] Incorrect results for a couple of real cpp_bin_float functions:ℝ→ℝ (#264)
What is the green line though?
The green line is +-|xf'(x)/f(x)|; the condition number of function evaluation.
From the graphs, it would appear that the green line is something like a one sigma indicator and two thirds of the values lie within the envelope.
Twice that + and – would seem to be a two-sigma indicator when 95% of values lie inside?
|
We've never intended to compete with mpfr/gmp: those guys have way more resources than we do. Besides they've been there and done that at very considerable effort already. That doesn't mean we can't be better though. |
You do for argument reduction unfortunately: if we assume the input is exact, then when you subtract NPI from the input and you can cancel at most B bits for a B-bit precision number, and therefore need an extra B bits in the value of NPI to fill in the missing data. Of course for inexact input this is a complete waste of time and we have no way of knowing what the missing bits should be. |
I'm not seeing how, what did you have in mind? |
I'm trying to decide whether that's a good way to think about it. In fact the process that generates these values is deterministic and can systematically bias in one way or another, but maybe it's not productive to think in this way. The way I think about is as follows: Given an input which is in error by half an ulp, what is the maximum ulp error computed via Taylor series? |
I had in mind the Taylor series expansion: having derivatives with ULP precision might help. Now I realize however, that there is nothing to use as an input function for the Taylor series, since the math functions are not there yet to use: they are being implemented. |
|
Is it realistic to expect only half an ulp? It's not what we see? What we see looks like a statistical distribution? Computation may be deterministic, but usually involve many instruction steps, each of which may be half a bit 'wrong' up or down. Some may have all 'wrong' ups and others all down, the worst outlier cases. The plots still look like a good way to tell where to focus efforts to improve - and even more important, perhaps, where not to bother to try. |
No, not half an ulp output, the assumption is half-ulp accurate input. |
|
Ah I see - yes half is expected |
|
Speaking of Kahan (12 posts up), are you using I have checked that it indeed helps with summations. |
|
@cosurgi : Kahan summation can be thought of as an emulation of twice the working precision. Most of the time it makes sense to just use twice the precision instead of using Kahan summation; e.g., it takes ~ IMO, the only goal should be error bound by the condition number under the assumption of half-ulp accurate inputs, making Kahan a bit overkill. Again, maybe that's a sensible solution for parts of certain algorithms, but I don't see it helping in general. |
jzmaddock
added a commit
to boostorg/math
that referenced
this issue
Aug 13, 2020
|
The PR mentioned above plucks a bit or low hanging fruit from erf, the error plot now looks like:
|
|
OK, a couple days ago the scanning of 2e7 evaluations was finished, using this code: https://gitlab.com/yade-dev/trunk/-/blob/master/py/high-precision/_RealHPDiagnostics.cpp#L445 This post is quite long, because I am posting more test cases: the worst arguments found. The loop was quite simple (linked above): whenever a larger ULP error was found, the offending arguments were added to the list. So it's not like I saved all which produced ULP error above threshold, just the largest found so far in the loop. I will be posting here only the tail of this list, because I suppose the start won't be that useful, but I can send you more if you want. All numbers have 207 bits of cpp_bin_float<62>. These are my results:
[*] - star means that in the first run I found a higher error than reproduced in the code below, but didn't store arguments back then. I could reproduce it, if you need these arguments. I would need to rerun with partially reverted yade commits 67786dfc77db71bcfd29961db0bf1007ea7f3ab0 and 3ffc9d82c6a878bf938efe3b48d8056fef25bb3e (that's just a note for future if needed later). Also in the code below I didn't add arguments for functions that had error < 4. I can provide them if you want. So this is the output of the test code provided below: The test code: |
|
This PR: #270 addresses a large chunk of this issue by improving the sin/cos/tan/asin/acos/atan functions. ULP plots now look like this:
Note that argument reduction in sin/cos/tan is still capped to 1/epsilon as beyond that the calculation fails without performing an expensive further-extended division. In addition to support arguments up to 2^N, we need > 2^N bit precision during argument reduction - and since the exponent type could be a 64 bit type - that would mean 2^64 bit temporaries for the reduction which would very likely take "forever" to process, even if the machine has enough memory :( @NAThompson : I did try making too large arguments a hard error with an exception being raised, but it broke the tanh-sinh tests - I had no idea that these were calling trig functions with such large arguments - even though the integral is no doubt zero in that region. |
Oh yeah; that's how tanh-sinh works. I'm not 100% convinced that the integrands will be zero there since we have so many domain transformations. But nonetheless it looks like this'll be a big improvement. |
|
How about using |
No that would make things far worse: currently we're at about 10eps accuracy near x = 1. What you can't see from the plot above is that very close to 1 we use an expansion in terms of 2F1 which is actually deadly accurate.... but leaving that aside, the result for |
|
Ah, I get it. The ulp errors are large because everything is very close to zero around 1. So the solution would be to use something else than Newton's method. I can't recall now anything that would help with this epsilon ill-conditioning here. Hmm.. something of higher order maybe... |
12 tasks
|
This should be addressed in the last release. |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
As requested I am opening another issue for problems with following functions, but only for cpp_bin_float real valued functions:
sin,cos,tanwith 7e7 ULP erroracoswith 10000 ULP errorerfcwith 20000 ULP errorlgammawith 70000 ULP error (edited)tgammawith 10000 ULP errorfmawith 2e5 ULP errorI prepared
std::vector<std::string> suspects = {……}so that only these arguments are tested. So I'm sorry that this test program is a little longer. I can push it somewhere via git if you want.So I do a domain check in the code below, but that's just the domain check via mathematical definition of the function. E.g. acos in (-1,1) range. But if there are some wider ill-conditioned
cpp_bin_floatranges then you tell me what are those :)The last couple tests are done on MPFR types, to show that they have no problems, even when I reduce the tolerance to 4 ULP.
This test is done on g++-10 and Boost 1.71
EDIT: Here's my (fixed lgamma domain) output:
The text was updated successfully, but these errors were encountered: