1 /* Copyright (C) 1997-2021 Free Software Foundation, Inc.
2 This file is part of the GNU C Library.
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.
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.
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/>. */
26 #define PRINT_ERRORS 0
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];
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
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
50 static const mp1 mp_exp1
= {
51 CONSTSZ (2, b7e15162
, 8aed2a6a
, bf715880
, 9cf4f3c7
, 62e7160f
, 38b4da56
,
52 a784d904
, 5190cfef
, 324e7738
, 926cfbe5
, f4bf8d8d
, 8c31d763
)
55 static const mp1 mp_log2
= {
56 CONSTSZ (0, b17217f7
, d1cf79ab
, c9e3b398
, 03f2f6af
, 40f34326
, 7298b62d
,
57 8a0d175b
, 8baafa2b
, e7b87620
, 6debac98
, 559552fb
, 4afa1b10
)
61 print_mpn_fp (const mp_limb_t
*x
, unsigned int dp
, unsigned int base
)
63 static const char hexdig
[16] = "0123456789abcdef";
67 memcpy (tx
, x
, sizeof (mp1
));
70 assert (x
[SZ
-1] < base
);
71 fputc (hexdig
[x
[SZ
- 1]], stdout
);
73 for (i
= 0; i
< dp
; i
++)
76 mpn_mul_1 (tx
, tx
, SZ
, base
);
77 assert (tx
[SZ
- 1] < base
);
78 fputc (hexdig
[tx
[SZ
- 1]], stdout
);
84 exp_mpn (mp1 ex
, mp1 x
)
89 mp_limb_t chk
__attribute__ ((unused
));
92 memset (xp
, 0, sizeof (mp1
));
93 memset (ex
, 0, sizeof (mp1
));
94 xp
[FRAC
/ mpbpl
] = (mp_limb_t
)1 << FRAC
% mpbpl
;
95 memset (tol
, 0, sizeof (mp1
));
96 tol
[(FRAC
- TOL
) / mpbpl
] = (mp_limb_t
)1 << (FRAC
- TOL
) % mpbpl
;
102 /* Calculate sum(x^n/n!) until the next term is sufficiently small. */
104 mpn_mul_n (tmp
, xp
, x
, SZ
);
105 assert(tmp
[SZ
* 2 - 1] == 0);
107 mpn_divmod_1 (xp
, tmp
+ FRAC
/ mpbpl
, SZ
, n
);
108 chk
= mpn_add_n (ex
, ex
, xp
, SZ
);
111 assert (n
< 80); /* Catch too-high TOL. */
113 while (n
< 10 || mpn_cmp (xp
, tol
, SZ
) >= 0);
118 exp2_mpn (mp1 ex
, mp1 x
)
121 mpn_mul_n (tmp
, x
, mp_log2
, SZ
);
122 assert(tmp
[SZ
* 2 - 1] == 0);
123 exp_mpn (ex
, tmp
+ FRAC
/ mpbpl
);
128 mpn_bitsize(const mp_limb_t
*SRC_PTR
, mp_size_t SIZE
)
131 for (i
= SIZE
- 1; i
> 0; --i
)
134 for (j
= mpbpl
- 1; j
>= 0; --j
)
135 if ((SRC_PTR
[i
] & (mp_limb_t
)1 << j
) != 0)
138 return i
* mpbpl
+ j
;
144 mp1 ex
, x
, xt
, e2
, e3
;
150 const double sf
= pow (2, mpbpl
);
152 /* assert(mpbpl == mp_bits_per_limb); */
153 assert(FRAC
/ mpbpl
* mpbpl
== FRAC
);
155 memset (maxerror
, 0, sizeof (mp1
));
156 memset (xt
, 0, sizeof (mp1
));
157 xt
[(FRAC
- N2
) / mpbpl
] = (mp_limb_t
)1 << (FRAC
- N2
) % mpbpl
;
159 for (i
= 0; i
< (1 << N2
); ++i
)
164 mpn_mul_1 (x
, xt
, SZ
, i
);
166 de2
= exp2 (i
/ (double) (1 << N2
));
167 for (j
= SZ
- 1; j
>= 0; --j
)
169 e2
[j
] = (mp_limb_t
) de2
;
170 de2
= (de2
- e2
[j
]) * sf
;
172 if (mpn_cmp (ex
, e2
, SZ
) >= 0)
173 mpn_sub_n (e3
, ex
, e2
, SZ
);
175 mpn_sub_n (e3
, e2
, ex
, SZ
);
177 e2s
= mpn_bitsize (e2
, SZ
);
178 e3s
= mpn_bitsize (e3
, SZ
);
179 if (e3s
>= 0 && e2s
- e3s
< 54)
182 printf ("%06x ", i
* (0x100000 / (1 << N2
)));
183 print_mpn_fp (ex
, (FRAC
/ 4) + 1, 16);
186 print_mpn_fp (e2
, (FRAC
/ 4) + 1, 16);
189 e2s
- e3s
< 54 ? e2s
- e3s
== 53 ? 'e' : 'F' : 'P');
190 print_mpn_fp (e3
, (FRAC
/ 4) + 1, 16);
193 errors
+= (e2s
- e3s
== 53);
194 failures
+= (e2s
- e3s
< 53);
196 if (e3s
>= maxerror_s
197 && mpn_cmp (e3
, maxerror
, SZ
) > 0)
199 memcpy (maxerror
, e3
, sizeof (mp1
));
204 /* Check exp_mpn against precomputed value of exp(1). */
205 memset (x
, 0, sizeof (mp1
));
206 x
[FRAC
/ mpbpl
] = (mp_limb_t
)1 << FRAC
% mpbpl
;
208 if (mpn_cmp (ex
, mp_exp1
, SZ
) >= 0)
209 mpn_sub_n (e3
, ex
, mp_exp1
, SZ
);
211 mpn_sub_n (e3
, mp_exp1
, ex
, SZ
);
213 printf ("%d failures; %d errors; error rate %0.2f%%\n", failures
, errors
,
214 errors
* 100.0 / (double) (1 << N2
));
215 fputs ("maximum error: ", stdout
);
216 print_mpn_fp (maxerror
, (FRAC
/ 4) + 1, 16);
218 fputs ("error in exp(1): ", stdout
);
219 print_mpn_fp (e3
, (FRAC
/ 4) + 1, 16);
222 return failures
== 0 ? 0 : 1;
226 #define TEST_FUNCTION do_test ()
227 #include "../test-skeleton.c"