1 // Auto-generated file. Do not edit!
2 // Template: src/f32-raddstoreexpminusmax/scalar-p5.c.in
3 // Generator: tools/xngen
4 //
5 // Copyright 2020 Google LLC
6 //
7 // This source code is licensed under the BSD-style license found in the
8 // LICENSE file in the root directory of this source tree.
9
10 #include <assert.h>
11
12 #include <xnnpack/common.h>
13 #include <xnnpack/raddstoreexpminusmax.h>
14
15 #include <fp16/bitcasts.h>
16
17
xnn_f32_raddstoreexpminusmax_ukernel__scalar_p5_x4(size_t elements,const float * input,float * output,float * sum,float vi_max)18 void xnn_f32_raddstoreexpminusmax_ukernel__scalar_p5_x4(
19 size_t elements,
20 const float* input,
21 float* output,
22 float* sum,
23 float vi_max)
24 {
25 assert(elements % sizeof(float) == 0);
26
27 const float vmagic_bias = 0x1.8000FEp23f;
28 // The smallest x for which expf(x) is normalized.
29 const float vdenorm_cutoff = -0x1.5D589Ep6f;
30 const float vlog2e = 0x1.715476p+0f;
31 // Last 7 bits are zeroes
32 const float vminus_ln2_hi = -0x1.62E400p-1f;
33 const float vminus_ln2_lo = -0x1.7F7D1Cp-20f;
34
35 const float vc1 = 0x1.FFFFF6p-1f;
36 const float vc2 = 0x1.FFFDC6p-2f;
37 const float vc3 = 0x1.555A80p-3f;
38 const float vc4 = 0x1.573A1Ap-5f;
39 const float vc5 = 0x1.0F9F9Cp-7f;
40
41 float vacc0 = 0.0f;
42 for (; elements >= 4 * sizeof(float); elements -= 4 * sizeof(float)) {
43 // Load 4 inputs at a time.
44 const float vi0 = input[0];
45 const float vi1 = input[1];
46 const float vi2 = input[2];
47 const float vi3 = input[3];
48 input += 4;
49
50 // Subtract maximum input x := i - i_max. This implies x <= 0.
51 const float vx0 = vi0 - vi_max;
52 const float vx1 = vi1 - vi_max;
53 const float vx2 = vi2 - vi_max;
54 const float vx3 = vi3 - vi_max;
55
56 // Compute reduced argument n := round(x / log(2)).
57 // We do it by adding a large number (magic bias) to the product x * (1/log(2)), which cause rounding of the result
58 // to an integer, then subtracing the large number back. The trick with adding large number is valid only within
59 // certain bounds (|x| <= 2**22), but thats ok, because inputs outside of [-87.336540, 0.0] underflow expf(x)
60 // anyway. We fixup the result for such inputs at the very end of the algorithm.
61 float vn0 = vx0 * vlog2e + vmagic_bias;
62 float vn1 = vx1 * vlog2e + vmagic_bias;
63 float vn2 = vx2 * vlog2e + vmagic_bias;
64 float vn3 = vx3 * vlog2e + vmagic_bias;
65
66 // Create a floating-point number s (scale) such that s == 2**n for inputs which don't cause underflow, i.e.
67 // -87.33642 <= x <= 0.0, and -126 <= n <= 0 accordingly.
68 const float vs0 = fp32_from_bits(fp32_to_bits(vn0) << 23);
69 const float vs1 = fp32_from_bits(fp32_to_bits(vn1) << 23);
70 const float vs2 = fp32_from_bits(fp32_to_bits(vn2) << 23);
71 const float vs3 = fp32_from_bits(fp32_to_bits(vn3) << 23);
72
73 // Subtract the large number back to get final n := round(x / log(2)).
74 vn0 -= vmagic_bias;
75 vn1 -= vmagic_bias;
76 vn2 -= vmagic_bias;
77 vn3 -= vmagic_bias;
78
79 // Compute reduced argument t := x - n * log(2).
80 // Use Cody-Waite range reduction method (note two constants to represent log(2)) to improve accuracy.
81 float vt0 = vn0 * vminus_ln2_hi + vx0;
82 float vt1 = vn1 * vminus_ln2_hi + vx1;
83 float vt2 = vn2 * vminus_ln2_hi + vx2;
84 float vt3 = vn3 * vminus_ln2_hi + vx3;
85
86 vt0 = vn0 * vminus_ln2_lo + vt0;
87 vt1 = vn1 * vminus_ln2_lo + vt1;
88 vt2 = vn2 * vminus_ln2_lo + vt2;
89 vt3 = vn3 * vminus_ln2_lo + vt3;
90
91 // Compute degree-5 polynomial approximation for exp(t) on [-log(2)/2, log(2)/2].
92 float vp0 = vc5 * vt0 + vc4;
93 float vp1 = vc5 * vt1 + vc4;
94 float vp2 = vc5 * vt2 + vc4;
95 float vp3 = vc5 * vt3 + vc4;
96
97 vp0 = vp0 * vt0 + vc3;
98 vp1 = vp1 * vt1 + vc3;
99 vp2 = vp2 * vt2 + vc3;
100 vp3 = vp3 * vt3 + vc3;
101
102 vp0 = vp0 * vt0 + vc2;
103 vp1 = vp1 * vt1 + vc2;
104 vp2 = vp2 * vt2 + vc2;
105 vp3 = vp3 * vt3 + vc2;
106
107 vp0 = vp0 * vt0 + vc1;
108 vp1 = vp1 * vt1 + vc1;
109 vp2 = vp2 * vt2 + vc1;
110 vp3 = vp3 * vt3 + vc1;
111
112 // Reconstruct the final f value:
113 // f = s * (1 + t * (c1 + t * (c2 + t * (c3 + t * (c4 + t * c5)))))
114 // = s + (t * s) * (c1 + t * (c2 + t * (c3 + t * (c4 + t * c5))))
115 // = s + (t * s) * p
116 vt0 *= vs0;
117 vt1 *= vs1;
118 vt2 *= vs2;
119 vt3 *= vs3;
120
121 float vf0 = vt0 * vp0 + vs0;
122 float vf1 = vt1 * vp1 + vs1;
123 float vf2 = vt2 * vp2 + vs2;
124 float vf3 = vt3 * vp3 + vs3;
125
126 // For inputs below denormal cutoff, replace output with +0.0f.
127 // Note that for NaN inputs, comparison result is false, and outputs are left unchanged.
128 if XNN_UNPREDICTABLE(vx0 < vdenorm_cutoff) {
129 vf0 = 0.0f;
130 }
131 if XNN_UNPREDICTABLE(vx1 < vdenorm_cutoff) {
132 vf1 = 0.0f;
133 }
134 if XNN_UNPREDICTABLE(vx2 < vdenorm_cutoff) {
135 vf2 = 0.0f;
136 }
137 if XNN_UNPREDICTABLE(vx3 < vdenorm_cutoff) {
138 vf3 = 0.0f;
139 }
140
141 // Store 4 outputs at a time.
142 output[0] = vf0;
143 output[1] = vf1;
144 output[2] = vf2;
145 output[3] = vf3;
146 output += 4;
147
148 // Accumulate computed exponents.
149 vacc0 += vf0;
150 vacc0 += vf1;
151 vacc0 += vf2;
152 vacc0 += vf3;
153 }
154
155 float vacc = vacc0;
156 for (; elements >= sizeof(float); elements -= sizeof(float)) {
157 // Load 1 input at a time.
158 const float vi = *input++;
159
160 // Subtract maximum input x := i - i_max. This implies x <= 0.
161 const float vx = vi - vi_max;
162
163 // Compute reduced argument n := round(x / log(2)).
164 // We do it by adding a large number (magic bias) to the product x * (1/log(2)), which cause rounding of the result
165 // to an integer, then subtracing the large number back. The trick with adding large number is valid only within
166 // certain bounds (|x| <= 2**22), but thats ok, because inputs outside of [-87.336540, 0.0] underflow expf(x)
167 // anyway. We fixup the result for such inputs at the very end of the algorithm.
168 float vn = vx * vlog2e + vmagic_bias;
169
170 // Create a floating-point number s (scale) such that s == 2**n for inputs which don't cause underflow, i.e.
171 // -87.33642 <= x <= 0.0, and -126 <= n <= 0 accordingly.
172 const float vs = fp32_from_bits(fp32_to_bits(vn) << 23);
173
174 // Subtract the large number back to get final n := round(x / log(2)).
175 vn -= vmagic_bias;
176
177 // Compute reduced argument t := x - n * log(2).
178 // Use Cody-Waite range reduction method (note two constants to represent log(2)) to improve accuracy.
179 float vt = vn * vminus_ln2_hi + vx;
180 vt = vn * vminus_ln2_lo + vt;
181
182 // Compute degree-5 polynomial approximation for exp(t) on [-log(2)/2, log(2)/2].
183 float vp = vc5 * vt + vc4;
184 vp = vp * vt + vc3;
185 vp = vp * vt + vc2;
186 vp = vp * vt + vc1;
187
188 // Reconstruct the final f value:
189 // f = s * (1 + t * (c1 + t * (c2 + t * (c3 + t * (c4 + t * c5)))))
190 // = s + (t * s) * (c1 + t * (c2 + t * (c3 + t * (c4 + t * c5))))
191 // = s + (t * s) * p
192 vt *= vs;
193 float vf = vt * vp + vs;
194
195 // For inputs below denormal cutoff, replace output with +0.0f.
196 // Note that for NaN inputs, comparison result is false, and outputs are left unchanged.
197 if XNN_UNPREDICTABLE(vx < vdenorm_cutoff) {
198 vf = 0.0f;
199 }
200
201 // Store 1 output at a time.
202 *output++ = vf;
203
204 // Accumulate computed exponents.
205 vacc += vf;
206 }
207 *sum = vacc;
208 }
209