summaryrefslogtreecommitdiffstats
path: root/lib/libm/common/atan2.c
diff options
context:
space:
mode:
Diffstat (limited to 'lib/libm/common/atan2.c')
-rw-r--r--lib/libm/common/atan2.c60
1 files changed, 30 insertions, 30 deletions
diff --git a/lib/libm/common/atan2.c b/lib/libm/common/atan2.c
index 958a154..b847a1d 100644
--- a/lib/libm/common/atan2.c
+++ b/lib/libm/common/atan2.c
@@ -38,21 +38,21 @@ static char sccsid[] = "@(#)atan2.c 8.1 (Berkeley) 6/4/93";
/* ATAN2(Y,X)
* RETURN ARG (X+iY)
* DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
- * CODED IN C BY K.C. NG, 1/8/85;
+ * CODED IN C BY K.C. NG, 1/8/85;
* REVISED BY K.C. NG on 2/7/85, 2/13/85, 3/7/85, 3/30/85, 6/29/85.
*
* Required system supported functions :
* copysign(x,y)
* scalb(x,y)
* logb(x)
- *
+ *
* Method :
* 1. Reduce y to positive by atan2(y,x)=-atan2(-y,x).
- * 2. Reduce x to positive by (if x and y are unexceptional):
+ * 2. Reduce x to positive by (if x and y are unexceptional):
* ARG (x+iy) = arctan(y/x) ... if x > 0,
* ARG (x+iy) = pi - arctan[y/(-x)] ... if x < 0,
- * 3. According to the integer k=4t+0.25 truncated , t=y/x, the argument
- * is further reduced to one of the following intervals and the
+ * 3. According to the integer k=4t+0.25 truncated , t=y/x, the argument
+ * is further reduced to one of the following intervals and the
* arctangent of y/x is evaluated by the corresponding formula:
*
* [0,7/16] atan(y/x) = t - t^3*(a1+t^2*(a2+...(a10+t^2*a11)...)
@@ -76,32 +76,32 @@ static char sccsid[] = "@(#)atan2.c 8.1 (Berkeley) 6/4/93";
* ARG( (anything but,0,NaN, and INF),+-INF ) is +-PI/2;
*
* Accuracy:
- * atan2(y,x) returns (PI/pi) * the exact ARG (x+iy) nearly rounded,
+ * atan2(y,x) returns (PI/pi) * the exact ARG (x+iy) nearly rounded,
* where
*
* in decimal:
- * pi = 3.141592653589793 23846264338327 .....
+ * pi = 3.141592653589793 23846264338327 .....
* 53 bits PI = 3.141592653589793 115997963 ..... ,
- * 56 bits PI = 3.141592653589793 227020265 ..... ,
+ * 56 bits PI = 3.141592653589793 227020265 ..... ,
*
* in hexadecimal:
* pi = 3.243F6A8885A308D313198A2E....
* 53 bits PI = 3.243F6A8885A30 = 2 * 1.921FB54442D18 error=.276ulps
* 56 bits PI = 3.243F6A8885A308 = 4 * .C90FDAA22168C2 error=.206ulps
- *
+ *
* In a test run with 356,000 random argument on [-1,1] * [-1,1] on a
* VAX, the maximum observed error was 1.41 ulps (units of the last place)
* compared with (PI/pi)*(the exact ARG(x+iy)).
*
* Note:
* We use machine PI (the true pi rounded) in place of the actual
- * value of pi for all the trig and inverse trig functions. In general,
- * if trig is one of sin, cos, tan, then computed trig(y) returns the
- * exact trig(y*pi/PI) nearly rounded; correspondingly, computed arctrig
- * returns the exact arctrig(y)*PI/pi nearly rounded. These guarantee the
+ * value of pi for all the trig and inverse trig functions. In general,
+ * if trig is one of sin, cos, tan, then computed trig(y) returns the
+ * exact trig(y*pi/PI) nearly rounded; correspondingly, computed arctrig
+ * returns the exact arctrig(y)*PI/pi nearly rounded. These guarantee the
* trig functions have period PI, and trig(arctrig(x)) returns x for
* all critical values x.
- *
+ *
* Constants:
* The hexadecimal values are the intended ones for the following constants.
* The decimal values may be used, provided that the compiler will convert
@@ -174,7 +174,7 @@ ic(a11, 1.6438029044759730479E-2 , -6, 1.0D52174A1BB54)
double atan2(y,x)
double y,x;
-{
+{
static const double zero=0, one=1, small=1.0E-9, big=1.0E18;
double t,z,signy,signx,hi,lo;
int k,m;
@@ -185,8 +185,8 @@ double y,x;
#endif /* !defined(vax)&&!defined(tahoe) */
/* copy down the sign of y and x */
- signy = copysign(one,y) ;
- signx = copysign(one,x) ;
+ signy = copysign(one,y) ;
+ signx = copysign(one,x) ;
/* if x is 1.0, goto begin */
if(x==1) { y=copysign(y,one); t=y; if(finite(t)) goto begin;}
@@ -196,10 +196,10 @@ double y,x;
/* when x = 0 */
if(x==zero) return(copysign(PIo2,signy));
-
+
/* when x is INF */
if(!finite(x))
- if(!finite(y))
+ if(!finite(y))
return(copysign((signx==one)?PIo4:3*PIo4,signy));
else
return(copysign((signx==one)?zero:PI,signy));
@@ -208,43 +208,43 @@ double y,x;
if(!finite(y)) return(copysign(PIo2,signy));
/* compute y/x */
- x=copysign(x,one);
- y=copysign(y,one);
- if((m=(k=logb(y))-logb(x)) > 60) t=big+big;
+ x=copysign(x,one);
+ y=copysign(y,one);
+ if((m=(k=logb(y))-logb(x)) > 60) t=big+big;
else if(m < -80 ) t=y/x;
else { t = y/x ; y = scalb(y,-k); x=scalb(x,-k); }
/* begin argument reduction */
begin:
- if (t < 2.4375) {
+ if (t < 2.4375) {
/* truncate 4(t+1/16) to integer for branching */
k = 4 * (t+0.0625);
switch (k) {
/* t is in [0,7/16] */
- case 0:
+ case 0:
case 1:
- if (t < small)
+ if (t < small)
{ big + small ; /* raise inexact flag */
return (copysign((signx>zero)?t:PI-t,signy)); }
hi = zero; lo = zero; break;
/* t is in [7/16,11/16] */
- case 2:
+ case 2:
hi = athfhi; lo = athflo;
z = x+x;
t = ( (y+y) - x ) / ( z + y ); break;
/* t is in [11/16,19/16] */
- case 3:
+ case 3:
case 4:
hi = PIo4; lo = zero;
t = ( y - x ) / ( x + y ); break;
/* t is in [19/16,39/16] */
- default:
+ default:
hi = at1fhi; lo = at1flo;
z = y-x; y=y+y+y; t = x+x;
t = ( (z+z)-x ) / ( t + y ); break;
@@ -252,7 +252,7 @@ begin:
}
/* end of if (t < 2.4375) */
- else
+ else
{
hi = PIo2; lo = zero;
@@ -260,7 +260,7 @@ begin:
if (t <= big) t = - x / y;
/* t is in [big, INF] */
- else
+ else
{ big+small; /* raise inexact flag */
t = zero; }
}
OpenPOWER on IntegriCloud