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