1 /* Copyright (C) 1997, 1998, 2000, 2006 Free Software Foundation, Inc.
2 This file is part of the GNU C Library.
3 Contributed by Geoffrey Keating <Geoff.Keating@anu.edu.au>, 1997.
5 The GNU C Library is free software; you can redistribute it and/or
6 modify it under the terms of the GNU Lesser General Public
7 License as published by the Free Software Foundation; either
8 version 2.1 of the License, or (at your option) any later version.
10 The GNU C Library is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 Lesser General Public License for more details.
15 You should have received a copy of the GNU Lesser General Public
16 License along with the GNU C Library; if not, write to the Free
17 Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
28 #define PRINT_ERRORS 0
34 #define mpbpl (CHAR_BIT * sizeof (mp_limb_t))
35 #define SZ (FRAC / mpbpl + 1)
36 typedef mp_limb_t mp1[SZ], mp2[SZ * 2];
38 /* This string has 101 hex digits. */
39 static const char exp1[102] = "2" /* point */
40 "b7e151628aed2a6abf7158809cf4f3c762e7160f38b4da56a7"
41 "84d9045190cfef324e7738926cfbe5f4bf8d8d8c31d763da07";
42 static const char exp_m1[102] = "0" /* point */
43 "5e2d58d8b3bcdf1abadec7829054f90dda9805aab56c773330"
44 "24b9d0a507daedb16400bf472b4215b8245b669d90d27a5aea";
46 static const char hexdig[] = "0123456789abcdef";
49 print_mpn_fp (const mp_limb_t *x, unsigned int dp, unsigned int base)
54 memcpy (tx, x, sizeof (mp1));
57 assert (x[SZ-1] < base);
58 fputc (hexdig[x[SZ - 1]], stdout);
60 for (i = 0; i < dp; i++)
63 mpn_mul_1 (tx, tx, SZ, base);
64 assert (tx[SZ - 1] < base);
65 fputc (hexdig[tx[SZ - 1]], stdout);
70 read_mpn_hex(mp_limb_t *x, const char *str)
74 memset (x, 0, sizeof (mp1));
75 for (i = -1; i < 100 && i < FRAC / 4; ++i)
76 x[(FRAC - i * 4 - 4) / mpbpl] |= ((mp_limb_t) (strchr (hexdig, str[i + 1])
78 << (FRAC - i * 4 - 4) % mpbpl);
81 static mp_limb_t *get_log2(void) __attribute__((const));
86 static int log2_m_inited = 0;
87 static const char log2[102] = "0" /* point */
88 "b17217f7d1cf79abc9e3b39803f2f6af40f343267298b62d8a"
89 "0d175b8baafa2be7b876206debac98559552fb4afa1b10ed2e";
93 read_mpn_hex (log2_m, log2);
101 exp_mpn (mp1 ex, mp1 x)
106 mp_limb_t chk, round;
109 memset (xp, 0, sizeof (mp1));
110 memset (ex, 0, sizeof (mp1));
111 xp[FRAC / mpbpl] = (mp_limb_t)1 << FRAC % mpbpl;
112 memset (tol, 0, sizeof (mp1));
113 tol[(FRAC - TOL) / mpbpl] = (mp_limb_t)1 << (FRAC - TOL) % mpbpl;
119 /* Calculate sum(x^n/n!) until the next term is sufficiently small. */
121 mpn_mul_n (tmp, xp, x, SZ);
122 assert(tmp[SZ * 2 - 1] == 0);
124 round = mpn_divmod_1 (xp, tmp + FRAC / mpbpl, SZ, n);
125 chk = mpn_add_n (ex, ex, xp, SZ);
128 assert (n < 80); /* Catch too-high TOL. */
130 while (n < 10 || mpn_cmp (xp, tol, SZ) >= 0);
135 exp2_mpn (mp1 ex, mp1 x)
138 mpn_mul_n (tmp, x, get_log2 (), SZ);
139 assert(tmp[SZ * 2 - 1] == 0);
140 exp_mpn (ex, tmp + FRAC / mpbpl);
145 mpn_bitsize(const mp_limb_t *SRC_PTR, mp_size_t SIZE)
148 for (i = SIZE - 1; i > 0; --i)
151 for (j = mpbpl - 1; j >= 0; --j)
152 if ((SRC_PTR[i] & (mp_limb_t)1 << j) != 0)
155 return i * mpbpl + j;
161 mp1 ex, x, xt, e2, e3;
167 const double sf = pow (2, mpbpl);
169 /* assert(mpbpl == mp_bits_per_limb); */
170 assert(FRAC / mpbpl * mpbpl == FRAC);
172 memset (maxerror, 0, sizeof (mp1));
173 memset (xt, 0, sizeof (mp1));
174 xt[(FRAC - N2) / mpbpl] = (mp_limb_t)1 << (FRAC - N2) % mpbpl;
176 for (i = 0; i < (1 << N2); ++i)
181 mpn_mul_1 (x, xt, SZ, i);
183 de2 = exp2 (i / (double) (1 << N2));
184 for (j = SZ - 1; j >= 0; --j)
186 e2[j] = (mp_limb_t) de2;
187 de2 = (de2 - e2[j]) * sf;
189 if (mpn_cmp (ex, e2, SZ) >= 0)
190 mpn_sub_n (e3, ex, e2, SZ);
192 mpn_sub_n (e3, e2, ex, SZ);
194 e2s = mpn_bitsize (e2, SZ);
195 e3s = mpn_bitsize (e3, SZ);
196 if (e3s >= 0 && e2s - e3s < 54)
199 printf ("%06x ", i * (0x100000 / (1 << N2)));
200 print_mpn_fp (ex, (FRAC / 4) + 1, 16);
203 print_mpn_fp (e2, (FRAC / 4) + 1, 16);
206 e2s - e3s < 54 ? e2s - e3s == 53 ? 'e' : 'F' : 'P');
207 print_mpn_fp (e3, (FRAC / 4) + 1, 16);
210 errors += (e2s - e3s == 53);
211 failures += (e2s - e3s < 53);
213 if (e3s >= maxerror_s
214 && mpn_cmp (e3, maxerror, SZ) > 0)
216 memcpy (maxerror, e3, sizeof (mp1));
221 /* Check exp_mpn against precomputed value of exp(1). */
222 memset (x, 0, sizeof (mp1));
223 x[FRAC / mpbpl] = (mp_limb_t)1 << FRAC % mpbpl;
225 read_mpn_hex (e2, exp1);
226 if (mpn_cmp (ex, e2, SZ) >= 0)
227 mpn_sub_n (e3, ex, e2, SZ);
229 mpn_sub_n (e3, e2, ex, SZ);
231 printf ("%d failures; %d errors; error rate %0.2f%%\n", failures, errors,
232 errors * 100.0 / (double) (1 << N2));
233 fputs ("maximum error: ", stdout);
234 print_mpn_fp (maxerror, (FRAC / 4) + 1, 16);
236 fputs ("error in exp(1): ", stdout);
237 print_mpn_fp (e3, (FRAC / 4) + 1, 16);
240 return failures == 0 ? 0 : 1;