Annotation of OpenXM_contrib/gmp/mpfr/atan.c, Revision 1.1
1.1 ! ohara 1: /* mpfr_atan -- arc-tangent of a floating-point number
! 2:
! 3: Copyright 2001 Free Software Foundation.
! 4:
! 5: This file is part of the MPFR Library, and was contributed by Mathieu Dutour.
! 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., 59 Temple Place - Suite 330, Boston,
! 20: MA 02111-1307, USA. */
! 21:
! 22: #include <stdio.h>
! 23: #include <stdlib.h>
! 24: #include "gmp.h"
! 25: #include "gmp-impl.h"
! 26: #include "mpfr.h"
! 27: #include "mpfr-impl.h"
! 28:
! 29: #define CST 2.27 /* CST=1+ln(2.4)/ln(2) */
! 30: #define CST2 1.45 /* CST2=1/ln(2) */
! 31:
! 32: static int mpfr_atan_aux _PROTO((mpfr_ptr, mpz_srcptr, int, int));
! 33:
! 34: #undef B
! 35: #define A
! 36: #define A1 1
! 37: #define A2 2
! 38: #define C
! 39: #define C1 3
! 40: #define C2 2
! 41: #define NO_FACTORIAL
! 42: #define GENERIC mpfr_atan_aux
! 43: #include "generic.c"
! 44: #undef C
! 45: #undef C1
! 46: #undef C2
! 47: #undef A
! 48: #undef A1
! 49: #undef A2
! 50: #undef NO_FACTORIAL
! 51: #undef GENERIC
! 52:
! 53: int
! 54: mpfr_atan (mpfr_ptr arctangent, mpfr_srcptr x, mp_rnd_t rnd_mode)
! 55: {
! 56: mpfr_t Pisur2;
! 57: mpfr_t xp;
! 58: mpfr_t arctgt;
! 59:
! 60: int comparaison, signe, suplement;
! 61:
! 62: mpfr_t t_arctan;
! 63: int i;
! 64: mpz_t ukz;
! 65: mpfr_t ukf;
! 66: mpfr_t sk,Ak;
! 67: mpz_t square;
! 68: mpfr_t tmp_arctan;
! 69: mpfr_t tmp, tmp2;
! 70: #ifdef DEBUG
! 71: mpfr_t tst;
! 72: #endif
! 73: int twopoweri;
! 74: int Prec;
! 75: int prec_x;
! 76: int prec_arctan;
! 77: int good = 0;
! 78: int realprec;
! 79: int estimated_delta;
! 80: /* calculation of the floor */
! 81: mp_exp_t exptol;
! 82:
! 83: int N0;
! 84: int logn;
! 85:
! 86: /* Trivial cases */
! 87: if (MPFR_IS_NAN(x))
! 88: {
! 89: MPFR_SET_NAN(arctangent);
! 90: MPFR_RET_NAN;
! 91: }
! 92:
! 93: if (MPFR_IS_INF(x))
! 94: {
! 95: MPFR_CLEAR_FLAGS(arctangent);
! 96: if (MPFR_SIGN(x) > 0) /* arctan(+inf) = Pi/2 */
! 97: mpfr_const_pi (arctangent, rnd_mode);
! 98: else /* arctan(-inf) = -Pi/2 */
! 99: {
! 100: if (rnd_mode == GMP_RNDU)
! 101: rnd_mode = GMP_RNDD;
! 102: else if (rnd_mode == GMP_RNDD)
! 103: rnd_mode = GMP_RNDU;
! 104: mpfr_const_pi (arctangent, rnd_mode);
! 105: }
! 106: MPFR_EXP(arctangent)--;
! 107: return 1; /* inexact */
! 108: }
! 109:
! 110: MPFR_CLEAR_FLAGS(arctangent);
! 111: if (MPFR_IS_ZERO(x))
! 112: {
! 113: mpfr_set_ui (arctangent, 0, GMP_RNDN);
! 114: return 0; /* exact result */
! 115: }
! 116:
! 117: signe = MPFR_SIGN(x);
! 118: prec_arctan = MPFR_PREC(arctangent);
! 119:
! 120: /* Set x_p=|x| */
! 121: mpfr_init2(xp, MPFR_PREC(x));
! 122: mpfr_set(xp, x, rnd_mode);
! 123: if (signe == -1)
! 124: MPFR_CHANGE_SIGN(xp);
! 125:
! 126: /* Other simple case arctang(-+1)=-+pi/4 */
! 127: comparaison=mpfr_cmp_ui(xp, 1);
! 128: if (comparaison == 0) {
! 129: mpfr_init2(Pisur2, prec_arctan);
! 130: mpfr_const_pi(Pisur2, rnd_mode);
! 131: mpfr_div_2ui(arctangent, Pisur2, 2, rnd_mode);
! 132: if (signe == -1)
! 133: MPFR_CHANGE_SIGN(arctangent);
! 134: mpfr_clear(Pisur2);
! 135: mpfr_clear(xp);
! 136: return 0; /* Result correct */
! 137: }
! 138: if (comparaison > 0)
! 139: suplement = 2;
! 140: else
! 141: suplement = 2-MPFR_EXP(xp);
! 142:
! 143: prec_x = _mpfr_ceil_log2 ((double) MPFR_PREC(x) / BITS_PER_MP_LIMB);
! 144: logn = _mpfr_ceil_log2 ((double) prec_x);
! 145: if (logn < 2) logn = 2;
! 146: realprec = prec_arctan + _mpfr_ceil_log2((double) prec_arctan) + 4;
! 147: mpz_init(ukz);
! 148: mpz_init(square);
! 149:
! 150:
! 151: while (!good){
! 152: N0 = _mpfr_ceil_log2((double) realprec + suplement + CST);
! 153: estimated_delta = 1 + suplement + _mpfr_ceil_log2((double) (3*N0-2));
! 154: Prec = realprec+estimated_delta;
! 155:
! 156: /* Initialisation */
! 157: mpfr_init2(sk,Prec);
! 158: mpfr_init2(ukf, Prec);
! 159: mpfr_init2(t_arctan, Prec);
! 160: mpfr_init2(tmp_arctan, Prec);
! 161: mpfr_init2(tmp, Prec);
! 162: mpfr_init2(tmp2, Prec);
! 163: mpfr_init2(Ak, Prec);
! 164: mpfr_init2(arctgt, Prec);
! 165:
! 166: #ifdef DEBUG
! 167: /* Tests */
! 168: mpfr_init2(tst, realprec);
! 169: #endif
! 170:
! 171: if (comparaison > 0)
! 172: {
! 173: mpfr_init2(Pisur2, Prec);
! 174: mpfr_const_pi(Pisur2, GMP_RNDN);
! 175: mpfr_div_2ui(Pisur2, Pisur2, 1, GMP_RNDN);
! 176: mpfr_ui_div(sk, 1, xp, GMP_RNDN);
! 177: }
! 178: else
! 179: mpfr_set(sk, xp, GMP_RNDN);
! 180:
! 181: /* Assignation */
! 182: mpfr_set_ui (tmp_arctan, 0, GMP_RNDN);
! 183: twopoweri = 1;
! 184: for(i = 0; i <= N0; i++){
! 185: mpfr_mul_2ui(tmp, sk, twopoweri, GMP_RNDN);
! 186: /* Calculation of trunc(tmp) --> mpz */
! 187: mpfr_trunc (ukf, tmp);
! 188: exptol = mpfr_get_z_exp (ukz, ukf);
! 189: if (exptol>0)
! 190: mpz_mul_2exp (ukz, ukz, exptol);
! 191: else
! 192: mpz_tdiv_q_2exp (ukz, ukz, (unsigned long int) (-exptol));
! 193:
! 194: /* Calculation of arctan(Ak) */
! 195: mpz_mul(square, ukz, ukz);
! 196: mpz_neg(square, square);
! 197: mpfr_atan_aux(t_arctan, square, 2*twopoweri, N0 - i);
! 198: mpfr_set_z(Ak, ukz, GMP_RNDN);
! 199: mpfr_div_2ui(Ak, Ak, twopoweri, GMP_RNDN);
! 200: mpfr_mul(t_arctan, t_arctan, Ak, GMP_RNDN);
! 201:
! 202: /* Addition and iteration */
! 203: mpfr_add(tmp_arctan, tmp_arctan, t_arctan, GMP_RNDN);
! 204: if (i<N0)
! 205: {
! 206: mpfr_sub(tmp, sk, Ak, GMP_RNDN);
! 207: mpfr_mul(tmp2, sk, Ak, GMP_RNDN);
! 208: mpfr_add_ui(tmp2, tmp2, 1, GMP_RNDN);
! 209: mpfr_div(sk, tmp, tmp2, GMP_RNDN);
! 210: twopoweri <<= 1;
! 211: }
! 212: }
! 213:
! 214: if (comparaison > 0)
! 215: {
! 216: mpfr_sub(arctgt, Pisur2, tmp_arctan, GMP_RNDN);
! 217: if (signe == -1)
! 218: MPFR_CHANGE_SIGN(arctgt);
! 219: }
! 220: else
! 221: {
! 222: mpfr_set(arctgt, tmp_arctan, GMP_RNDN);
! 223: if (signe == -1)
! 224: MPFR_CHANGE_SIGN(arctgt);
! 225: }
! 226:
! 227: #ifdef DEBUG
! 228: mpfr_set(tst, arctgt, rnd_mode);
! 229: #endif
! 230:
! 231: if (mpfr_can_round(arctgt, realprec, GMP_RNDN, rnd_mode, MPFR_PREC(arctangent)))
! 232: {
! 233: mpfr_set(arctangent, arctgt, rnd_mode);
! 234: good = 1;
! 235: realprec += 1;
! 236: }
! 237: else
! 238: {
! 239: realprec += _mpfr_ceil_log2 ((double) realprec);
! 240: }
! 241:
! 242: mpfr_clear(sk);
! 243: mpfr_clear(ukf);
! 244: mpfr_clear(t_arctan);
! 245: mpfr_clear(tmp_arctan);
! 246: mpfr_clear(tmp);
! 247: mpfr_clear(tmp2);
! 248: mpfr_clear(Ak);
! 249: mpfr_clear(arctgt);
! 250:
! 251: #ifdef DEBUG
! 252: mpfr_clear(tst);
! 253: #endif
! 254: if (comparaison > 0)
! 255: mpfr_clear(Pisur2);
! 256: }
! 257:
! 258: mpfr_clear(xp);
! 259: mpz_clear(ukz);
! 260: mpz_clear(square);
! 261:
! 262: return 1; /* inexact result */
! 263: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to atan.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpfr