OpenXM_contrib/gmp/mpfr/log.c - diff

Return to log.c CVS log

Up to [local] / OpenXM_contrib / gmp / mpfr

Diff for /OpenXM_contrib/gmp/mpfr/Attic/log.c between version 1.1 and 1.1.1.2

version 1.1, 2000/09/09 14:12:19

version 1.1.1.2, 2003/08/25 16:06:07

Line 1

/* mpfr_log -- natural logarithm of a floating-point number

This file is part of the MPFR Library.

The MPFR Library is free software; you can redistribute it and/or modify

it under the terms of the GNU Library General Public License as published by

it under the terms of the GNU Lesser General Public License as published by

the Free Software Foundation; either version 2 of the License, or (at your

the Free Software Foundation; either version 2.1 of the License, or (at your

option) any later version.

The MPFR Library is distributed in the hope that it will be useful, but

WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY

or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public

or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public

License for more details.

You should have received a copy of the GNU Library General Public License

You should have received a copy of the GNU Lesser General Public License

along with the MPFR Library; see the file COPYING.LIB. If not, write to

the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,

MA 02111-1307, USA. */

#include <stdio.h>

#include <math.h>

#include "gmp.h"

#include "gmp-impl.h"

#include "mpfr.h"

#include "mpfr-impl.h"

/* The computation of log(a) is done using the formula :

if we want p bits of the result,

Line 40 MA 02111-1307, USA. */

Line 39 MA 02111-1307, USA. */

so the relative error 4/s^2 is < 4/2^p i.e. 4 ulps.

#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;

/* #define DEBUG */

int

#if __STDC__

mpfr_log (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t rnd_mode)

mpfr_log(mpfr_ptr r, mpfr_srcptr a, unsigned char rnd_mode)

#else

mpfr_log(r, a, rnd_mode)

mpfr_ptr r;

mpfr_srcptr a;

unsigned char rnd_mode;

#endif

{

int p, m, q, bool, size, cancel;

int m, bool, size, cancel, inexact = 0;

mp_prec_t p, q;

mpfr_t cst, rapport, agm, tmp1, tmp2, s, mm;

mp_limb_t *cstp, *rapportp, *agmp, *tmp1p, *tmp2p, *sp, *mmp;

double ref;

TMP_DECL(marker);

/* If a is NaN or a is negative or null, the result is NaN */

/* If a is NaN, the result is NaN */

if (FLAG_NAN(a) || (SIGN(a)<=0))

if (MPFR_IS_NAN(a))

{ SET_NAN(r); return 1; }

{

MPFR_SET_NAN(r);

MPFR_RET_NAN;

}

MPFR_CLEAR_NAN(r);

/* check for infinity before zero */

if (MPFR_IS_INF(a))

{

if (MPFR_SIGN(a) < 0) /* log(-Inf) = NaN */

{

MPFR_SET_NAN(r);

MPFR_RET_NAN;

}

else /* log(+Inf) = +Inf */

{

MPFR_SET_INF(r);

MPFR_SET_POS(r);

MPFR_RET(0);

}

/* Now we can clear the flags without damage even if r == a */

MPFR_CLEAR_INF(r);

if (MPFR_IS_ZERO(a))

{

MPFR_SET_INF(r);

MPFR_SET_NEG(r);

MPFR_RET(0); /* log(0) is an exact -infinity */

}

/* If a is negative, the result is NaN */

if (MPFR_SIGN(a) < 0)

{

MPFR_SET_NAN(r);

MPFR_RET_NAN;

}

/* If a is 1, the result is 0 */

if (mpfr_cmp_ui_2exp(a,1,0)==0){

if (mpfr_cmp_ui (a, 1) == 0)

SET_ZERO(r);

{

return 0; /* only case where the result is exact */

MPFR_SET_ZERO(r);

}

MPFR_SET_POS(r);

MPFR_RET(0); /* only "normal" case where the result is exact */

}

q=PREC(r);

q=MPFR_PREC(r);

ref=mpfr_get_d(a)-1.0;

ref = mpfr_get_d1 (a) - 1.0;

if (ref<0)

ref=-ref;

p=q+4;

/* use initial precision about q+lg(q)+5 */

p=q+5; m=q; while (m) { p++; m >>= 1; }

/* adjust to entire limb */

if (p%BITS_PER_MP_LIMB) p += BITS_PER_MP_LIMB - (p%BITS_PER_MP_LIMB);

Line 85 mpfr_log(r, a, rnd_mode)

Line 117 mpfr_log(r, a, rnd_mode)

while (bool==1) {

#ifdef DEBUG

printf("a="); mpfr_print_raw(a); putchar('\n');

printf("a="); mpfr_print_binary(a); putchar('\n');

printf("p=%d\n", p);

#endif

/* Calculus of m (depends on p) */

m=(int) ceil(((double) p)/2.0) -EXP(a)+1;

m = (p + 1) / 2 - MPFR_EXP(a) + 1;

/* All the mpfr_t needed have a precision of p */

TMP_MARK(marker);

size=(p-1)/BITS_PER_MP_LIMB+1;

MON_INIT(cstp, cst, p, size);

MPFR_INIT(cstp, cst, p, size);

MON_INIT(rapportp, rapport, p, size);

MPFR_INIT(rapportp, rapport, p, size);

MON_INIT(agmp, agm, p, size);

MPFR_INIT(agmp, agm, p, size);

MON_INIT(tmp1p, tmp1, p, size);

MPFR_INIT(tmp1p, tmp1, p, size);

MON_INIT(tmp2p, tmp2, p, size);

MPFR_INIT(tmp2p, tmp2, p, size);

MON_INIT(sp, s, p, size);

MPFR_INIT(sp, s, p, size);

MON_INIT(mmp, mm, p, size);

MPFR_INIT(mmp, mm, p, size);

mpfr_set_si(mm,m,GMP_RNDN); /* I have m, supposed exact */

mpfr_set_si (mm, m, GMP_RNDN); /* I have m, supposed exact */

mpfr_set_si(tmp1,1,GMP_RNDN); /* I have 1, exact */

mpfr_set_si (tmp1, 1, GMP_RNDN); /* I have 1, exact */

mpfr_set_si(tmp2,4,GMP_RNDN); /* I have 4, exact */

mpfr_set_si (tmp2, 4, GMP_RNDN); /* I have 4, exact */

mpfr_mul_2exp(s,a,m,GMP_RNDN); /* I compute s=a*2^m, err <= 1 ulp */

mpfr_mul_2si (s, a, m, GMP_RNDN); /* I compute s=a*2^m, err <= 1 ulp */

mpfr_div(rapport,tmp2,s,GMP_RNDN); /* I compute 4/s, err <= 2 ulps */

mpfr_div (rapport, tmp2, s, GMP_RNDN);/* I compute 4/s, err <= 2 ulps */

mpfr_agm(agm,tmp1,rapport,GMP_RNDN); /* AG(1,4/s), err<=3 ulps */

mpfr_agm (agm, tmp1, rapport, GMP_RNDN); /* AG(1,4/s), err<=3 ulps */

mpfr_mul_2exp(tmp1,agm,1,GMP_RNDN); /* 2*AG(1,4/s), still err<=3 ulps */

mpfr_mul_2ui (tmp1, agm, 1, GMP_RNDN); /* 2*AG(1,4/s), still err<=3 ulps */

mpfr_pi(cst, GMP_RNDN); /* I compute pi, err<=1ulp */

mpfr_const_pi (cst, GMP_RNDN); /* compute pi, err<=1ulp */

mpfr_div(tmp2,cst,tmp1,GMP_RNDN); /* pi/2*AG(1,4/s), err<=5ulps */

mpfr_div (tmp2, cst, tmp1, GMP_RNDN); /* pi/2*AG(1,4/s), err<=5ulps */

mpfr_log2(cst,GMP_RNDN); /* I compute log(2), err<=1ulp */

mpfr_const_log2 (cst, GMP_RNDN); /* compute log(2), err<=1ulp */

mpfr_mul(tmp1,cst,mm,GMP_RNDN); /* I compute m*log(2), err<=2ulps */

cancel = EXP(tmp2);

cancel = MPFR_EXP(tmp2);

mpfr_sub(cst,tmp2,tmp1,GMP_RNDN); /* log(a), err<=7ulps+cancel */

cancel -= EXP(cst);

cancel -= MPFR_EXP(cst);

#ifdef DEBUG

printf("cancelled bits=%d\n", cancel);

printf("canceled bits=%d\n", cancel);

printf("approx="); mpfr_print_raw(cst); putchar('\n');

printf("approx="); mpfr_print_binary(cst); putchar('\n');

#endif

if (cancel<0) cancel=0;

Line 126 mpfr_log(r, a, rnd_mode)

Line 158 mpfr_log(r, a, rnd_mode)

/* we have 7 ulps of error from the above roundings,

4 ulps from the 4/s^2 second order term,

plus the cancelled bits */

plus the canceled bits */

if (mpfr_can_round(cst,p-cancel-4,GMP_RNDN,rnd_mode,q)==1) {

if (mpfr_can_round (cst, p - cancel - 4, GMP_RNDN, rnd_mode, q) == 1) {

mpfr_set(r,cst,rnd_mode);

inexact = mpfr_set (r, cst, rnd_mode);

#ifdef DEBUG

printf("result="); mpfr_print_raw(r); putchar('\n');

printf("result="); mpfr_print_binary(r); putchar('\n');

#endif

bool=0;

}

/* else we increase the precision */

else {

p += BITS_PER_MP_LIMB+cancel;

p += BITS_PER_MP_LIMB + cancel;

TMP_FREE(marker);

}

/* We clean */

TMP_FREE(marker);

}

return 1; /* result is inexact */

return inexact; /* result is inexact */

}

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>