Annotation of OpenXM_contrib/gmp/mpfr/const_log2.c, Revision 1.1.1.1
1.1 ohara 1: /* mpfr_const_log2 -- compute natural logarithm of 2
2:
3: Copyright 1999, 2001 Free Software Foundation, Inc.
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 Lesser General Public License as published by
9: the Free Software Foundation; either version 2.1 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 Lesser General Public
15: License for more details.
16:
17: You should have received a copy of the GNU Lesser 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 "gmp.h"
24: #include "gmp-impl.h"
25: #include "longlong.h"
26: #include "mpfr.h"
27: #include "mpfr-impl.h"
28:
29: mpfr_t __mpfr_const_log2; /* stored value of log(2) */
30: mp_prec_t __mpfr_const_log2_prec=0; /* precision of stored value */
31: mp_rnd_t __mpfr_const_log2_rnd; /* rounding mode of stored value */
32:
33: static int mpfr_aux_log2 _PROTO ((mpfr_ptr, mpz_srcptr, int, int));
34: static int mpfr_const_aux_log2 _PROTO ((mpfr_ptr, mp_rnd_t));
35:
36: #define A
37: #define A1 1
38: #define A2 1
39: #undef B
40: #define C
41: #define C1 2
42: #define C2 1
43: #define NO_FACTORIAL
44: #undef R_IS_RATIONAL
45: #define GENERIC mpfr_aux_log2
46: #include "generic.c"
47: #undef A
48: #undef A1
49: #undef A2
50: #undef NO_FACTORIAL
51: #undef GENERIC
52: #undef C
53: #undef C1
54: #undef C2
55:
56: static int
57: mpfr_const_aux_log2 (mpfr_ptr mylog, mp_rnd_t rnd_mode)
58: {
59: mp_prec_t prec;
60: mpfr_t tmp1, tmp2, result,tmp3;
61: mpz_t cst;
62: int good = 0;
63: int logn;
64: mp_prec_t prec_i_want = MPFR_PREC(mylog);
65: mp_prec_t prec_x;
66:
67: mpz_init(cst);
68: logn = _mpfr_ceil_log2 ((double) MPFR_PREC(mylog));
69: prec_x = prec_i_want + logn;
70: while (!good){
71: prec = _mpfr_ceil_log2 ((double) prec_x);
72: mpfr_init2(tmp1, prec_x);
73: mpfr_init2(result, prec_x);
74: mpfr_init2(tmp2, prec_x);
75: mpfr_init2(tmp3, prec_x);
76: mpz_set_ui(cst, 1);
77: mpfr_aux_log2(tmp1, cst, 4, prec-2);
78: mpfr_div_2ui(tmp1, tmp1, 4, GMP_RNDD);
79: mpfr_mul_ui(tmp1, tmp1, 15, GMP_RNDD);
80:
81: mpz_set_ui(cst, 3);
82: mpfr_aux_log2(tmp2, cst, 7, prec-2);
83: mpfr_div_2ui(tmp2, tmp2, 7, GMP_RNDD);
84: mpfr_mul_ui(tmp2, tmp2, 5*3, GMP_RNDD);
85: mpfr_sub(result, tmp1, tmp2, GMP_RNDD);
86:
87: mpz_set_ui(cst, 13);
88: mpfr_aux_log2(tmp3, cst, 8, prec-2);
89: mpfr_div_2ui(tmp3, tmp3, 8, GMP_RNDD);
90: mpfr_mul_ui(tmp3, tmp3, 3*13, GMP_RNDD);
91: mpfr_sub(result, result, tmp3, GMP_RNDD);
92:
93: mpfr_clear(tmp1);
94: mpfr_clear(tmp2);
95: mpfr_clear(tmp3);
96: if (mpfr_can_round(result, prec_x, GMP_RNDD, rnd_mode, prec_i_want)){
97: mpfr_set(mylog, result, rnd_mode);
98: good = 1;
99: } else
100: {
101: prec_x += logn;
102: }
103: mpfr_clear(result);
104: }
105: mpz_clear(cst);
106: return 0;
107: }
108:
109: /* Cross-over point from nai"ve Taylor series to binary splitting,
110: obtained experimentally on a Pentium II. Optimal value for
111: target machine should be determined by tuneup. */
112: #define LOG2_THRESHOLD 25000
113:
114: /* set x to log(2) rounded to precision MPFR_PREC(x) with direction rnd_mode
115:
116: use formula log(2) = sum(1/k/2^k, k=1..infinity)
117:
118: whence 2^N*log(2) = S(N) + R(N)
119:
120: where S(N) = sum(2^(N-k)/k, k=1..N-1)
121: and R(N) = sum(1/k/2^(k-N), k=N..infinity) < 2/N
122:
123: Let S'(N) = sum(floor(2^(N-k)/k), k=1..N-1)
124:
125: Then 2^N*log(2)-S'(N) <= N-1+2/N <= N for N>=2.
126: */
127: void
128: mpfr_const_log2 (mpfr_ptr x, mp_rnd_t rnd_mode)
129: {
130: mp_prec_t N, k, precx;
131: mpz_t s, t, u;
132:
133: precx = MPFR_PREC(x);
134: MPFR_CLEAR_FLAGS(x);
135:
136: /* has stored value enough precision ? */
137: if (precx <= __mpfr_const_log2_prec)
138: {
139: if ((rnd_mode == __mpfr_const_log2_rnd) ||
140: mpfr_can_round (__mpfr_const_log2, __mpfr_const_log2_prec - 1,
141: __mpfr_const_log2_rnd, rnd_mode, precx))
142: {
143: mpfr_set (x, __mpfr_const_log2, rnd_mode);
144: return;
145: }
146: }
147:
148: /* need to recompute */
149: if (precx < LOG2_THRESHOLD) /* use nai"ve Taylor series evaluation */
150: {
151: /* the following was checked by exhaustive search to give a correct
152: result for all 4 rounding modes up to precx = 13500 */
153: N = precx + 2 * _mpfr_ceil_log2 ((double) precx) + 1;
154:
155: mpz_init (s); /* set to zero */
156: mpz_init (u);
157: mpz_init_set_ui (t, 1);
158:
159: /* use log(2) = sum((6*k-1)/(2*k^2-k)/2^(2*k+1), k=1..infinity) */
160: mpz_mul_2exp (t, t, N-1);
161: for (k=1; k<=N/2; k++)
162: {
163: mpz_div_2exp (t, t, 2);
164: mpz_mul_ui (u, t, 6*k-1);
165: mpz_fdiv_q_ui (u, u, k*(2*k-1));
166: mpz_add (s, s, u);
167: }
168:
169: mpfr_set_z (x, s, rnd_mode);
170: MPFR_EXP(x) -= N;
171: mpz_clear (s);
172: mpz_clear (t);
173: mpz_clear (u);
174: }
175: else /* use binary splitting method */
176: mpfr_const_aux_log2(x, rnd_mode);
177:
178: /* store computed value */
179: if (__mpfr_const_log2_prec == 0)
180: mpfr_init2 (__mpfr_const_log2, precx);
181: else
182: mpfr_set_prec (__mpfr_const_log2, precx);
183:
184: mpfr_set (__mpfr_const_log2, x, rnd_mode);
185: __mpfr_const_log2_prec = precx;
186: __mpfr_const_log2_rnd = rnd_mode;
187: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to const_log2.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpfr