.
[glibc/history.git] / math / atest-exp2.c
blob2c028a2b8386852dd7c5ab6693dc280b8563ba2c
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
18 02111-1307 USA. */
20 #include <stdio.h>
21 #include <math.h>
22 #include <gmp.h>
23 #include <string.h>
24 #include <limits.h>
25 #include <assert.h>
26 #include <stdlib.h>
28 #define PRINT_ERRORS 0
30 #define TOL 80
31 #define N2 18
32 #define FRAC (32*4)
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";
48 static void
49 print_mpn_fp (const mp_limb_t *x, unsigned int dp, unsigned int base)
51 unsigned int i;
52 mp1 tx;
54 memcpy (tx, x, sizeof (mp1));
55 if (base == 16)
56 fputs ("0x", stdout);
57 assert (x[SZ-1] < base);
58 fputc (hexdig[x[SZ - 1]], stdout);
59 fputc ('.', stdout);
60 for (i = 0; i < dp; i++)
62 tx[SZ - 1] = 0;
63 mpn_mul_1 (tx, tx, SZ, base);
64 assert (tx[SZ - 1] < base);
65 fputc (hexdig[tx[SZ - 1]], stdout);
69 static void
70 read_mpn_hex(mp_limb_t *x, const char *str)
72 int i;
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])
77 - hexdig)
78 << (FRAC - i * 4 - 4) % mpbpl);
81 static mp_limb_t *get_log2(void) __attribute__((const));
82 static mp_limb_t *
83 get_log2(void)
85 static mp1 log2_m;
86 static int log2_m_inited = 0;
87 static const char log2[102] = "0" /* point */
88 "b17217f7d1cf79abc9e3b39803f2f6af40f343267298b62d8a"
89 "0d175b8baafa2be7b876206debac98559552fb4afa1b10ed2e";
91 if (!log2_m_inited)
93 read_mpn_hex (log2_m, log2);
94 log2_m_inited = 1;
96 return log2_m;
99 /* Compute e^x. */
100 static void
101 exp_mpn (mp1 ex, mp1 x)
103 unsigned int n;
104 mp1 xp;
105 mp2 tmp;
106 mp_limb_t chk, round;
107 mp1 tol;
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;
115 n = 0;
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);
123 if (n > 0)
124 round = mpn_divmod_1 (xp, tmp + FRAC / mpbpl, SZ, n);
125 chk = mpn_add_n (ex, ex, xp, SZ);
126 assert (chk == 0);
127 ++n;
128 assert (n < 80); /* Catch too-high TOL. */
130 while (n < 10 || mpn_cmp (xp, tol, SZ) >= 0);
133 /* Calculate 2^x. */
134 static void
135 exp2_mpn (mp1 ex, mp1 x)
137 mp2 tmp;
138 mpn_mul_n (tmp, x, get_log2 (), SZ);
139 assert(tmp[SZ * 2 - 1] == 0);
140 exp_mpn (ex, tmp + FRAC / mpbpl);
144 static int
145 mpn_bitsize(const mp_limb_t *SRC_PTR, mp_size_t SIZE)
147 int i, j;
148 for (i = SIZE - 1; i > 0; --i)
149 if (SRC_PTR[i] != 0)
150 break;
151 for (j = mpbpl - 1; j >= 0; --j)
152 if ((SRC_PTR[i] & (mp_limb_t)1 << j) != 0)
153 break;
155 return i * mpbpl + j;
159 main (void)
161 mp1 ex, x, xt, e2, e3;
162 int i;
163 int errors = 0;
164 int failures = 0;
165 mp1 maxerror;
166 int maxerror_s = 0;
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)
178 int e2s, e3s, j;
179 double de2;
181 mpn_mul_1 (x, xt, SZ, i);
182 exp2_mpn (ex, x);
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);
191 else
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)
198 #if PRINT_ERRORS
199 printf ("%06x ", i * (0x100000 / (1 << N2)));
200 print_mpn_fp (ex, (FRAC / 4) + 1, 16);
201 putchar ('\n');
202 fputs (" ",stdout);
203 print_mpn_fp (e2, (FRAC / 4) + 1, 16);
204 putchar ('\n');
205 printf (" %c ",
206 e2s - e3s < 54 ? e2s - e3s == 53 ? 'e' : 'F' : 'P');
207 print_mpn_fp (e3, (FRAC / 4) + 1, 16);
208 putchar ('\n');
209 #endif
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));
217 maxerror_s = e3s;
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;
224 exp_mpn (ex, x);
225 read_mpn_hex (e2, exp1);
226 if (mpn_cmp (ex, e2, SZ) >= 0)
227 mpn_sub_n (e3, ex, e2, SZ);
228 else
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);
235 putchar ('\n');
236 fputs ("error in exp(1): ", stdout);
237 print_mpn_fp (e3, (FRAC / 4) + 1, 16);
238 putchar ('\n');
240 return failures == 0 ? 0 : 1;