mirror of
https://sourceware.org/git/glibc.git
synced 2024-12-26 12:41:05 +00:00
650ef4bd79
http://sourceware.org/ml/libc-alpha/2013-08/msg00084.html Another batch of ieee854 macros and union replacement. These four files also have bugs fixed with this patch. The fact that the two doubles in an IBM long double may have different signs means that negation and absolute value operations can't just twiddle one sign bit as you can with ieee864 style extended double. fmodl, remainderl, erfl and erfcl all had errors of this type. erfl also returned +1 for large magnitude negative input where it should return -1. The hypotl error is innocuous since the value adjusted twice is only used as a flag. The e_hypotl.c tests for large "a" and small "b" are mutually exclusive because we've already exited when x/y > 2**120. That allows some further small simplifications. [BZ #15734], [BZ #15735] * sysdeps/ieee754/ldbl-128ibm/e_fmodl.c (__ieee754_fmodl): Rewrite all uses of ieee875 long double macros and unions. Simplify test for 0.0L. Correct |x|<|y| and |x|=|y| test. Use ldbl_extract_mantissa value for ix,iy exponents. Properly normalize after ldbl_extract_mantissa, and don't add hidden bit already handled. Don't treat low word of ieee854 mantissa like low word of IBM long double and mask off bit when testing for zero. * sysdeps/ieee754/ldbl-128ibm/e_hypotl.c (__ieee754_hypotl): Rewrite all uses of ieee875 long double macros and unions. Simplify tests for 0.0L and inf. Correct double adjustment of k. Delete dead code adjusting ha,hb. Simplify code setting kld. Delete two600 and two1022, instead use their values. Recognise that tests for large "a" and small "b" are mutually exclusive. Rename vars. Comment. * sysdeps/ieee754/ldbl-128ibm/e_remainderl.c (__ieee754_remainderl): Rewrite all uses of ieee875 long double macros and unions. Simplify test for 0.0L and nan. Correct negation. * sysdeps/ieee754/ldbl-128ibm/s_erfl.c (__erfl): Rewrite all uses of ieee875 long double macros and unions. Correct output for large magnitude x. Correct absolute value calculation. (__erfcl): Likewise. * math/libm-test.inc: Add tests for errors discovered in IBM long double versions of fmodl, remainderl, erfl and erfcl.
133 lines
3.4 KiB
C
133 lines
3.4 KiB
C
/* @(#)e_hypotl.c 5.1 93/09/24 */
|
|
/*
|
|
* ====================================================
|
|
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
|
*
|
|
* Developed at SunPro, a Sun Microsystems, Inc. business.
|
|
* Permission to use, copy, modify, and distribute this
|
|
* software is freely granted, provided that this notice
|
|
* is preserved.
|
|
* ====================================================
|
|
*/
|
|
|
|
/* __ieee754_hypotl(x,y)
|
|
*
|
|
* Method :
|
|
* If (assume round-to-nearest) z=x*x+y*y
|
|
* has error less than sqrtl(2)/2 ulp, than
|
|
* sqrtl(z) has error less than 1 ulp (exercise).
|
|
*
|
|
* So, compute sqrtl(x*x+y*y) with some care as
|
|
* follows to get the error below 1 ulp:
|
|
*
|
|
* Assume x>y>0;
|
|
* (if possible, set rounding to round-to-nearest)
|
|
* 1. if x > 2y use
|
|
* x1*x1+(y*y+(x2*(x+x1))) for x*x+y*y
|
|
* where x1 = x with lower 53 bits cleared, x2 = x-x1; else
|
|
* 2. if x <= 2y use
|
|
* t1*y1+((x-y)*(x-y)+(t1*y2+t2*y))
|
|
* where t1 = 2x with lower 53 bits cleared, t2 = 2x-t1,
|
|
* y1= y with lower 53 bits chopped, y2 = y-y1.
|
|
*
|
|
* NOTE: scaling may be necessary if some argument is too
|
|
* large or too tiny
|
|
*
|
|
* Special cases:
|
|
* hypotl(x,y) is INF if x or y is +INF or -INF; else
|
|
* hypotl(x,y) is NAN if x or y is NAN.
|
|
*
|
|
* Accuracy:
|
|
* hypotl(x,y) returns sqrtl(x^2+y^2) with error less
|
|
* than 1 ulps (units in the last place)
|
|
*/
|
|
|
|
#include <math.h>
|
|
#include <math_private.h>
|
|
|
|
long double
|
|
__ieee754_hypotl(long double x, long double y)
|
|
{
|
|
long double a,b,a1,a2,b1,b2,w,kld;
|
|
int64_t j,k,ha,hb;
|
|
double xhi, yhi, hi, lo;
|
|
|
|
xhi = ldbl_high (x);
|
|
EXTRACT_WORDS64 (ha, xhi);
|
|
yhi = ldbl_high (y);
|
|
EXTRACT_WORDS64 (hb, yhi);
|
|
ha &= 0x7fffffffffffffffLL;
|
|
hb &= 0x7fffffffffffffffLL;
|
|
if(hb > ha) {a=y;b=x;j=ha; ha=hb;hb=j;} else {a=x;b=y;}
|
|
a = fabsl(a); /* a <- |a| */
|
|
b = fabsl(b); /* b <- |b| */
|
|
if((ha-hb)>0x0780000000000000LL) {return a+b;} /* x/y > 2**120 */
|
|
k=0;
|
|
kld = 1.0L;
|
|
if(ha > 0x5f30000000000000LL) { /* a>2**500 */
|
|
if(ha >= 0x7ff0000000000000LL) { /* Inf or NaN */
|
|
w = a+b; /* for sNaN */
|
|
if(ha == 0x7ff0000000000000LL)
|
|
w = a;
|
|
if(hb == 0x7ff0000000000000LL)
|
|
w = b;
|
|
return w;
|
|
}
|
|
/* scale a and b by 2**-600 */
|
|
a *= 0x1p-600L;
|
|
b *= 0x1p-600L;
|
|
k = 600;
|
|
kld = 0x1p+600L;
|
|
}
|
|
else if(hb < 0x23d0000000000000LL) { /* b < 2**-450 */
|
|
if(hb <= 0x000fffffffffffffLL) { /* subnormal b or 0 */
|
|
if(hb==0) return a;
|
|
a *= 0x1p+1022L;
|
|
b *= 0x1p+1022L;
|
|
k = -1022;
|
|
kld = 0x1p-1022L;
|
|
} else { /* scale a and b by 2^600 */
|
|
a *= 0x1p+600L;
|
|
b *= 0x1p+600L;
|
|
k = -600;
|
|
kld = 0x1p-600L;
|
|
}
|
|
}
|
|
/* medium size a and b */
|
|
w = a-b;
|
|
if (w>b) {
|
|
ldbl_unpack (a, &hi, &lo);
|
|
a1 = hi;
|
|
a2 = lo;
|
|
/* a*a + b*b
|
|
= (a1+a2)*a + b*b
|
|
= a1*a + a2*a + b*b
|
|
= a1*(a1+a2) + a2*a + b*b
|
|
= a1*a1 + a1*a2 + a2*a + b*b
|
|
= a1*a1 + a2*(a+a1) + b*b */
|
|
w = __ieee754_sqrtl(a1*a1-(b*(-b)-a2*(a+a1)));
|
|
} else {
|
|
a = a+a;
|
|
ldbl_unpack (b, &hi, &lo);
|
|
b1 = hi;
|
|
b2 = lo;
|
|
ldbl_unpack (a, &hi, &lo);
|
|
a1 = hi;
|
|
a2 = lo;
|
|
/* a*a + b*b
|
|
= a*a + (a-b)*(a-b) - (a-b)*(a-b) + b*b
|
|
= a*a + w*w - (a*a - 2*a*b + b*b) + b*b
|
|
= w*w + 2*a*b
|
|
= w*w + (a1+a2)*b
|
|
= w*w + a1*b + a2*b
|
|
= w*w + a1*(b1+b2) + a2*b
|
|
= w*w + a1*b1 + a1*b2 + a2*b */
|
|
w = __ieee754_sqrtl(a1*b1-(w*(-w)-(a1*b2+a2*b)));
|
|
}
|
|
if(k!=0)
|
|
return w*kld;
|
|
else
|
|
return w;
|
|
}
|
|
strong_alias (__ieee754_hypotl, __hypotl_finite)
|