diff options
Diffstat (limited to 'mpfr/exp3.c')
| -rw-r--r-- | mpfr/exp3.c | 94 |
1 files changed, 25 insertions, 69 deletions
diff --git a/mpfr/exp3.c b/mpfr/exp3.c index 175ef143d..8441850b9 100644 --- a/mpfr/exp3.c +++ b/mpfr/exp3.c @@ -1,20 +1,20 @@ /* mpfr_exp -- exponential of a floating-point number -Copyright (C) 1999 Free Software Foundation. +Copyright (C) 1999, 2001 Free Software Foundation, Inc. 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 -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +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. */ @@ -25,21 +25,11 @@ MA 02111-1307, USA. */ #include "mpfr.h" #include "mpfr-impl.h" -/* #define DEBUG */ - -int mpfr_exp_rational _PROTO ((mpfr_ptr, mpz_srcptr, int, int)); +static int mpfr_exp_rational _PROTO ((mpfr_ptr, mpz_srcptr, int, int)); int mpfr_exp3 _PROTO ((mpfr_ptr, mpfr_srcptr, mp_rnd_t)); -int -#if __STDC__ +static int mpfr_exp_rational (mpfr_ptr y, mpz_srcptr p, int r, int m) -#else -mpfr_exp_rational (y, p, r, m) - mpfr_ptr y; - mpz_srcptr p; - int r; - int m; -#endif { int n,i,k,j,l; mpz_t* P,*S; @@ -129,14 +119,7 @@ mpfr_exp_rational (y, p, r, m) #define shift (BITS_PER_MP_LIMB/2) int -#if __STDC__ mpfr_exp3 (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) -#else -mpfr_exp3 (y, x, rnd_mode) - mpfr_ptr y; - mpfr_srcptr x; - mp_rnd_t rnd_mode; -#endif { mpfr_t t; mpfr_t x_copy; @@ -152,27 +135,10 @@ mpfr_exp3 (y, x, rnd_mode) int good = 0; int realprec = 0; int iter; - int logn; - - /* commencons par 0 */ - if (MPFR_IS_NAN(x)) { MPFR_SET_NAN(y); return 1; } - - if (MPFR_IS_INF(x)) - { - if (MPFR_SIGN(x) > 0) - { MPFR_SET_INF(y); if (MPFR_SIGN(y) == -1) { MPFR_CHANGE_SIGN(y); } } - else - { MPFR_SET_ZERO(y); if (MPFR_SIGN(y) == -1) { MPFR_CHANGE_SIGN(y); } } - /* TODO: conflits entre infinis et zeros ? */ - } - - if (!MPFR_NOTZERO(x)) { - mpfr_set_ui(y, 1, GMP_RNDN); - return 0; - } + int logn, inexact = 0; - /* Decomposer x */ - /* on commence par ecrire x = 1.xxxxxxxxxxxxx + /* decompose x */ + /* we first write x = 1.xxxxxxxxxxxxx ----- k bits -- */ prec_x = _mpfr_ceil_log2 ((double) (MPFR_PREC(x)) / BITS_PER_MP_LIMB); if (prec_x < 0) prec_x = 0; @@ -181,7 +147,7 @@ mpfr_exp3 (y, x, rnd_mode) ttt = MPFR_EXP(x); mpfr_init2(x_copy,MPFR_PREC(x)); mpfr_set(x_copy,x,GMP_RNDD); - /* on fait le shift pour que le nombre soit inferieur a 1 */ + /* we shift to get a number less than 1 */ if (ttt > 0) { shift_x = ttt; @@ -195,52 +161,42 @@ mpfr_exp3 (y, x, rnd_mode) k = _mpfr_ceil_log2 ((double) Prec / BITS_PER_MP_LIMB); /* now we have to extract */ - mpfr_init2(t, Prec); - mpfr_init2(tmp, Prec); + mpfr_init2 (t, Prec); + mpfr_init2 (tmp, Prec); mpfr_set_ui(tmp,1,GMP_RNDN); twopoweri = BITS_PER_MP_LIMB; if (k <= prec_x) iter = k; else iter= prec_x; for(i = 0; i <= iter; i++){ mpfr_extract (uk, x_copy, i); -#ifdef DEBUG - mpz_out_str(stderr,2, uk); - fprintf(stderr, "---\n"); - fprintf(stderr, "---%d\n", twopoweri - ttt); -#endif if (i) mpfr_exp_rational (t, uk, twopoweri - ttt, k - i + 1); else { - /* cas particulier : on est oblige de faire les calculs avec x/2^. - puis elever au carre (plus rapide) */ + /* particular case: we have to compute with x/2^., then + do squarings (this is faster) */ mpfr_exp_rational (t, uk, shift + twopoweri - ttt, k+1); for (loop= 0 ; loop < shift; loop++) mpfr_mul(t,t,t,GMP_RNDD); } mpfr_mul(tmp,tmp,t,GMP_RNDD); -#ifdef DEBUG - fprintf(stderr, "fin\n"); - mpfr_out_str(stderr, 2, MPFR_PREC(y), t, GMP_RNDD); - fprintf(stderr, "\n ii --- ii \n"); -#endif twopoweri <<= 1; } + mpfr_clear (t); for (loop= 0 ; loop < shift_x; loop++) - mpfr_mul(tmp,tmp,tmp,GMP_RNDD); - mpfr_clear(t); - if (mpfr_can_round(tmp, realprec, GMP_RNDD, rnd_mode, MPFR_PREC(y))){ - mpfr_set(y,tmp,rnd_mode); - mpfr_clear(tmp); - good = 1; - } else { - mpfr_clear(tmp); + mpfr_mul (tmp, tmp, tmp, GMP_RNDD); + if (mpfr_can_round (tmp, realprec, GMP_RNDD, rnd_mode, MPFR_PREC(y))) + { + inexact = mpfr_set (y, tmp, rnd_mode); + good = 1; + } + else realprec += 3*logn; - } + mpfr_clear (tmp); } mpz_clear (uk); mpfr_clear(x_copy); - return 0; + return inexact; } |