]>
sourceware.org Git - glibc.git/blob - math/atest-exp.c
1 /* Copyright (C) 1997-2024 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/>. */
25 #define PRINT_ERRORS 0
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];
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";
42 print_mpn_hex (const mp_limb_t
*x
, unsigned size
)
46 const unsigned final
= (size
* 4 > SZ
* mpbpl
) ? SZ
* mpbpl
/ 4 : size
;
48 memset (value
, '0', size
);
50 for (i
= 0; i
< final
; i
++)
51 value
[size
- 1 - i
] = hexdig
[x
[i
* 4 / mpbpl
] >> (i
* 4) % mpbpl
& 0xf];
54 fputs (value
, stdout
);
58 exp_mpn (mp1 ex
, mp1 x
)
63 mp_limb_t chk
__attribute__ ((unused
));
66 memset (xp
, 0, sizeof (mp1
));
67 memset (ex
, 0, sizeof (mp1
));
68 xp
[FRAC
/ mpbpl
] = (mp_limb_t
)1 << FRAC
% mpbpl
;
69 memset (tol
,0, sizeof (mp1
));
70 tol
[(FRAC
- TOL
) / mpbpl
] = (mp_limb_t
)1 << (FRAC
- TOL
) % mpbpl
;
76 /* Calculate sum(x^n/n!) until the next term is sufficiently small. */
78 mpn_mul_n (tmp
, xp
, x
, SZ
);
79 assert (tmp
[SZ
* 2 - 1] == 0);
81 mpn_divmod_1 (xp
, tmp
+ FRAC
/ mpbpl
, SZ
, n
);
82 chk
= mpn_add_n (ex
, ex
, xp
, SZ
);
85 assert (n
< 80); /* Catch too-high TOL. */
87 while (n
< 10 || mpn_cmp (xp
, tol
, SZ
) >= 0);
91 mpn_bitsize(const mp_limb_t
*SRC_PTR
, mp_size_t SIZE
)
94 for (i
= SIZE
- 1; i
> 0; i
--)
97 for (j
= mpbpl
- 1; j
>= 0; j
--)
98 if ((SRC_PTR
[i
] & (mp_limb_t
)1 << j
) != 0)
101 return i
* mpbpl
+ j
;
107 mp1 ex
, x
, xt
, e2
, e3
;
113 const double sf
= pow (2, mpbpl
);
115 /* assert (mpbpl == mp_bits_per_limb); */
116 assert (FRAC
/ mpbpl
* mpbpl
== FRAC
);
118 memset (maxerror
, 0, sizeof (mp1
));
119 memset (xt
, 0, sizeof (mp1
));
120 xt
[(FRAC
- N2
) / mpbpl
] = (mp_limb_t
)1 << (FRAC
- N2
) % mpbpl
;
122 for (i
= 0; i
< 1 << N2
; i
++)
127 mpn_mul_1 (x
,xt
,SZ
,i
);
129 de2
= exp (i
/ (double) (1 << N2
));
130 for (j
= SZ
-1; j
>= 0; j
--)
132 e2
[j
] = (mp_limb_t
) de2
;
133 de2
= (de2
- e2
[j
]) * sf
;
135 if (mpn_cmp (ex
,e2
,SZ
) >= 0)
136 mpn_sub_n (e3
,ex
,e2
,SZ
);
138 mpn_sub_n (e3
,e2
,ex
,SZ
);
140 e2s
= mpn_bitsize (e2
,SZ
);
141 e3s
= mpn_bitsize (e3
,SZ
);
142 if (e3s
>= 0 && e2s
- e3s
< 54)
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);
150 e2s
- e3s
< 54 ? e2s
- e3s
== 53 ? 'e' : 'F' : 'P');
151 print_mpn_hex (e3
, (FRAC
/ 4) + 1);
154 errors
+= (e2s
- e3s
== 53);
155 failures
+= (e2s
- e3s
< 53);
157 if (e3s
>= maxerror_s
158 && mpn_cmp (e3
, maxerror
, SZ
) > 0)
160 memcpy (maxerror
, e3
, sizeof (mp1
));
165 /* Check exp_mpn against precomputed value of exp(1). */
166 memset (x
, '\0', sizeof (mp1
));
167 x
[FRAC
/ mpbpl
] = (mp_limb_t
)1 << FRAC
% mpbpl
;
170 memset (e2
, '\0', sizeof (mp1
));
171 for (i
= -1; i
< 100 && i
< FRAC
/ 4; i
++)
172 e2
[(FRAC
- i
* 4 - 4) / mpbpl
] |= ((mp_limb_t
) (strchr (hexdig
,
175 << (FRAC
- i
* 4 - 4) % mpbpl
);
177 if (mpn_cmp (ex
, e2
, SZ
) >= 0)
178 mpn_sub_n (e3
, ex
, e2
, SZ
);
180 mpn_sub_n (e3
, e2
, ex
, SZ
);
182 printf ("%d failures; %d errors; error rate %0.2f%%\n", failures
, errors
,
183 errors
* 100.0 / (double) (1 << N2
));
184 fputs ("maximum error: ", stdout
);
185 print_mpn_hex (maxerror
, (FRAC
/ 4) + 1);
186 fputs ("\nerror in exp(1): ", stdout
);
187 print_mpn_hex (e3
, (FRAC
/ 4) + 1);
190 return failures
== 0 ? 0 : 1;
194 #define TEST_FUNCTION do_test ()
195 #include "../test-skeleton.c"
This page took 0.39647 seconds and 6 git commands to generate.