← Back to branch summary

~ubuntu-branches/ubuntu/precise/eglibc/precise-201308281639

« back to all changes in this revision

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

  • Committer: Package Import Robot
  • Author(s): Matthias Klose
  • Date: 2012-02-08 01:58:09 UTC
  • mfrom: (1.5.3) (288.1.12 precise)
  • Revision ID: package-import@ubuntu.com-20120208015809-ulscst7uteq3e22z
Tags: 2.15~pre6-0ubuntu10
Merge from Debian (r5151, 2.13-26).

Show diffs side-by-side

added added

removed removed

Lines of Context:
sysdeps/ieee754/dbl-64/s_atan.c
1
1
/*
2
2
 * IBM Accurate Mathematical Library
3
3
 * written by International Business Machines Corp.
4
 
 * Copyright (C) 2001 Free Software Foundation
 
4
 * Copyright (C) 2001, 2011 Free Software Foundation
5
5
 *
6
6
 * This program is free software; you can redistribute it and/or modify
7
7
 * it under the terms of the GNU Lesser General Public License as published by
37
37
/*                                                                      */
38
38
/************************************************************************/
39
39
 
40
 
#include "dla.h"
 
40
#include <dla.h>
41
41
#include "mpa.h"
42
42
#include "MathLib.h"
43
43
#include "uatan.tbl"
46
46
 
47
47
void __mpatan(mp_no *,mp_no *,int);          /* see definition in mpatan.c */
48
48
static double atanMp(double,const int[]);
49
 
double __signArctan(double,double);
 
49
 
 
50
  /* Fix the sign of y and return */
 
51
static double  __signArctan(double x,double y){
 
52
  return __copysign(y, x);
 
53
}
 
54
 
 
55
 
50
56
/* An ultimate atan() routine. Given an IEEE double machine number x,    */
51
57
/* routine computes the correctly rounded (to nearest) value of atan(x). */
52
58
double atan(double x) {
53
59
 
54
60
 
55
 
  double cor,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,u,u2,u3,
56
 
         v,vv,w,ww,y,yy,z,zz;
 
61
  double cor,s1,ss1,s2,ss2,t1,t2,t3,t7,t8,t9,t10,u,u2,u3,
 
62
         v,vv,w,ww,y,yy,z,zz;
 
63
#ifndef DLA_FMS
 
64
  double t4,t5,t6;
 
65
#endif
57
66
#if 0
58
67
  double y1,y2;
59
68
#endif
78
87
  if (u<C) {
79
88
    if (u<B) {
80
89
      if (u<A) {                                           /* u < A */
81
 
         return x; }
 
90
         return x; }
82
91
      else {                                               /* A <= u < B */
83
 
        v=x*x;  yy=x*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
84
 
        if ((y=x+(yy-U1*x)) == x+(yy+U1*x))  return y;
85
 
 
86
 
        EMULV(x,x,v,vv,t1,t2,t3,t4,t5)                       /* v+vv=x^2 */
87
 
        s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d))));
88
 
        ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2)
89
 
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
90
 
        ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2)
91
 
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
92
 
        ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2)
93
 
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
94
 
        ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2)
95
 
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
96
 
        MUL2(x,ZERO,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
97
 
        ADD2(x,ZERO,s2,ss2,s1,ss1,t1,t2)
98
 
        if ((y=s1+(ss1-U5*s1)) == s1+(ss1+U5*s1))  return y;
99
 
 
100
 
        return atanMp(x,pr);
 
92
        v=x*x;  yy=x*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
 
93
        if ((y=x+(yy-U1*x)) == x+(yy+U1*x))  return y;
 
94
 
 
95
        EMULV(x,x,v,vv,t1,t2,t3,t4,t5)                       /* v+vv=x^2 */
 
96
        s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d))));
 
97
        ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2)
 
98
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 
99
        ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2)
 
100
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 
101
        ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2)
 
102
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 
103
        ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2)
 
104
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 
105
        MUL2(x,ZERO,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
 
106
        ADD2(x,ZERO,s2,ss2,s1,ss1,t1,t2)
 
107
        if ((y=s1+(ss1-U5*s1)) == s1+(ss1+U5*s1))  return y;
 
108
 
 
109
        return atanMp(x,pr);
101
110
      } }
102
111
    else {  /* B <= u < C */
103
112
      i=(TWO52+TWO8*u)-TWO52;  i-=16;
104
113
      z=u-cij[i][0].d;
105
114
      yy=z*(cij[i][2].d+z*(cij[i][3].d+z*(cij[i][4].d+
106
 
                        z*(cij[i][5].d+z* cij[i][6].d))));
 
115
                        z*(cij[i][5].d+z* cij[i][6].d))));
107
116
      t1=cij[i][1].d;
108
117
      if (i<112) {
109
 
        if (i<48)  u2=U21;    /* u < 1/4        */
110
 
        else       u2=U22; }  /* 1/4 <= u < 1/2 */
 
118
        if (i<48)  u2=U21;    /* u < 1/4        */
 
119
        else       u2=U22; }  /* 1/4 <= u < 1/2 */
111
120
      else {
112
 
        if (i<176) u2=U23;    /* 1/2 <= u < 3/4 */
113
 
        else       u2=U24; }  /* 3/4 <= u <= 1  */
 
121
        if (i<176) u2=U23;    /* 1/2 <= u < 3/4 */
 
122
        else       u2=U24; }  /* 3/4 <= u <= 1  */
114
123
      if ((y=t1+(yy-u2*t1)) == t1+(yy+u2*t1))  return __signArctan(x,y);
115
124
 
116
125
      z=u-hij[i][0].d;
117
126
      s1=z*(hij[i][11].d+z*(hij[i][12].d+z*(hij[i][13].d+
118
 
         z*(hij[i][14].d+z* hij[i][15].d))));
 
127
         z*(hij[i][14].d+z* hij[i][15].d))));
119
128
      ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
120
129
      MUL2(z,ZERO,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
121
130
      ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2)
138
147
      i=(TWO52+TWO8*w)-TWO52;  i-=16;
139
148
      z=(w-cij[i][0].d)+ww;
140
149
      yy=HPI1-z*(cij[i][2].d+z*(cij[i][3].d+z*(cij[i][4].d+
141
 
                             z*(cij[i][5].d+z* cij[i][6].d))));
 
150
                             z*(cij[i][5].d+z* cij[i][6].d))));
142
151
      t1=HPI-cij[i][1].d;
143
152
      if (i<112)  u3=U31;  /* w <  1/2 */
144
153
      else        u3=U32;  /* w >= 1/2 */
148
157
      t1=w-hij[i][0].d;
149
158
      EADD(t1,ww,z,zz)
150
159
      s1=z*(hij[i][11].d+z*(hij[i][12].d+z*(hij[i][13].d+
151
 
         z*(hij[i][14].d+z* hij[i][15].d))));
 
160
         z*(hij[i][14].d+z* hij[i][15].d))));
152
161
      ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
153
162
      MUL2(z,zz,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
154
163
      ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2)
165
174
    }
166
175
    else {
167
176
      if (u<E) { /* D <= u < E */
168
 
        w=ONE/u;   v=w*w;
169
 
        EMULV(w,u,t1,t2,t3,t4,t5,t6,t7)
170
 
        yy=w*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
171
 
        ww=w*((ONE-t1)-t2);
172
 
        ESUB(HPI,w,t3,cor)
173
 
        yy=((HPI1+cor)-ww)-yy;
174
 
        if ((y=t3+(yy-U4)) == t3+(yy+U4))  return __signArctan(x,y);
 
177
        w=ONE/u;   v=w*w;
 
178
        EMULV(w,u,t1,t2,t3,t4,t5,t6,t7)
 
179
        yy=w*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
 
180
        ww=w*((ONE-t1)-t2);
 
181
        ESUB(HPI,w,t3,cor)
 
182
        yy=((HPI1+cor)-ww)-yy;
 
183
        if ((y=t3+(yy-U4)) == t3+(yy+U4))  return __signArctan(x,y);
175
184
 
176
 
        DIV2(ONE,ZERO,u,ZERO,w,ww,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
177
 
        MUL2(w,ww,w,ww,v,vv,t1,t2,t3,t4,t5,t6,t7,t8)
178
 
        s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d))));
179
 
        ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2)
180
 
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
181
 
        ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2)
182
 
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
183
 
        ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2)
184
 
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
185
 
        ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2)
186
 
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
187
 
        MUL2(w,ww,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
188
 
        ADD2(w,ww,s2,ss2,s1,ss1,t1,t2)
189
 
        SUB2(HPI,HPI1,s1,ss1,s2,ss2,t1,t2)
190
 
        if ((y=s2+(ss2-U8)) == s2+(ss2+U8))  return __signArctan(x,y);
 
185
        DIV2(ONE,ZERO,u,ZERO,w,ww,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
 
186
        MUL2(w,ww,w,ww,v,vv,t1,t2,t3,t4,t5,t6,t7,t8)
 
187
        s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d))));
 
188
        ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2)
 
189
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 
190
        ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2)
 
191
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 
192
        ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2)
 
193
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 
194
        ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2)
 
195
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 
196
        MUL2(w,ww,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
 
197
        ADD2(w,ww,s2,ss2,s1,ss1,t1,t2)
 
198
        SUB2(HPI,HPI1,s1,ss1,s2,ss2,t1,t2)
 
199
        if ((y=s2+(ss2-U8)) == s2+(ss2+U8))  return __signArctan(x,y);
191
200
 
192
201
      return atanMp(x,pr);
193
202
      }
194
203
      else {
195
 
        /* u >= E */
196
 
        if (x>0) return  HPI;
197
 
        else     return MHPI; }
 
204
        /* u >= E */
 
205
        if (x>0) return  HPI;
 
206
        else     return MHPI; }
198
207
    }
199
208
  }
200
209
 
201
210
}
202
211
 
203
 
 
204
 
  /* Fix the sign of y and return */
205
 
double  __signArctan(double x,double y){
206
 
 
207
 
    if (x<ZERO) return -y;
208
 
    else        return  y;
209
 
}
210
 
 
211
212
 /* Final stages. Compute atan(x) by multiple precision arithmetic */
212
213
static double atanMp(double x,const int pr[]){
213
214
  mp_no mpx,mpy,mpy2,mperr,mpt1,mpy1;