| 1 |
/* -*- Mode: C; tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*- |
| 2 |
* |
| 3 |
* ***** BEGIN LICENSE BLOCK ***** |
| 4 |
* Version: MPL 1.1/GPL 2.0/LGPL 2.1 |
| 5 |
* |
| 6 |
* The contents of this file are subject to the Mozilla Public License Version |
| 7 |
* 1.1 (the "License"); you may not use this file except in compliance with |
| 8 |
* the License. You may obtain a copy of the License at |
| 9 |
* http://www.mozilla.org/MPL/ |
| 10 |
* |
| 11 |
* Software distributed under the License is distributed on an "AS IS" basis, |
| 12 |
* WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License |
| 13 |
* for the specific language governing rights and limitations under the |
| 14 |
* License. |
| 15 |
* |
| 16 |
* The Original Code is Mozilla Communicator client code, released |
| 17 |
* March 31, 1998. |
| 18 |
* |
| 19 |
* The Initial Developer of the Original Code is |
| 20 |
* Sun Microsystems, Inc. |
| 21 |
* Portions created by the Initial Developer are Copyright (C) 1998 |
| 22 |
* the Initial Developer. All Rights Reserved. |
| 23 |
* |
| 24 |
* Contributor(s): |
| 25 |
* |
| 26 |
* Alternatively, the contents of this file may be used under the terms of |
| 27 |
* either of the GNU General Public License Version 2 or later (the "GPL"), |
| 28 |
* or the GNU Lesser General Public License Version 2.1 or later (the "LGPL"), |
| 29 |
* in which case the provisions of the GPL or the LGPL are applicable instead |
| 30 |
* of those above. If you wish to allow use of your version of this file only |
| 31 |
* under the terms of either the GPL or the LGPL, and not to allow others to |
| 32 |
* use your version of this file under the terms of the MPL, indicate your |
| 33 |
* decision by deleting the provisions above and replace them with the notice |
| 34 |
* and other provisions required by the GPL or the LGPL. If you do not delete |
| 35 |
* the provisions above, a recipient may use your version of this file under |
| 36 |
* the terms of any one of the MPL, the GPL or the LGPL. |
| 37 |
* |
| 38 |
* ***** END LICENSE BLOCK ***** */ |
| 39 |
|
| 40 |
/* @(#)e_atan2.c 1.3 95/01/18 */ |
| 41 |
/* |
| 42 |
* ==================================================== |
| 43 |
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
| 44 |
* |
| 45 |
* Developed at SunSoft, a Sun Microsystems, Inc. business. |
| 46 |
* Permission to use, copy, modify, and distribute this |
| 47 |
* software is freely granted, provided that this notice |
| 48 |
* is preserved. |
| 49 |
* ==================================================== |
| 50 |
* |
| 51 |
*/ |
| 52 |
|
| 53 |
/* __ieee754_atan2(y,x) |
| 54 |
* Method : |
| 55 |
* 1. Reduce y to positive by atan2(y,x)=-atan2(-y,x). |
| 56 |
* 2. Reduce x to positive by (if x and y are unexceptional): |
| 57 |
* ARG (x+iy) = arctan(y/x) ... if x > 0, |
| 58 |
* ARG (x+iy) = pi - arctan[y/(-x)] ... if x < 0, |
| 59 |
* |
| 60 |
* Special cases: |
| 61 |
* |
| 62 |
* ATAN2((anything), NaN ) is NaN; |
| 63 |
* ATAN2(NAN , (anything) ) is NaN; |
| 64 |
* ATAN2(+-0, +(anything but NaN)) is +-0 ; |
| 65 |
* ATAN2(+-0, -(anything but NaN)) is +-pi ; |
| 66 |
* ATAN2(+-(anything but 0 and NaN), 0) is +-pi/2; |
| 67 |
* ATAN2(+-(anything but INF and NaN), +INF) is +-0 ; |
| 68 |
* ATAN2(+-(anything but INF and NaN), -INF) is +-pi; |
| 69 |
* ATAN2(+-INF,+INF ) is +-pi/4 ; |
| 70 |
* ATAN2(+-INF,-INF ) is +-3pi/4; |
| 71 |
* ATAN2(+-INF, (anything but,0,NaN, and INF)) is +-pi/2; |
| 72 |
* |
| 73 |
* Constants: |
| 74 |
* The hexadecimal values are the intended ones for the following |
| 75 |
* constants. The decimal values may be used, provided that the |
| 76 |
* compiler will convert from decimal to binary accurately enough |
| 77 |
* to produce the hexadecimal values shown. |
| 78 |
*/ |
| 79 |
|
| 80 |
#include "fdlibm.h" |
| 81 |
|
| 82 |
#ifdef __STDC__ |
| 83 |
static const double |
| 84 |
#else |
| 85 |
static double |
| 86 |
#endif |
| 87 |
tiny = 1.0e-300, |
| 88 |
zero = 0.0, |
| 89 |
pi_o_4 = 7.8539816339744827900E-01, /* 0x3FE921FB, 0x54442D18 */ |
| 90 |
pi_o_2 = 1.5707963267948965580E+00, /* 0x3FF921FB, 0x54442D18 */ |
| 91 |
pi = 3.1415926535897931160E+00, /* 0x400921FB, 0x54442D18 */ |
| 92 |
pi_lo = 1.2246467991473531772E-16; /* 0x3CA1A626, 0x33145C07 */ |
| 93 |
|
| 94 |
#ifdef __STDC__ |
| 95 |
double __ieee754_atan2(double y, double x) |
| 96 |
#else |
| 97 |
double __ieee754_atan2(y,x) |
| 98 |
double y,x; |
| 99 |
#endif |
| 100 |
{ |
| 101 |
fd_twoints ux, uy, uz; |
| 102 |
double z; |
| 103 |
int k,m,hx,hy,ix,iy; |
| 104 |
unsigned lx,ly; |
| 105 |
|
| 106 |
ux.d = x; uy.d = y; |
| 107 |
hx = __HI(ux); ix = hx&0x7fffffff; |
| 108 |
lx = __LO(ux); |
| 109 |
hy = __HI(uy); iy = hy&0x7fffffff; |
| 110 |
ly = __LO(uy); |
| 111 |
if(((ix|((lx|-(int)lx)>>31))>0x7ff00000)|| |
| 112 |
((iy|((ly|-(int)ly)>>31))>0x7ff00000)) /* x or y is NaN */ |
| 113 |
return x+y; |
| 114 |
if(((hx-0x3ff00000)|lx)==0) return fd_atan(y); /* x=1.0 */ |
| 115 |
m = ((hy>>31)&1)|((hx>>30)&2); /* 2*sign(x)+sign(y) */ |
| 116 |
|
| 117 |
/* when y = 0 */ |
| 118 |
if((iy|ly)==0) { |
| 119 |
switch(m) { |
| 120 |
case 0: |
| 121 |
case 1: return y; /* atan(+-0,+anything)=+-0 */ |
| 122 |
case 2: return pi+tiny;/* atan(+0,-anything) = pi */ |
| 123 |
case 3: return -pi-tiny;/* atan(-0,-anything) =-pi */ |
| 124 |
} |
| 125 |
} |
| 126 |
/* when x = 0 */ |
| 127 |
if((ix|lx)==0) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny; |
| 128 |
|
| 129 |
/* when x is INF */ |
| 130 |
if(ix==0x7ff00000) { |
| 131 |
if(iy==0x7ff00000) { |
| 132 |
switch(m) { |
| 133 |
case 0: return pi_o_4+tiny;/* atan(+INF,+INF) */ |
| 134 |
case 1: return -pi_o_4-tiny;/* atan(-INF,+INF) */ |
| 135 |
case 2: return 3.0*pi_o_4+tiny;/*atan(+INF,-INF)*/ |
| 136 |
case 3: return -3.0*pi_o_4-tiny;/*atan(-INF,-INF)*/ |
| 137 |
} |
| 138 |
} else { |
| 139 |
switch(m) { |
| 140 |
case 0: return zero ; /* atan(+...,+INF) */ |
| 141 |
case 1: return -zero ; /* atan(-...,+INF) */ |
| 142 |
case 2: return pi+tiny ; /* atan(+...,-INF) */ |
| 143 |
case 3: return -pi-tiny ; /* atan(-...,-INF) */ |
| 144 |
} |
| 145 |
} |
| 146 |
} |
| 147 |
/* when y is INF */ |
| 148 |
if(iy==0x7ff00000) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny; |
| 149 |
|
| 150 |
/* compute y/x */ |
| 151 |
k = (iy-ix)>>20; |
| 152 |
if(k > 60) z=pi_o_2+0.5*pi_lo; /* |y/x| > 2**60 */ |
| 153 |
else if(hx<0&&k<-60) z=0.0; /* |y|/x < -2**60 */ |
| 154 |
else z=fd_atan(fd_fabs(y/x)); /* safe to do y/x */ |
| 155 |
switch (m) { |
| 156 |
case 0: return z ; /* atan(+,+) */ |
| 157 |
case 1: uz.d = z; |
| 158 |
__HI(uz) ^= 0x80000000; |
| 159 |
z = uz.d; |
| 160 |
return z ; /* atan(-,+) */ |
| 161 |
case 2: return pi-(z-pi_lo);/* atan(+,-) */ |
| 162 |
default: /* case 3 */ |
| 163 |
return (z-pi_lo)-pi;/* atan(-,-) */ |
| 164 |
} |
| 165 |
} |
Parent Directory
| 
Revision Log