summaryrefslogtreecommitdiff
path: root/log2.c
diff options
context:
space:
mode:
authorzimmerma <zimmerma@280ebfd0-de03-0410-8827-d642c229c3f4>2001-11-30 13:26:03 +0000
committerzimmerma <zimmerma@280ebfd0-de03-0410-8827-d642c229c3f4>2001-11-30 13:26:03 +0000
commit8456e81ad0e5afa3e233d8459a82a6ab44534aeb (patch)
tree8993a410e651bff05c409f6e81ccf31e734e01e8 /log2.c
parent8d6244b0cdfb304290f19fabea7d91de590adf7f (diff)
downloadmpfr-8456e81ad0e5afa3e233d8459a82a6ab44534aeb.tar.gz
rewritten part with Taylor series
git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@1600 280ebfd0-de03-0410-8827-d642c229c3f4
Diffstat (limited to 'log2.c')
-rw-r--r--log2.c34
1 files changed, 17 insertions, 17 deletions
diff --git a/log2.c b/log2.c
index fcfc76920..774eb072c 100644
--- a/log2.c
+++ b/log2.c
@@ -127,37 +127,37 @@ mpfr_const_aux_log2 (mpfr_ptr mylog, mp_rnd_t rnd_mode)
void
mpfr_const_log2 (mpfr_ptr x, mp_rnd_t rnd_mode)
{
- int N, oldN, k, precx; mpz_t s, t, u;
+ int N, k, precx; mpz_t s, t, u;
precx = MPFR_PREC(x);
MPFR_CLEAR_FLAGS(x);
/* has stored value enough precision ? */
- if (precx <= __mpfr_const_log2_prec) {
- if (rnd_mode==__mpfr_const_log2_rnd || mpfr_can_round(__mpfr_const_log2,
- __mpfr_const_log2_prec, __mpfr_const_log2_rnd, rnd_mode, precx))
- {
- mpfr_set(x, __mpfr_const_log2, rnd_mode); return;
- }
- }
-
- /* need to recompute */
- if (precx < LOG2_THRESHOLD) /* use nai"ve Taylor series evaluation */
+ if (precx <= __mpfr_const_log2_prec)
{
- N=2;
- do
+ if ((rnd_mode == __mpfr_const_log2_rnd) ||
+ mpfr_can_round (__mpfr_const_log2, __mpfr_const_log2_prec - 1,
+ __mpfr_const_log2_rnd, rnd_mode, precx))
{
- oldN = N;
- N = precx + _mpfr_ceil_log2 ((double) N);
+ mpfr_set (x, __mpfr_const_log2, rnd_mode);
+ return;
}
- while (N != oldN);
+ }
+
+ /* need to recompute */
+ if (precx < LOG2_THRESHOLD) /* use nai"ve Taylor series evaluation */
+ {
+ /* the following was checked by exhaustive search to give a correct
+ result for all 4 rounding modes up to precx = 13500 */
+ N = precx + 2 * _mpfr_ceil_log2 ((double) precx) + 1;
+
mpz_init (s); /* set to zero */
mpz_init (u);
mpz_init_set_ui (t, 1);
/* use log(2) = sum((6*k-1)/(2*k^2-k)/2^(2*k+1), k=1..infinity) */
mpz_mul_2exp (t, t, N-1);
- for (k=1; k<N/2; k++)
+ for (k=1; k<=N/2; k++)
{
mpz_div_2exp (t, t, 2);
mpz_mul_ui (u, t, 6*k-1);