1 /*
2  * SPDX-License-Identifier: BSD-3-Clause
3  *
4  * Copyright © 2023 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 __STDC_WANT_IEC_60559_BFP_EXT__
37 #include <math.h>
38 #include <stdio.h>
39 #include <float.h>
40 #include <stdlib.h>
41 #include <stdbool.h>
42 #include <fenv.h>
43 
44 #if FLT_MANT_DIG == 24
45 
46 static int roundings[] = {
47 #ifdef FE_TONEAREST
48     FE_TONEAREST,
49 #else
50     -1,
51 #endif
52 #ifdef FE_UPWARD
53     FE_UPWARD,
54 #endif
55 #ifdef FE_DOWNWARD
56     FE_DOWNWARD,
57 #endif
58 #ifdef FE_TOWARDZERO
59     FE_TOWARDZERO,
60 #endif
61 };
62 
63 static int rounding_map[] = {
64 #ifdef FE_TONEAREST
65     0,
66 #endif
67 #ifdef FE_UPWARD
68     1,
69 #endif
70 #ifdef FE_DOWNWARD
71     2,
72 #endif
73 #ifdef FE_TOWARDZERO
74     3,
75 #endif
76 };
77 
78 static const char *rounding_names[] = {
79 #ifdef FE_TONEAREST
80     "FE_TONEAREST",
81 #else
82     "default",
83 #endif
84 #ifdef FE_UPWARD
85     "FE_UPWARD",
86 #endif
87 #ifdef FE_DOWNWARD
88     "FE_DOWNWARD",
89 #endif
90 #ifdef FE_TOWARDZERO
91     "FE_TOWARDZERO",
92 #endif
93 };
94 
95 #define NROUND  (sizeof(roundings)/sizeof(roundings[0]))
96 
97 #ifdef PICOLIBC_DOUBLE_NOEXCEPT
98 /*
99  * Assume that a lack of exceptions for doubles also means a lack of
100  * support for non-default rounding modes for doubles.
101  */
102 #define SKIP_ROUNDINGS_FLOAT
103 #endif
104 
105 #ifdef SKIP_ROUNDINGS_FLOAT
106 #define NROUND_FLOAT    1
107 #else
108 #define NROUND_FLOAT NROUND
109 #endif
110 
111 #define nanf __builtin_nanf("")
112 #define inff INFINITY
113 #define nan __builtin_nan("")
114 #define inf ((double) INFINITY)
115 #define nanl __builtin_nanl("")
116 #define infl ((long double) INFINITY)
117 
118 struct fmaf_vec {
119     float       x, y, z;
120     float       r[4];
121 };
122 
123 struct fma_vec {
124     double      x, y, z;
125     double      r[4];
126 };
127 
128 #if LDBL_MANT_DIG != 0
129 struct fmal_vec {
130     long double      x, y, z;
131     long double      r[4];
132 };
133 #endif
134 
135 #include "fma_vec.h"
136 
137 #define NUM_FMAF_VEC (sizeof(fmaf_vec)/sizeof(fmaf_vec[0]))
138 #define NUM_FMA_VEC (sizeof(fma_vec)/sizeof(fma_vec[0]))
139 #define NUM_FMAL_VEC (sizeof(fmal_vec)/sizeof(fmal_vec[0]))
140 
141 #endif
142 
143 #ifdef HAVE_FLOAT_FMA_VEC
144 
145 static bool
equalf(float x,float y)146 equalf(float x, float y)
147 {
148     if (isnan(x) && isnan(y))
149         return true;
150     return x == y;
151 }
152 
153 #if defined(__PICOLIBC__) && !defined(TINY_STDIO)
154 static
strfromf(char * str,size_t n,const char * format,float f)155 int strfromf(char *str, size_t n, const char *format, float f)
156 {
157     return snprintf(str, n, format, (double) f);
158 }
159 #endif
160 
161 static int
test_fmaf(void)162 test_fmaf(void)
163 {
164     char        xs[20], ys[20], zs[20], rs[20], ss[20];
165     int ret = 0;
166     unsigned int ro;
167     int defround = fegetround();
168     unsigned int t;
169 
170     printf("float\n");
171     for (ro = 0; ro < NROUND_FLOAT; ro++) {
172         if (roundings[ro] == -1)
173             fesetround(defround);
174         else
175             fesetround(roundings[ro]);
176         for (t = 0; t < NUM_FMAF_VEC; t++) {
177             float x = fmaf_vec[t].x;
178             float y = fmaf_vec[t].y;
179             float z = fmaf_vec[t].z;
180             volatile float r = fmaf(x, y, z);
181             float s = fmaf_vec[t].r[rounding_map[ro]];
182             if (!equalf(r, s)) {
183                 strfromf(xs, sizeof(xs), "%a", x);
184                 strfromf(ys, sizeof(xs), "%a", y);
185                 strfromf(zs, sizeof(xs), "%a", z);
186                 strfromf(rs, sizeof(xs), "%a", r);
187                 strfromf(ss, sizeof(xs), "%a", s);
188                 printf("%u: round %s %s %s %s -> got %s want %s\n",
189                        t, rounding_names[ro],
190                        xs, ys, zs, rs, ss);
191                 ret = 1;
192             }
193         }
194     }
195     fesetround(defround);
196     return ret;
197 }
198 #else
199 #define test_fmaf() 0
200 #endif
201 
202 #ifdef HAVE_DOUBLE_FMA_VEC
203 
204 #ifdef PICOLIBC_DOUBLE_NOEXCEPT
205 /*
206  * Assume that a lack of exceptions for doubles also means a lack of
207  * support for non-default rounding modes for doubles.
208  */
209 #define SKIP_ROUNDINGS_DOUBLE
210 #endif
211 
212 #ifdef SKIP_ROUNDINGS_DOUBLE
213 #define NROUND_DOUBLE    1
214 #else
215 #define NROUND_DOUBLE NROUND
216 #endif
217 
218 static bool
equal(double x,double y)219 equal(double x, double y)
220 {
221     if (isnan(x) && isnan(y))
222         return true;
223     return x == y;
224 }
225 
226 static int
test_fma(void)227 test_fma(void)
228 {
229     int ret = 0;
230     unsigned int ro;
231     int defround = fegetround();
232     unsigned int t;
233 
234     printf("double\n");
235     for (ro = 0; ro < NROUND_DOUBLE; ro++) {
236         if (roundings[ro] == -1)
237             fesetround(defround);
238         else
239             fesetround(roundings[ro]);
240         for (t = 0; t < NUM_FMA_VEC; t++) {
241             double x = fma_vec[t].x;
242             double y = fma_vec[t].y;
243             double z = fma_vec[t].z;
244             volatile double r = fma(x, y, z);
245             double s = fma_vec[t].r[rounding_map[ro]];
246             if (!equal(r, s)) {
247                 printf("%u: round %s %a * %a + %a -> got %a want %a\n",
248                        t, rounding_names[ro],
249                        x, y, z, r, s);
250                 ret = 1;
251             }
252         }
253     }
254     fesetround(defround);
255     return ret;
256 }
257 #else
258 #define test_fma() 0
259 #endif
260 
261 #if defined (_TEST_LONG_DOUBLE) && defined(HAVE_LONG_DOUBLE_FMA_VEC)
262 
263 #ifdef PICOLIBC_LONG_DOUBLE_NOEXCEPT
264 /*
265  * Assume that a lack of exceptions for doubles also means a lack of
266  * support for non-default rounding modes for doubles.
267  */
268 #define SKIP_ROUNDINGS_LONG_DOUBLE
269 #endif
270 
271 #ifdef SKIP_ROUNDINGS_LONG_DOUBLE
272 #define NROUND_LONG_DOUBLE    1
273 #else
274 #define NROUND_LONG_DOUBLE NROUND
275 #endif
276 
277 static bool
equall(long double x,long double y)278 equall(long double x, long double y)
279 {
280     if (isnan(x) && isnan(y))
281         return true;
282     return x == y;
283 }
284 
285 static int
test_fmal(void)286 test_fmal(void)
287 {
288     int ret = 0;
289     unsigned int ro;
290     int defround = fegetround();
291     unsigned int t;
292 
293     printf("long double\n");
294     for (ro = 0; ro < NROUND_LONG_DOUBLE; ro++) {
295         if (roundings[ro] == -1)
296             fesetround(defround);
297         else
298             fesetround(roundings[ro]);
299         for (t = 0; t < NUM_FMAL_VEC; t++) {
300             if (t < 30)
301                 continue;
302             long double x = fmal_vec[t].x;
303             long double y = fmal_vec[t].y;
304             long double z = fmal_vec[t].z;
305             long double r = fmal(x, y, z);
306             long double s = fmal_vec[t].r[rounding_map[ro]];
307             if (!equall(r, s)) {
308                 printf("%u: round %s %La * %La + %La -> got %La want %La\n",
309                        t, rounding_names[ro],
310                        x, y, z, r, s);
311                 ret = 1;
312             }
313         }
314     }
315     fesetround(defround);
316     return ret;
317 }
318 #else
319 #define test_fmal() 0
320 #endif
321 
322 int
main(void)323 main(void)
324 {
325 #ifdef __arc__
326     volatile float x = 0x1.000002p-2f;
327     volatile float y = 0x1.000002p-126f;
328     volatile float z = 0x0.400002p-126f;
329     if (x * y != z) {
330         printf("ARC soft float bug, skipping FMA tests\n");
331         return 77;
332     }
333 #endif
334     (void) rounding_names;
335     (void) roundings;
336     int ret = 0;
337     ret |= test_fmaf();
338     ret |= test_fma();
339 #if defined(__m68k__) && __LDBL_MIN_EXP__ == -16382
340     volatile long double big = 0x1p+16383l;
341     volatile long double small = 0x1p-16446l;
342     if (big * small != 0x1p-63) {
343         printf("m68k qemu long double denorm bug, skipping long double tests\n");
344         return ret;
345     }
346 #endif
347     ret |= test_fmal();
348     return ret;
349 }
350