[PATCH] implementation for expm1l()
Bruce Evans
brde at optusnet.com.au
Tue Sep 25 18:39:47 UTC 2012
On Tue, 25 Sep 2012, Steve Kargl wrote:
> The patch following my .sig implements expm1l(). An example
> of testing 10M values uniformly distributed in each interval
> gives
>
> amd64:
> % ./testl -n 10
> Interval tested: [-11355.000000:-0.250000]
> Max ULP: 0.50325
> x Max ULP: -5.21316172131617162e+00, -0x1.4da4710f73f83p+2
> Interval tested: [-0.250000:0.250000]
> Max ULP: 0.53926
> x Max ULP: -1.15709936570993657e-01, -0x1.d9f2a996507ec308p-4
> Interval tested: [0.250000:11356.000000]
> Max ULP: 0.50584
> x Max ULP: 8.79177827642782764e+03, 0x1.12be39e8fde623b2p+13
>
> sparc64 (only very limited testing due to speed of flame):
> % ./testl -n 10
> Interval tested: [-11355.000000:-0.287682]
> Max ULP: 0.50165
> x Max ULP: -5.34961332997700764996445801944829055e+00,
> -0x1.566010969fcd47a521247b4128p+2
> Interval tested: [-0.287682:0.223144]
> Max ULP: 0.54843
> x Max ULP: -2.86792265208762197924244858751372070e-01,
> -0x1.25acdf1f45027985d7c04fa34168p-2
> Interval tested: [0.223144:11356.000000]
> Max ULP: 0.50577
> x Max ULP: 7.83174792710817081424055325173740399e+03,
> 0x1.e97bf7826a562b03aafc33340393p+12
>
> Now, the problem! On i386 in the interval [-0.287682:0.223144],
> I see up to 15.5 ULP for the ld80 patch. Note, if I don't
ulps or ULPs
> use ENTERI and RETURNI in the ld80 code, I see ULP that are
> consistent with the amd64 results. This of course leads back
> to the 53- versus 64-bit precision settings on i386.
I thought I fixed it, and you have all my patches, but...
It can't possibly work better without ENTERI. Maybe a bug in the
test code.
You should test the internal values before they are added together,
and expect an error of < ~1/256 ulps (or ~0.5 ulps in 72-bit precision).
I thought I did that. Maybe I only tested amd64. Now that I look at
it again, I see that it wasn't worth testing on i386 since it can't
possibly work there...
> I have tried re-arranging the order of evaluation of the polynomial
> and some of the intermediate results with little improvement. I've
> also increased the order of the polynomial. Everything I've tried
> seems to fail on i386. So, anyone got a suggestion?
>
> PS: the relevant section of code is in the 'if (T1 < x && x < T2)'
> block of code.
>
> --
> Steve
>
> Index: ld80/s_expl.c
> ===================================================================
> --- ld80/s_expl.c (revision 240889)
> +++ ld80/s_expl.c (working copy)
> @@ -301,3 +301,149 @@
> ...
> +#if 0
> +static const long double
> +T1 = -0x1.269621134db92786p-2L, /* log(1-1/4) = -2.876820724517809275e-01L */
> +T2 = 0x1.c8ff7c79a9a21ac4p-3L; /* log(1+1/4) = 2.231435513142097558e-01L */
> +#else
> +static const long double
> +T1 = -0.25L,
> +T2 = 0.25L;
> +#endif
Apparently these thresholds aren't critical (since fuzzy ones work on amd64).
The fuzzy ones should be declared double.
> +
> +/*
> + * Remes polynomial coefficients on the interval [T1:T2] with an error
> + * estimate of log2(|expm1(x)-p(x)|) ~ -74. These are the coefficients
It's expm1(x)/x that must be approximated by p(x)/x, not expm1(x) by
p(x). This combined with the nonstandard comment style confused me for
a while. (Plain expm1(x) is just as easy to approximate as exp(x) and
doesn't require any long double coefficients.)
> + * rounded to long double precision where B5 through B14 can be used with
> + * a double type.
> + */
> +static const long double
> +B3 = 0x1.5555555555555556p-3L, /* 1.66666666666666666671e-01 */
> +B4 = 0x1.5555555555555556p-5L; /* 4.16666666666666666678e-02 */
These 2 aren't actually rounded to long double precision on i386, so the
point of having them in more precision than the others is completely lost.
gcc rounds them to nearest, so the result is the same as:
static const non_long double
B3 = 0x1.5555555555555800p-3, /* some wrong value */
B4 = 0x1.5555555555555800p-5; /* some wrong value */
This is the only obvious problem.
For a quick check that this is the only problem, try compiling with clang.
I think it is incompatible with gcc here and rounds to 64 bits.
LD80C must be used for most long double constants in ld80 (unless the
constants can be double).
I generated the following polynomial. Using only 2 full-precision coeffs
limit the accuracy, so B14 doesn't improve accuracy significantly.
@ P0 = 1.00000000000000000000e0L, /* 0x8000000000000000.0p-63L */
@ P1 = 5.00000000000000000000e-1L, /* 0x8000000000000000.0p-64L */
P0 and P1 are implicit and/or should be manually reduced to doubles.
The above was automatically generated in a standard format that is too
general for them.
@ P2 = 1.66666666666666666671e-1L, /* 0xaaaaaaaaaaaaaaab.0p-66L */
@ P2hi = 1.6666666666666666e-1, /* 0x15555555555555.0p-55 */
@ P2lo = 9.2563760475949941e-18, /* 0x15580000000000.0p-109 */
@ P3 = 4.16666666666666670913e-2L, /* 0xaaaaaaaaaaaaab28.0p-68L */
@ P3hi = 4.1666666666666664e-2, /* 0x15555555555555.0p-57 */
@ P3lo = 2.7376104855258987e-18, /* 0x19400000000000.0p-111 */
This is still automatically generated for the old pre-LD80C() format in
which P2 is written as (P2hi + P2lo) and some manual editing was needed
to put P2lo in the volatile section. The hex values are perhaps not
quite write for manyally copying them into LD80().
@ P4 = 8.3333333333333332e-3, /* 0x11111111111111.0p-59 */
@ P5 = 1.3888888888888057e-3, /* 0x16c16c16c16a97.0p-62 */
@ P6 = 1.9841269841266616e-4, /* 0x1a01a01a019b74.0p-65 */
@ P7 = 2.4801587307106348e-5, /* 0x1a01a01a1a7b21.0p-68 */
@ P8 = 2.7557319244569047e-6, /* 0x171de3a5a0efd0.0p-71 */
@ P9 = 2.7557302872739145e-7, /* 0x127e4eff5f5ccd.0p-74 */
@ P10 = 2.5052063395588939e-8, /* 0x1ae6423d7e31ae.0p-78 */
@ P11 = 2.0899093104590682e-9, /* 0x11f3c259569a3b.0p-81 */
@ P12 = 1.6096692806293200e-10, /* 0x161f8531c9dc24.0p-85 */
The others are in a standard format.
Note that my P12 is your B13 (different by 1 for dividing by x in
expm1(x)/x).
@
@ 1.147e-22 |"'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''|
@ |: |
@ |: |
@ |:: |
@ |:: "" |
@ |:: _x |
@ : " x " _ |
@ : : : : " : |
@ : : : : : __ |
@ : : : " """_ " " :: |
@ : : " : " x : : : |
@ :,,:,,:,,,,,:,,,,,,,x,,,,,,"x___,,,,,,,,,,,,,_,,,,,,:,,,,,_,,:,,
@ _ : : x ""x : x : : |
@ | " : : x " : : : _ |
@ | " : " " : : : |
@ | _ x x " : : : |
@ | x x " " : |
@ | x " x _ : :|
@ | x" "x _ : :|
@ | " :|
@ | " :|
@ -9.54e-23 |............................................................._"
@ -0.25 0.25
@ adj cp = 64 range ~[-1.1479e-22, 1.1501e-22] 2**-72.881
> Index: ld128/s_expl.c
> ===================================================================
> --- ld128/s_expl.c (revision 240889)
> +++ ld128/s_expl.c (working copy)
> @@ -259,3 +259,132 @@
> return (t * twopkp10000 * twom10000);
> }
> }
> +
> +static const long double
> +T1 = -2.87682072451780927439219005993827491e-01L, /* log(1-1/4) */
> +T2 = 2.23143551314209755766295090309834524e-01L; /* log(1+1/4) */
Now I see that these aren't critical. But isn't T1 = -0.25 not far
enough negative? Probably T2 being larger than necessary compensates.
These should probably be double in all cases.
> +
> +/*
> + * Remes polynomial coefficients on the interval [T1:T2] with an error
> + * estimate of log2(|expm1(x)-p(x)|) ~ -118.
> + */
> +static const long double
> +B3 = 1.6666666666666666666666666666666666467e-01L,
> ...
> +B20 = 4.0010880199748106995049690682517608541e-19L;
I get 123 bits with the same number of terms (P19 corresponds to B20).
P0 = 1.0000000000000000000000000000000000e0L, /* 0x10000000000000000000000000000.0p-112L */
P1 = 5.0000000000000000000000000000000000e-1L, /* 0x10000000000000000000000000000.0p-113L */
P2 = 1.6666666666666666666666666666666661e-1L, /* 0x15555555555555555555555555553.0p-115L */
P3 = 4.1666666666666666666666666666665461e-2L, /* 0x1555555555555555555555555548d.0p-117L */
P4 = 8.3333333333333333333333333333702524e-3L, /* 0x111111111111111111111111170ea.0p-119L */
P5 = 1.3888888888888888888888888894097559e-3L, /* 0x16c16c16c16c16c16c16c16ebae13.0p-122L */
P6 = 1.9841269841269841269841269028212328e-4L, /* 0x1a01a01a01a01a01a01a005649497.0p-125L */
P7 = 2.4801587301587301587301496436041811e-5L, /* 0x1a01a01a01a01a01a0198e708bd2f.0p-128L */
P8 = 2.7557319223985890652565668610907848e-6L, /* 0x171de3a556c7338faae2dbd716389.0p-131L */
P9 = 2.7557319223985890653389692780995245e-7L, /* 0x127e4fb7789f5c72f94f52e6db855.0p-134L */
P10 = 2.5052108385441718730021176147351305e-8L, /* 0x1ae64567f544e38c6b59f0b86d695.0p-138L */
P11 = 2.0876756987868094553775599688331605e-9L, /* 0x11eed8eff8d896a3590b97ce89c87.0p-141L */
@ P12 = 1.6059043836821743738564862857627589e-10L, /* 0x16124613a86d3b71f680b5698fe4d.0p-145L */
@ P13 = 1.1470745597743859829054288214163328e-11L, /* 0x193974a8c09eca572389f431d83ef.0p-149L */
@ P14 = 7.6471637316522595640322188859803229e-13L, /* 0x1ae7f3e7312feabddd5c9a9a658ec.0p-153L */
@ P15 = 4.7794773056701410114955405123070665e-14L, /* 0x1ae7f3e4adbc63a0143e13d8c8b13.0p-157L */
@ P16 = 2.8114572392723202593376369659548583e-15L, /* 0x1952c76de98f8b161d8d0e78ecd0a.0p-161L */
@ P17 = 1.5619481870532875016483016145454221e-16L, /* 0x1682925a721faed3a1323d09bba02.0p-165L */
@ P18 = 8.2232942744013145419986831227112038e-18L, /* 0x12f62d18a8dc6089b1fa968c8c030.0p-169L */
@ P19 = 3.9923442242776910307469064711618566e-19L, /* 0x1d75533495953f6181159d629556a.0p-174L */
@
@ 7.3e-38 "'''''x'''''''_''''''''''''''''''''''''''''''''''''''''''''''''_
@ : : ": x " :
@ : :: :: "" _ : :
@ : " :" : : _ : x: x :
@ |:: : : : : " " : : :: ::|
@ |:: : : : x : : : : : ::|
@ |:: : : _ : _ : : : : : ::|
@ |:: : : : : : x _ : : : ::|
@ |": x : : : x : : : x : _ :x|
@ | : : : : : " x x : : : : : : |
@ ,,:,,:,,:,,,:,,,:,,,,,:,",,,,,,,","__,,,,,,,_,:,,,,:,,,:,,:,,:,,
@ | : : : : : : : _ x : : : : |
@ | : : _ : x x : : x : " : : |
@ | : : : : : : x : : : : : |
@ | :: : : : : : x : : :: |
@ | :: : " : : : : _ : :: |
@ | :: : : : : _ : : : :: |
@ | :: : : : " " : : :: |
@ | _: :: _ : __ :: :" |
@ | : :: : : :x x |
@ | " x: :_ x |
@ -7.32e-38 |........."......._............................................|
@ -0.288 0.224
@ adj cp = 113 range ~[-7.4093e-38, 7.3114e-38] 2**-123.34
Next I will apply coefficient reduction to double. This is a much large
rptimizations for ld128 (for the full optimizations, lower terms of
the poly (most of the ones using double coeffs) should also use the
arg reduced to double). It will take a while to search for the best
number of double coeffs.
Bruce
More information about the freebsd-numerics
mailing list