Date: Wed, 26 Sep 2012 08:16:05 +1000 (EST) From: Bruce Evans <brde@optusnet.com.au> To: Steve Kargl <sgk@troutmask.apl.washington.edu> Cc: numerics@FreeBSD.org Subject: Re: [PATCH] implementation for expm1l() Message-ID: <20120926071841.T3491@besplex.bde.org> In-Reply-To: <20120925202228.GA31099@troutmask.apl.washington.edu> References: <20120925154606.GA28919@troutmask.apl.washington.edu> <20120926021900.R2081@besplex.bde.org> <20120925202228.GA31099@troutmask.apl.washington.edu>
next in thread | previous in thread | raw e-mail | index | archive | help
On Tue, 25 Sep 2012, Steve Kargl wrote: > On Wed, Sep 26, 2012 at 04:39:36AM +1000, Bruce Evans wrote: >> On Tue, 25 Sep 2012, Steve Kargl wrote: > ... >>> 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. > > The fuzzy ones are left over from an attempt to use a 2nd lookup > table on the range [-0.25:0.25]. This has the nice property that > T2 - T1 = 0.5 and dividing this up into 2**n segments gives a > delete of the form 0x1.0p-7. > > The log(1+-1/4) expressions are adopted from Tang's paper where he > gives an error analysis to establish an upper bound on the max > ULP. Do the intervals have anything to do with this polynomial? As I understand it (incompletely), we just use a special approximation for |x| < ~0.25 where the one for exp(x) is not good enough, and otherwise we use (exp(x) - 1) manually inlined. Anyway, I regenerated the polynomial for Tang's range. It is almost 1 bit more accurate. >>> + >>> +/* >>> + * 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.) > > The actual function that I approximated was (expm1(x)-1-x-x**2/2)/x**3. > Tang's algorithm assumes that one compute x+x**2/2 with extra precision. > He does this by decomposing x into high and low parts. I suppose I > should update the comment to reflect the function I approximate. I should have done that too, but it complicates understanding and calculating the relative error. The result is near x, so the relative error is obtained by dividing by nearly x. This happens naturally for (exp(x)-1)/x. For your function (fixed) of (exp(x)-1-x-x**2/2)/x**3, the first term is 1/6 and I think the error must be divided by 1/6 -- minimizing the error gives the same result, but this division is needed to give the actual error. What I actually do gives an approximation with first to terms exactly 1 and 0.5, so it gives the same result. I check this manually. My algorithm might prefer to change these constants, so I use special code to force them for several functions that want them for the first 2 coeffs. I do this for exmp1 for the first coeff only, although I do it for the second coeff for most functions where it applies. >> I generated the following polynomial. Using only 2 full-precision coeffs >> limit the accuracy, so B14 doesn't improve accuracy significantly. > > (coefficients deleted) Now using Tang's interval: @ P0 = 1.00000000000000000000e0L, /* 0x8000000000000000.0p-63L */ @ P1 = 5.00000000000000000000e-1L, /* 0x8000000000000000.0p-64L */ @ P2 = 1.66666666666666666712e-1L, /* 0xaaaaaaaaaaaaaaae.0p-66L */ The hex formats are designed for showing the mantissa bits unshifted, but bits after the point would work equally well and would obfuscate the exponent less. E.g., 1.0 shouldn't be expressed as 0x1.00...p0L since this shifts the bits, but it can be 0x0.8000000000000000p1L, which can be read as 0.5*2**1. @ P2hi = 1.6666666666666666e-1, /* 0x15555555555555.0p-55 */ @ P2lo = 9.2970336290632005e-18, /* 0x15700000000000.0p-109 */ @ P3 = 4.16666666666666666678e-2L, /* 0xaaaaaaaaaaaaaaab.0p-68L */ @ P3hi = 4.1666666666666664e-2, /* 0x15555555555555.0p-57 */ @ P3lo = 2.3140940118987485e-18, /* 0x15580000000000.0p-111 */ @ P4 = 8.3333333333333211e-3, /* 0x1111111111110a.0p-59 */ @ P5 = 1.3888888888888794e-3, /* 0x16c16c16c16beb.0p-62 */ @ P6 = 1.9841269841376342e-4, /* 0x1a01a01a023996.0p-65 */ @ P7 = 2.4801587302997418e-5, /* 0x1a01a01a07f9db.0p-68 */ @ P8 = 2.7557318820808094e-6, /* 0x171de39faa2d52.0p-71 */ @ P9 = 2.7557312008157410e-7, /* 0x127e4f663a69b6.0p-74 */ @ P10 = 2.5052802881898768e-8, /* 0x1ae67646e4deb3.0p-78 */ @ P11 = 2.0892043051363171e-9, /* 0x11f235775d156c.0p-81 */ @ P12 = 1.5615878424451485e-10, /* 0x15765948fa3f05.0p-85 */ @ @ 6.857e-23 |"''''''"''''''''''''''''''''''''''''''''''''''''x'''''''''''''" @ |: x " _ : @ |: _ : "_ : : @ |: : : " _ : _ : @ |: : : _ : : : : @ |:: : _ _ : : x :" : @ |:: : : : : : : : : : @ |:: : : : : _ : : : : @ |:: _ : : " : : " : : @ |:x : : _ : : : : : : @ : : : : : : : x : : : @ :`:```:````"``````:``````x`````````""````````_``````:````:```::` @ : : : : : : " " : : : ::| @ : : : : : : _ _ : : : _:| @ : : : : " _ : : _ ::| @ : : : : : : _ _ " _ : ::| @ : : " x : : : : ::| @ : : : : : x _ _ _ : : :| @ : : : : " _ : : :| @ : ": : "_" " : _ :| @ _ : " _ " :| @ -6.64e-23 |..."........._........................................_......_| @ -0.288 0.224 @ adj cp = 64 range ~[-6.6356e-23, 6.8973e-23] 2**-73.618 When choosing endpoints, I just use decimal and round up the last place. Experience shows that using a slightly wider than necessary interval doesn't allow using fewer coeffs. Even rounding down is usually harmless, since the endpoints shouldn't be poles so everything is continuous and the approximation doesn't diverge very far for a small rounding error. Now for ld128: # P0 = 1.0000000000000000000000000000000000e0L, /* 0x10000000000000000000000000000.0p-112L */ Now 1.0 is well expressed as 0x1.0000000000000000000000000000.0p0L. Similarly for d64 (the number of mantissa bits is 4n+1). # 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.3888888888888888888888888894888658e-3L, /* 0x16c16c16c16c16c16c16c16f2191f.0p-122L */ # P6 = 1.9841269841269841269841269038522581e-4L, /* 0x1a01a01a01a01a01a01a005a77f70.0p-125L */ # P7 = 2.4801587301587301587301479770823380e-5L, /* 0x1a01a01a01a01a01a01979505fa4a.0p-128L */ # P8 = 2.7557319223985890652565291171436921e-6L, /* 0x171de3a556c7338faae15d1232a5c.0p-131L */ # P9 = 2.7557319223985890653527586372574888e-7L, /* 0x127e4fb7789f5c72fb0453795a335.0p-134L */ # P10 = 2.5052108385441718734295029309253960e-8L, /* 0x1ae64567f544e38cc000f34e989ac.0p-138L */ # P11 = 2.0876756987868093985249237137054359e-9L, /* 0x11eed8eff8d896802862ae8d8dbd9.0p-141L */ # P12 = 1.6059043836821721373405224365590722e-10L, /* 0x16124613a86d32cb06578b4f2c75d.0p-145L */ # P13 = 1.1470745597745053e-11, /* 0x193974a8c0a1ad.0p-89 */ # P14 = 7.6471637317124909e-13, /* 0x1ae7f3e73218ed.0p-93 */ # P15 = 4.7794773046283608e-14, /* 0x1ae7f3e4948c15.0p-97 */ # P16 = 2.8114571580397259e-15, /* 0x1952c761a306e5.0p-101 */ # P17 = 1.5619480336357893e-16, /* 0x168292355a0a63.0p-105 */ # P18 = 8.2237292848800199e-18, /* 0x12f66ed54ff940.0p-109 */ # P19 = 3.9981067521176347e-19, /* 0x1d8035cafa761f.0p-114 */ # # 1.048e-37 _'''''_'''''''''''''''''''''''''''''''''''''''''''''''''''''''"| # : x : :| # : : :: _ :| # : : _ : _ " :| # : : : : _x x : ": :| # : :: : _ :: : : x_ : : ::| # |: : : : : : : x x : : " " ::| # |: :: : : _ : : : : : : ::| # |: :: : : : x : x " " : : : ::| # |: ": : _ : : _ : : : : x ::| # |: " : : : : : x"x " : : : x : x:| # `_``````:```:```:```:``````:```"```""``````_```:````:```:`:``:`: # | : : : : x _ : : " : : : : # | : : " x __ : " : : : : : # | _ : : : : : : : : : : # | : : : : " : : : : : : # | : _ : : x : _ :: : # | : : _ _ " : : :: : # | : : ": x: : # | :: x : :: : # | :: " : : # -1.01e-37 |........"_................................................._.." # -0.288 0.224 # adj cp = 113 range ~[-1.0660e-37, 1.0660e-37] 2**-122.82 ld128 is much more resistant to reducing to double precision -- it needs more than half of the coeffs to have full precision. Double precision for P11 onwards is unusable. Double precision for P12 onwards loses a couple of bits. There are a couple to spare, but I prefer not to reduce if the result loses many bits compared with infinite precision. Bruce
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20120926071841.T3491>