| /* |
| * IBM Accurate Mathematical Library |
| * written by International Business Machines Corp. |
| * Copyright (C) 2001-2014 Free Software Foundation, Inc. |
| * |
| * This program is free software; you can redistribute it and/or modify |
| * it under the terms of the GNU Lesser General Public License as published by |
| * the Free Software Foundation; either version 2.1 of the License, or |
| * (at your option) any later version. |
| * |
| * This program is distributed in the hope that it will be useful, |
| * but WITHOUT ANY WARRANTY; without even the implied warranty of |
| * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
| * GNU Lesser General Public License for more details. |
| * |
| * You should have received a copy of the GNU Lesser General Public License |
| * along with this program; if not, see <http://www.gnu.org/licenses/>. |
| */ |
| /************************************************************************/ |
| /* MODULE_NAME: atnat2.c */ |
| /* */ |
| /* FUNCTIONS: uatan2 */ |
| /* atan2Mp */ |
| /* signArctan2 */ |
| /* normalized */ |
| /* */ |
| /* FILES NEEDED: dla.h endian.h mpa.h mydefs.h atnat2.h */ |
| /* mpatan.c mpatan2.c mpsqrt.c */ |
| /* uatan.tbl */ |
| /* */ |
| /* An ultimate atan2() routine. Given two IEEE double machine numbers y,*/ |
| /* x it computes the correctly rounded (to nearest) value of atan2(y,x).*/ |
| /* */ |
| /* Assumption: Machine arithmetic operations are performed in */ |
| /* round to nearest mode of IEEE 754 standard. */ |
| /* */ |
| /************************************************************************/ |
| |
| #include <dla.h> |
| #include "mpa.h" |
| #include "MathLib.h" |
| #include "uatan.tbl" |
| #include "atnat2.h" |
| #include <math_private.h> |
| #include <stap-probe.h> |
| |
| #ifndef SECTION |
| # define SECTION |
| #endif |
| |
| /************************************************************************/ |
| /* An ultimate atan2 routine. Given two IEEE double machine numbers y,x */ |
| /* it computes the correctly rounded (to nearest) value of atan2(y,x). */ |
| /* Assumption: Machine arithmetic operations are performed in */ |
| /* round to nearest mode of IEEE 754 standard. */ |
| /************************************************************************/ |
| static double atan2Mp (double, double, const int[]); |
| /* Fix the sign and return after stage 1 or stage 2 */ |
| static double |
| signArctan2 (double y, double z) |
| { |
| return __copysign (z, y); |
| } |
| |
| static double normalized (double, double, double, double); |
| void __mpatan2 (mp_no *, mp_no *, mp_no *, int); |
| |
| double |
| SECTION |
| __ieee754_atan2 (double y, double x) |
| { |
| int i, de, ux, dx, uy, dy; |
| static const int pr[MM] = { 6, 8, 10, 20, 32 }; |
| double ax, ay, u, du, u9, ua, v, vv, dv, t1, t2, t3, t7, t8, |
| z, zz, cor, s1, ss1, s2, ss2; |
| #ifndef DLA_FMS |
| double t4, t5, t6; |
| #endif |
| number num; |
| |
| static const int ep = 59768832, /* 57*16**5 */ |
| em = -59768832; /* -57*16**5 */ |
| |
| /* x=NaN or y=NaN */ |
| num.d = x; |
| ux = num.i[HIGH_HALF]; |
| dx = num.i[LOW_HALF]; |
| if ((ux & 0x7ff00000) == 0x7ff00000) |
| { |
| if (((ux & 0x000fffff) | dx) != 0x00000000) |
| return x + x; |
| } |
| num.d = y; |
| uy = num.i[HIGH_HALF]; |
| dy = num.i[LOW_HALF]; |
| if ((uy & 0x7ff00000) == 0x7ff00000) |
| { |
| if (((uy & 0x000fffff) | dy) != 0x00000000) |
| return y + y; |
| } |
| |
| /* y=+-0 */ |
| if (uy == 0x00000000) |
| { |
| if (dy == 0x00000000) |
| { |
| if ((ux & 0x80000000) == 0x00000000) |
| return 0; |
| else |
| return opi.d; |
| } |
| } |
| else if (uy == 0x80000000) |
| { |
| if (dy == 0x00000000) |
| { |
| if ((ux & 0x80000000) == 0x00000000) |
| return -0.0; |
| else |
| return mopi.d; |
| } |
| } |
| |
| /* x=+-0 */ |
| if (x == 0) |
| { |
| if ((uy & 0x80000000) == 0x00000000) |
| return hpi.d; |
| else |
| return mhpi.d; |
| } |
| |
| /* x=+-INF */ |
| if (ux == 0x7ff00000) |
| { |
| if (dx == 0x00000000) |
| { |
| if (uy == 0x7ff00000) |
| { |
| if (dy == 0x00000000) |
| return qpi.d; |
| } |
| else if (uy == 0xfff00000) |
| { |
| if (dy == 0x00000000) |
| return mqpi.d; |
| } |
| else |
| { |
| if ((uy & 0x80000000) == 0x00000000) |
| return 0; |
| else |
| return -0.0; |
| } |
| } |
| } |
| else if (ux == 0xfff00000) |
| { |
| if (dx == 0x00000000) |
| { |
| if (uy == 0x7ff00000) |
| { |
| if (dy == 0x00000000) |
| return tqpi.d; |
| } |
| else if (uy == 0xfff00000) |
| { |
| if (dy == 0x00000000) |
| return mtqpi.d; |
| } |
| else |
| { |
| if ((uy & 0x80000000) == 0x00000000) |
| return opi.d; |
| else |
| return mopi.d; |
| } |
| } |
| } |
| |
| /* y=+-INF */ |
| if (uy == 0x7ff00000) |
| { |
| if (dy == 0x00000000) |
| return hpi.d; |
| } |
| else if (uy == 0xfff00000) |
| { |
| if (dy == 0x00000000) |
| return mhpi.d; |
| } |
| |
| /* either x/y or y/x is very close to zero */ |
| ax = (x < 0) ? -x : x; |
| ay = (y < 0) ? -y : y; |
| de = (uy & 0x7ff00000) - (ux & 0x7ff00000); |
| if (de >= ep) |
| { |
| return ((y > 0) ? hpi.d : mhpi.d); |
| } |
| else if (de <= em) |
| { |
| if (x > 0) |
| { |
| if ((z = ay / ax) < TWOM1022) |
| return normalized (ax, ay, y, z); |
| else |
| return signArctan2 (y, z); |
| } |
| else |
| { |
| return ((y > 0) ? opi.d : mopi.d); |
| } |
| } |
| |
| /* if either x or y is extremely close to zero, scale abs(x), abs(y). */ |
| if (ax < twom500.d || ay < twom500.d) |
| { |
| ax *= two500.d; |
| ay *= two500.d; |
| } |
| |
| /* Likewise for large x and y. */ |
| if (ax > two500.d || ay > two500.d) |
| { |
| ax *= twom500.d; |
| ay *= twom500.d; |
| } |
| |
| /* x,y which are neither special nor extreme */ |
| if (ay < ax) |
| { |
| u = ay / ax; |
| EMULV (ax, u, v, vv, t1, t2, t3, t4, t5); |
| du = ((ay - v) - vv) / ax; |
| } |
| else |
| { |
| u = ax / ay; |
| EMULV (ay, u, v, vv, t1, t2, t3, t4, t5); |
| du = ((ax - v) - vv) / ay; |
| } |
| |
| if (x > 0) |
| { |
| /* (i) x>0, abs(y)< abs(x): atan(ay/ax) */ |
| if (ay < ax) |
| { |
| if (u < inv16.d) |
| { |
| v = u * u; |
| |
| zz = du + u * v * (d3.d |
| + v * (d5.d |
| + v * (d7.d |
| + v * (d9.d |
| + v * (d11.d |
| + v * d13.d))))); |
| |
| if ((z = u + (zz - u1.d * u)) == u + (zz + u1.d * u)) |
| return signArctan2 (y, z); |
| |
| MUL2 (u, du, u, du, v, vv, t1, t2, t3, t4, t5, t6, t7, t8); |
| s1 = v * (f11.d + v * (f13.d |
| + v * (f15.d + v * (f17.d + v * f19.d)))); |
| ADD2 (f9.d, ff9.d, s1, 0, s2, ss2, t1, t2); |
| MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
| ADD2 (f7.d, ff7.d, s1, ss1, s2, ss2, t1, t2); |
| MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
| ADD2 (f5.d, ff5.d, s1, ss1, s2, ss2, t1, t2); |
| MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
| ADD2 (f3.d, ff3.d, s1, ss1, s2, ss2, t1, t2); |
| MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
| MUL2 (u, du, s1, ss1, s2, ss2, t1, t2, t3, t4, t5, t6, t7, t8); |
| ADD2 (u, du, s2, ss2, s1, ss1, t1, t2); |
| |
| if ((z = s1 + (ss1 - u5.d * s1)) == s1 + (ss1 + u5.d * s1)) |
| return signArctan2 (y, z); |
| |
| return atan2Mp (x, y, pr); |
| } |
| |
| i = (TWO52 + TWO8 * u) - TWO52; |
| i -= 16; |
| t3 = u - cij[i][0].d; |
| EADD (t3, du, v, dv); |
| t1 = cij[i][1].d; |
| t2 = cij[i][2].d; |
| zz = v * t2 + (dv * t2 |
| + v * v * (cij[i][3].d |
| + v * (cij[i][4].d |
| + v * (cij[i][5].d |
| + v * cij[i][6].d)))); |
| if (i < 112) |
| { |
| if (i < 48) |
| u9 = u91.d; /* u < 1/4 */ |
| else |
| u9 = u92.d; |
| } /* 1/4 <= u < 1/2 */ |
| else |
| { |
| if (i < 176) |
| u9 = u93.d; /* 1/2 <= u < 3/4 */ |
| else |
| u9 = u94.d; |
| } /* 3/4 <= u <= 1 */ |
| if ((z = t1 + (zz - u9 * t1)) == t1 + (zz + u9 * t1)) |
| return signArctan2 (y, z); |
| |
| t1 = u - hij[i][0].d; |
| EADD (t1, du, v, vv); |
| s1 = v * (hij[i][11].d |
| + v * (hij[i][12].d |
| + v * (hij[i][13].d |
| + v * (hij[i][14].d |
| + v * hij[i][15].d)))); |
| ADD2 (hij[i][9].d, hij[i][10].d, s1, 0, s2, ss2, t1, t2); |
| MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
| ADD2 (hij[i][7].d, hij[i][8].d, s1, ss1, s2, ss2, t1, t2); |
| MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
| ADD2 (hij[i][5].d, hij[i][6].d, s1, ss1, s2, ss2, t1, t2); |
| MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
| ADD2 (hij[i][3].d, hij[i][4].d, s1, ss1, s2, ss2, t1, t2); |
| MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
| ADD2 (hij[i][1].d, hij[i][2].d, s1, ss1, s2, ss2, t1, t2); |
| |
| if ((z = s2 + (ss2 - ub.d * s2)) == s2 + (ss2 + ub.d * s2)) |
| return signArctan2 (y, z); |
| return atan2Mp (x, y, pr); |
| } |
| |
| /* (ii) x>0, abs(x)<=abs(y): pi/2-atan(ax/ay) */ |
| if (u < inv16.d) |
| { |
| v = u * u; |
| zz = u * v * (d3.d |
| + v * (d5.d |
| + v * (d7.d |
| + v * (d9.d |
| + v * (d11.d |
| + v * d13.d))))); |
| ESUB (hpi.d, u, t2, cor); |
| t3 = ((hpi1.d + cor) - du) - zz; |
| if ((z = t2 + (t3 - u2.d)) == t2 + (t3 + u2.d)) |
| return signArctan2 (y, z); |
| |
| MUL2 (u, du, u, du, v, vv, t1, t2, t3, t4, t5, t6, t7, t8); |
| s1 = v * (f11.d |
| + v * (f13.d |
| + v * (f15.d + v * (f17.d + v * f19.d)))); |
| ADD2 (f9.d, ff9.d, s1, 0, s2, ss2, t1, t2); |
| MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
| ADD2 (f7.d, ff7.d, s1, ss1, s2, ss2, t1, t2); |
| MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
| ADD2 (f5.d, ff5.d, s1, ss1, s2, ss2, t1, t2); |
| MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
| ADD2 (f3.d, ff3.d, s1, ss1, s2, ss2, t1, t2); |
| MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
| MUL2 (u, du, s1, ss1, s2, ss2, t1, t2, t3, t4, t5, t6, t7, t8); |
| ADD2 (u, du, s2, ss2, s1, ss1, t1, t2); |
| SUB2 (hpi.d, hpi1.d, s1, ss1, s2, ss2, t1, t2); |
| |
| if ((z = s2 + (ss2 - u6.d)) == s2 + (ss2 + u6.d)) |
| return signArctan2 (y, z); |
| return atan2Mp (x, y, pr); |
| } |
| |
| i = (TWO52 + TWO8 * u) - TWO52; |
| i -= 16; |
| v = (u - cij[i][0].d) + du; |
| |
| zz = hpi1.d - v * (cij[i][2].d |
| + v * (cij[i][3].d |
| + v * (cij[i][4].d |
| + v * (cij[i][5].d |
| + v * cij[i][6].d)))); |
| t1 = hpi.d - cij[i][1].d; |
| if (i < 112) |
| ua = ua1.d; /* w < 1/2 */ |
| else |
| ua = ua2.d; /* w >= 1/2 */ |
| if ((z = t1 + (zz - ua)) == t1 + (zz + ua)) |
| return signArctan2 (y, z); |
| |
| t1 = u - hij[i][0].d; |
| EADD (t1, du, v, vv); |
| |
| s1 = v * (hij[i][11].d |
| + v * (hij[i][12].d |
| + v * (hij[i][13].d |
| + v * (hij[i][14].d |
| + v * hij[i][15].d)))); |
| |
| ADD2 (hij[i][9].d, hij[i][10].d, s1, 0, s2, ss2, t1, t2); |
| MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
| ADD2 (hij[i][7].d, hij[i][8].d, s1, ss1, s2, ss2, t1, t2); |
| MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
| ADD2 (hij[i][5].d, hij[i][6].d, s1, ss1, s2, ss2, t1, t2); |
| MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
| ADD2 (hij[i][3].d, hij[i][4].d, s1, ss1, s2, ss2, t1, t2); |
| MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
| ADD2 (hij[i][1].d, hij[i][2].d, s1, ss1, s2, ss2, t1, t2); |
| SUB2 (hpi.d, hpi1.d, s2, ss2, s1, ss1, t1, t2); |
| |
| if ((z = s1 + (ss1 - uc.d)) == s1 + (ss1 + uc.d)) |
| return signArctan2 (y, z); |
| return atan2Mp (x, y, pr); |
| } |
| |
| /* (iii) x<0, abs(x)< abs(y): pi/2+atan(ax/ay) */ |
| if (ax < ay) |
| { |
| if (u < inv16.d) |
| { |
| v = u * u; |
| zz = u * v * (d3.d |
| + v * (d5.d |
| + v * (d7.d |
| + v * (d9.d |
| + v * (d11.d + v * d13.d))))); |
| EADD (hpi.d, u, t2, cor); |
| t3 = ((hpi1.d + cor) + du) + zz; |
| if ((z = t2 + (t3 - u3.d)) == t2 + (t3 + u3.d)) |
| return signArctan2 (y, z); |
| |
| MUL2 (u, du, u, du, v, vv, t1, t2, t3, t4, t5, t6, t7, t8); |
| s1 = v * (f11.d |
| + v * (f13.d + v * (f15.d + v * (f17.d + v * f19.d)))); |
| ADD2 (f9.d, ff9.d, s1, 0, s2, ss2, t1, t2); |
| MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
| ADD2 (f7.d, ff7.d, s1, ss1, s2, ss2, t1, t2); |
| MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
| ADD2 (f5.d, ff5.d, s1, ss1, s2, ss2, t1, t2); |
| MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
| ADD2 (f3.d, ff3.d, s1, ss1, s2, ss2, t1, t2); |
| MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
| MUL2 (u, du, s1, ss1, s2, ss2, t1, t2, t3, t4, t5, t6, t7, t8); |
| ADD2 (u, du, s2, ss2, s1, ss1, t1, t2); |
| ADD2 (hpi.d, hpi1.d, s1, ss1, s2, ss2, t1, t2); |
| |
| if ((z = s2 + (ss2 - u7.d)) == s2 + (ss2 + u7.d)) |
| return signArctan2 (y, z); |
| return atan2Mp (x, y, pr); |
| } |
| |
| i = (TWO52 + TWO8 * u) - TWO52; |
| i -= 16; |
| v = (u - cij[i][0].d) + du; |
| zz = hpi1.d + v * (cij[i][2].d |
| + v * (cij[i][3].d |
| + v * (cij[i][4].d |
| + v * (cij[i][5].d |
| + v * cij[i][6].d)))); |
| t1 = hpi.d + cij[i][1].d; |
| if (i < 112) |
| ua = ua1.d; /* w < 1/2 */ |
| else |
| ua = ua2.d; /* w >= 1/2 */ |
| if ((z = t1 + (zz - ua)) == t1 + (zz + ua)) |
| return signArctan2 (y, z); |
| |
| t1 = u - hij[i][0].d; |
| EADD (t1, du, v, vv); |
| s1 = v * (hij[i][11].d |
| + v * (hij[i][12].d |
| + v * (hij[i][13].d |
| + v * (hij[i][14].d |
| + v * hij[i][15].d)))); |
| ADD2 (hij[i][9].d, hij[i][10].d, s1, 0, s2, ss2, t1, t2); |
| MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
| ADD2 (hij[i][7].d, hij[i][8].d, s1, ss1, s2, ss2, t1, t2); |
| MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
| ADD2 (hij[i][5].d, hij[i][6].d, s1, ss1, s2, ss2, t1, t2); |
| MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
| ADD2 (hij[i][3].d, hij[i][4].d, s1, ss1, s2, ss2, t1, t2); |
| MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
| ADD2 (hij[i][1].d, hij[i][2].d, s1, ss1, s2, ss2, t1, t2); |
| ADD2 (hpi.d, hpi1.d, s2, ss2, s1, ss1, t1, t2); |
| |
| if ((z = s1 + (ss1 - uc.d)) == s1 + (ss1 + uc.d)) |
| return signArctan2 (y, z); |
| return atan2Mp (x, y, pr); |
| } |
| |
| /* (iv) x<0, abs(y)<=abs(x): pi-atan(ax/ay) */ |
| if (u < inv16.d) |
| { |
| v = u * u; |
| zz = u * v * (d3.d |
| + v * (d5.d |
| + v * (d7.d |
| + v * (d9.d + v * (d11.d + v * d13.d))))); |
| ESUB (opi.d, u, t2, cor); |
| t3 = ((opi1.d + cor) - du) - zz; |
| if ((z = t2 + (t3 - u4.d)) == t2 + (t3 + u4.d)) |
| return signArctan2 (y, z); |
| |
| MUL2 (u, du, u, du, v, vv, t1, t2, t3, t4, t5, t6, t7, t8); |
| s1 = v * (f11.d + v * (f13.d + v * (f15.d + v * (f17.d + v * f19.d)))); |
| ADD2 (f9.d, ff9.d, s1, 0, s2, ss2, t1, t2); |
| MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
| ADD2 (f7.d, ff7.d, s1, ss1, s2, ss2, t1, t2); |
| MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
| ADD2 (f5.d, ff5.d, s1, ss1, s2, ss2, t1, t2); |
| MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
| ADD2 (f3.d, ff3.d, s1, ss1, s2, ss2, t1, t2); |
| MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
| MUL2 (u, du, s1, ss1, s2, ss2, t1, t2, t3, t4, t5, t6, t7, t8); |
| ADD2 (u, du, s2, ss2, s1, ss1, t1, t2); |
| SUB2 (opi.d, opi1.d, s1, ss1, s2, ss2, t1, t2); |
| |
| if ((z = s2 + (ss2 - u8.d)) == s2 + (ss2 + u8.d)) |
| return signArctan2 (y, z); |
| return atan2Mp (x, y, pr); |
| } |
| |
| i = (TWO52 + TWO8 * u) - TWO52; |
| i -= 16; |
| v = (u - cij[i][0].d) + du; |
| zz = opi1.d - v * (cij[i][2].d |
| + v * (cij[i][3].d |
| + v * (cij[i][4].d |
| + v * (cij[i][5].d + v * cij[i][6].d)))); |
| t1 = opi.d - cij[i][1].d; |
| if (i < 112) |
| ua = ua1.d; /* w < 1/2 */ |
| else |
| ua = ua2.d; /* w >= 1/2 */ |
| if ((z = t1 + (zz - ua)) == t1 + (zz + ua)) |
| return signArctan2 (y, z); |
| |
| t1 = u - hij[i][0].d; |
| |
| EADD (t1, du, v, vv); |
| |
| s1 = v * (hij[i][11].d |
| + v * (hij[i][12].d |
| + v * (hij[i][13].d |
| + v * (hij[i][14].d + v * hij[i][15].d)))); |
| |
| ADD2 (hij[i][9].d, hij[i][10].d, s1, 0, s2, ss2, t1, t2); |
| MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
| ADD2 (hij[i][7].d, hij[i][8].d, s1, ss1, s2, ss2, t1, t2); |
| MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
| ADD2 (hij[i][5].d, hij[i][6].d, s1, ss1, s2, ss2, t1, t2); |
| MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
| ADD2 (hij[i][3].d, hij[i][4].d, s1, ss1, s2, ss2, t1, t2); |
| MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
| ADD2 (hij[i][1].d, hij[i][2].d, s1, ss1, s2, ss2, t1, t2); |
| SUB2 (opi.d, opi1.d, s2, ss2, s1, ss1, t1, t2); |
| |
| if ((z = s1 + (ss1 - uc.d)) == s1 + (ss1 + uc.d)) |
| return signArctan2 (y, z); |
| return atan2Mp (x, y, pr); |
| } |
| |
| #ifndef __ieee754_atan2 |
| strong_alias (__ieee754_atan2, __atan2_finite) |
| #endif |
| |
| /* Treat the Denormalized case */ |
| static double |
| SECTION |
| normalized (double ax, double ay, double y, double z) |
| { |
| int p; |
| mp_no mpx, mpy, mpz, mperr, mpz2, mpt1; |
| p = 6; |
| __dbl_mp (ax, &mpx, p); |
| __dbl_mp (ay, &mpy, p); |
| __dvd (&mpy, &mpx, &mpz, p); |
| __dbl_mp (ue.d, &mpt1, p); |
| __mul (&mpz, &mpt1, &mperr, p); |
| __sub (&mpz, &mperr, &mpz2, p); |
| __mp_dbl (&mpz2, &z, p); |
| return signArctan2 (y, z); |
| } |
| |
| /* Stage 3: Perform a multi-Precision computation */ |
| static double |
| SECTION |
| atan2Mp (double x, double y, const int pr[]) |
| { |
| double z1, z2; |
| int i, p; |
| mp_no mpx, mpy, mpz, mpz1, mpz2, mperr, mpt1; |
| for (i = 0; i < MM; i++) |
| { |
| p = pr[i]; |
| __dbl_mp (x, &mpx, p); |
| __dbl_mp (y, &mpy, p); |
| __mpatan2 (&mpy, &mpx, &mpz, p); |
| __dbl_mp (ud[i].d, &mpt1, p); |
| __mul (&mpz, &mpt1, &mperr, p); |
| __add (&mpz, &mperr, &mpz1, p); |
| __sub (&mpz, &mperr, &mpz2, p); |
| __mp_dbl (&mpz1, &z1, p); |
| __mp_dbl (&mpz2, &z2, p); |
| if (z1 == z2) |
| { |
| LIBC_PROBE (slowatan2, 4, &p, &x, &y, &z1); |
| return z1; |
| } |
| } |
| LIBC_PROBE (slowatan2_inexact, 4, &p, &x, &y, &z1); |
| return z1; /*if impossible to do exact computing */ |
| } |
