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 | |
25 | #define PRINT_ERRORS 0 |
26 | |
27 | #define TOL 80 |
28 | #define N2 17 |
29 | #define FRAC (32*4) |
30 | |
31 | #define mpbpl (CHAR_BIT * sizeof (mp_limb_t)) |
32 | #define SZ (FRAC / mpbpl + 1) |
33 | typedef mp_limb_t mp1[SZ], mp2[SZ * 2]; |
34 | |
35 | /* This string has 101 hex digits. */ |
36 | static const char exp1[102] = "2" /* point */ |
37 | "b7e151628aed2a6abf7158809cf4f3c762e7160f38b4da56a7" |
38 | "84d9045190cfef324e7738926cfbe5f4bf8d8d8c31d763da07" ; |
39 | static const char hexdig[] = "0123456789abcdef" ; |
40 | |
41 | static void |
42 | print_mpn_hex (const mp_limb_t *x, unsigned size) |
43 | { |
44 | char value[size + 1]; |
45 | unsigned i; |
46 | const unsigned final = (size * 4 > SZ * mpbpl) ? SZ * mpbpl / 4 : size; |
47 | |
48 | memset (s: value, c: '0', n: size); |
49 | |
50 | for (i = 0; i < final ; i++) |
51 | value[size - 1 - i] = hexdig[x[i * 4 / mpbpl] >> (i * 4) % mpbpl & 0xf]; |
52 | |
53 | value[size] = '\0'; |
54 | fputs (s: value, stdout); |
55 | } |
56 | |
57 | static void |
58 | exp_mpn (mp1 ex, mp1 x) |
59 | { |
60 | unsigned n; |
61 | mp1 xp; |
62 | mp2 tmp; |
63 | mp_limb_t chk __attribute__ ((unused)); |
64 | mp1 tol; |
65 | |
66 | memset (s: xp, c: 0, n: sizeof (mp1)); |
67 | memset (s: ex, c: 0, n: sizeof (mp1)); |
68 | xp[FRAC / mpbpl] = (mp_limb_t)1 << FRAC % mpbpl; |
69 | memset (s: tol,c: 0, n: sizeof (mp1)); |
70 | tol[(FRAC - TOL) / mpbpl] = (mp_limb_t)1 << (FRAC - TOL) % mpbpl; |
71 | |
72 | n = 0; |
73 | |
74 | do |
75 | { |
76 | /* Calculate sum(x^n/n!) until the next term is sufficiently small. */ |
77 | |
78 | mpn_mul_n (tmp, xp, x, SZ); |
79 | assert (tmp[SZ * 2 - 1] == 0); |
80 | if (n > 0) |
81 | mpn_divmod_1 (xp, tmp + FRAC / mpbpl, SZ, n); |
82 | chk = mpn_add_n (ex, ex, xp, SZ); |
83 | assert (chk == 0); |
84 | n++; |
85 | assert (n < 80); /* Catch too-high TOL. */ |
86 | } |
87 | while (n < 10 || mpn_cmp (xp, tol, SZ) >= 0); |
88 | } |
89 | |
90 | static int |
91 | mpn_bitsize(const mp_limb_t *SRC_PTR, mp_size_t SIZE) |
92 | { |
93 | int i, j; |
94 | for (i = SIZE - 1; i > 0; i--) |
95 | if (SRC_PTR[i] != 0) |
96 | break; |
97 | for (j = mpbpl - 1; j >= 0; j--) |
98 | if ((SRC_PTR[i] & (mp_limb_t)1 << j) != 0) |
99 | break; |
100 | |
101 | return i * mpbpl + j; |
102 | } |
103 | |
104 | static int |
105 | do_test (void) |
106 | { |
107 | mp1 ex, x, xt, e2, e3; |
108 | int i; |
109 | int errors = 0; |
110 | int failures = 0; |
111 | mp1 maxerror; |
112 | int maxerror_s = 0; |
113 | const double sf = pow (x: 2, mpbpl); |
114 | |
115 | /* assert (mpbpl == mp_bits_per_limb); */ |
116 | assert (FRAC / mpbpl * mpbpl == FRAC); |
117 | |
118 | memset (s: maxerror, c: 0, n: sizeof (mp1)); |
119 | memset (s: xt, c: 0, n: sizeof (mp1)); |
120 | xt[(FRAC - N2) / mpbpl] = (mp_limb_t)1 << (FRAC - N2) % mpbpl; |
121 | |
122 | for (i = 0; i < 1 << N2; i++) |
123 | { |
124 | int e2s, e3s, j; |
125 | double de2; |
126 | |
127 | mpn_mul_1 (x,xt,SZ,i); |
128 | exp_mpn (ex, x); |
129 | de2 = exp (x: i / (double) (1 << N2)); |
130 | for (j = SZ-1; j >= 0; j--) |
131 | { |
132 | e2[j] = (mp_limb_t) de2; |
133 | de2 = (de2 - e2[j]) * sf; |
134 | } |
135 | if (mpn_cmp (ex,e2,SZ) >= 0) |
136 | mpn_sub_n (e3,ex,e2,SZ); |
137 | else |
138 | mpn_sub_n (e3,e2,ex,SZ); |
139 | |
140 | e2s = mpn_bitsize (SRC_PTR: e2,SZ); |
141 | e3s = mpn_bitsize (SRC_PTR: e3,SZ); |
142 | if (e3s >= 0 && e2s - e3s < 54) |
143 | { |
144 | #if PRINT_ERRORS |
145 | printf ("%06x " , i * (0x100000 / (1 << N2))); |
146 | print_mpn_hex (ex, (FRAC / 4) + 1); |
147 | fputs ("\n " ,stdout); |
148 | print_mpn_hex (e2, (FRAC / 4) + 1); |
149 | printf ("\n %c " , |
150 | e2s - e3s < 54 ? e2s - e3s == 53 ? 'e' : 'F' : 'P'); |
151 | print_mpn_hex (e3, (FRAC / 4) + 1); |
152 | putchar ('\n'); |
153 | #endif |
154 | errors += (e2s - e3s == 53); |
155 | failures += (e2s - e3s < 53); |
156 | } |
157 | if (e3s >= maxerror_s |
158 | && mpn_cmp (e3, maxerror, SZ) > 0) |
159 | { |
160 | memcpy (dest: maxerror, src: e3, n: sizeof (mp1)); |
161 | maxerror_s = e3s; |
162 | } |
163 | } |
164 | |
165 | /* Check exp_mpn against precomputed value of exp(1). */ |
166 | memset (s: x, c: '\0', n: sizeof (mp1)); |
167 | x[FRAC / mpbpl] = (mp_limb_t)1 << FRAC % mpbpl; |
168 | exp_mpn (ex, x); |
169 | |
170 | memset (s: e2, c: '\0', n: sizeof (mp1)); |
171 | for (i = -1; i < 100 && i < FRAC / 4; i++) |
172 | e2[(FRAC - i * 4 - 4) / mpbpl] |= ((mp_limb_t) (strchr (s: hexdig, |
173 | c: exp1[i + 1]) |
174 | - hexdig) |
175 | << (FRAC - i * 4 - 4) % mpbpl); |
176 | |
177 | if (mpn_cmp (ex, e2, SZ) >= 0) |
178 | mpn_sub_n (e3, ex, e2, SZ); |
179 | else |
180 | mpn_sub_n (e3, e2, ex, SZ); |
181 | |
182 | printf (format: "%d failures; %d errors; error rate %0.2f%%\n" , failures, errors, |
183 | errors * 100.0 / (double) (1 << N2)); |
184 | fputs (s: "maximum error: " , stdout); |
185 | print_mpn_hex (x: maxerror, size: (FRAC / 4) + 1); |
186 | fputs (s: "\nerror in exp(1): " , stdout); |
187 | print_mpn_hex (x: e3, size: (FRAC / 4) + 1); |
188 | putchar (c: '\n'); |
189 | |
190 | return failures == 0 ? 0 : 1; |
191 | } |
192 | |
193 | #define TIMEOUT 200 |
194 | #define TEST_FUNCTION do_test () |
195 | #include "../test-skeleton.c" |
196 | |