Correct the range of the log1p polynomial - #1432
Open
ibmibmibm wants to merge 1 commit into
Open
Conversation
The log1p function calculated a Taylor polynomial of log(1 + x) for all values of x in the range (-1, 0.5]. The test had an upper limit at x = 0.5, but it had no lower limit. The coefficient tables were not long enough for this range. The function gave a result with too few correct digits, and it reported no error. At x = -0.99, the relative error was 34 percent for decimal32_t. It was 25 percent for decimal64_t, and 16 percent for decimal128_t. All six types had this problem, because the fast types hold a copy of the same values. The tables now hold the coefficients 2/(2k+1) of the series 2 * atanh(w), with w = x / (2 + x). The counts are 7, 16 and 34 terms, in place of 12, 20 and 36. For |x| not more than 1/2, the value of |w| is not more than 1/3. A test of abs(x) gives the lower limit that the code did not have. For all larger values of |x|, the code calculates log(x + 1). For |x| less than 10^-digits10, the function now returns x. The sum of the terms after x is less than one half of the last digit of x. This branch also keeps w * w away from the exponents where fma gives an incorrect result for decimal32_t and decimal64_t. The exponents are 1e-21 to 1e-31 for decimal32_t, and 1e-42 to 1e-62 for decimal64_t. A test with MPFR at 1610 points of the branch gives these results: decimal32_t 3.37e+05 eps -> 1.02 eps, 551 ns -> 478 ns decimal64_t 2.52e+14 eps -> 0.88 eps, 3299 ns -> 2613 ns decimal128_t 1.61e+32 eps -> 0.91 eps, 5291 ns -> 5291 ns The largest error is now at x = -0.64. This point is in the branch that calculates log(x + 1), thus the value is the error of the log function. The atanh, asinh and acosh functions call log1p. They get the correction with no change. For decimal128_t, the error of atanh at x = 0.49 decreases from 5.13e+19 eps to 0.23 eps. The error of acosh at x = 1.08 decreases from 6.45e+19 eps to 0.41 eps. The old tests for these four functions used decimal32_t and compared with float. A test of this form cannot find an error that is smaller than the epsilon of float. The old tests also did not use values of x that are less than -0.375. Each test now has a table of control values with 36 digits, for all six types. The table for log1p has 32 points, and four of them are near zero. Each new test fails with the old code. Fixes boostorg#1430.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
The log1p function calculated a Taylor polynomial of log(1 + x) for all values of x in the range (-1, 0.5]. The test had an upper limit at x = 0.5, but it had no lower limit. The coefficient tables were not long enough for this range.
The function gave a result with too few correct digits, and it reported no error. At x = -0.99, the relative error was 34 percent for decimal32_t. It was 25 percent for decimal64_t, and 16 percent for decimal128_t. All six types had this problem, because the fast types hold a copy of the same values.
The tables now hold the coefficients 2/(2k+1) of the series 2 * atanh(w), with w = x / (2 + x). The counts are 7, 16 and 34 terms, in place of 12, 20 and 36. For |x| not more than 1/2, the value of |w| is not more than 1/3. A test of abs(x) gives the lower limit that the code did not have. For all larger values of |x|, the code calculates log(x + 1).
For |x| less than 10^-digits10, the function now returns x. The sum of the terms after x is less than one half of the last digit of x. This branch also keeps w * w away from the exponents where fma gives an incorrect result for decimal32_t and decimal64_t. The exponents are 1e-21 to 1e-31 for decimal32_t, and 1e-42 to 1e-62 for decimal64_t.
A test with MPFR at 1610 points of the branch gives these results:
decimal32_t 3.37e+05 eps -> 1.02 eps, 551 ns -> 478 ns
decimal64_t 2.52e+14 eps -> 0.88 eps, 3299 ns -> 2613 ns
decimal128_t 1.61e+32 eps -> 0.91 eps, 5291 ns -> 5291 ns
The largest error is now at x = -0.64. This point is in the branch that calculates log(x + 1), thus the value is the error of the log function.
The atanh, asinh and acosh functions call log1p. They get the correction with no change. For decimal128_t, the error of atanh at x = 0.49 decreases from 5.13e+19 eps to 0.23 eps. The error of acosh at x = 1.08 decreases from 6.45e+19 eps to 0.41 eps.
The old tests for these four functions used decimal32_t and compared with float. A test of this form cannot find an error that is smaller than the epsilon of float. The old tests also did not use values of x that are less than -0.375. Each test now has a table of control values with 36 digits, for all six types. The table for log1p has 32 points, and four of them are near zero. Each new test fails with the old code.
Fixes #1430.