← Back to branch summary

~ubuntu-branches/ubuntu/natty/eglibc/natty-security

« back to all changes in this revision

Viewing changes to sysdeps/ieee754/dbl-64/e_atan2.c

  • Committer: Bazaar Package Importer
  • Author(s): Aurelien Jarno
  • Date: 2009-05-05 09:54:14 UTC
  • Revision ID: james.westby@ubuntu.com-20090505095414-c45qsg9ixjheohru
Import upstream version 2.9

Show diffs side-by-side

added added

removed removed

Lines of Context:
sysdeps/ieee754/dbl-64/e_atan2.c
 
1
/*
 
2
 * IBM Accurate Mathematical Library
 
3
 * written by International Business Machines Corp.
 
4
 * Copyright (C) 2001 Free Software Foundation
 
5
 *
 
6
 * This program is free software; you can redistribute it and/or modify
 
7
 * it under the terms of the GNU Lesser General Public License as published by
 
8
 * the Free Software Foundation; either version 2.1 of the License, or
 
9
 * (at your option) any later version.
 
10
 *
 
11
 * This program is distributed in the hope that it will be useful,
 
12
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 
13
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
14
 * GNU Lesser General Public License for more details.
 
15
 *
 
16
 * You should have received a copy of the GNU Lesser General Public License
 
17
 * along with this program; if not, write to the Free Software
 
18
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
 
19
 */
 
20
/************************************************************************/
 
21
/*  MODULE_NAME: atnat2.c                                               */
 
22
/*                                                                      */
 
23
/*  FUNCTIONS: uatan2                                                   */
 
24
/*             atan2Mp                                                  */
 
25
/*             signArctan2                                              */
 
26
/*             normalized                                               */
 
27
/*                                                                      */
 
28
/*  FILES NEEDED: dla.h endian.h mpa.h mydefs.h atnat2.h                */
 
29
/*                mpatan.c mpatan2.c mpsqrt.c                           */
 
30
/*                uatan.tbl                                             */
 
31
/*                                                                      */
 
32
/* An ultimate atan2() routine. Given two IEEE double machine numbers y,*/
 
33
/* x it computes the correctly rounded (to nearest) value of atan2(y,x).*/
 
34
/*                                                                      */
 
35
/* Assumption: Machine arithmetic operations are performed in           */
 
36
/* round to nearest mode of IEEE 754 standard.                          */
 
37
/*                                                                      */
 
38
/************************************************************************/
 
39
 
 
40
#include "dla.h"
 
41
#include "mpa.h"
 
42
#include "MathLib.h"
 
43
#include "uatan.tbl"
 
44
#include "atnat2.h"
 
45
#include "math_private.h"
 
46
 
 
47
/************************************************************************/
 
48
/* An ultimate atan2 routine. Given two IEEE double machine numbers y,x */
 
49
/* it computes the correctly rounded (to nearest) value of atan2(y,x).  */
 
50
/* Assumption: Machine arithmetic operations are performed in           */
 
51
/* round to nearest mode of IEEE 754 standard.                          */
 
52
/************************************************************************/
 
53
static double atan2Mp(double ,double ,const int[]);
 
54
static double signArctan2(double ,double);
 
55
static double normalized(double ,double,double ,double);
 
56
void __mpatan2(mp_no *,mp_no *,mp_no *,int);
 
57
 
 
58
double __ieee754_atan2(double y,double x) {
 
59
 
 
60
  int i,de,ux,dx,uy,dy;
 
61
#if 0
 
62
  int p;
 
63
#endif
 
64
  static const int pr[MM]={6,8,10,20,32};
 
65
  double ax,ay,u,du,u9,ua,v,vv,dv,t1,t2,t3,t4,t5,t6,t7,t8,
 
66
         z,zz,cor,s1,ss1,s2,ss2;
 
67
#if 0
 
68
  double z1,z2;
 
69
#endif
 
70
  number num;
 
71
#if 0
 
72
  mp_no mperr,mpt1,mpx,mpy,mpz,mpz1,mpz2;
 
73
#endif
 
74
 
 
75
  static const int ep= 59768832,   /*  57*16**5   */
 
76
                   em=-59768832;   /* -57*16**5   */
 
77
 
 
78
  /* x=NaN or y=NaN */
 
79
  num.d = x;  ux = num.i[HIGH_HALF];  dx = num.i[LOW_HALF];
 
80
  if   ((ux&0x7ff00000)    ==0x7ff00000) {
 
81
    if (((ux&0x000fffff)|dx)!=0x00000000) return x+x; }
 
82
  num.d = y;  uy = num.i[HIGH_HALF];  dy = num.i[LOW_HALF];
 
83
  if   ((uy&0x7ff00000)    ==0x7ff00000) {
 
84
    if (((uy&0x000fffff)|dy)!=0x00000000) return y+y; }
 
85
 
 
86
  /* y=+-0 */
 
87
  if      (uy==0x00000000) {
 
88
    if    (dy==0x00000000) {
 
89
      if  ((ux&0x80000000)==0x00000000)  return ZERO;
 
90
      else                               return opi.d; } }
 
91
  else if (uy==0x80000000) {
 
92
    if    (dy==0x00000000) {
 
93
      if  ((ux&0x80000000)==0x00000000)  return MZERO;
 
94
      else                               return mopi.d;} }
 
95
 
 
96
  /* x=+-0 */
 
97
  if (x==ZERO) {
 
98
    if ((uy&0x80000000)==0x00000000)     return hpi.d;
 
99
    else                                 return mhpi.d; }
 
100
 
 
101
  /* x=+-INF */
 
102
  if          (ux==0x7ff00000) {
 
103
    if        (dx==0x00000000) {
 
104
      if      (uy==0x7ff00000) {
 
105
        if    (dy==0x00000000)  return qpi.d; }
 
106
      else if (uy==0xfff00000) {
 
107
        if    (dy==0x00000000)  return mqpi.d; }
 
108
      else {
 
109
        if    ((uy&0x80000000)==0x00000000)  return ZERO;
 
110
        else                                 return MZERO; }
 
111
    }
 
112
  }
 
113
  else if     (ux==0xfff00000) {
 
114
    if        (dx==0x00000000) {
 
115
      if      (uy==0x7ff00000) {
 
116
        if    (dy==0x00000000)  return tqpi.d; }
 
117
      else if (uy==0xfff00000) {
 
118
        if    (dy==0x00000000)  return mtqpi.d; }
 
119
      else                     {
 
120
        if    ((uy&0x80000000)==0x00000000)  return opi.d;
 
121
        else                                 return mopi.d; }
 
122
    }
 
123
  }
 
124
 
 
125
  /* y=+-INF */
 
126
  if      (uy==0x7ff00000) {
 
127
    if    (dy==0x00000000)  return hpi.d; }
 
128
  else if (uy==0xfff00000) {
 
129
    if    (dy==0x00000000)  return mhpi.d; }
 
130
 
 
131
  /* either x/y or y/x is very close to zero */
 
132
  ax = (x<ZERO) ? -x : x;    ay = (y<ZERO) ? -y : y;
 
133
  de = (uy & 0x7ff00000) - (ux & 0x7ff00000);
 
134
  if      (de>=ep)  { return ((y>ZERO) ? hpi.d : mhpi.d); }
 
135
  else if (de<=em)  {
 
136
    if    (x>ZERO)  {
 
137
      if  ((z=ay/ax)<TWOM1022)  return normalized(ax,ay,y,z);
 
138
      else                      return signArctan2(y,z); }
 
139
    else            { return ((y>ZERO) ? opi.d : mopi.d); } }
 
140
 
 
141
  /* if either x or y is extremely close to zero, scale abs(x), abs(y). */
 
142
  if (ax<twom500.d || ay<twom500.d) { ax*=two500.d;  ay*=two500.d; }
 
143
 
 
144
  /* x,y which are neither special nor extreme */
 
145
  if (ay<ax) {
 
146
    u=ay/ax;
 
147
    EMULV(ax,u,v,vv,t1,t2,t3,t4,t5)
 
148
    du=((ay-v)-vv)/ax; }
 
149
  else {
 
150
    u=ax/ay;
 
151
    EMULV(ay,u,v,vv,t1,t2,t3,t4,t5)
 
152
    du=((ax-v)-vv)/ay; }
 
153
 
 
154
  if (x>ZERO) {
 
155
 
 
156
    /* (i)   x>0, abs(y)< abs(x):  atan(ay/ax) */
 
157
    if (ay<ax) {
 
158
      if (u<inv16.d) {
 
159
        v=u*u;  zz=du+u*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
 
160
        if ((z=u+(zz-u1.d*u)) == u+(zz+u1.d*u))  return signArctan2(y,z);
 
161
 
 
162
        MUL2(u,du,u,du,v,vv,t1,t2,t3,t4,t5,t6,t7,t8)
 
163
        s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d))));
 
164
        ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2)
 
165
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 
166
        ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2)
 
167
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 
168
        ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2)
 
169
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 
170
        ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2)
 
171
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 
172
        MUL2(u,du,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
 
173
        ADD2(u,du,s2,ss2,s1,ss1,t1,t2)
 
174
        if ((z=s1+(ss1-u5.d*s1)) == s1+(ss1+u5.d*s1))  return signArctan2(y,z);
 
175
        return atan2Mp(x,y,pr);
 
176
      }
 
177
      else {
 
178
        i=(TWO52+TWO8*u)-TWO52;  i-=16;
 
179
        t3=u-cij[i][0].d;
 
180
        EADD(t3,du,v,dv)
 
181
        t1=cij[i][1].d;  t2=cij[i][2].d;
 
182
        zz=v*t2+(dv*t2+v*v*(cij[i][3].d+v*(cij[i][4].d+
 
183
                         v*(cij[i][5].d+v* cij[i][6].d))));
 
184
        if (i<112) {
 
185
          if (i<48)  u9=u91.d;    /* u < 1/4        */
 
186
          else       u9=u92.d; }  /* 1/4 <= u < 1/2 */
 
187
        else {
 
188
          if (i<176) u9=u93.d;    /* 1/2 <= u < 3/4 */
 
189
          else       u9=u94.d; }  /* 3/4 <= u <= 1  */
 
190
        if ((z=t1+(zz-u9*t1)) == t1+(zz+u9*t1))  return signArctan2(y,z);
 
191
 
 
192
        t1=u-hij[i][0].d;
 
193
        EADD(t1,du,v,vv)
 
194
        s1=v*(hij[i][11].d+v*(hij[i][12].d+v*(hij[i][13].d+
 
195
           v*(hij[i][14].d+v* hij[i][15].d))));
 
196
        ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
 
197
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 
198
        ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2)
 
199
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 
200
        ADD2(hij[i][5].d,hij[i][6].d,s1,ss1,s2,ss2,t1,t2)
 
201
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 
202
        ADD2(hij[i][3].d,hij[i][4].d,s1,ss1,s2,ss2,t1,t2)
 
203
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 
204
        ADD2(hij[i][1].d,hij[i][2].d,s1,ss1,s2,ss2,t1,t2)
 
205
        if ((z=s2+(ss2-ub.d*s2)) == s2+(ss2+ub.d*s2))  return signArctan2(y,z);
 
206
        return atan2Mp(x,y,pr);
 
207
      }
 
208
    }
 
209
 
 
210
    /* (ii)  x>0, abs(x)<=abs(y):  pi/2-atan(ax/ay) */
 
211
    else {
 
212
      if (u<inv16.d) {
 
213
        v=u*u;
 
214
        zz=u*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
 
215
        ESUB(hpi.d,u,t2,cor)
 
216
        t3=((hpi1.d+cor)-du)-zz;
 
217
        if ((z=t2+(t3-u2.d)) == t2+(t3+u2.d))  return signArctan2(y,z);
 
218
 
 
219
        MUL2(u,du,u,du,v,vv,t1,t2,t3,t4,t5,t6,t7,t8)
 
220
        s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d))));
 
221
        ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2)
 
222
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 
223
        ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2)
 
224
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 
225
        ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2)
 
226
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 
227
        ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2)
 
228
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 
229
        MUL2(u,du,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
 
230
        ADD2(u,du,s2,ss2,s1,ss1,t1,t2)
 
231
        SUB2(hpi.d,hpi1.d,s1,ss1,s2,ss2,t1,t2)
 
232
        if ((z=s2+(ss2-u6.d)) == s2+(ss2+u6.d))  return signArctan2(y,z);
 
233
        return atan2Mp(x,y,pr);
 
234
      }
 
235
      else {
 
236
        i=(TWO52+TWO8*u)-TWO52;  i-=16;
 
237
        v=(u-cij[i][0].d)+du;
 
238
        zz=hpi1.d-v*(cij[i][2].d+v*(cij[i][3].d+v*(cij[i][4].d+
 
239
                                 v*(cij[i][5].d+v* cij[i][6].d))));
 
240
        t1=hpi.d-cij[i][1].d;
 
241
        if (i<112)  ua=ua1.d;  /* w <  1/2 */
 
242
        else        ua=ua2.d;  /* w >= 1/2 */
 
243
        if ((z=t1+(zz-ua)) == t1+(zz+ua))  return signArctan2(y,z);
 
244
 
 
245
        t1=u-hij[i][0].d;
 
246
        EADD(t1,du,v,vv)
 
247
        s1=v*(hij[i][11].d+v*(hij[i][12].d+v*(hij[i][13].d+
 
248
           v*(hij[i][14].d+v* hij[i][15].d))));
 
249
        ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
 
250
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 
251
        ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2)
 
252
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 
253
        ADD2(hij[i][5].d,hij[i][6].d,s1,ss1,s2,ss2,t1,t2)
 
254
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 
255
        ADD2(hij[i][3].d,hij[i][4].d,s1,ss1,s2,ss2,t1,t2)
 
256
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 
257
        ADD2(hij[i][1].d,hij[i][2].d,s1,ss1,s2,ss2,t1,t2)
 
258
        SUB2(hpi.d,hpi1.d,s2,ss2,s1,ss1,t1,t2)
 
259
        if ((z=s1+(ss1-uc.d)) == s1+(ss1+uc.d))  return signArctan2(y,z);
 
260
        return atan2Mp(x,y,pr);
 
261
      }
 
262
    }
 
263
  }
 
264
  else {
 
265
 
 
266
    /* (iii) x<0, abs(x)< abs(y):  pi/2+atan(ax/ay) */
 
267
    if (ax<ay) {
 
268
      if (u<inv16.d) {
 
269
        v=u*u;
 
270
        zz=u*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
 
271
        EADD(hpi.d,u,t2,cor)
 
272
        t3=((hpi1.d+cor)+du)+zz;
 
273
        if ((z=t2+(t3-u3.d)) == t2+(t3+u3.d))  return signArctan2(y,z);
 
274
 
 
275
        MUL2(u,du,u,du,v,vv,t1,t2,t3,t4,t5,t6,t7,t8)
 
276
        s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d))));
 
277
        ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2)
 
278
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 
279
        ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2)
 
280
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 
281
        ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2)
 
282
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 
283
        ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2)
 
284
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 
285
        MUL2(u,du,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
 
286
        ADD2(u,du,s2,ss2,s1,ss1,t1,t2)
 
287
        ADD2(hpi.d,hpi1.d,s1,ss1,s2,ss2,t1,t2)
 
288
        if ((z=s2+(ss2-u7.d)) == s2+(ss2+u7.d))  return signArctan2(y,z);
 
289
        return atan2Mp(x,y,pr);
 
290
      }
 
291
      else {
 
292
        i=(TWO52+TWO8*u)-TWO52;  i-=16;
 
293
        v=(u-cij[i][0].d)+du;
 
294
        zz=hpi1.d+v*(cij[i][2].d+v*(cij[i][3].d+v*(cij[i][4].d+
 
295
                                 v*(cij[i][5].d+v* cij[i][6].d))));
 
296
        t1=hpi.d+cij[i][1].d;
 
297
        if (i<112)  ua=ua1.d;  /* w <  1/2 */
 
298
        else        ua=ua2.d;  /* w >= 1/2 */
 
299
        if ((z=t1+(zz-ua)) == t1+(zz+ua))  return signArctan2(y,z);
 
300
 
 
301
        t1=u-hij[i][0].d;
 
302
        EADD(t1,du,v,vv)
 
303
        s1=v*(hij[i][11].d+v*(hij[i][12].d+v*(hij[i][13].d+
 
304
           v*(hij[i][14].d+v* hij[i][15].d))));
 
305
        ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
 
306
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 
307
        ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2)
 
308
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 
309
        ADD2(hij[i][5].d,hij[i][6].d,s1,ss1,s2,ss2,t1,t2)
 
310
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 
311
        ADD2(hij[i][3].d,hij[i][4].d,s1,ss1,s2,ss2,t1,t2)
 
312
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 
313
        ADD2(hij[i][1].d,hij[i][2].d,s1,ss1,s2,ss2,t1,t2)
 
314
        ADD2(hpi.d,hpi1.d,s2,ss2,s1,ss1,t1,t2)
 
315
        if ((z=s1+(ss1-uc.d)) == s1+(ss1+uc.d))  return signArctan2(y,z);
 
316
        return atan2Mp(x,y,pr);
 
317
      }
 
318
    }
 
319
 
 
320
    /* (iv)  x<0, abs(y)<=abs(x):  pi-atan(ax/ay) */
 
321
    else {
 
322
      if (u<inv16.d) {
 
323
        v=u*u;
 
324
        zz=u*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
 
325
        ESUB(opi.d,u,t2,cor)
 
326
        t3=((opi1.d+cor)-du)-zz;
 
327
        if ((z=t2+(t3-u4.d)) == t2+(t3+u4.d))  return signArctan2(y,z);
 
328
 
 
329
        MUL2(u,du,u,du,v,vv,t1,t2,t3,t4,t5,t6,t7,t8)
 
330
        s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d))));
 
331
        ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2)
 
332
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 
333
        ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2)
 
334
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 
335
        ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2)
 
336
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 
337
        ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2)
 
338
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 
339
        MUL2(u,du,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
 
340
        ADD2(u,du,s2,ss2,s1,ss1,t1,t2)
 
341
        SUB2(opi.d,opi1.d,s1,ss1,s2,ss2,t1,t2)
 
342
        if ((z=s2+(ss2-u8.d)) == s2+(ss2+u8.d))  return signArctan2(y,z);
 
343
        return atan2Mp(x,y,pr);
 
344
      }
 
345
      else {
 
346
        i=(TWO52+TWO8*u)-TWO52;  i-=16;
 
347
        v=(u-cij[i][0].d)+du;
 
348
        zz=opi1.d-v*(cij[i][2].d+v*(cij[i][3].d+v*(cij[i][4].d+
 
349
                                 v*(cij[i][5].d+v* cij[i][6].d))));
 
350
        t1=opi.d-cij[i][1].d;
 
351
        if (i<112)  ua=ua1.d;  /* w <  1/2 */
 
352
        else        ua=ua2.d;  /* w >= 1/2 */
 
353
        if ((z=t1+(zz-ua)) == t1+(zz+ua))  return signArctan2(y,z);
 
354
 
 
355
        t1=u-hij[i][0].d;
 
356
        EADD(t1,du,v,vv)
 
357
        s1=v*(hij[i][11].d+v*(hij[i][12].d+v*(hij[i][13].d+
 
358
           v*(hij[i][14].d+v* hij[i][15].d))));
 
359
        ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
 
360
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 
361
        ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2)
 
362
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 
363
        ADD2(hij[i][5].d,hij[i][6].d,s1,ss1,s2,ss2,t1,t2)
 
364
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 
365
        ADD2(hij[i][3].d,hij[i][4].d,s1,ss1,s2,ss2,t1,t2)
 
366
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 
367
        ADD2(hij[i][1].d,hij[i][2].d,s1,ss1,s2,ss2,t1,t2)
 
368
        SUB2(opi.d,opi1.d,s2,ss2,s1,ss1,t1,t2)
 
369
        if ((z=s1+(ss1-uc.d)) == s1+(ss1+uc.d))  return signArctan2(y,z);
 
370
        return atan2Mp(x,y,pr);
 
371
      }
 
372
    }
 
373
  }
 
374
}
 
375
  /* Treat the Denormalized case */
 
376
static double  normalized(double ax,double ay,double y, double z)
 
377
    { int p;
 
378
      mp_no mpx,mpy,mpz,mperr,mpz2,mpt1;
 
379
  p=6;
 
380
  __dbl_mp(ax,&mpx,p);  __dbl_mp(ay,&mpy,p);  __dvd(&mpy,&mpx,&mpz,p);
 
381
  __dbl_mp(ue.d,&mpt1,p);   __mul(&mpz,&mpt1,&mperr,p);
 
382
  __sub(&mpz,&mperr,&mpz2,p);  __mp_dbl(&mpz2,&z,p);
 
383
  return signArctan2(y,z);
 
384
}
 
385
  /* Fix the sign and return after stage 1 or stage 2 */
 
386
static double signArctan2(double y,double z)
 
387
{
 
388
  return ((y<ZERO) ? -z : z);
 
389
}
 
390
  /* Stage 3: Perform a multi-Precision computation */
 
391
static double  atan2Mp(double x,double y,const int pr[])
 
392
{
 
393
  double z1,z2;
 
394
  int i,p;
 
395
  mp_no mpx,mpy,mpz,mpz1,mpz2,mperr,mpt1;
 
396
  for (i=0; i<MM; i++) {
 
397
    p = pr[i];
 
398
    __dbl_mp(x,&mpx,p);  __dbl_mp(y,&mpy,p);
 
399
    __mpatan2(&mpy,&mpx,&mpz,p);
 
400
    __dbl_mp(ud[i].d,&mpt1,p);   __mul(&mpz,&mpt1,&mperr,p);
 
401
    __add(&mpz,&mperr,&mpz1,p);  __sub(&mpz,&mperr,&mpz2,p);
 
402
    __mp_dbl(&mpz1,&z1,p);       __mp_dbl(&mpz2,&z2,p);
 
403
    if (z1==z2)   return z1;
 
404
  }
 
405
  return z1; /*if unpossible to do exact computing */
 
406
}