1 | /* Single-precision e^x function. |
---|---|

2 | Copyright (C) 2017-2019 Free Software Foundation, Inc. |

3 | This file is part of the GNU C Library. |

4 | |

5 | The GNU C Library is free software; you can redistribute it and/or |

6 | modify it under the terms of the GNU Lesser General Public |

7 | License as published by the Free Software Foundation; either |

8 | version 2.1 of the License, or (at your option) any later version. |

9 | |

10 | The GNU C Library is distributed in the hope that it will be useful, |

11 | but WITHOUT ANY WARRANTY; without even the implied warranty of |

12 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |

13 | Lesser General Public License for more details. |

14 | |

15 | You should have received a copy of the GNU Lesser General Public |

16 | License along with the GNU C Library; if not, see |

17 | <http://www.gnu.org/licenses/>. */ |

18 | |

19 | #ifdef __expf |

20 | # undef libm_hidden_proto |

21 | # define libm_hidden_proto(ignored) |

22 | #endif |

23 | |

24 | #include <math.h> |

25 | #include <math-narrow-eval.h> |

26 | #include <stdint.h> |

27 | #include <shlib-compat.h> |

28 | #include <libm-alias-float.h> |

29 | #include "math_config.h" |

30 | |

31 | /* |

32 | EXP2F_TABLE_BITS = 5 |

33 | EXP2F_POLY_ORDER = 3 |

34 | |

35 | ULP error: 0.502 (nearest rounding.) |

36 | Relative error: 1.69 * 2^-34 in [-ln2/64, ln2/64] (before rounding.) |

37 | Wrong count: 170635 (all nearest rounding wrong results with fma.) |

38 | Non-nearest ULP error: 1 (rounded ULP error) |

39 | */ |

40 | |

41 | #define N (1 << EXP2F_TABLE_BITS) |

42 | #define InvLn2N __exp2f_data.invln2_scaled |

43 | #define T __exp2f_data.tab |

44 | #define C __exp2f_data.poly_scaled |

45 | |

46 | static inline uint32_t |

47 | top12 (float x) |

48 | { |

49 | return asuint (x) >> 20; |

50 | } |

51 | |

52 | float |

53 | __expf (float x) |

54 | { |

55 | uint32_t abstop; |

56 | uint64_t ki, t; |

57 | /* double_t for better performance on targets with FLT_EVAL_METHOD==2. */ |

58 | double_t kd, xd, z, r, r2, y, s; |

59 | |

60 | xd = (double_t) x; |

61 | abstop = top12 (x) & 0x7ff; |

62 | if (__glibc_unlikely (abstop >= top12 (88.0f))) |

63 | { |

64 | /* |x| >= 88 or x is nan. */ |

65 | if (asuint (x) == asuint (-INFINITY)) |

66 | return 0.0f; |

67 | if (abstop >= top12 (INFINITY)) |

68 | return x + x; |

69 | if (x > 0x1.62e42ep6f) /* x > log(0x1p128) ~= 88.72 */ |

70 | return __math_oflowf (0); |

71 | if (x < -0x1.9fe368p6f) /* x < log(0x1p-150) ~= -103.97 */ |

72 | return __math_uflowf (0); |

73 | #if WANT_ERRNO_UFLOW |

74 | if (x < -0x1.9d1d9ep6f) /* x < log(0x1p-149) ~= -103.28 */ |

75 | return __math_may_uflowf (0); |

76 | #endif |

77 | } |

78 | |

79 | /* x*N/Ln2 = k + r with r in [-1/2, 1/2] and int k. */ |

80 | z = InvLn2N * xd; |

81 | |

82 | /* Round and convert z to int, the result is in [-150*N, 128*N] and |

83 | ideally ties-to-even rule is used, otherwise the magnitude of r |

84 | can be bigger which gives larger approximation error. */ |

85 | #if TOINT_INTRINSICS |

86 | kd = roundtoint (z); |

87 | ki = converttoint (z); |

88 | #else |

89 | # define SHIFT __exp2f_data.shift |

90 | kd = math_narrow_eval ((double) (z + SHIFT)); /* Needs to be double. */ |

91 | ki = asuint64 (kd); |

92 | kd -= SHIFT; |

93 | #endif |

94 | r = z - kd; |

95 | |

96 | /* exp(x) = 2^(k/N) * 2^(r/N) ~= s * (C0*r^3 + C1*r^2 + C2*r + 1) */ |

97 | t = T[ki % N]; |

98 | t += ki << (52 - EXP2F_TABLE_BITS); |

99 | s = asdouble (t); |

100 | z = C[0] * r + C[1]; |

101 | r2 = r * r; |

102 | y = C[2] * r + 1; |

103 | y = z * r2 + y; |

104 | y = y * s; |

105 | return (float) y; |

106 | } |

107 | |

108 | #ifndef __expf |

109 | hidden_def (__expf) |

110 | strong_alias (__expf, __ieee754_expf) |

111 | strong_alias (__expf, __expf_finite) |

112 | versioned_symbol (libm, __expf, expf, GLIBC_2_27); |

113 | libm_alias_float_other (__exp, exp) |

114 | #endif |

115 |