1
/* mpfr_cos -- cosine of a floating-point number
3
Copyright 2001, 2002 Free Software Foundation.
5
This file is part of the MPFR Library.
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.
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.
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., 59 Temple Place - Suite 330, Boston,
20
MA 02111-1307, USA. */
26
#include "mpfr-impl.h"
28
static int mpfr_cos2_aux _PROTO ((mpfr_ptr, mpfr_srcptr));
31
mpfr_cos (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode)
33
int K0, K, precy, m, k, l, inexact;
36
if (MPFR_IS_NAN(x) || MPFR_IS_INF(x))
44
mpfr_set_ui (y, 1, GMP_RNDN);
50
K0 = _mpfr_isqrt(precy / 2);
51
/* we need at least K + log2(precy/K) extra bits */
52
m = precy + 3 * K0 + 3;
59
mpfr_mul (r, x, x, GMP_RNDU); /* err <= 1 ulp */
61
/* we need that |r| < 1 for mpfr_cos2_aux, i.e. up(x^2)/2^(2K) < 1 */
62
K = K0 + MAX(MPFR_EXP(r), 0);
64
mpfr_div_2ui (r, r, 2 * K, GMP_RNDN); /* r = (x/2^K)^2, err <= 1 ulp */
66
/* s <- 1 - r/2! + ... + (-1)^l r^l/(2l)! */
67
l = mpfr_cos2_aux (s, r);
69
for (k = 0; k < K; k++)
71
mpfr_mul (s, s, s, GMP_RNDU); /* err <= 2*olderr */
72
mpfr_mul_2ui (s, s, 1, GMP_RNDU); /* err <= 4*olderr */
73
mpfr_sub_ui (s, s, 1, GMP_RNDN);
76
/* absolute error on s is bounded by (2l+1/3)*2^(2K-m) */
77
for (k = 2 * K, l = 2 * l + 1; l > 1; k++, l = (l + 1) >> 1);
78
/* now the error is bounded by 2^(k-m) = 2^(EXP(s)-err) */
80
l = mpfr_can_round (s, MPFR_EXP(s) + m - k, GMP_RNDN, rnd_mode, precy);
84
m += BITS_PER_MP_LIMB;
91
inexact = mpfr_set (y, s, rnd_mode);
99
/* s <- 1 - r/2! + r^2/4! + ... + (-1)^l r^l/(2l)! + ...
101
Returns the index l0 of the last term (-1)^l r^l/(2l)!.
102
The absolute error on s is at most 2 * l0 * 2^(-m).
105
mpfr_cos2_aux (mpfr_ptr s, mpfr_srcptr r)
107
unsigned int l, b = 2;
108
long int prec_t, m = MPFR_PREC(s);
111
MPFR_ASSERTN (MPFR_EXP(r) <= 0);
113
mpfr_set_ui (t, 1, GMP_RNDN);
114
mpfr_set_ui(s, 1, GMP_RNDN);
116
for (l = 1; MPFR_EXP(t) + m >= 0; l++)
118
mpfr_mul (t, t, r, GMP_RNDU); /* err <= (3l-1) ulp */
119
mpfr_div_ui (t, t, (2*l-1)*(2*l), GMP_RNDU); /* err <= 3l ulp */
121
mpfr_add (s, s, t, GMP_RNDD);
123
mpfr_sub (s, s, t, GMP_RNDD);
124
MPFR_ASSERTN (MPFR_EXP(s) == 0); /* check 1/2 <= s < 1 */
125
/* err(s) <= l * 2^(-m) */
126
if (3 * l > (1 << b))
128
/* now 3l <= 2^b, we want 3l*ulp(t) <= 2^(-m)
129
i.e. b+EXP(t)-PREC(t) <= -m */
130
prec_t = m + MPFR_EXP(t) + b;
131
if (prec_t >= MPFR_PREC_MIN)
132
mpfr_round_prec (t, GMP_RNDN, prec_t);