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 /*
45  * fall-through case statement annotations
46  */
47 #if __cplusplus >= 201703L || __STDC_VERSION__ > 201710L
48 /* Standard C++17/C23 attribute */
49 #define __TEST_PICOLIBC_FALLTHROUGH [[fallthrough]]
50 #elif __has_attribute(fallthrough)
51 /* Non-standard but supported by at least gcc and clang */
52 #define __TEST_PICOLIBC_FALLTHROUGH __attribute__((fallthrough))
53 #else
54 #define __TEST_PICOLIBC_FALLTHROUGH do { } while(0)
55 #endif
56 
57 #ifdef _TEST_LONG_DOUBLE
58 
59 static long double max_error;
60 
61 static bool
within_error(long double expect,long double result,long double error)62 within_error(long double expect, long double result, long double error)
63 {
64     long double difference;
65     long double e = 1.0L;
66 
67     if (isnan(expect) && isnan(result))
68         return true;
69 
70     if (expect == result)
71         return true;
72 
73     if (expect != 0)
74         e = scalbnl(1.0L, -ilogbl(expect));
75 
76     difference = fabsl(expect - result) * e;
77 
78     if (difference > max_error)
79         max_error = difference;
80 
81     return difference <= error;
82 }
83 
84 static int
check_long_double(const char * name,int i,long double prec,long double expect,long double result)85 check_long_double(const char *name, int i, long double prec, long double expect, long double result)
86 {
87     if (!within_error(expect, result, prec)) {
88         long double diff = fabsl(expect - result);
89 #ifdef __PICOLIBC__
90 #ifdef _WANT_IO_LONG_DOUBLE
91         printf("%s test %d got %La expect %La diff %La\n", name, i, result, expect, diff);
92 #else
93         printf("%s test %d got %a expect %a diff %a\n", name, i, (double) result, (double) expect, (double) diff);
94 #endif
95 #else
96 //        printf("%s test %d got %.33Lg expect %.33Lg diff %.33Lg\n", name, i, result, expect, diff);
97         printf("%s test %d got %La expect %La diff %La\n", name, i, result, expect, diff);
98 #endif
99         return 1;
100     }
101     return 0;
102 }
103 
104 static int
check_long_long(const char * name,int i,long long expect,long long result)105 check_long_long(const char *name, int i, long long expect, long long result)
106 {
107     if (expect != result) {
108         long long diff = expect - result;
109         printf("%s test %d got %lld expect %lld diff %lld\n", name, i, result, expect, diff);
110         return 1;
111     }
112     return 0;
113 }
114 
115 typedef const struct {
116     const char *name;
117     int (*test)(void);
118 } long_double_test_t;
119 
120 typedef const struct {
121     int line;
122     long double x;
123     long double y;
124 } long_double_test_f_f_t;
125 
126 typedef const struct {
127     int line;
128     long double x0;
129     long double x1;
130     long double y;
131 } long_double_test_f_ff_t;
132 
133 typedef const struct {
134     int line;
135     long double x0;
136     long double x1;
137     long double x2;
138     long double y;
139 } long_double_test_f_fff_t;
140 
141 typedef const struct {
142     int line;
143     long double x0;
144     int x1;
145     long double y;
146 } long_double_test_f_fi_t;
147 
148 typedef const struct {
149     int line;
150     long double x;
151     long long y;
152 } long_double_test_i_f_t;
153 
154 /*
155  * sqrtl is "exact", but can have up to one bit of error as it might
156  * not have enough guard bits to correctly perform rounding, leading
157  * to some answers getting rounded to an even value instead of the
158  * (more accurate) odd value
159  */
160 #define FMA_PREC 0
161 #if LDBL_MANT_DIG == 64
162 #define DEFAULT_PREC 0x1p-55L
163 #define SQRTL_PREC 0x1.0p-63L
164 #define FULL_LONG_DOUBLE
165 #elif LDBL_MANT_DIG == 113
166 #define FULL_LONG_DOUBLE
167 #define DEFAULT_PREC 0x1p-105L
168 #define SQRTL_PREC 0x1.0p-112L
169 #elif LDBL_MANT_DIG == 106
170 #define DEFAULT_PREC 0x1p-97L
171 #define SQRTL_PREC 0x1.0p-105L
172 #define PART_LONG_DOUBLE
173 #elif LDBL_MANT_DIG == 53
174 #define DEFAULT_PREC 0x1p-48L
175 #define SQRTL_PREC 0x1.0p-52L
176 #else
177 #error unsupported long double
178 #endif
179 
180 #define HYPOTL_PREC     SQRTL_PREC
181 #define CBRTL_PREC      SQRTL_PREC
182 #define CEILL_PREC      0
183 #define FLOORL_PREC     0
184 #define LOGBL_PREC      0
185 #define RINTL_PREC      0
186 #define FMINL_PREC      0
187 #define FMAXL_PREC      0
188 #define SCALBNL_PREC    0
189 #define SCALBL_PREC     0
190 #define LDEXPL_PREC     0
191 #define COPYSIGNL_PREC  0
192 #define NEARBYINTL_PREC 0
193 #define ROUNDL_PREC     0
194 #define TRUNCL_PREC     0
195 
196 #include "long_double_vec.h"
197 
198 #if !defined(__PICOLIBC__) || (defined(_WANT_IO_LONG_DOUBLE) && (defined(TINY_STDIO) || defined(FLOATING_POINT)))
199 #define TEST_IO_LONG_DOUBLE
200 #endif
201 
202 #if defined(__PICOLIBC__) && defined(__m68k__) && !defined(TINY_STDIO)
203 #undef TEST_IO_LONG_DOUBLE
204 #endif
205 
206 #if defined(__PICOLIBC__) && !defined(TINY_STDIO) && __LDBL_MANT_DIG__ != 64
207 #undef TEST_IO_LONG_DOUBLE
208 #endif
209 
210 #ifdef TEST_IO_LONG_DOUBLE
211 static long double vals[] = {
212     1.0L,
213     0x1.8p0L,
214     3.141592653589793238462643383279502884197169L,
215     0.0L,
216     1.0L/0.0L,
217     0.0L/0.0L,
218 };
219 
220 #define NVALS   (sizeof(vals)/sizeof(vals[0]))
221 
222 static long double
naive_strtold(const char * buf)223 naive_strtold(const char *buf)
224 {
225     long double v = 0.0L;
226     long exp = 0;
227     long double frac_mul;
228     long exp_sign = 1;
229     long double base = 10.0L;
230     char c;
231     enum {
232         LDOUBLE_INT,
233         LDOUBLE_FRAC,
234         LDOUBLE_EXP,
235     } state = LDOUBLE_INT;
236 
237     if (strncmp(buf, "0x", 2) == 0) {
238         base = 16.0L;
239         buf += 2;
240     }
241     frac_mul = 1.0L / base;
242     while ((c = *buf++)) {
243         int digit;
244         switch (c) {
245         case '.':
246             if (state == LDOUBLE_INT) {
247                 state = LDOUBLE_FRAC;
248                 continue;
249             }
250             return -(long double)INFINITY;
251         case 'p':
252             if (base == 16.0L) {
253                 if (state == LDOUBLE_INT || state == LDOUBLE_FRAC) {
254                     state = LDOUBLE_EXP;
255                     continue;
256                 }
257             }
258             return -(long double)INFINITY;
259         case '-':
260             if (state == LDOUBLE_EXP) {
261                 exp_sign = -1;
262                 continue;
263             }
264             return -(long double)INFINITY;
265         case '+':
266             if (state == LDOUBLE_EXP) {
267                 exp_sign = 1;
268                 continue;
269             }
270             return -(long double)INFINITY;
271         case '0': case '1': case '2': case '3': case '4':
272         case '5': case '6': case '7': case '8': case '9':
273             digit = c - '0';
274             break;
275         case 'E':
276             if (base == 10.0L) {
277                 if (state == LDOUBLE_INT || state == LDOUBLE_FRAC) {
278                     state = LDOUBLE_EXP;
279                     continue;
280                 }
281                 return -(long double)INFINITY;
282             }
283             if (base == 10.0L) {
284                 if (state == LDOUBLE_INT || state == LDOUBLE_FRAC) {
285                     state = LDOUBLE_EXP;
286                     continue;
287                 }
288                 return -(long double)INFINITY;
289             }
290             __TEST_PICOLIBC_FALLTHROUGH;
291         case 'A': case 'B': case 'C':
292         case 'D': case 'F':
293             if (state == LDOUBLE_INT || state == LDOUBLE_FRAC) {
294                 digit = c - 'A' + 10;
295                 break;
296             }
297             return -(long double)INFINITY;
298         case 'e':
299             if (base == 10.0L) {
300                 if (state == LDOUBLE_INT || state == LDOUBLE_FRAC) {
301                     state = LDOUBLE_EXP;
302                     continue;
303                 }
304                 return -(long double)INFINITY;
305             }
306             __TEST_PICOLIBC_FALLTHROUGH;
307         case 'a': case 'b': case 'c':
308         case 'd': case 'f':
309             if (state == LDOUBLE_INT || state == LDOUBLE_FRAC) {
310                 digit = c - 'a' + 10;
311                 break;
312             }
313             return -(long double)INFINITY;
314         default:
315             return -(long double)INFINITY;
316         }
317         switch (state) {
318         case LDOUBLE_INT:
319             v = v * base + digit;
320             break;
321         case LDOUBLE_FRAC:
322             v = v + digit * frac_mul;
323             frac_mul *= 1.0L / base;
324             break;
325         case LDOUBLE_EXP:
326             exp = exp * 10 + digit;
327             break;
328         }
329     }
330     if (base == 10.0L) {
331         long etop = exp / 2;
332         long ebot = exp - etop;
333         long double epow_top = powl(10.0L, etop * exp_sign);
334         long double epow_bot = powl(10.0L, ebot * exp_sign);
335         long double vpow = v * epow_top;
336         long double r = vpow * epow_bot;
337         return r;
338     } else
339         return ldexpl(v, exp * exp_sign);
340 }
341 
342 static const char *formats[] = { "%La", "%.30Le", };
343 
344 #define NFMTS (sizeof (formats)/sizeof(formats[0]))
345 
346 static bool
close(long double have,long double want,long double max_error)347 close(long double have, long double want, long double max_error)
348 {
349     if (have == want)
350         return true;
351 
352     if (max_error == 0.0L)
353         return false;
354 
355     if (want == 0.0L)
356         return fabsl(have) <= max_error;
357     return fabsl((have - want) / want) <= max_error;
358 }
359 
360 static const int test_exp[] = {
361     __LDBL_MIN_EXP__ - __LDBL_MANT_DIG__ - 1,
362     __LDBL_MIN_EXP__ - __LDBL_MANT_DIG__,
363     __LDBL_MIN_EXP__ - __LDBL_MANT_DIG__ + 1,
364     __LDBL_MIN_EXP__ - __LDBL_MANT_DIG__ + 2,
365     __LDBL_MIN_EXP__ - __LDBL_MANT_DIG__ + 3,
366     __LDBL_MIN_EXP__ - 3,
367     __LDBL_MIN_EXP__ - 2,
368     __LDBL_MIN_EXP__ - 1,
369     __LDBL_MIN_EXP__,
370     __LDBL_MIN_EXP__ + 1,
371     __LDBL_MIN_EXP__ + 2,
372     __LDBL_MIN_EXP__ + 3,
373     -3,
374     -2,
375     -1,
376     0,
377     1,
378     2,
379     3,
380     __LDBL_MAX_EXP__ - 3,
381     __LDBL_MAX_EXP__ - 2,
382     __LDBL_MAX_EXP__ - 1,
383     __LDBL_MAX_EXP__,
384     __LDBL_MAX_EXP__ + 1,
385 };
386 
387 #define NEXPS (sizeof (test_exp)/ sizeof(test_exp[0]))
388 
389 /*
390  * For 64-bit values, we may have exact conversions. Otherwise, allow
391  * some error
392  */
393 #ifdef _IO_FLOAT_EXACT
394 # if __SIZEOF_LONG_DOUBLE__ == 8
395 #  define MAX_DECIMAL_ERROR       0
396 # else
397 #  define MAX_DECIMAL_ERROR     1e-10L
398 # endif
399 #else
400 # if __SIZEOF_LONG_DOUBLE__ == 8
401 #  define MAX_DECIMAL_ERROR       1e-5L
402 # else
403 #  define MAX_DECIMAL_ERROR       1e-10L
404 # endif
405 #endif
406 
407 static int
test_io(void)408 test_io(void)
409 {
410     unsigned e;
411     int result = 0;
412     char buf[80];
413     unsigned i, j;
414     long double max_error, max_error_naive;
415     char *end;
416 
417     for (e = 0; e < NEXPS; e++)
418     {
419         for (i = 0; i < NVALS; i++) {
420 
421             long double v, r;
422             v = ldexpl(vals[i], test_exp[e]);
423 
424             for (j = 0; j < NFMTS; j++) {
425 
426                 if (j == 0) {
427                     max_error = 0;
428                     max_error_naive = 0;
429                 } else {
430                     max_error = MAX_DECIMAL_ERROR;
431                     max_error_naive = 1e-6L;
432                 }
433 
434                 sprintf(buf, formats[j], v);
435                 if (isinf(v)) {
436                     if (strcmp(buf, "inf") != 0) {
437                         printf("test_io i %d val %La exp %d: is %s should be inf\n", i, vals[i], test_exp[e], buf);
438                         result++;
439                     }
440                 } else if (isnan(v)) {
441                     if (strcmp(buf, "nan") != 0) {
442                         printf("test_io is %s should be nan\n", buf);
443                         result++;
444                     }
445                 } else {
446                     r = naive_strtold(buf);
447                     if (!close(r, v, max_error_naive)) {
448                         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);
449                         result++;
450                     }
451                 }
452                 sscanf(buf, "%Lf", &r);
453                 if (!close(r, v, max_error) && !(isnan(v) && isnan(r))) {
454                     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);
455                     result++;
456                 }
457                 r = strtold(buf, &end);
458                 if ((!close(r, v, max_error) && !(isnan(v) && isnan(r)))|| end != buf + strlen(buf)) {
459                     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);
460                     result++;
461                 }
462             }
463         }
464     }
465     return result;
466 }
467 #endif
468 
main(void)469 int main(void)
470 {
471     int result = 0;
472     unsigned int i;
473 
474     printf("LDBL_MANT_DIG %d\n", LDBL_MANT_DIG);
475 #ifdef __m68k__
476     volatile long double zero = 0.0L;
477     volatile long double one = 1.0L;
478     volatile long double check = nextafterl(zero, one);
479     if (check + check == zero) {
480         printf("m68k emulating long double with double, skipping\n");
481         return 77;
482     }
483 #endif
484 #ifdef TEST_IO_LONG_DOUBLE
485     result += test_io();
486 #endif
487     for (i = 0; i < sizeof(long_double_tests) / sizeof(long_double_tests[0]); i++) {
488         result += long_double_tests[i].test();
489     }
490     return result != 0;
491 }
492 
493 #else
main(void)494 int main(void)
495 {
496     printf("no long double support\n");
497     return 0;
498 }
499 #endif
500