mirror of
https://github.com/opnsense/src.git
synced 2026-03-23 19:23:10 -04:00
6f1b8a0792
independent of the precision in most cases. This is mainly to simplify checking for errors. r176266 did this for e_pow[f].c using a less refined expression that often didn't work. r176276 fixes an error in the log message for r176266. The main refinement is to always expand to long double precision. See old log messages (especially these 2) and the comment on the macro for more general details. Specific details: - using nan_mix() consistently for the new and old pow*() functions was the only thing needed to make my consistency test for powl() vs pow() pass on amd64. - catrig[fl].c already had all the refinements, but open-coded. - e_atan2[fl].c, e_fmod[fl].c and s_remquo[fl] only had primitive NaN mixing. - e_hypot[fl].c already had a different refined version of r176266. Refine this further. nan_mix() is not directly usable here since we want to clear the sign bit. - e_remainder[f].c already had an earlier version of r176266. - s_ccosh[f].c,/s_csinh[f].c already had a version equivalent to r176266. Refine this further. nan_mix() is not directly usable here since the expression has to handle some non-NaN cases. - s_csqrt.[fl]: the mixing was special and mostly wrong. Partially fix the special version. - s_ctanh[f].c already had a version of r176266.
131 lines
3.3 KiB
C
131 lines
3.3 KiB
C
|
|
/* @(#)e_hypot.c 1.3 95/01/18 */
|
|
/*
|
|
* ====================================================
|
|
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
|
*
|
|
* Developed at SunSoft, a Sun Microsystems, Inc. business.
|
|
* Permission to use, copy, modify, and distribute this
|
|
* software is freely granted, provided that this notice
|
|
* is preserved.
|
|
* ====================================================
|
|
*/
|
|
|
|
#include <sys/cdefs.h>
|
|
__FBSDID("$FreeBSD$");
|
|
|
|
/* __ieee754_hypot(x,y)
|
|
*
|
|
* Method :
|
|
* If (assume round-to-nearest) z=x*x+y*y
|
|
* has error less than sqrt(2)/2 ulp, than
|
|
* sqrt(z) has error less than 1 ulp (exercise).
|
|
*
|
|
* So, compute sqrt(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 32 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 32 bits cleared, t2 = 2x-t1,
|
|
* y1= y with lower 32 bits chopped, y2 = y-y1.
|
|
*
|
|
* NOTE: scaling may be necessary if some argument is too
|
|
* large or too tiny
|
|
*
|
|
* Special cases:
|
|
* hypot(x,y) is INF if x or y is +INF or -INF; else
|
|
* hypot(x,y) is NAN if x or y is NAN.
|
|
*
|
|
* Accuracy:
|
|
* hypot(x,y) returns sqrt(x^2+y^2) with error less
|
|
* than 1 ulps (units in the last place)
|
|
*/
|
|
|
|
#include <float.h>
|
|
|
|
#include "math.h"
|
|
#include "math_private.h"
|
|
|
|
double
|
|
__ieee754_hypot(double x, double y)
|
|
{
|
|
double a,b,t1,t2,y1,y2,w;
|
|
int32_t j,k,ha,hb;
|
|
|
|
GET_HIGH_WORD(ha,x);
|
|
ha &= 0x7fffffff;
|
|
GET_HIGH_WORD(hb,y);
|
|
hb &= 0x7fffffff;
|
|
if(hb > ha) {a=y;b=x;j=ha; ha=hb;hb=j;} else {a=x;b=y;}
|
|
a = fabs(a);
|
|
b = fabs(b);
|
|
if((ha-hb)>0x3c00000) {return a+b;} /* x/y > 2**60 */
|
|
k=0;
|
|
if(ha > 0x5f300000) { /* a>2**500 */
|
|
if(ha >= 0x7ff00000) { /* Inf or NaN */
|
|
u_int32_t low;
|
|
/* Use original arg order iff result is NaN; quieten sNaNs. */
|
|
w = fabsl(x+0.0L)-fabs(y+0);
|
|
GET_LOW_WORD(low,a);
|
|
if(((ha&0xfffff)|low)==0) w = a;
|
|
GET_LOW_WORD(low,b);
|
|
if(((hb^0x7ff00000)|low)==0) w = b;
|
|
return w;
|
|
}
|
|
/* scale a and b by 2**-600 */
|
|
ha -= 0x25800000; hb -= 0x25800000; k += 600;
|
|
SET_HIGH_WORD(a,ha);
|
|
SET_HIGH_WORD(b,hb);
|
|
}
|
|
if(hb < 0x20b00000) { /* b < 2**-500 */
|
|
if(hb <= 0x000fffff) { /* subnormal b or 0 */
|
|
u_int32_t low;
|
|
GET_LOW_WORD(low,b);
|
|
if((hb|low)==0) return a;
|
|
t1=0;
|
|
SET_HIGH_WORD(t1,0x7fd00000); /* t1=2^1022 */
|
|
b *= t1;
|
|
a *= t1;
|
|
k -= 1022;
|
|
} else { /* scale a and b by 2^600 */
|
|
ha += 0x25800000; /* a *= 2^600 */
|
|
hb += 0x25800000; /* b *= 2^600 */
|
|
k -= 600;
|
|
SET_HIGH_WORD(a,ha);
|
|
SET_HIGH_WORD(b,hb);
|
|
}
|
|
}
|
|
/* medium size a and b */
|
|
w = a-b;
|
|
if (w>b) {
|
|
t1 = 0;
|
|
SET_HIGH_WORD(t1,ha);
|
|
t2 = a-t1;
|
|
w = sqrt(t1*t1-(b*(-b)-t2*(a+t1)));
|
|
} else {
|
|
a = a+a;
|
|
y1 = 0;
|
|
SET_HIGH_WORD(y1,hb);
|
|
y2 = b - y1;
|
|
t1 = 0;
|
|
SET_HIGH_WORD(t1,ha+0x00100000);
|
|
t2 = a - t1;
|
|
w = sqrt(t1*y1-(w*(-w)-(t1*y2+t2*b)));
|
|
}
|
|
if(k!=0) {
|
|
u_int32_t high;
|
|
t1 = 1.0;
|
|
GET_HIGH_WORD(high,t1);
|
|
SET_HIGH_WORD(t1,high+(k<<20));
|
|
return t1*w;
|
|
} else return w;
|
|
}
|
|
|
|
#if LDBL_MANT_DIG == 53
|
|
__weak_reference(hypot, hypotl);
|
|
#endif
|