← Back to branch summary

~ubuntu-branches/ubuntu/intrepid/parted/intrepid

« back to all changes in this revision

Viewing changes to gnulib/lib/tanl.c

  • Committer: Bazaar Package Importer
  • Author(s): Colin Watson
  • Date: 2008-06-24 14:31:05 UTC
  • mfrom: (1.1.3 upstream)
  • Revision ID: james.westby@ubuntu.com-20080624143105-rd7yw67a9qnvh51i
Tags: 1.8.8.git.2008.03.24-7ubuntu1
https://launchpad.net/bugs/237568
* Resynchronise with Debian (LP: #237568). Remaining changes:
  - swap-uuid.dpatch: Create UUIDs on new swap partitions.
  - gptsync.dpatch: On Intel Mac systems, write a synced MBR rather than a
    protective MBR.
  - Add -fno-stack-protector on sparc.
  - sparc-new-label.dpatch: Fix sparc disk label generation. This is
    required for LDOM and parallel installations with Solaris 10.
  - loop-partitions.dpatch: Loop devices can only have one partition, so
    don't generate device names such as "/dev/loop0p1".
  - unpartitioned-disks.dpatch: Don't try to call BLKPG ioctls on
    unpartitionable disks (only implemented for loop devices at the
    moment), as they will always fail.
  - When building with gcc-4.3, add -Wno-array-bounds to CFLAGS.
  - Cell partition tables are misdetected as pc98, so disable pc98 support
    on powerpc.
  - array-bounds.dpatch: Backport patch from git to allow building with
    gcc-4.3.

Show diffs side-by-side

added added

removed removed

Lines of Context:
gnulib/lib/tanl.c
 
1
/* s_tanl.c -- long double version of s_tan.c.
 
2
 * Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz.
 
3
 */
 
4
 
 
5
/* @(#)s_tan.c 5.1 93/09/24 */
 
6
/*
 
7
 * ====================================================
 
8
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 
9
 *
 
10
 * Developed at SunPro, a Sun Microsystems, Inc. business.
 
11
 * Permission to use, copy, modify, and distribute this
 
12
 * software is freely granted, provided that this notice
 
13
 * is preserved.
 
14
 * ====================================================
 
15
 */
 
16
 
 
17
#include <config.h>
 
18
 
 
19
/* Specification.  */
 
20
#include <math.h>
 
21
 
 
22
/* tanl(x)
 
23
 * Return tangent function of x.
 
24
 *
 
25
 * kernel function:
 
26
 *      __kernel_tanl           ... tangent function on [-pi/4,pi/4]
 
27
 *      __ieee754_rem_pio2l     ... argument reduction routine
 
28
 *
 
29
 * Method.
 
30
 *      Let S,C and T denote the sin, cos and tan respectively on
 
31
 *      [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
 
32
 *      in [-pi/4 , +pi/4], and let n = k mod 4.
 
33
 *      We have
 
34
 *
 
35
 *          n        sin(x)      cos(x)        tan(x)
 
36
 *     ----------------------------------------------------------
 
37
 *          0          S           C             T
 
38
 *          1          C          -S            -1/T
 
39
 *          2         -S          -C             T
 
40
 *          3         -C           S            -1/T
 
41
 *     ----------------------------------------------------------
 
42
 *
 
43
 * Special cases:
 
44
 *      Let trig be any of sin, cos, or tan.
 
45
 *      trig(+-INF)  is NaN, with signals;
 
46
 *      trig(NaN)    is that NaN;
 
47
 *
 
48
 * Accuracy:
 
49
 *      TRIG(x) returns trig(x) nearly rounded
 
50
 */
 
51
 
 
52
#include "trigl.h"
 
53
#ifdef HAVE_SINL
 
54
#ifdef HAVE_COSL
 
55
#include "trigl.c"
 
56
#endif
 
57
#endif
 
58
#include "isnanl.h"
 
59
 
 
60
/*
 
61
 * ====================================================
 
62
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 
63
 *
 
64
 * Developed at SunPro, a Sun Microsystems, Inc. business.
 
65
 * Permission to use, copy, modify, and distribute this
 
66
 * software is freely granted, provided that this notice
 
67
 * is preserved.
 
68
 * ====================================================
 
69
 */
 
70
 
 
71
/*
 
72
  Long double expansions contributed by
 
73
  Stephen L. Moshier <moshier@na-net.ornl.gov>
 
74
*/
 
75
 
 
76
/* __kernel_tanl( x, y, k )
 
77
 * kernel tan function on [-pi/4, pi/4], pi/4 ~ 0.7854
 
78
 * Input x is assumed to be bounded by ~pi/4 in magnitude.
 
79
 * Input y is the tail of x.
 
80
 * Input k indicates whether tan (if k=1) or
 
81
 * -1/tan (if k= -1) is returned.
 
82
 *
 
83
 * Algorithm
 
84
 *      1. Since tan(-x) = -tan(x), we need only to consider positive x.
 
85
 *      2. if x < 2^-57, return x with inexact if x!=0.
 
86
 *      3. tan(x) is approximated by a rational form x + x^3 / 3 + x^5 R(x^2)
 
87
 *          on [0,0.67433].
 
88
 *
 
89
 *         Note: tan(x+y) = tan(x) + tan'(x)*y
 
90
 *                        ~ tan(x) + (1+x*x)*y
 
91
 *         Therefore, for better accuracy in computing tan(x+y), let
 
92
 *              r = x^3 * R(x^2)
 
93
 *         then
 
94
 *              tan(x+y) = x + (x^3 / 3 + (x^2 *(r+y)+y))
 
95
 *
 
96
 *      4. For x in [0.67433,pi/4],  let y = pi/4 - x, then
 
97
 *              tan(x) = tan(pi/4-y) = (1-tan(y))/(1+tan(y))
 
98
 *                     = 1 - 2*(tan(y) - (tan(y)^2)/(1+tan(y)))
 
99
 */
 
100
 
 
101
 
 
102
static const long double
 
103
  pio4hi = 7.8539816339744830961566084581987569936977E-1L,
 
104
  pio4lo = 2.1679525325309452561992610065108379921906E-35L,
 
105
 
 
106
  /* tan x = x + x^3 / 3 + x^5 T(x^2)/U(x^2)
 
107
     0 <= x <= 0.6743316650390625
 
108
     Peak relative error 8.0e-36  */
 
109
 TH =  3.333333333333333333333333333333333333333E-1L,
 
110
 T0 = -1.813014711743583437742363284336855889393E7L,
 
111
 T1 =  1.320767960008972224312740075083259247618E6L,
 
112
 T2 = -2.626775478255838182468651821863299023956E4L,
 
113
 T3 =  1.764573356488504935415411383687150199315E2L,
 
114
 T4 = -3.333267763822178690794678978979803526092E-1L,
 
115
 
 
116
 U0 = -1.359761033807687578306772463253710042010E8L,
 
117
 U1 =  6.494370630656893175666729313065113194784E7L,
 
118
 U2 = -4.180787672237927475505536849168729386782E6L,
 
119
 U3 =  8.031643765106170040139966622980914621521E4L,
 
120
 U4 = -5.323131271912475695157127875560667378597E2L;
 
121
  /* 1.000000000000000000000000000000000000000E0 */
 
122
 
 
123
 
 
124
long double
 
125
kernel_tanl (long double x, long double y, int iy)
 
126
{
 
127
  long double z, r, v, w, s, u, u1;
 
128
  int flag, sign;
 
129
 
 
130
  sign = 1;
 
131
  if (x < 0)
 
132
    {
 
133
      x = -x;
 
134
      y = -y;
 
135
      sign = -1;
 
136
    }
 
137
 
 
138
  if (x < 0.000000000000000006938893903907228377647697925567626953125L) /* x < 2**-57 */
 
139
    {
 
140
      if ((int) x == 0)
 
141
        {                       /* generate inexact */
 
142
          if (iy == -1 && x == 0.0)
 
143
            return 1.0L / fabs (x);
 
144
          else
 
145
            return (iy == 1) ? x : -1.0L / x;
 
146
        }
 
147
    }
 
148
  if (x >= 0.6743316650390625) /* |x| >= 0.6743316650390625 */
 
149
    {
 
150
      flag = 1;
 
151
 
 
152
      z = pio4hi - x;
 
153
      w = pio4lo - y;
 
154
      x = z + w;
 
155
      y = 0.0;
 
156
    }
 
157
  z = x * x;
 
158
  r = T0 + z * (T1 + z * (T2 + z * (T3 + z * T4)));
 
159
  v = U0 + z * (U1 + z * (U2 + z * (U3 + z * (U4 + z))));
 
160
  r = r / v;
 
161
 
 
162
  s = z * x;
 
163
  r = y + z * (s * r + y);
 
164
  r += TH * s;
 
165
  w = x + r;
 
166
  if (flag)
 
167
    {
 
168
      v = (long double) iy;
 
169
      w = (v - 2.0 * (x - (w * w / (w + v) - r)));
 
170
      if (sign < 0)
 
171
        w = -w;
 
172
      return w;
 
173
    }
 
174
  if (iy == 1)
 
175
    return w;
 
176
  else
 
177
    {                           /* if allow error up to 2 ulp,
 
178
                                   simply return -1.0/(x+r) here */
 
179
      /*  compute -1.0/(x+r) accurately */
 
180
      u1 = (double) w;
 
181
      v = r - (u1 - x);
 
182
      z = -1.0 / w;
 
183
      u = (double) z;
 
184
      s = 1.0 + u * u1;
 
185
      return u + z * (s + u * v);
 
186
    }
 
187
}
 
188
 
 
189
long double
 
190
tanl (long double x)
 
191
{
 
192
  long double y[2], z = 0.0L;
 
193
  int n;
 
194
 
 
195
  /* tanl(NaN) is NaN */
 
196
  if (isnanl (x))
 
197
    return x;
 
198
 
 
199
  /* |x| ~< pi/4 */
 
200
  if (x >= -0.7853981633974483096156608458198757210492 &&
 
201
      x <= 0.7853981633974483096156608458198757210492)
 
202
    return kernel_tanl (x, z, 1);
 
203
 
 
204
  /* tanl(Inf) is NaN, tanl(0) is 0 */
 
205
  else if (x + x == x)
 
206
    return x - x;               /* NaN */
 
207
 
 
208
  /* argument reduction needed */
 
209
  else
 
210
    {
 
211
      n = ieee754_rem_pio2l (x, y);
 
212
      /* 1 -- n even, -1 -- n odd */
 
213
      return kernel_tanl (y[0], y[1], 1 - ((n & 1) << 1));
 
214
    }
 
215
}
 
216
 
 
217
#if 0
 
218
int
 
219
main (void)
 
220
{
 
221
  printf ("%.16Lg\n", tanl(0.7853981633974483096156608458198757210492));
 
222
  printf ("%.16Lg\n", tanl(-0.7853981633974483096156608458198757210492));
 
223
  printf ("%.16Lg\n", tanl(0.7853981633974483096156608458198757210492 *3));
 
224
  printf ("%.16Lg\n", tanl(-0.7853981633974483096156608458198757210492 *31));
 
225
  printf ("%.16Lg\n", tanl(0.7853981633974483096156608458198757210492 / 2));
 
226
  printf ("%.16Lg\n", tanl(0.7853981633974483096156608458198757210492 * 3/2));
 
227
  printf ("%.16Lg\n", tanl(0.7853981633974483096156608458198757210492 * 5/2));
 
228
}
 
229
#endif