ftz 0.0.1
Fast, reproducible floating-point arithmetic
Loading...
Searching...
No Matches
tanh.h
Go to the documentation of this file.
1#pragma once
2#include "ftz/config.h"
3#include "native/attributes.h"
4#include "ftz/math/float.h"
5
6namespace ftz { namespace detail { namespace math {
7 // Fixed odd graph x*P(x*x). Every active stage is normal or exact zero:
8 // the tiny bypass prevents a square near underflow, and retained interval
9 // bounds keep every Horner stage away from zero. Neither FTZ policy relies
10 // on hardware's ambiguous minimum-normal rounding strip. RNE and real fused
11 // FMA are required; no implicit contraction/reassociation is admitted.
12 template <bool Hardware = FTZ_FP32_HARDWARE_FTZ != 0>
13 [[nodiscard]] native_constexpr native_inline native_const unsigned int tanh_words(unsigned int word) {
14 unsigned int magnitude = word & 0x7fffffffu;
15 unsigned int sign = word & 0x80000000u;
16 if (magnitude > 0x7f800000u) {
17 return 0x7fc00000u;
18 }
19 if (magnitude < 0x00800000u) {
20 return sign;
21 }
22 // tanh(x) rounds to x throughout this small interval. Preserve raw -0
23 // above and normal bits here without evaluating an underflowing product.
24 if (magnitude <= 0x39800000u) {
25 return word;
26 }
27 // tanh(10) differs from one by less than half a binary32 ULP below one.
28 // This also handles either infinity before any arithmetic.
29 if (magnitude >= 0x41200000u) {
30 return sign | 0x3f800000u;
31 }
32 float x = fp32_decode(magnitude);
33 float z = fp32_mul<true,Hardware>(x, x);
34 float t, h;
35 if (magnitude <= 0x3f800000u) {
36 t = fp32_fma<true,Hardware>(z, fp32_decode(0x40000000u), fp32_decode(0xbf800000u));
37 h = fp32_decode(0x34facb37u);
38 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0xb63a0d2du));
39 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0x37813497u));
40 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0xb8bfb3f8u));
41 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0x3a0e6d24u));
42 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0xbb535f6cu));
43 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0x3c9d20e4u));
44 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0xbded544du));
45 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0x3f5c6e3eu));
46 }
47 else if (magnitude <= 0x40000000u) {
48 t = fp32_fma<true,Hardware>(z, fp32_decode(0x3f000000u), fp32_decode(0xbfa00000u));
49 h = fp32_decode(0x38752140u);
50 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0xb9183513u));
51 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0x398d9ee0u));
52 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0xba2fdf3au));
53 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0x3ae1130du));
54 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0xbb8bc302u));
55 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0x3c2d6773u));
56 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0xbcd7a178u));
57 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0x3d86d45du));
58 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0xbe2e2df9u));
59 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0x3f14c222u));
60 }
61 else if (magnitude <= 0x40400000u) {
62 t = fp32_fma<true,Hardware>(z, fp32_decode(0x3e800000u), fp32_decode(0xbfd00000u));
63 h = fp32_decode(0xb947bd66u);
64 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0x39dfe5a4u));
65 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0xba4a3864u));
66 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0x3ae2b84bu));
67 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0xbb813c08u));
68 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0x3c1129ceu));
69 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0xbca3c0b8u));
70 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0x3d3bcdf4u));
71 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0xbde4f8bcu));
72 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0x3ec662fcu));
73 }
74 else if (magnitude <= 0x40800000u) {
75 t = fp32_fma<true,Hardware>(z, fp32_decode(0x3e800000u), fp32_decode(0xc0480000u));
76 h = fp32_decode(0xb58f6d10u);
77 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0x36863471u));
78 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0xb758ec9du));
79 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0x384b3fb7u));
80 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0xb9404721u));
81 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0x3a356f7bu));
82 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0xbb2d3d49u));
83 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0x3c2a7e14u));
84 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0xbd36d397u));
85 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0x3e9091d7u));
86 }
87 else if (magnitude <= 0x40c00000u) {
88 t = fp32_fma<true,Hardware>(z, fp32_decode(0x3d800000u), fp32_decode(0xbfd00000u));
89 h = fp32_decode(0xb9405fbfu);
90 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0x39ab4d5fu));
91 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0xb9c08feeu));
92 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0x3a2bf5fau));
93 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0xbaa65b2au));
94 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0x3b1584c1u));
95 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0xbb86bc79u));
96 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0x3bf71905u));
97 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0xbc67da1fu));
98 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0x3ce35ff2u));
99 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0xbd76f5e5u));
100 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0x3e48ced7u));
101 }
102 else if (magnitude <= 0x41000000u) {
103 t = fp32_fma<true,Hardware>(z, fp32_decode(0x3d800000u), fp32_decode(0xc0480000u));
104 h = fp32_decode(0xb59171b1u);
105 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0x3671d749u));
106 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0xb72757e1u));
107 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0x380d4deeu));
108 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0xb8f4646bu));
109 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0x39d46b54u));
110 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0xbabdbc02u));
111 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0x3bb1ed8au));
112 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0xbcb95c19u));
113 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0x3e10d0b5u));
114 }
115 else {
116 t = fp32_fma<true,Hardware>(z, fp32_decode(0x3d000000u), fp32_decode(0xc0240000u));
117 h = fp32_decode(0xb811c415u);
118 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0x38c8e73cu));
119 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0xb980a2b8u));
120 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0x3a37296bu));
121 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0xbb06690eu));
122 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0x3bcea86au));
123 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0xbcb08499u));
124 h = fp32_fma<true,Hardware>(h, t, fp32_decode(0x3de229ecu));
125 }
126 // Keep the mathematical range despite the last multiply's rounding error.
127 unsigned int rounded_magnitude = fp32_encode(fp32_mul<true,Hardware>(x, h));
128 return (rounded_magnitude > 0x3f800000u ? 0x3f800000u : rounded_magnitude) | sign;
129 }
130 template <bool Hardware = FTZ_FP32_HARDWARE_FTZ != 0>
131 [[nodiscard]] native_constexpr native_inline native_const float tanh(float value) {
132 return fp32_decode(tanh_words<Hardware>(fp32_encode(value)));
133 }
134}}}
135
float.h
FP32 bit casts, precise arithmetic and signed flush-to-zero normalization.