← Back to branch summary

~ubuntu-branches/debian/sid/genius/sid

« back to all changes in this revision

Viewing changes to mpfr/sin.c

  • Committer: Bazaar Package Importer
  • Author(s): Daniel Holbach
  • Date: 2006-08-21 12:57:45 UTC
  • Revision ID: james.westby@ubuntu.com-20060821125745-sl9ks8v7fq324bdf
Tags: upstream-0.7.6.1
Import upstream version 0.7.6.1

Show diffs side-by-side

added added

removed removed

Lines of Context:
mpfr/sin.c
 
1
/* mpfr_sin -- sine of a floating-point number
 
2
 
 
3
Copyright 2001, 2002, 2003, 2004, 2005 Free Software Foundation, Inc.
 
4
 
 
5
This file is part of the MPFR Library.
 
6
 
 
7
The MPFR Library is free software; you can redistribute it and/or modify
 
8
it under the terms of the GNU Lesser General Public License as published by
 
9
the Free Software Foundation; either version 2.1 of the License, or (at your
 
10
option) any later version.
 
11
 
 
12
The MPFR Library is distributed in the hope that it will be useful, but
 
13
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
 
14
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
 
15
License for more details.
 
16
 
 
17
You should have received a copy of the GNU Lesser General Public License
 
18
along with the MPFR Library; see the file COPYING.LIB.  If not, write to
 
19
the Free Software Foundation, Inc., 51 Franklin Place, Fifth Floor, Boston,
 
20
MA 02110-1301, USA. */
 
21
 
 
22
#define MPFR_NEED_LONGLONG_H
 
23
#include "mpfr-impl.h"
 
24
 
 
25
/* determine the sign of sin(x) using argument reduction.
 
26
   Assumes x is not an exact multiple of Pi (this excludes x=0). */
 
27
static int
 
28
mpfr_sin_sign (mpfr_srcptr x)
 
29
{
 
30
  mpfr_t c, k;
 
31
  mp_exp_t K;
 
32
  int sign;
 
33
  mp_prec_t m;
 
34
  mpfr_srcptr y;
 
35
  MPFR_ZIV_DECL (loop);
 
36
 
 
37
  K = MPFR_GET_EXP(x);
 
38
  if (K < 0)  /* Trivial case if ABS(x) < 1 */
 
39
    return MPFR_SIGN (x);
 
40
 
 
41
  m = K + BITS_PER_MP_LIMB;
 
42
  mpfr_init2 (c, m);
 
43
  mpfr_init2 (k, m);
 
44
 
 
45
  MPFR_ZIV_INIT (loop, m);
 
46
  for (;;)
 
47
    {
 
48
      /* first determine round(x/Pi): does not have to be exact since
 
49
         the result is an integer */
 
50
      mpfr_const_pi (c, GMP_RNDN); /* err <= 1/2*ulp(c) = 2^(1-m) */
 
51
      /* we need that k is not-to-badly rounded to an integer,
 
52
         i.e. ulp(k) <= 1, so m >= EXP(k). */
 
53
      mpfr_div (k, x, c, GMP_RNDN);
 
54
      mpfr_round (k, k);
 
55
 
 
56
      sign = 1;
 
57
 
 
58
      if (!MPFR_IS_ZERO (k)) /* subtract k*approx(Pi) */
 
59
        {
 
60
          /* determine parity of k for sign */
 
61
          if (MPFR_GET_EXP (k) <= 0 || (mpfr_uexp_t) MPFR_GET_EXP (k) <= m)
 
62
            {
 
63
              mp_size_t j = BITS_PER_MP_LIMB * MPFR_LIMB_SIZE(k)
 
64
                - MPFR_GET_EXP(k);
 
65
              mp_size_t l = j / BITS_PER_MP_LIMB;
 
66
              /* parity bit is j-th bit starting from least significant bits */
 
67
              if ((MPFR_MANT(k)[l] >> (j % BITS_PER_MP_LIMB)) & 1)
 
68
                sign = -1; /* k is odd */
 
69
            }
 
70
          K = MPFR_GET_EXP (k); /* k is an integer, thus K >= 1, k < 2^K */
 
71
          mpfr_mul (k, k, c, GMP_RNDN); /* err <= oldk*err(c) + 1/2*ulp(k)
 
72
                                               <= 2^(K+2-m) */
 
73
          mpfr_sub (k, x, k, GMP_RNDN);
 
74
          /* assuming |k| <= Pi, err <= 2^(1-m)+2^(K+2-m) < 2^(K+3-m) */
 
75
          MPFR_ASSERTN (MPFR_IS_ZERO (k) || MPFR_GET_EXP (k) <= 2);
 
76
          y = k;
 
77
        }
 
78
      else
 
79
        {
 
80
          K = 1;
 
81
          y = x;
 
82
        }
 
83
      /* sign of sign(y) is uncertain if |y| <= err < 2^(K+3-m),
 
84
         thus EXP(y) < K+4-m */
 
85
      if (MPFR_LIKELY (!MPFR_IS_ZERO (y)
 
86
                       && MPFR_GET_EXP (y) >= K + 4 - (mp_exp_t) m))
 
87
        break;
 
88
      MPFR_ZIV_NEXT (loop, m);
 
89
      mpfr_set_prec (c, m);
 
90
      mpfr_set_prec (k, m);
 
91
    }
 
92
 
 
93
  if (MPFR_IS_NEG (y))
 
94
    sign = -sign;
 
95
 
 
96
  mpfr_clear (k);
 
97
  mpfr_clear (c);
 
98
 
 
99
  return sign;
 
100
}
 
101
 
 
102
int
 
103
mpfr_sin (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode)
 
104
{
 
105
  mpfr_t c;
 
106
  mp_exp_t e;
 
107
  mp_prec_t precy, m;
 
108
  int inexact, sign;
 
109
  MPFR_ZIV_DECL (loop);
 
110
 
 
111
  MPFR_LOG_FUNC (("x[%#R]=%R rnd=%d", x, x, rnd_mode),
 
112
                  ("y[%#R]=%R inexact=%d", y, y, inexact));
 
113
 
 
114
  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
 
115
    {
 
116
      if (MPFR_IS_NAN (x) || MPFR_IS_INF (x))
 
117
        {
 
118
          MPFR_SET_NAN (y);
 
119
          MPFR_RET_NAN;
 
120
 
 
121
        }
 
122
      else /* x is zero */
 
123
        {
 
124
          MPFR_ASSERTD (MPFR_IS_ZERO (x));
 
125
          MPFR_SET_ZERO (y);
 
126
          MPFR_SET_SAME_SIGN (y, x);
 
127
          MPFR_RET (0);
 
128
        }
 
129
    }
 
130
 
 
131
  /* sin(x) = x - x^3/6 + ... so the error is < 2^(3*EXP(x)-2) */
 
132
  MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, -2*MPFR_GET_EXP (x)+2,0,rnd_mode,{});
 
133
 
 
134
  /* Compute initial precision */
 
135
  precy = MPFR_PREC (y);
 
136
  m = precy + MPFR_INT_CEIL_LOG2 (precy) + 13;
 
137
  e = MPFR_GET_EXP (x);
 
138
  m += (e < 0) ? -2*e : e;
 
139
 
 
140
  sign = mpfr_sin_sign (x);
 
141
  mpfr_init2 (c, m);
 
142
 
 
143
  MPFR_ZIV_INIT (loop, m);
 
144
  for (;;)
 
145
    {
 
146
      mpfr_cos (c, x, GMP_RNDZ);    /* can't be exact */
 
147
      mpfr_nexttoinf (c);           /* now c = cos(x) rounded away */
 
148
      mpfr_mul (c, c, c, GMP_RNDU); /* away */
 
149
      mpfr_ui_sub (c, 1, c, GMP_RNDZ);
 
150
      mpfr_sqrt (c, c, GMP_RNDZ);
 
151
      if (MPFR_IS_NEG_SIGN(sign))
 
152
        MPFR_CHANGE_SIGN(c);
 
153
 
 
154
      /* Warning c may be 0 ! */
 
155
      if (MPFR_UNLIKELY (MPFR_IS_ZERO (c)))
 
156
        {
 
157
          /* Huge cancellation: increase prec a lot! */
 
158
          m = MAX (m, MPFR_PREC (x));
 
159
          m = 2*m;
 
160
        }
 
161
      else
 
162
        {
 
163
          /* the absolute error on c is at most 2^(3-m-EXP(c)) */
 
164
          e = 2 * MPFR_GET_EXP (c) + m - 3;
 
165
          if (mpfr_can_round (c, e, GMP_RNDN, GMP_RNDZ,
 
166
                              precy + (rnd_mode == GMP_RNDN)))
 
167
            /* WARNING: even if we know c <= sin(x), don't give GMP_RNDZ
 
168
               as 3rd argument to mpfr_can_round, since if c is exactly
 
169
               representable to the target precision (inexact = 0 below),
 
170
               we would have to add one ulp when rounding away from 0. */
 
171
            break;
 
172
 
 
173
          /* check for huge cancellation (Near 0) */
 
174
          if (e < (mp_exp_t) MPFR_PREC (y))
 
175
            m += MPFR_PREC (y) - e;
 
176
          /* Check if near 1 */
 
177
          if (MPFR_GET_EXP (c) == 1)
 
178
            m += m;
 
179
        }
 
180
 
 
181
      /* Else generic increase */
 
182
      MPFR_ZIV_NEXT (loop, m);
 
183
      mpfr_set_prec (c, m);
 
184
    }
 
185
  MPFR_ZIV_FREE (loop);
 
186
 
 
187
  inexact = mpfr_set (y, c, rnd_mode);
 
188
  /* inexact cannot be 0, since this would mean that c was representable
 
189
     within the target precision, but in that case mpfr_can_round will fail */
 
190
 
 
191
  mpfr_clear (c);
 
192
 
 
193
  return inexact; /* inexact */
 
194
}