diff options
author | Danny Smith <dannysmith@users.sourceforge.net> | 2002-09-02 07:00:37 +0400 |
---|---|---|
committer | Danny Smith <dannysmith@users.sourceforge.net> | 2002-09-02 07:00:37 +0400 |
commit | 169618f29fbf81921ed941f5512810822b575339 (patch) | |
tree | af9efb1c40e75d28bc0a4cc29d51672d56cacfc1 /winsup/mingw/mingwex | |
parent | c8bef4002688ea8a5d919c312f25c899f5552625 (diff) |
* mingwex/math/hypotl.c: Replace with version based on cephes
library.
Diffstat (limited to 'winsup/mingw/mingwex')
-rw-r--r-- | winsup/mingw/mingwex/math/hypotl.c | 101 |
1 files changed, 58 insertions, 43 deletions
diff --git a/winsup/mingw/mingwex/math/hypotl.c b/winsup/mingw/mingwex/math/hypotl.c index 2bcfc8638..2a25b82c3 100644 --- a/winsup/mingw/mingwex/math/hypotl.c +++ b/winsup/mingw/mingwex/math/hypotl.c @@ -1,58 +1,73 @@ #include <math.h> +#include <float.h> +#include <errno.h> -typedef union -{ - long double value; - struct - { - unsigned mantissa[2]; - unsigned sign_exponent : 16; - unsigned empty : 16; - } parts; -} ieee_long_double_shape_type; +/* +This implementation is based largely on Cephes library +function cabsl (cmplxl.c), which bears the following notice: +Cephes Math Library Release 2.1: January, 1989 +Copyright 1984, 1987, 1989 by Stephen L. Moshier +Direct inquiries to 30 Frost Street, Cambridge, MA 02140 +*/ +/* + Modified for use in libmingwex.a + 02 Sept 2002 Danny Smith <dannysmith@users.sourceforege.net> + Calls to ldexpl replaced by logbl and calls to frexpl replaced + by scalbnl to avoid duplicated range checks. +*/ -/* Get int from the exponent of a long double. */ -static __inline__ void -get_ld_exp (unsigned exp,long double d) -{ - ieee_long_double_shape_type u; - u.value = d; - exp = u.parts.sign_exponent; -} - -/* Set exponent of a long double from an int. */ -static __inline__ void -set_ld_exp (long double d,unsigned exp) -{ - ieee_long_double_shape_type u; - u.value = d; - u.parts.sign_exponent = exp; - d = u.value; -} +extern long double __INFL; +#define PRECL 32 long double hypotl (long double x, long double y) { - unsigned exx = 0U; - unsigned eyy = 0U; - unsigned scale; + int exx; + int eyy; + int scale; long double xx =fabsl(x); long double yy =fabsl(y); if (!isfinite(xx) || !isfinite(yy)) return xx + yy; /* Return INF or NAN. */ - /* Scale to avoid overflow.*/ - get_ld_exp (exx, xx); - get_ld_exp (eyy, yy); - scale = (exx > eyy ? exx : eyy); - if (scale == 0) - return 0.0L; - set_ld_exp (xx, exx - scale); - set_ld_exp (yy, eyy - scale); - xx = sqrtl(xx * xx + yy * yy); - get_ld_exp (exx,xx); - set_ld_exp (xx, exx + scale); - return xx; + if (xx == 0.0L) + return yy; + if (yy == 0.0L) + return xx; + + /* Get exponents */ + exx = logbl (xx); + eyy = logbl (yy); + + /* Check if large differences in scale */ + scale = exx - eyy; + if ( scale > PRECL) + return xx; + if ( scale < -PRECL) + return yy; + + /* Exponent of approximate geometric mean (x 2) */ + scale = (exx + eyy) >> 1; + + /* Rescale: Geometric mean is now about 2 */ + x = scalbnl(xx, -scale); + y = scalbnl(yy, -scale); + + xx = sqrtl(x * x + y * y); + + /* Check for overflow and underflow */ + exx = logbl(xx); + exx += scale; + if (exx > LDBL_MAX_EXP) + { + errno = ERANGE; + return __INFL; + } + if (exx < LDBL_MIN_EXP) + return 0.0L; + + /* Undo scaling */ + return (scalbnl (xx, scale)); } |