1 /*
2 Copyright (C) 2002 by Red Hat, Incorporated. All rights reserved.
3
4 Permission to use, copy, modify, and distribute this software
5 is freely granted, provided that this notice is preserved.
6 */
7 /*
8 FUNCTION
9 <<fma>>, <<fmaf>>---floating multiply add
10 INDEX
11 fma
12 INDEX
13 fmaf
14
15 SYNOPSIS
16 #include <math.h>
17 double fma(double <[x]>, double <[y]>, double <[z]>);
18 float fmaf(float <[x]>, float <[y]>, float <[z]>);
19
20 DESCRIPTION
21 The <<fma>> functions compute (<[x]> * <[y]>) + <[z]>, rounded as one ternary
22 operation: they compute the value (as if) to infinite precision and round once
23 to the result format, according to the rounding mode characterized by the value
24 of FLT_ROUNDS. That is, they are supposed to do this: see below.
25
26 RETURNS
27 The <<fma>> functions return (<[x]> * <[y]>) + <[z]>, rounded as one ternary
28 operation.
29
30 BUGS
31 This implementation does not provide the function that it should, purely
32 returning "(<[x]> * <[y]>) + <[z]>;" with no attempt at all to provide the
33 simulated infinite precision intermediates which are required. DO NOT USE THEM.
34
35 If double has enough more precision than float, then <<fmaf>> should provide
36 the expected numeric results, as it does use double for the calculation. But
37 since this is not the case for all platforms, this manual cannot determine
38 if it is so for your case.
39
40 PORTABILITY
41 ANSI C, POSIX.
42
43 */
44
45 #include "fdlibm.h"
46
47 #if !_HAVE_FAST_FMA
48
49 #ifdef _NEED_FLOAT64
50
51 #if __FLT_EVAL_METHOD__ == 2 && defined(_HAVE_LONG_DOUBLE)
52
53 __float64
fma64(__float64 x,__float64 y,__float64 z)54 fma64(__float64 x, __float64 y, __float64 z)
55 {
56 return (__float64) fmal((long double) x, (long double) y, (long double) z);
57 }
58
59 #else
60
61 typedef __float64 FLOAT_T;
62
63 #define FMA fma64
64 #define NEXTAFTER nextafter64
65 #define LDEXP ldexp64
66 #define FREXP frexp64
67 #define SCALBN scalbn64
68 #define ILOGB ilogb64
69 #define COPYSIGN copysign64
70
71 #define SPLIT ((FLOAT_T) 0x1p26 + (FLOAT_T) 1.0)
72 #define FLOAT_MANT_DIG _FLOAT64_MANT_DIG
73 #define FLOAT_MAX_EXP _FLOAT64_MAX_EXP
74 #define FLOAT_MIN _FLOAT64_MIN
75
76 static inline int
odd_mant(FLOAT_T x)77 odd_mant(FLOAT_T x)
78 {
79 return asuint64(x) & 1;
80 }
81
82 static unsigned int
EXPONENT(FLOAT_T x)83 EXPONENT(FLOAT_T x)
84 {
85 return _exponent64(asuint64(x));
86 }
87
88 #ifdef PICOLIBC_FLOAT64_NOEXCEPT
89 #define feraiseexcept(x) ((void) (x))
90 #endif
91
92 #include "fma_inc.h"
93
94 #endif
95
96 _MATH_ALIAS_d_ddd(fma)
97
98 #endif /* _NEED_FLOAT64 */
99
100 #endif /* !_HAVE_FAST_FMA */
101
