path: root/src/math/hypot.c
AgeCommit message (Collapse)AuthorLines
2015-11-21math: explicitly promote expressions to excess-precision typesRich Felker-2/+2
a conforming compiler for an arch with excess precision floating point (FLT_EVAL_METHOD!=0; presently i386 is the only such arch supported) computes all intermediate results in the types float_t and double_t rather than the nominal type of the expression. some incorrect compilers, however, only keep excess precision in registers, and convert down to the nominal type when spilling intermediate results to memory, yielding unpredictable results that depend on the compiler's choices of what/when to spill. in particular, this happens on old gcc versions with -ffloat-store, which we need in order to work around bugs where the compiler wrongly keeps explicitly-dropped excess precision. by explicitly converting to double_t where expressions are expected be be evaluated in double_t precision, we can avoid depending on the compiler to get types correct when spilling; the nominal and intermediate precision now match. this commit should not change the code generated by correct compilers, or by old ones on non-i386 archs where double_t is defined as double. this fixes a serious bug in argument reduction observed on i386 with gcc 4.2: for values of x outside the unit circle, sin(x) was producing results outside the interval [-1,1]. changes made in commit 0ce946cf808274c2d6e5419b139e130c8ad4bd30 were likely responsible for breaking compatibility with this and other old gcc versions. patch by Szabolcs Nagy.
2013-09-05math: rewrite hypotSzabolcs Nagy-115/+59
method: if there is a large difference between the scale of x and y then the larger magnitude dominates, otherwise reduce x,y so the argument of sqrt (x*x+y*y) does not overflow or underflow and calculate the argument precisely using exact multiplication. If the argument has less error than 1/sqrt(2) ~ 0.7 ulp, then the result has less error than 1 ulp in nearest rounding mode. the original fdlibm method was the same, except it used bit hacks instead of dekker-veltkamp algorithm, which is problematic for long double where different representations are supported. (the new hypot and hypotl code should be smaller and faster on 32bit cpu archs with fast fpu), the new code behaves differently in non-nearest rounding, but the error should be still less than 2ulps. ld80 and ld128 are supported
2012-11-13math: simplify hypot and hypotf using scalbnSzabolcs Nagy-7/+2
this also fixes overflow/underflow raising and excess precision issues (as those are handled well in scalbn)
2012-03-13first commit of the new libm!Rich Felker-0/+128
thanks to the hard work of Szabolcs Nagy (nsz), identifying the best (from correctness and license standpoint) implementations from freebsd and openbsd and cleaning them up! musl should now fully support c99 float and long double math functions, and has near-complete complex math support. tgmath should also work (fully on gcc-compatible compilers, and mostly on any c99 compiler). based largely on commit 0376d44a890fea261506f1fc63833e7a686dca19 from nsz's libm git repo, with some additions (dummy versions of a few missing long double complex functions, etc.) by me. various cleanups still need to be made, including re-adding (if they're correct) some asm functions that were dropped.