diff options
author | jeandel <jeandel@280ebfd0-de03-0410-8827-d642c229c3f4> | 2000-07-13 08:54:06 +0000 |
---|---|---|
committer | jeandel <jeandel@280ebfd0-de03-0410-8827-d642c229c3f4> | 2000-07-13 08:54:06 +0000 |
commit | 2887b99b09bbb669cc58f8314c57e182f6688976 (patch) | |
tree | b30e319f8f6ee8cabb8195c02d1bdfa172f39a58 /log2.c | |
parent | 7145862827ac6f63c6bdcc58328dc0a0c152db03 (diff) | |
download | mpfr-2887b99b09bbb669cc58f8314c57e182f6688976.tar.gz |
New algorithm, new wrapper
git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@680 280ebfd0-de03-0410-8827-d642c229c3f4
Diffstat (limited to 'log2.c')
-rw-r--r-- | log2.c | 88 |
1 files changed, 87 insertions, 1 deletions
@@ -29,6 +29,88 @@ MA 02111-1307, USA. */ mpfr_t __mpfr_const_log2; /* stored value of log(2) with rnd_mode=GMP_RNDZ */ int __mpfr_const_log2_prec=0; /* precision of stored value */ + +#define A +#define A1 1 +#define A2 1 +#undef B +#define C +#define C1 2 +#define C2 1 +#define NO_FACTORIAL +#undef R_IS_RATIONAL +#define GENERIC mpfr_aux_log2 +#include "generic.c" +#undef A +#undef A1 +#undef A2 +#undef NO_FACTORIAL +#undef GENERIC +#undef C +#undef C1 +#undef C2 + +int +#if __STDC__ +mpfr_const_aux_log2(mpfr_ptr mylog, mp_rnd_t rnd_mode) +#else +mpfr_const_aux_log2(mylog, rnd_mode) mpfr_ptr +mylog; mp_rnd_t rnd_mode; +#endif +{ + int prec; + mpfr_t tmp1, tmp2, result,tmp3; + mpz_t cst; + int good = 0; + int logn; + int prec_i_want = PREC(mylog); + int prec_x; + mpz_init(cst); + logn = (int) ceil(log + ((double) PREC(mylog)) + /log(2.0)); + prec_x = prec_i_want + logn; + while (!good){ + prec = (int) ceil(log + ((double) (prec_x)) + /log(2.0)); + mpfr_init2(tmp1, prec_x); + mpfr_init2(result, prec_x); + mpfr_init2(tmp2, prec_x); + mpfr_init2(tmp3, prec_x); + mpz_set_ui(cst, 1); + mpfr_aux_log2(tmp1, cst, 4, prec-2); + mpfr_div_2exp(tmp1, tmp1, 4,GMP_RNDD); + mpfr_mul_ui(tmp1, tmp1, 15,GMP_RNDD); + + mpz_set_ui(cst, 3); + mpfr_aux_log2(tmp2, cst, 7, prec-2); + mpfr_div_2exp(tmp2, tmp2, 7,GMP_RNDD); + mpfr_mul_ui(tmp2, tmp2, 5*3,GMP_RNDD); + mpfr_sub(result, tmp1, tmp2, GMP_RNDD); + + mpz_set_ui(cst, 13); + mpfr_aux_log2(tmp3, cst, 8, prec-2); + mpfr_div_2exp(tmp3, tmp3, 8,GMP_RNDD); + mpfr_mul_ui(tmp3, tmp3, 3*13,GMP_RNDD); + mpfr_sub(result, result, tmp3, GMP_RNDD); + + mpfr_clear(tmp1); + mpfr_clear(tmp2); + if (mpfr_can_round(result, prec_x, GMP_RNDD, rnd_mode, prec_i_want)){ + mpfr_set(mylog, result, rnd_mode); + mpfr_clear(result); + good = 1; + } else + { + mpfr_clear(result); + prec_x += logn; + } + } + mpz_clear(cst); + return 0; +} + /* set x to log(2) rounded to precision PREC(x) with direction rnd_mode use formula log(2) = sum(1/k/2^k, k=1..infinity) @@ -63,6 +145,7 @@ mpfr_const_log2(x, rnd_mode) mpfr_ptr x; mp_rnd_t rnd_mode; } } + if (precx < 30000){ /* need to recompute */ N=2; do { @@ -92,7 +175,10 @@ mpfr_const_log2(x, rnd_mode) mpfr_ptr x; mp_rnd_t rnd_mode; #endif mpfr_set_z(x, s, rnd_mode); EXP(x) -= N; - + } else + { + mpfr_const_aux_log2(x, rnd_mode); + } /* stored computed value */ if (__mpfr_const_log2_prec==0) mpfr_init2(__mpfr_const_log2, precx); else mpfr_set_prec(__mpfr_const_log2, precx); |