Date: Tue, 23 Oct 2012 03:09:25 +1100 (EST) From: Bruce Evans <brde@optusnet.com.au> To: Warner Losh <imp@bsdimp.com> Cc: svn-src-head@FreeBSD.org, svn-src-all@FreeBSD.org, src-committers@FreeBSD.org, Warner Losh <imp@FreeBSD.org>, Bruce Evans <brde@optusnet.com.au> Subject: Re: svn commit: r241755 - head/lib/msun/src Message-ID: <20121023011459.G2088@besplex.bde.org> In-Reply-To: <80A9038D-C84B-47D3-A137-F281449F4D88@bsdimp.com> References: <201210192246.q9JMkm4R092929@svn.freebsd.org> <20121020150917.I1095@besplex.bde.org> <18177777-6EE0-4103-98B0-272EFF98FE96@bsdimp.com> <20121022213348.T1373@besplex.bde.org> <80A9038D-C84B-47D3-A137-F281449F4D88@bsdimp.com>
next in thread | previous in thread | raw e-mail | index | archive | help
On Mon, 22 Oct 2012, Warner Losh wrote: > On Oct 22, 2012, at 5:14 AM, Bruce Evans wrote: > >> On Sun, 21 Oct 2012, Warner Losh wrote: >> >>> Feel free to fix them however. I added the comments because the algorithms weren't quite the same... If you have a better way, feel free to back my stuff out on the way to it. >> >> But the algorithms are identical to a fault. Inside the functions, all >> lines except 1 in each correspond exactly, and the exception is a style >> bug. Only about 30 lines in each are not lexically identical. The >> non-lexical differences are for things like different magic numbers. > > Except that's not true. For expf, only the first two terms of the 5 are computed in Remes expansion. That's why I bothered to document the difference. For logf, there are at least two additional terms in the intermediate form. No, they both use a Remez-type polynomial approximation. The only differences are that the approximation for logf needs _fewer_ terms, and in fact has 3 fewer (except in the Cygnus version it had the same number of terms with the lower ones useless at best), and that the grouping of the terms is different for efficiency. I say "Remez-type" because it is unclear if anyone actually uses pure Remez. I use my own idea of what it is, without having looked at documentation of the original algorithm. I probably have different tweaks than normal. fdlibm's approximations often seem to be too inaccurate to have been generated by full Remez. The ones here in the Cygnus version were generated by starting with the double precision coeffs and blindly rounding them to float precision. This is a bad non-Remez method (but pure Remez doesn't understand limited precision AFAIK). In practice it gave 6 extra bits using Lg1-Lg7, where my coeffs generated for float precision give 10 extra bits using only Lg1-Lg4. The extra coeffs tend to be actively harmful (and IIRC they were here), because when Remez generates the lower coeffs it tunes them for the lower coeffs actually working, all in infinite precision, and after rounding to finite precision the delicate tuning for this just doesn't work. The comments actually say "special Remes". Who knows what that is? :-) Even running same algorithm, there is a lot of noise in the lower bits for all except the first few coeffs, so it is hard to duplicate coeffs produced by other special versions. > The differences here are: > /* logf */ > t1= w*(Lg2+w*Lg4); > t2= z*(Lg1+w*Lg3); > /* log */ > t1= w*(Lg2+w*(Lg4+w*Lg6)); > t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7))); > > which to even the most casual observer are clearly different. Too casual :-). Both use a moderately efficient grouping of terms for efficiency. Both are more optimal than the simple Horner's method grouping which would look similar even to a too-casual observer: R = s*s*(Lg1+s*(Lg2+s*(...+Lg[4or7])...)) Neither is highly optimized, but the logf one is is better, mainly because a power of 2 number of terms is easiest to get right and a small number of terms is especially easy to get right. In fact, the logf one is just the log one with higher terms omitted. Simply omitting higher terms is usually too simple, but here it happens to improve the grouping. The grouping is a general polynomial evaluation technique and isn't part of Remez. The details of it aren't documented at all. In fact, the pseudo-code doesn't even mention z completely: it says R(z) ~= Lg1*s**2 + ... +Lg7*s**14 (using superscripts for the powers) where the code says t2 = z*(Lg1+w*(Lg3+...)... R = t2+t1. Many other details involve unmentioned parts of general polynymonal evaluation technique (mainly, add terms from high to low for accuracy, but violate this rule as much as possible for efficiency, and avoid underflow by filtering out small args early and by not calculating high powers directly). >> The old fdlibm comments aren't too careful about keeping magic numbers >> out of the algorithm description so that the algorithm description is >> as general as possible. The precise magic numbers are often critical >> to the details of the implementation of the algorithm but not really >> to the algorithm itself. > > The current code isn't careful at all about magic numbers. Which ones are magic hex constants for IEEE stuff, and which ones are hex numbers for the exact representation of important constants? None for the hex FP constants. Hex FP constants are for the programmer's convenience, but I are just style bugs where I use put them in e_expf.c. Decimal FP constants work with any C compiler. They are only pseudocode in comments in fdlibm, and no one every used their values from there AFAIK, since C compilers were common when fdlibm was released. The style bug is to mix C99 hex FP constants in old fdlibm (Cygnus) code. Though only in comments, it looks strange since the active code still laboriously spells out (float) casts to get float constants. Even C90 has float constants via the F suffix. But I recently decided that many of the FP constants are style bugs. When you just want to add 2 as in e_log*.c you should spell 2 as 2 and not as 2.0 for the double case, 2.0F or (float)2.0 for the float case, and maybe 2.0L or (long double)2.0 for the long double case. C's promotion rules work correctly, so the 2 in (2+f) is promoted to the same type as f in all cases. By depending in this, the code becomes independent of the precision. A more interesting case is multiplying by 0.5 or (float)0.5 or 0.5F or ... This can be written as division byi integer 2, since the compiler is permitted to turn this into multiplication by 0.5 with the correct type. gcc has done this for many uears, so there is no need to manually optimize this. > Anyway, it is clear you guys don't want them in there so I've reverted them. I'm totally baffled by this request because these algorithms are similar not identical, but since you guys are the active maintainers, I'll accede to your wishes after expressing my utter bafflement at the justifications. Well, we (I) wrote the changes in the polynomial evaluation for both e_expf.c and e_logf.c, and was careful to a fault not to change the algorithm, so I know that it is more identical in these files than in others. Bruce
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20121023011459.G2088>