Annotation of OpenXM_contrib/gmp/mpfr/log.c, Revision 1.1
1.1 ! maekawa 1: /* mpfr_log -- natural logarithm of a floating-point number
! 2:
! 3: Copyright (C) 1999 PolKA project, Inria Lorraine and Loria
! 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 Library General Public License as published by
! 9: the Free Software Foundation; either version 2 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 Library General Public
! 15: License for more details.
! 16:
! 17: You should have received a copy of the GNU Library 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 <math.h>
! 24: #include "gmp.h"
! 25: #include "gmp-impl.h"
! 26: #include "mpfr.h"
! 27:
! 28:
! 29: /* The computation of log(a) is done using the formula :
! 30: if we want p bits of the result,
! 31: pi
! 32: log(a) ~ ------------ - m log 2
! 33: 2 AG(1,4/s)
! 34:
! 35: where s = x 2^m > 2^(p/2)
! 36:
! 37: More precisely, if F(x) = int(1/sqrt(1-(1-x^2)*sin(t)^2), t=0..PI/2),
! 38: then for s>=1.26 we have log(s) < F(4/s) < log(s)*(1+4/s^2)
! 39: from which we deduce pi/2/AG(1,4/s)*(1-4/s^2) < log(s) < pi/2/AG(1,4/s)
! 40: so the relative error 4/s^2 is < 4/2^p i.e. 4 ulps.
! 41: */
! 42:
! 43:
! 44: #define MON_INIT(xp, x, p, s) xp = (mp_ptr) TMP_ALLOC(s*BYTES_PER_MP_LIMB); x -> _mp_prec = p; x -> _mp_d = xp; x -> _mp_size = s; x -> _mp_exp = 0;
! 45:
! 46: /* #define DEBUG */
! 47:
! 48: int
! 49: #if __STDC__
! 50: mpfr_log(mpfr_ptr r, mpfr_srcptr a, unsigned char rnd_mode)
! 51: #else
! 52: mpfr_log(r, a, rnd_mode)
! 53: mpfr_ptr r;
! 54: mpfr_srcptr a;
! 55: unsigned char rnd_mode;
! 56: #endif
! 57: {
! 58: int p, m, q, bool, size, cancel;
! 59: mpfr_t cst, rapport, agm, tmp1, tmp2, s, mm;
! 60: mp_limb_t *cstp, *rapportp, *agmp, *tmp1p, *tmp2p, *sp, *mmp;
! 61: double ref;
! 62: TMP_DECL(marker);
! 63:
! 64: /* If a is NaN or a is negative or null, the result is NaN */
! 65: if (FLAG_NAN(a) || (SIGN(a)<=0))
! 66: { SET_NAN(r); return 1; }
! 67:
! 68: /* If a is 1, the result is 0 */
! 69: if (mpfr_cmp_ui_2exp(a,1,0)==0){
! 70: SET_ZERO(r);
! 71: return 0; /* only case where the result is exact */
! 72: }
! 73:
! 74: q=PREC(r);
! 75:
! 76: ref=mpfr_get_d(a)-1.0;
! 77: if (ref<0)
! 78: ref=-ref;
! 79:
! 80: p=q+4;
! 81: /* adjust to entire limb */
! 82: if (p%BITS_PER_MP_LIMB) p += BITS_PER_MP_LIMB - (p%BITS_PER_MP_LIMB);
! 83:
! 84: bool=1;
! 85:
! 86: while (bool==1) {
! 87: #ifdef DEBUG
! 88: printf("a="); mpfr_print_raw(a); putchar('\n');
! 89: printf("p=%d\n", p);
! 90: #endif
! 91: /* Calculus of m (depends on p) */
! 92: m=(int) ceil(((double) p)/2.0) -EXP(a)+1;
! 93:
! 94: /* All the mpfr_t needed have a precision of p */
! 95: TMP_MARK(marker);
! 96: size=(p-1)/BITS_PER_MP_LIMB+1;
! 97: MON_INIT(cstp, cst, p, size);
! 98: MON_INIT(rapportp, rapport, p, size);
! 99: MON_INIT(agmp, agm, p, size);
! 100: MON_INIT(tmp1p, tmp1, p, size);
! 101: MON_INIT(tmp2p, tmp2, p, size);
! 102: MON_INIT(sp, s, p, size);
! 103: MON_INIT(mmp, mm, p, size);
! 104:
! 105: mpfr_set_si(mm,m,GMP_RNDN); /* I have m, supposed exact */
! 106: mpfr_set_si(tmp1,1,GMP_RNDN); /* I have 1, exact */
! 107: mpfr_set_si(tmp2,4,GMP_RNDN); /* I have 4, exact */
! 108: mpfr_mul_2exp(s,a,m,GMP_RNDN); /* I compute s=a*2^m, err <= 1 ulp */
! 109: mpfr_div(rapport,tmp2,s,GMP_RNDN); /* I compute 4/s, err <= 2 ulps */
! 110: mpfr_agm(agm,tmp1,rapport,GMP_RNDN); /* AG(1,4/s), err<=3 ulps */
! 111: mpfr_mul_2exp(tmp1,agm,1,GMP_RNDN); /* 2*AG(1,4/s), still err<=3 ulps */
! 112: mpfr_pi(cst, GMP_RNDN); /* I compute pi, err<=1ulp */
! 113: mpfr_div(tmp2,cst,tmp1,GMP_RNDN); /* pi/2*AG(1,4/s), err<=5ulps */
! 114: mpfr_log2(cst,GMP_RNDN); /* I compute log(2), err<=1ulp */
! 115: mpfr_mul(tmp1,cst,mm,GMP_RNDN); /* I compute m*log(2), err<=2ulps */
! 116: cancel = EXP(tmp2);
! 117: mpfr_sub(cst,tmp2,tmp1,GMP_RNDN); /* log(a), err<=7ulps+cancel */
! 118: cancel -= EXP(cst);
! 119: #ifdef DEBUG
! 120: printf("cancelled bits=%d\n", cancel);
! 121: printf("approx="); mpfr_print_raw(cst); putchar('\n');
! 122: #endif
! 123: if (cancel<0) cancel=0;
! 124:
! 125: /* If we can round the result, we set it and go out of the loop */
! 126:
! 127: /* we have 7 ulps of error from the above roundings,
! 128: 4 ulps from the 4/s^2 second order term,
! 129: plus the cancelled bits */
! 130: if (mpfr_can_round(cst,p-cancel-4,GMP_RNDN,rnd_mode,q)==1) {
! 131: mpfr_set(r,cst,rnd_mode);
! 132: #ifdef DEBUG
! 133: printf("result="); mpfr_print_raw(r); putchar('\n');
! 134: #endif
! 135: bool=0;
! 136: }
! 137: /* else we increase the precision */
! 138: else {
! 139: p += BITS_PER_MP_LIMB+cancel;
! 140: TMP_FREE(marker);
! 141: }
! 142:
! 143: /* We clean */
! 144: TMP_FREE(marker);
! 145:
! 146: }
! 147: return 1; /* result is inexact */
! 148: }
! 149:
! 150:
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