1 /*
2  * SPDX-License-Identifier: BSD-3-Clause
3  *
4  * Copyright © 2022 Keith Packard
5  *
6  * Redistribution and use in source and binary forms, with or without
7  * modification, are permitted provided that the following conditions
8  * are met:
9  *
10  * 1. Redistributions of source code must retain the above copyright
11  *    notice, this list of conditions and the following disclaimer.
12  *
13  * 2. Redistributions in binary form must reproduce the above
14  *    copyright notice, this list of conditions and the following
15  *    disclaimer in the documentation and/or other materials provided
16  *    with the distribution.
17  *
18  * 3. Neither the name of the copyright holder nor the names of its
19  *    contributors may be used to endorse or promote products derived
20  *    from this software without specific prior written permission.
21  *
22  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
23  * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
24  * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
25  * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
26  * COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
27  * INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
28  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
29  * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
31  * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
32  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
33  * OF THE POSSIBILITY OF SUCH DAMAGE.
34  */
35 
36 #define _GNU_SOURCE
37 #include <math.h>
38 #include <stdbool.h>
39 #include <stdio.h>
40 #include <float.h>
41 #include <stdlib.h>
42 #include <string.h>
43 
44 #ifdef _TEST_LONG_DOUBLE
45 
46 static long double max_error;
47 
48 bool
within_error(long double expect,long double result,long double error)49 within_error(long double expect, long double result, long double error)
50 {
51     long double difference;
52     long double e = 1.0L;
53 
54     if (isnan(expect) && isnan(result))
55         return true;
56 
57     if (expect == result)
58         return true;
59 
60     if (expect != 0)
61         e = scalbnl(1.0L, -ilogbl(expect));
62 
63     difference = fabsl(expect - result) * e;
64 
65     if (difference > max_error)
66         max_error = difference;
67 
68     return difference <= error;
69 }
70 
71 int
check_long_double(const char * name,int i,long double prec,long double expect,long double result)72 check_long_double(const char *name, int i, long double prec, long double expect, long double result)
73 {
74     if (!within_error(expect, result, prec)) {
75         long double diff = fabsl(expect - result);
76 #ifdef __PICOLIBC__
77 #ifdef _WANT_IO_LONG_DOUBLE
78         printf("%s test %d got %La expect %La diff %La\n", name, i, result, expect, diff);
79 #else
80         printf("%s test %d got %a expect %a diff %a\n", name, i, (double) result, (double) expect, (double) diff);
81 #endif
82 #else
83 //        printf("%s test %d got %.33Lg expect %.33Lg diff %.33Lg\n", name, i, result, expect, diff);
84         printf("%s test %d got %La expect %La diff %La\n", name, i, result, expect, diff);
85 #endif
86         return 1;
87     }
88     return 0;
89 }
90 
91 int
check_long_long(const char * name,int i,long long expect,long long result)92 check_long_long(const char *name, int i, long long expect, long long result)
93 {
94     if (expect != result) {
95         long long diff = expect - result;
96         printf("%s test %d got %lld expect %lld diff %lld\n", name, i, result, expect, diff);
97         return 1;
98     }
99     return 0;
100 }
101 
102 typedef const struct {
103     const char *name;
104     int (*test)(void);
105 } long_double_test_t;
106 
107 typedef const struct {
108     int line;
109     long double x;
110     long double y;
111 } long_double_test_f_f_t;
112 
113 typedef const struct {
114     int line;
115     long double x0;
116     long double x1;
117     long double y;
118 } long_double_test_f_ff_t;
119 
120 typedef const struct {
121     int line;
122     long double x0;
123     long double x1;
124     long double x2;
125     long double y;
126 } long_double_test_f_fff_t;
127 
128 typedef const struct {
129     int line;
130     long double x0;
131     int x1;
132     long double y;
133 } long_double_test_f_fi_t;
134 
135 typedef const struct {
136     int line;
137     long double x;
138     long long y;
139 } long_double_test_i_f_t;
140 
141 /*
142  * sqrtl is "exact", but can have up to one bit of error as it might
143  * not have enough guard bits to correctly perform rounding, leading
144  * to some answers getting rounded to an even value instead of the
145  * (more accurate) odd value
146  */
147 #define FMA_PREC 0
148 #if LDBL_MANT_DIG == 64
149 #define DEFAULT_PREC 0x1p-55L
150 #define SQRTL_PREC 0x1.0p-63L
151 #define FULL_LONG_DOUBLE
152 #elif LDBL_MANT_DIG == 113
153 #define FULL_LONG_DOUBLE
154 #define DEFAULT_PREC 0x1p-105L
155 #define SQRTL_PREC 0x1.0p-112L
156 #elif LDBL_MANT_DIG == 106
157 #define DEFAULT_PREC 0x1p-97L
158 #define SQRTL_PREC 0x1.0p-105L
159 #define PART_LONG_DOUBLE
160 #elif LDBL_MANT_DIG == 53
161 #define DEFAULT_PREC 0x1p-48L
162 #define SQRTL_PREC 0x1.0p-52L
163 #else
164 #error unsupported long double
165 #endif
166 
167 #define HYPOTL_PREC     SQRTL_PREC
168 #define CBRTL_PREC      SQRTL_PREC
169 #define CEILL_PREC      0
170 #define FLOORL_PREC     0
171 #define LOGBL_PREC      0
172 #define RINTL_PREC      0
173 #define FMINL_PREC      0
174 #define FMAXL_PREC      0
175 #define SCALBNL_PREC    0
176 #define SCALBL_PREC     0
177 #define LDEXPL_PREC     0
178 #define COPYSIGNL_PREC  0
179 #define NEARBYINTL_PREC 0
180 #define ROUNDL_PREC     0
181 #define TRUNCL_PREC     0
182 
183 #include "long_double_vec.h"
184 
185 #if !defined(__PICOLIBC__) || (defined(_WANT_IO_LONG_DOUBLE) && (defined(TINY_STDIO) || defined(FLOATING_POINT)))
186 #define TEST_IO_LONG_DOUBLE
187 #endif
188 
189 #if defined(__PICOLIBC__) && defined(__m68k__) && !defined(TINY_STDIO)
190 #undef TEST_IO_LONG_DOUBLE
191 #endif
192 
193 #if defined(__PICOLIBC__) && !defined(TINY_STDIO) && __LDBL_MANT_DIG__ != 64
194 #undef TEST_IO_LONG_DOUBLE
195 #endif
196 
197 #ifdef TEST_IO_LONG_DOUBLE
198 static long double vals[] = {
199     1.0L,
200     0x1.8p0L,
201     3.141592653589793238462643383279502884197169L,
202     0.0L,
203     1.0L/0.0L,
204     0.0L/0.0L,
205 };
206 
207 #define NVALS   (sizeof(vals)/sizeof(vals[0]))
208 
209 static long double
naive_strtold(const char * buf)210 naive_strtold(const char *buf)
211 {
212     long double v = 0.0L;
213     long exp = 0;
214     long double frac_mul;
215     long exp_sign = 1;
216     long double base = 10.0L;
217     char c;
218     enum {
219         LDOUBLE_INT,
220         LDOUBLE_FRAC,
221         LDOUBLE_EXP,
222     } state = LDOUBLE_INT;
223 
224     if (strncmp(buf, "0x", 2) == 0) {
225         base = 16.0L;
226         buf += 2;
227     }
228     frac_mul = 1.0L / base;
229     while ((c = *buf++)) {
230         int digit;
231         switch (c) {
232         case '.':
233             if (state == LDOUBLE_INT) {
234                 state = LDOUBLE_FRAC;
235                 continue;
236             }
237             return -(long double)INFINITY;
238         case 'p':
239             if (base == 16.0L) {
240                 if (state == LDOUBLE_INT || state == LDOUBLE_FRAC) {
241                     state = LDOUBLE_EXP;
242                     continue;
243                 }
244             }
245             return -(long double)INFINITY;
246         case '-':
247             if (state == LDOUBLE_EXP) {
248                 exp_sign = -1;
249                 continue;
250             }
251             return -(long double)INFINITY;
252         case '+':
253             if (state == LDOUBLE_EXP) {
254                 exp_sign = 1;
255                 continue;
256             }
257             return -(long double)INFINITY;
258         case '0': case '1': case '2': case '3': case '4':
259         case '5': case '6': case '7': case '8': case '9':
260             digit = c - '0';
261             break;
262         case 'E':
263             if (base == 10.0L) {
264                 if (state == LDOUBLE_INT || state == LDOUBLE_FRAC) {
265                     state = LDOUBLE_EXP;
266                     continue;
267                 }
268                 return -(long double)INFINITY;
269             }
270             /* FALLTHROUGH */
271             if (base == 10.0L) {
272                 if (state == LDOUBLE_INT || state == LDOUBLE_FRAC) {
273                     state = LDOUBLE_EXP;
274                     continue;
275                 }
276                 return -(long double)INFINITY;
277             }
278             /* FALLTHROUGH */
279         case 'A': case 'B': case 'C':
280         case 'D': case 'F':
281             if (state == LDOUBLE_INT || state == LDOUBLE_FRAC) {
282                 digit = c - 'A' + 10;
283                 break;
284             }
285             return -(long double)INFINITY;
286         case 'e':
287             if (base == 10.0L) {
288                 if (state == LDOUBLE_INT || state == LDOUBLE_FRAC) {
289                     state = LDOUBLE_EXP;
290                     continue;
291                 }
292                 return -(long double)INFINITY;
293             }
294             /* FALLTHROUGH */
295         case 'a': case 'b': case 'c':
296         case 'd': case 'f':
297             if (state == LDOUBLE_INT || state == LDOUBLE_FRAC) {
298                 digit = c - 'a' + 10;
299                 break;
300             }
301             return -(long double)INFINITY;
302         default:
303             return -(long double)INFINITY;
304         }
305         switch (state) {
306         case LDOUBLE_INT:
307             v = v * base + digit;
308             break;
309         case LDOUBLE_FRAC:
310             v = v + digit * frac_mul;
311             frac_mul *= 1.0L / base;
312             break;
313         case LDOUBLE_EXP:
314             exp = exp * 10 + digit;
315             break;
316         }
317     }
318     if (base == 10.0L) {
319         long etop = exp / 2;
320         long ebot = exp - etop;
321         long double epow_top = powl(10.0L, etop * exp_sign);
322         long double epow_bot = powl(10.0L, ebot * exp_sign);
323         long double vpow = v * epow_top;
324         long double r = vpow * epow_bot;
325         return r;
326     } else
327         return ldexpl(v, exp * exp_sign);
328 }
329 
330 static const char *formats[] = { "%La", "%.30Le", };
331 
332 #define NFMTS (sizeof (formats)/sizeof(formats[0]))
333 
334 static bool
close(long double have,long double want,long double max_error)335 close(long double have, long double want, long double max_error)
336 {
337     if (have == want)
338         return true;
339 
340     if (max_error == 0.0L)
341         return false;
342 
343     if (want == 0.0L)
344         return fabsl(have) <= max_error;
345     return fabsl((have - want) / want) <= max_error;
346 }
347 
348 static const int test_exp[] = {
349     __LDBL_MIN_EXP__ - __LDBL_MANT_DIG__ - 1,
350     __LDBL_MIN_EXP__ - __LDBL_MANT_DIG__,
351     __LDBL_MIN_EXP__ - __LDBL_MANT_DIG__ + 1,
352     __LDBL_MIN_EXP__ - __LDBL_MANT_DIG__ + 2,
353     __LDBL_MIN_EXP__ - __LDBL_MANT_DIG__ + 3,
354     __LDBL_MIN_EXP__ - 3,
355     __LDBL_MIN_EXP__ - 2,
356     __LDBL_MIN_EXP__ - 1,
357     __LDBL_MIN_EXP__,
358     __LDBL_MIN_EXP__ + 1,
359     __LDBL_MIN_EXP__ + 2,
360     __LDBL_MIN_EXP__ + 3,
361     -3,
362     -2,
363     -1,
364     0,
365     1,
366     2,
367     3,
368     __LDBL_MAX_EXP__ - 3,
369     __LDBL_MAX_EXP__ - 2,
370     __LDBL_MAX_EXP__ - 1,
371     __LDBL_MAX_EXP__,
372     __LDBL_MAX_EXP__ + 1,
373 };
374 
375 #define NEXPS (sizeof (test_exp)/ sizeof(test_exp[0]))
376 
377 /*
378  * For 64-bit values, we may have exact conversions. Otherwise, allow
379  * some error
380  */
381 #ifdef _IO_FLOAT_EXACT
382 # if __SIZEOF_LONG_DOUBLE__ == 8
383 #  define MAX_DECIMAL_ERROR       0
384 # else
385 #  define MAX_DECIMAL_ERROR     1e-10L
386 # endif
387 #else
388 # if __SIZEOF_LONG_DOUBLE__ == 8
389 #  define MAX_DECIMAL_ERROR       1e-5L
390 # else
391 #  define MAX_DECIMAL_ERROR       1e-10L
392 # endif
393 #endif
394 
395 static int
test_io(void)396 test_io(void)
397 {
398     unsigned e;
399     int result = 0;
400     char buf[80];
401     unsigned i, j;
402     long double max_error, max_error_naive;
403     char *end;
404 
405     for (e = 0; e < NEXPS; e++)
406     {
407         for (i = 0; i < NVALS; i++) {
408 
409             long double v, r;
410             v = ldexpl(vals[i], test_exp[e]);
411 
412             for (j = 0; j < NFMTS; j++) {
413 
414                 if (j == 0) {
415                     max_error = 0;
416                     max_error_naive = 0;
417                 } else {
418                     max_error = MAX_DECIMAL_ERROR;
419                     max_error_naive = 1e-6L;
420                 }
421 
422                 sprintf(buf, formats[j], v);
423                 if (isinf(v)) {
424                     if (strcmp(buf, "inf") != 0) {
425                         printf("test_io i %d val %La exp %d: is %s should be inf\n", i, vals[i], test_exp[e], buf);
426                         result++;
427                     }
428                 } else if (isnan(v)) {
429                     if (strcmp(buf, "nan") != 0) {
430                         printf("test_io is %s should be nan\n", buf);
431                         result++;
432                     }
433                 } else {
434                     r = naive_strtold(buf);
435                     if (!close(r, v, max_error_naive)) {
436                         printf("test_io naive i %d val %La exp %d: \"%s\", is %La should be %La\n", i, vals[i], test_exp[e], buf, r, v);
437                         result++;
438                     }
439                 }
440                 sscanf(buf, "%Lf", &r);
441                 if (!close(r, v, max_error) && !(isnan(v) && isnan(r))) {
442                     printf("test_io scanf i %d val %La exp %d: \"%s\", is %La should be %La\n", i, vals[i], test_exp[e], buf, r, v);
443                     result++;
444                 }
445                 r = strtold(buf, &end);
446                 if ((!close(r, v, max_error) && !(isnan(v) && isnan(r)))|| end != buf + strlen(buf)) {
447                     printf("test_io strtold i %d val %La exp %d: \"%s\", is %La should be %La\n", i, vals[i], test_exp[e], buf, r, v);
448                     result++;
449                 }
450             }
451         }
452     }
453     return result;
454 }
455 #endif
456 
main(void)457 int main(void)
458 {
459     int result = 0;
460     unsigned int i;
461 
462     printf("LDBL_MANT_DIG %d\n", LDBL_MANT_DIG);
463 #ifdef __m68k__
464     volatile long double zero = 0.0L;
465     volatile long double one = 1.0L;
466     volatile long double check = nextafterl(zero, one);
467     if (check + check == zero) {
468         printf("m68k emulating long double with double, skipping\n");
469         return 77;
470     }
471 #endif
472 #ifdef TEST_IO_LONG_DOUBLE
473     result += test_io();
474 #endif
475     for (i = 0; i < sizeof(long_double_tests) / sizeof(long_double_tests[0]); i++) {
476         result += long_double_tests[i].test();
477     }
478     return result != 0;
479 }
480 
481 #else
main(void)482 int main(void)
483 {
484     printf("no long double support\n");
485     return 0;
486 }
487 #endif
488