Annotation of OpenXM_contrib/gmp/mpfr/log.c, Revision 1.1.1.2
1.1 maekawa 1: /* mpfr_log -- natural logarithm of a floating-point number
2:
1.1.1.2 ! ohara 3: Copyright 1999, 2000, 2001, 2002 Free Software Foundation.
1.1 maekawa 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
1.1.1.2 ! ohara 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
1.1 maekawa 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
1.1.1.2 ! ohara 14: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
1.1 maekawa 15: License for more details.
16:
1.1.1.2 ! ohara 17: You should have received a copy of the GNU Lesser General Public License
1.1 maekawa 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 "gmp.h"
24: #include "gmp-impl.h"
25: #include "mpfr.h"
1.1.1.2 ! ohara 26: #include "mpfr-impl.h"
1.1 maekawa 27:
28: /* The computation of log(a) is done using the formula :
29: if we want p bits of the result,
30: pi
31: log(a) ~ ------------ - m log 2
32: 2 AG(1,4/s)
33:
34: where s = x 2^m > 2^(p/2)
35:
36: More precisely, if F(x) = int(1/sqrt(1-(1-x^2)*sin(t)^2), t=0..PI/2),
37: then for s>=1.26 we have log(s) < F(4/s) < log(s)*(1+4/s^2)
38: from which we deduce pi/2/AG(1,4/s)*(1-4/s^2) < log(s) < pi/2/AG(1,4/s)
39: so the relative error 4/s^2 is < 4/2^p i.e. 4 ulps.
40: */
41:
42: /* #define DEBUG */
43:
44: int
1.1.1.2 ! ohara 45: mpfr_log (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t rnd_mode)
1.1 maekawa 46: {
1.1.1.2 ! ohara 47: int m, bool, size, cancel, inexact = 0;
! 48: mp_prec_t p, q;
1.1 maekawa 49: mpfr_t cst, rapport, agm, tmp1, tmp2, s, mm;
50: mp_limb_t *cstp, *rapportp, *agmp, *tmp1p, *tmp2p, *sp, *mmp;
51: double ref;
52: TMP_DECL(marker);
53:
1.1.1.2 ! ohara 54: /* If a is NaN, the result is NaN */
! 55: if (MPFR_IS_NAN(a))
! 56: {
! 57: MPFR_SET_NAN(r);
! 58: MPFR_RET_NAN;
! 59: }
! 60:
! 61: MPFR_CLEAR_NAN(r);
! 62:
! 63: /* check for infinity before zero */
! 64: if (MPFR_IS_INF(a))
! 65: {
! 66: if (MPFR_SIGN(a) < 0) /* log(-Inf) = NaN */
! 67: {
! 68: MPFR_SET_NAN(r);
! 69: MPFR_RET_NAN;
! 70: }
! 71: else /* log(+Inf) = +Inf */
! 72: {
! 73: MPFR_SET_INF(r);
! 74: MPFR_SET_POS(r);
! 75: MPFR_RET(0);
! 76: }
! 77: }
! 78:
! 79: /* Now we can clear the flags without damage even if r == a */
! 80: MPFR_CLEAR_INF(r);
! 81:
! 82: if (MPFR_IS_ZERO(a))
! 83: {
! 84: MPFR_SET_INF(r);
! 85: MPFR_SET_NEG(r);
! 86: MPFR_RET(0); /* log(0) is an exact -infinity */
! 87: }
! 88:
! 89: /* If a is negative, the result is NaN */
! 90: if (MPFR_SIGN(a) < 0)
! 91: {
! 92: MPFR_SET_NAN(r);
! 93: MPFR_RET_NAN;
! 94: }
1.1 maekawa 95:
96: /* If a is 1, the result is 0 */
1.1.1.2 ! ohara 97: if (mpfr_cmp_ui (a, 1) == 0)
! 98: {
! 99: MPFR_SET_ZERO(r);
! 100: MPFR_SET_POS(r);
! 101: MPFR_RET(0); /* only "normal" case where the result is exact */
! 102: }
1.1 maekawa 103:
1.1.1.2 ! ohara 104: q=MPFR_PREC(r);
1.1 maekawa 105:
1.1.1.2 ! ohara 106: ref = mpfr_get_d1 (a) - 1.0;
1.1 maekawa 107: if (ref<0)
108: ref=-ref;
109:
1.1.1.2 ! ohara 110: /* use initial precision about q+lg(q)+5 */
! 111: p=q+5; m=q; while (m) { p++; m >>= 1; }
! 112:
1.1 maekawa 113: /* adjust to entire limb */
114: if (p%BITS_PER_MP_LIMB) p += BITS_PER_MP_LIMB - (p%BITS_PER_MP_LIMB);
115:
116: bool=1;
117:
118: while (bool==1) {
119: #ifdef DEBUG
1.1.1.2 ! ohara 120: printf("a="); mpfr_print_binary(a); putchar('\n');
1.1 maekawa 121: printf("p=%d\n", p);
122: #endif
123: /* Calculus of m (depends on p) */
1.1.1.2 ! ohara 124: m = (p + 1) / 2 - MPFR_EXP(a) + 1;
1.1 maekawa 125:
126: /* All the mpfr_t needed have a precision of p */
127: TMP_MARK(marker);
128: size=(p-1)/BITS_PER_MP_LIMB+1;
1.1.1.2 ! ohara 129: MPFR_INIT(cstp, cst, p, size);
! 130: MPFR_INIT(rapportp, rapport, p, size);
! 131: MPFR_INIT(agmp, agm, p, size);
! 132: MPFR_INIT(tmp1p, tmp1, p, size);
! 133: MPFR_INIT(tmp2p, tmp2, p, size);
! 134: MPFR_INIT(sp, s, p, size);
! 135: MPFR_INIT(mmp, mm, p, size);
! 136:
! 137: mpfr_set_si (mm, m, GMP_RNDN); /* I have m, supposed exact */
! 138: mpfr_set_si (tmp1, 1, GMP_RNDN); /* I have 1, exact */
! 139: mpfr_set_si (tmp2, 4, GMP_RNDN); /* I have 4, exact */
! 140: mpfr_mul_2si (s, a, m, GMP_RNDN); /* I compute s=a*2^m, err <= 1 ulp */
! 141: mpfr_div (rapport, tmp2, s, GMP_RNDN);/* I compute 4/s, err <= 2 ulps */
! 142: mpfr_agm (agm, tmp1, rapport, GMP_RNDN); /* AG(1,4/s), err<=3 ulps */
! 143: mpfr_mul_2ui (tmp1, agm, 1, GMP_RNDN); /* 2*AG(1,4/s), still err<=3 ulps */
! 144: mpfr_const_pi (cst, GMP_RNDN); /* compute pi, err<=1ulp */
! 145: mpfr_div (tmp2, cst, tmp1, GMP_RNDN); /* pi/2*AG(1,4/s), err<=5ulps */
! 146: mpfr_const_log2 (cst, GMP_RNDN); /* compute log(2), err<=1ulp */
1.1 maekawa 147: mpfr_mul(tmp1,cst,mm,GMP_RNDN); /* I compute m*log(2), err<=2ulps */
1.1.1.2 ! ohara 148: cancel = MPFR_EXP(tmp2);
1.1 maekawa 149: mpfr_sub(cst,tmp2,tmp1,GMP_RNDN); /* log(a), err<=7ulps+cancel */
1.1.1.2 ! ohara 150: cancel -= MPFR_EXP(cst);
1.1 maekawa 151: #ifdef DEBUG
1.1.1.2 ! ohara 152: printf("canceled bits=%d\n", cancel);
! 153: printf("approx="); mpfr_print_binary(cst); putchar('\n');
1.1 maekawa 154: #endif
155: if (cancel<0) cancel=0;
156:
157: /* If we can round the result, we set it and go out of the loop */
158:
159: /* we have 7 ulps of error from the above roundings,
160: 4 ulps from the 4/s^2 second order term,
1.1.1.2 ! ohara 161: plus the canceled bits */
! 162: if (mpfr_can_round (cst, p - cancel - 4, GMP_RNDN, rnd_mode, q) == 1) {
! 163: inexact = mpfr_set (r, cst, rnd_mode);
1.1 maekawa 164: #ifdef DEBUG
1.1.1.2 ! ohara 165: printf("result="); mpfr_print_binary(r); putchar('\n');
1.1 maekawa 166: #endif
167: bool=0;
168: }
169: /* else we increase the precision */
170: else {
1.1.1.2 ! ohara 171: p += BITS_PER_MP_LIMB + cancel;
1.1 maekawa 172: }
173:
174: /* We clean */
175: TMP_FREE(marker);
176:
177: }
1.1.1.2 ! ohara 178: return inexact; /* result is inexact */
1.1 maekawa 179: }
180:
181:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to log.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpfr