← Back to branch summary

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

« back to all changes in this revision

Viewing changes to sysdeps/ieee754/ldbl-128ibm/e_acosl.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/ldbl-128ibm/e_acosl.c
 
1
/*
 
2
 * ====================================================
 
3
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 
4
 *
 
5
 * Developed at SunPro, a Sun Microsystems, Inc. business.
 
6
 * Permission to use, copy, modify, and distribute this
 
7
 * software is freely granted, provided that this notice
 
8
 * is preserved.
 
9
 * ====================================================
 
10
 */
 
11
 
 
12
/*
 
13
   Long double expansions are
 
14
   Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
 
15
   and are incorporated herein by permission of the author.  The author
 
16
   reserves the right to distribute this material elsewhere under different
 
17
   copying permissions.  These modifications are distributed here under
 
18
   the following terms:
 
19
 
 
20
    This library is free software; you can redistribute it and/or
 
21
    modify it under the terms of the GNU Lesser General Public
 
22
    License as published by the Free Software Foundation; either
 
23
    version 2.1 of the License, or (at your option) any later version.
 
24
 
 
25
    This library is distributed in the hope that it will be useful,
 
26
    but WITHOUT ANY WARRANTY; without even the implied warranty of
 
27
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 
28
    Lesser General Public License for more details.
 
29
 
 
30
    You should have received a copy of the GNU Lesser General Public
 
31
    License along with this library; if not, write to the Free Software
 
32
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA */
 
33
 
 
34
/* __ieee754_acosl(x)
 
35
 * Method :
 
36
 *      acos(x)  = pi/2 - asin(x)
 
37
 *      acos(-x) = pi/2 + asin(x)
 
38
 * For |x| <= 0.375
 
39
 *      acos(x) = pi/2 - asin(x)
 
40
 * Between .375 and .5 the approximation is
 
41
 *      acos(0.4375 + x) = acos(0.4375) + x P(x) / Q(x)
 
42
 * Between .5 and .625 the approximation is
 
43
 *      acos(0.5625 + x) = acos(0.5625) + x rS(x) / sS(x)
 
44
 * For x > 0.625,
 
45
 *      acos(x) = 2 asin(sqrt((1-x)/2))
 
46
 *      computed with an extended precision square root in the leading term.
 
47
 * For x < -0.625
 
48
 *      acos(x) = pi - 2 asin(sqrt((1-|x|)/2))
 
49
 *
 
50
 * Special cases:
 
51
 *      if x is NaN, return x itself;
 
52
 *      if |x|>1, return NaN with invalid signal.
 
53
 *
 
54
 * Functions needed: __ieee754_sqrtl.
 
55
 */
 
56
 
 
57
#include "math.h"
 
58
#include "math_private.h"
 
59
 
 
60
#ifdef __STDC__
 
61
static const long double
 
62
#else
 
63
static long double
 
64
#endif
 
65
  one = 1.0L,
 
66
  pio2_hi = 1.5707963267948966192313216916397514420986L,
 
67
  pio2_lo = 4.3359050650618905123985220130216759843812E-35L,
 
68
 
 
69
  /* acos(0.5625 + x) = acos(0.5625) + x rS(x) / sS(x)
 
70
     -0.0625 <= x <= 0.0625
 
71
     peak relative error 3.3e-35  */
 
72
 
 
73
  rS0 =  5.619049346208901520945464704848780243887E0L,
 
74
  rS1 = -4.460504162777731472539175700169871920352E1L,
 
75
  rS2 =  1.317669505315409261479577040530751477488E2L,
 
76
  rS3 = -1.626532582423661989632442410808596009227E2L,
 
77
  rS4 =  3.144806644195158614904369445440583873264E1L,
 
78
  rS5 =  9.806674443470740708765165604769099559553E1L,
 
79
  rS6 = -5.708468492052010816555762842394927806920E1L,
 
80
  rS7 = -1.396540499232262112248553357962639431922E1L,
 
81
  rS8 =  1.126243289311910363001762058295832610344E1L,
 
82
  rS9 =  4.956179821329901954211277873774472383512E-1L,
 
83
  rS10 = -3.313227657082367169241333738391762525780E-1L,
 
84
 
 
85
  sS0 = -4.645814742084009935700221277307007679325E0L,
 
86
  sS1 =  3.879074822457694323970438316317961918430E1L,
 
87
  sS2 = -1.221986588013474694623973554726201001066E2L,
 
88
  sS3 =  1.658821150347718105012079876756201905822E2L,
 
89
  sS4 = -4.804379630977558197953176474426239748977E1L,
 
90
  sS5 = -1.004296417397316948114344573811562952793E2L,
 
91
  sS6 =  7.530281592861320234941101403870010111138E1L,
 
92
  sS7 =  1.270735595411673647119592092304357226607E1L,
 
93
  sS8 = -1.815144839646376500705105967064792930282E1L,
 
94
  sS9 = -7.821597334910963922204235247786840828217E-2L,
 
95
  /* 1.000000000000000000000000000000000000000E0 */
 
96
 
 
97
  acosr5625 = 9.7338991014954640492751132535550279812151E-1L,
 
98
  pimacosr5625 = 2.1682027434402468335351320579240000860757E0L,
 
99
 
 
100
  /* acos(0.4375 + x) = acos(0.4375) + x rS(x) / sS(x)
 
101
     -0.0625 <= x <= 0.0625
 
102
     peak relative error 2.1e-35  */
 
103
 
 
104
  P0 =  2.177690192235413635229046633751390484892E0L,
 
105
  P1 = -2.848698225706605746657192566166142909573E1L,
 
106
  P2 =  1.040076477655245590871244795403659880304E2L,
 
107
  P3 = -1.400087608918906358323551402881238180553E2L,
 
108
  P4 =  2.221047917671449176051896400503615543757E1L,
 
109
  P5 =  9.643714856395587663736110523917499638702E1L,
 
110
  P6 = -5.158406639829833829027457284942389079196E1L,
 
111
  P7 = -1.578651828337585944715290382181219741813E1L,
 
112
  P8 =  1.093632715903802870546857764647931045906E1L,
 
113
  P9 =  5.448925479898460003048760932274085300103E-1L,
 
114
  P10 = -3.315886001095605268470690485170092986337E-1L,
 
115
  Q0 = -1.958219113487162405143608843774587557016E0L,
 
116
  Q1 =  2.614577866876185080678907676023269360520E1L,
 
117
  Q2 = -9.990858606464150981009763389881793660938E1L,
 
118
  Q3 =  1.443958741356995763628660823395334281596E2L,
 
119
  Q4 = -3.206441012484232867657763518369723873129E1L,
 
120
  Q5 = -1.048560885341833443564920145642588991492E2L,
 
121
  Q6 =  6.745883931909770880159915641984874746358E1L,
 
122
  Q7 =  1.806809656342804436118449982647641392951E1L,
 
123
  Q8 = -1.770150690652438294290020775359580915464E1L,
 
124
  Q9 = -5.659156469628629327045433069052560211164E-1L,
 
125
  /* 1.000000000000000000000000000000000000000E0 */
 
126
 
 
127
  acosr4375 = 1.1179797320499710475919903296900511518755E0L,
 
128
  pimacosr4375 = 2.0236129215398221908706530535894517323217E0L,
 
129
 
 
130
  /* asin(x) = x + x^3 pS(x^2) / qS(x^2)
 
131
     0 <= x <= 0.5
 
132
     peak relative error 1.9e-35  */
 
133
  pS0 = -8.358099012470680544198472400254596543711E2L,
 
134
  pS1 =  3.674973957689619490312782828051860366493E3L,
 
135
  pS2 = -6.730729094812979665807581609853656623219E3L,
 
136
  pS3 =  6.643843795209060298375552684423454077633E3L,
 
137
  pS4 = -3.817341990928606692235481812252049415993E3L,
 
138
  pS5 =  1.284635388402653715636722822195716476156E3L,
 
139
  pS6 = -2.410736125231549204856567737329112037867E2L,
 
140
  pS7 =  2.219191969382402856557594215833622156220E1L,
 
141
  pS8 = -7.249056260830627156600112195061001036533E-1L,
 
142
  pS9 =  1.055923570937755300061509030361395604448E-3L,
 
143
 
 
144
  qS0 = -5.014859407482408326519083440151745519205E3L,
 
145
  qS1 =  2.430653047950480068881028451580393430537E4L,
 
146
  qS2 = -4.997904737193653607449250593976069726962E4L,
 
147
  qS3 =  5.675712336110456923807959930107347511086E4L,
 
148
  qS4 = -3.881523118339661268482937768522572588022E4L,
 
149
  qS5 =  1.634202194895541569749717032234510811216E4L,
 
150
  qS6 = -4.151452662440709301601820849901296953752E3L,
 
151
  qS7 =  5.956050864057192019085175976175695342168E2L,
 
152
  qS8 = -4.175375777334867025769346564600396877176E1L;
 
153
  /* 1.000000000000000000000000000000000000000E0 */
 
154
 
 
155
#ifdef __STDC__
 
156
long double
 
157
__ieee754_acosl (long double x)
 
158
#else
 
159
long double
 
160
__ieee754_acosl (x)
 
161
     long double x;
 
162
#endif
 
163
{
 
164
  long double z, r, w, p, q, s, t, f2;
 
165
  int32_t ix, sign;
 
166
  ieee854_long_double_shape_type u;
 
167
 
 
168
  u.value = x;
 
169
  sign = u.parts32.w0;
 
170
  ix = sign & 0x7fffffff;
 
171
  u.parts32.w0 = ix;            /* |x| */
 
172
  if (ix >= 0x3ff00000)         /* |x| >= 1 */
 
173
    {
 
174
      if (ix == 0x3ff00000
 
175
          && (u.parts32.w1 | (u.parts32.w2&0x7fffffff) | u.parts32.w3) == 0)
 
176
        {                       /* |x| == 1 */
 
177
          if ((sign & 0x80000000) == 0)
 
178
            return 0.0;         /* acos(1) = 0  */
 
179
          else
 
180
            return (2.0 * pio2_hi) + (2.0 * pio2_lo);   /* acos(-1)= pi */
 
181
        }
 
182
      return (x - x) / (x - x); /* acos(|x| > 1) is NaN */
 
183
    }
 
184
  else if (ix < 0x3fe00000)     /* |x| < 0.5 */
 
185
    {
 
186
      if (ix < 0x3c600000)      /* |x| < 2**-57 */
 
187
        return pio2_hi + pio2_lo;
 
188
      if (ix < 0x3fde0000)      /* |x| < .4375 */
 
189
        {
 
190
          /* Arcsine of x.  */
 
191
          z = x * x;
 
192
          p = (((((((((pS9 * z
 
193
                       + pS8) * z
 
194
                      + pS7) * z
 
195
                     + pS6) * z
 
196
                    + pS5) * z
 
197
                   + pS4) * z
 
198
                  + pS3) * z
 
199
                 + pS2) * z
 
200
                + pS1) * z
 
201
               + pS0) * z;
 
202
          q = (((((((( z
 
203
                       + qS8) * z
 
204
                     + qS7) * z
 
205
                    + qS6) * z
 
206
                   + qS5) * z
 
207
                  + qS4) * z
 
208
                 + qS3) * z
 
209
                + qS2) * z
 
210
               + qS1) * z
 
211
            + qS0;
 
212
          r = x + x * p / q;
 
213
          z = pio2_hi - (r - pio2_lo);
 
214
          return z;
 
215
        }
 
216
      /* .4375 <= |x| < .5 */
 
217
      t = u.value - 0.4375L;
 
218
      p = ((((((((((P10 * t
 
219
                    + P9) * t
 
220
                   + P8) * t
 
221
                  + P7) * t
 
222
                 + P6) * t
 
223
                + P5) * t
 
224
               + P4) * t
 
225
              + P3) * t
 
226
             + P2) * t
 
227
            + P1) * t
 
228
           + P0) * t;
 
229
 
 
230
      q = (((((((((t
 
231
                   + Q9) * t
 
232
                  + Q8) * t
 
233
                 + Q7) * t
 
234
                + Q6) * t
 
235
               + Q5) * t
 
236
              + Q4) * t
 
237
             + Q3) * t
 
238
            + Q2) * t
 
239
           + Q1) * t
 
240
        + Q0;
 
241
      r = p / q;
 
242
      if (sign & 0x80000000)
 
243
        r = pimacosr4375 - r;
 
244
      else
 
245
        r = acosr4375 + r;
 
246
      return r;
 
247
    }
 
248
  else if (ix < 0x3fe40000)     /* |x| < 0.625 */
 
249
    {
 
250
      t = u.value - 0.5625L;
 
251
      p = ((((((((((rS10 * t
 
252
                    + rS9) * t
 
253
                   + rS8) * t
 
254
                  + rS7) * t
 
255
                 + rS6) * t
 
256
                + rS5) * t
 
257
               + rS4) * t
 
258
              + rS3) * t
 
259
             + rS2) * t
 
260
            + rS1) * t
 
261
           + rS0) * t;
 
262
 
 
263
      q = (((((((((t
 
264
                   + sS9) * t
 
265
                  + sS8) * t
 
266
                 + sS7) * t
 
267
                + sS6) * t
 
268
               + sS5) * t
 
269
              + sS4) * t
 
270
             + sS3) * t
 
271
            + sS2) * t
 
272
           + sS1) * t
 
273
        + sS0;
 
274
      if (sign & 0x80000000)
 
275
        r = pimacosr5625 - p / q;
 
276
      else
 
277
        r = acosr5625 + p / q;
 
278
      return r;
 
279
    }
 
280
  else
 
281
    {                           /* |x| >= .625 */
 
282
      z = (one - u.value) * 0.5;
 
283
      s = __ieee754_sqrtl (z);
 
284
      /* Compute an extended precision square root from
 
285
         the Newton iteration  s -> 0.5 * (s + z / s).
 
286
         The change w from s to the improved value is
 
287
            w = 0.5 * (s + z / s) - s  = (s^2 + z)/2s - s = (z - s^2)/2s.
 
288
          Express s = f1 + f2 where f1 * f1 is exactly representable.
 
289
          w = (z - s^2)/2s = (z - f1^2 - 2 f1 f2 - f2^2)/2s .
 
290
          s + w has extended precision.  */
 
291
      u.value = s;
 
292
      u.parts32.w2 = 0;
 
293
      u.parts32.w3 = 0;
 
294
      f2 = s - u.value;
 
295
      w = z - u.value * u.value;
 
296
      w = w - 2.0 * u.value * f2;
 
297
      w = w - f2 * f2;
 
298
      w = w / (2.0 * s);
 
299
      /* Arcsine of s.  */
 
300
      p = (((((((((pS9 * z
 
301
                   + pS8) * z
 
302
                  + pS7) * z
 
303
                 + pS6) * z
 
304
                + pS5) * z
 
305
               + pS4) * z
 
306
              + pS3) * z
 
307
             + pS2) * z
 
308
            + pS1) * z
 
309
           + pS0) * z;
 
310
      q = (((((((( z
 
311
                   + qS8) * z
 
312
                 + qS7) * z
 
313
                + qS6) * z
 
314
               + qS5) * z
 
315
              + qS4) * z
 
316
             + qS3) * z
 
317
            + qS2) * z
 
318
           + qS1) * z
 
319
        + qS0;
 
320
      r = s + (w + s * p / q);
 
321
 
 
322
      if (sign & 0x80000000)
 
323
        w = pio2_hi + (pio2_lo - r);
 
324
      else
 
325
        w = r;
 
326
      return 2.0 * w;
 
327
    }
 
328
}