1 | /* Copyright (C) 1997-2022 Free Software Foundation, Inc. |
2 | This file is part of the GNU C Library. |
3 | |
4 | The GNU C Library is free software; you can redistribute it and/or |
5 | modify it under the terms of the GNU Lesser General Public |
6 | License as published by the Free Software Foundation; either |
7 | version 2.1 of the License, or (at your option) any later version. |
8 | |
9 | The GNU C Library is distributed in the hope that it will be useful, |
10 | but WITHOUT ANY WARRANTY; without even the implied warranty of |
11 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
12 | Lesser General Public License for more details. |
13 | |
14 | You should have received a copy of the GNU Lesser General Public |
15 | License along with the GNU C Library; if not, see |
16 | <https://www.gnu.org/licenses/>. */ |
17 | |
18 | #include <stdio.h> |
19 | #include <math.h> |
20 | #include <gmp.h> |
21 | #include <string.h> |
22 | #include <limits.h> |
23 | #include <assert.h> |
24 | #include <stdlib.h> |
25 | |
26 | #define PRINT_ERRORS 0 |
27 | |
28 | #define TOL 80 |
29 | #define N2 18 |
30 | #define FRAC (32*4) |
31 | |
32 | #define mpbpl (CHAR_BIT * sizeof (mp_limb_t)) |
33 | #define SZ (FRAC / mpbpl + 1) |
34 | typedef mp_limb_t mp1[SZ], mp2[SZ * 2]; |
35 | |
36 | #if BITS_PER_MP_LIMB == 64 |
37 | # define LIMB64(L, H) 0x ## H ## L |
38 | #elif BITS_PER_MP_LIMB == 32 |
39 | # define LIMB64(L, H) 0x ## L, 0x ## H |
40 | #else |
41 | # error |
42 | #endif |
43 | |
44 | /* Once upon a time these constants were generated to 400 bits. |
45 | We only need FRAC bits (128) at present, but we retain 384 bits |
46 | in the text Just In Case. */ |
47 | #define CONSTSZ(INT, F1, F2, F3, F4, F5, F6, F7, F8, F9, Fa, Fb, Fc) \ |
48 | LIMB64(F4, F3), LIMB64(F2, F1), INT |
49 | |
50 | static const mp1 mp_exp1 = { |
51 | CONSTSZ (2, b7e15162, 8aed2a6a, bf715880, 9cf4f3c7, 62e7160f, 38b4da56, |
52 | a784d904, 5190cfef, 324e7738, 926cfbe5, f4bf8d8d, 8c31d763) |
53 | }; |
54 | |
55 | static const mp1 mp_log2 = { |
56 | CONSTSZ (0, b17217f7, d1cf79ab, c9e3b398, 03f2f6af, 40f34326, 7298b62d, |
57 | 8a0d175b, 8baafa2b, e7b87620, 6debac98, 559552fb, 4afa1b10) |
58 | }; |
59 | |
60 | static void |
61 | print_mpn_fp (const mp_limb_t *x, unsigned int dp, unsigned int base) |
62 | { |
63 | static const char hexdig[16] = "0123456789abcdef" ; |
64 | unsigned int i; |
65 | mp1 tx; |
66 | |
67 | memcpy (dest: tx, src: x, n: sizeof (mp1)); |
68 | if (base == 16) |
69 | fputs (s: "0x" , stdout); |
70 | assert (x[SZ-1] < base); |
71 | fputc (c: hexdig[x[SZ - 1]], stdout); |
72 | fputc (c: '.', stdout); |
73 | for (i = 0; i < dp; i++) |
74 | { |
75 | tx[SZ - 1] = 0; |
76 | mpn_mul_1 (tx, tx, SZ, base); |
77 | assert (tx[SZ - 1] < base); |
78 | fputc (c: hexdig[tx[SZ - 1]], stdout); |
79 | } |
80 | } |
81 | |
82 | /* Compute e^x. */ |
83 | static void |
84 | exp_mpn (mp1 ex, mp1 x) |
85 | { |
86 | unsigned int n; |
87 | mp1 xp; |
88 | mp2 tmp; |
89 | mp_limb_t chk __attribute__ ((unused)); |
90 | mp1 tol; |
91 | |
92 | memset (s: xp, c: 0, n: sizeof (mp1)); |
93 | memset (s: ex, c: 0, n: sizeof (mp1)); |
94 | xp[FRAC / mpbpl] = (mp_limb_t)1 << FRAC % mpbpl; |
95 | memset (s: tol, c: 0, n: sizeof (mp1)); |
96 | tol[(FRAC - TOL) / mpbpl] = (mp_limb_t)1 << (FRAC - TOL) % mpbpl; |
97 | |
98 | n = 0; |
99 | |
100 | do |
101 | { |
102 | /* Calculate sum(x^n/n!) until the next term is sufficiently small. */ |
103 | |
104 | mpn_mul_n (tmp, xp, x, SZ); |
105 | assert(tmp[SZ * 2 - 1] == 0); |
106 | if (n > 0) |
107 | mpn_divmod_1 (xp, tmp + FRAC / mpbpl, SZ, n); |
108 | chk = mpn_add_n (ex, ex, xp, SZ); |
109 | assert (chk == 0); |
110 | ++n; |
111 | assert (n < 80); /* Catch too-high TOL. */ |
112 | } |
113 | while (n < 10 || mpn_cmp (xp, tol, SZ) >= 0); |
114 | } |
115 | |
116 | /* Calculate 2^x. */ |
117 | static void |
118 | exp2_mpn (mp1 ex, mp1 x) |
119 | { |
120 | mp2 tmp; |
121 | mpn_mul_n (tmp, x, mp_log2, SZ); |
122 | assert(tmp[SZ * 2 - 1] == 0); |
123 | exp_mpn (ex, x: tmp + FRAC / mpbpl); |
124 | } |
125 | |
126 | |
127 | static int |
128 | mpn_bitsize(const mp_limb_t *SRC_PTR, mp_size_t SIZE) |
129 | { |
130 | int i, j; |
131 | for (i = SIZE - 1; i > 0; --i) |
132 | if (SRC_PTR[i] != 0) |
133 | break; |
134 | for (j = mpbpl - 1; j >= 0; --j) |
135 | if ((SRC_PTR[i] & (mp_limb_t)1 << j) != 0) |
136 | break; |
137 | |
138 | return i * mpbpl + j; |
139 | } |
140 | |
141 | static int |
142 | do_test (void) |
143 | { |
144 | mp1 ex, x, xt, e2, e3; |
145 | int i; |
146 | int errors = 0; |
147 | int failures = 0; |
148 | mp1 maxerror; |
149 | int maxerror_s = 0; |
150 | const double sf = pow (x: 2, mpbpl); |
151 | |
152 | /* assert(mpbpl == mp_bits_per_limb); */ |
153 | assert(FRAC / mpbpl * mpbpl == FRAC); |
154 | |
155 | memset (s: maxerror, c: 0, n: sizeof (mp1)); |
156 | memset (s: xt, c: 0, n: sizeof (mp1)); |
157 | xt[(FRAC - N2) / mpbpl] = (mp_limb_t)1 << (FRAC - N2) % mpbpl; |
158 | |
159 | for (i = 0; i < (1 << N2); ++i) |
160 | { |
161 | int e2s, e3s, j; |
162 | double de2; |
163 | |
164 | mpn_mul_1 (x, xt, SZ, i); |
165 | exp2_mpn (ex, x); |
166 | de2 = exp2 (x: i / (double) (1 << N2)); |
167 | for (j = SZ - 1; j >= 0; --j) |
168 | { |
169 | e2[j] = (mp_limb_t) de2; |
170 | de2 = (de2 - e2[j]) * sf; |
171 | } |
172 | if (mpn_cmp (ex, e2, SZ) >= 0) |
173 | mpn_sub_n (e3, ex, e2, SZ); |
174 | else |
175 | mpn_sub_n (e3, e2, ex, SZ); |
176 | |
177 | e2s = mpn_bitsize (SRC_PTR: e2, SZ); |
178 | e3s = mpn_bitsize (SRC_PTR: e3, SZ); |
179 | if (e3s >= 0 && e2s - e3s < 54) |
180 | { |
181 | #if PRINT_ERRORS |
182 | printf ("%06x " , i * (0x100000 / (1 << N2))); |
183 | print_mpn_fp (ex, (FRAC / 4) + 1, 16); |
184 | putchar ('\n'); |
185 | fputs (" " ,stdout); |
186 | print_mpn_fp (e2, (FRAC / 4) + 1, 16); |
187 | putchar ('\n'); |
188 | printf (" %c " , |
189 | e2s - e3s < 54 ? e2s - e3s == 53 ? 'e' : 'F' : 'P'); |
190 | print_mpn_fp (e3, (FRAC / 4) + 1, 16); |
191 | putchar ('\n'); |
192 | #endif |
193 | errors += (e2s - e3s == 53); |
194 | failures += (e2s - e3s < 53); |
195 | } |
196 | if (e3s >= maxerror_s |
197 | && mpn_cmp (e3, maxerror, SZ) > 0) |
198 | { |
199 | memcpy (dest: maxerror, src: e3, n: sizeof (mp1)); |
200 | maxerror_s = e3s; |
201 | } |
202 | } |
203 | |
204 | /* Check exp_mpn against precomputed value of exp(1). */ |
205 | memset (s: x, c: 0, n: sizeof (mp1)); |
206 | x[FRAC / mpbpl] = (mp_limb_t)1 << FRAC % mpbpl; |
207 | exp_mpn (ex, x); |
208 | if (mpn_cmp (ex, mp_exp1, SZ) >= 0) |
209 | mpn_sub_n (e3, ex, mp_exp1, SZ); |
210 | else |
211 | mpn_sub_n (e3, mp_exp1, ex, SZ); |
212 | |
213 | printf (format: "%d failures; %d errors; error rate %0.2f%%\n" , failures, errors, |
214 | errors * 100.0 / (double) (1 << N2)); |
215 | fputs (s: "maximum error: " , stdout); |
216 | print_mpn_fp (x: maxerror, dp: (FRAC / 4) + 1, base: 16); |
217 | putchar (c: '\n'); |
218 | fputs (s: "error in exp(1): " , stdout); |
219 | print_mpn_fp (x: e3, dp: (FRAC / 4) + 1, base: 16); |
220 | putchar (c: '\n'); |
221 | |
222 | return failures == 0 ? 0 : 1; |
223 | } |
224 | |
225 | #define TIMEOUT 300 |
226 | #define TEST_FUNCTION do_test () |
227 | #include "../test-skeleton.c" |
228 | |