1 // Auto-generated file. Do not edit!
2 // Template: src/f32-raddstoreexpminusmax/wasmsimd-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 <wasm_simd128.h>
13
14 #include <xnnpack/common.h>
15 #include <xnnpack/raddstoreexpminusmax.h>
16
17
xnn_f32_raddstoreexpminusmax_ukernel__wasmsimd_rr2_p5_x12(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)])18 void xnn_f32_raddstoreexpminusmax_ukernel__wasmsimd_rr2_p5_x12(
19 size_t elements,
20 const float* input,
21 const float* max,
22 float* output,
23 float* sum,
24 const union xnn_f32_expminus_params params[restrict XNN_MIN_ELEMENTS(1)]) XNN_OOB_READS
25 {
26 assert(elements % sizeof(float) == 0);
27
28 const v128_t vi_max = wasm_v128_load32_splat(max);
29 const v128_t vlog2e = wasm_v128_load64_splat(params->wasmsimd_rr2_p5.log2e);
30 const v128_t vmagic_bias = wasm_v128_load64_splat(params->wasmsimd_rr2_p5.magic_bias);
31 const v128_t vminus_ln2_hi = wasm_v128_load64_splat(params->wasmsimd_rr2_p5.minus_ln2_hi);
32 const v128_t vminus_ln2_lo = wasm_v128_load64_splat(params->wasmsimd_rr2_p5.minus_ln2_lo);
33 const v128_t vc5 = wasm_v128_load64_splat(params->wasmsimd_rr2_p5.c5);
34 const v128_t vc4 = wasm_v128_load64_splat(params->wasmsimd_rr2_p5.c4);
35 const v128_t vc3 = wasm_v128_load64_splat(params->wasmsimd_rr2_p5.c3);
36 const v128_t vc2 = wasm_v128_load64_splat(params->wasmsimd_rr2_p5.c2);
37 const v128_t vc1 = wasm_v128_load64_splat(params->wasmsimd_rr2_p5.c1);
38 const v128_t vdenorm_cutoff = wasm_v128_load64_splat(params->wasmsimd_rr2_p5.denorm_cutoff);
39
40 v128_t vacc0 = wasm_f32x4_const_splat(0.0f);
41 for (; elements >= 12 * sizeof(float); elements -= 12 * sizeof(float)) {
42 // Load 12 (3x4) inputs at a time.
43 const v128_t vi0123 = wasm_v128_load(input);
44 const v128_t vi4567 = wasm_v128_load(input + 4);
45 const v128_t vi89AB = wasm_v128_load(input + 8);
46 input += 12;
47
48 // Subtract maximum input x := i - i_max. This implies x <= 0.
49 const v128_t vx0123 = wasm_f32x4_sub(vi0123, vi_max);
50 const v128_t vx4567 = wasm_f32x4_sub(vi4567, vi_max);
51 const v128_t vx89AB = wasm_f32x4_sub(vi89AB, vi_max);
52
53 // Compute reduced argument elements := round(x / log(2)).
54 v128_t vn0123 = wasm_f32x4_add(vmagic_bias, wasm_f32x4_mul(vx0123, vlog2e));
55 v128_t vn4567 = wasm_f32x4_add(vmagic_bias, wasm_f32x4_mul(vx4567, vlog2e));
56 v128_t vn89AB = wasm_f32x4_add(vmagic_bias, wasm_f32x4_mul(vx89AB, vlog2e));
57
58 // Create a floating-point number s (scale) such that s == 2**elements for inputs which don't cause underflow, i.e.
59 // -87.33642 <= x <= 0.0, and -126 <= elements <= 0 accordingly.
60 const v128_t vs0123 = wasm_i32x4_shl(vn0123, 23);
61 const v128_t vs4567 = wasm_i32x4_shl(vn4567, 23);
62 const v128_t vs89AB = wasm_i32x4_shl(vn89AB, 23);
63
64 // Subtract the large number back to get final elements := round(x / log(2)).
65 vn0123 = wasm_f32x4_sub(vn0123, vmagic_bias);
66 vn4567 = wasm_f32x4_sub(vn4567, vmagic_bias);
67 vn89AB = wasm_f32x4_sub(vn89AB, vmagic_bias);
68
69 // Compute reduced argument t := x - elements * log(2).
70 // Use Cody-Waite range reduction method (note two constants to represent log(2)) to improve accuracy.
71 v128_t vt0123 = wasm_f32x4_add(vx0123, wasm_f32x4_mul(vn0123, vminus_ln2_hi));
72 v128_t vt4567 = wasm_f32x4_add(vx4567, wasm_f32x4_mul(vn4567, vminus_ln2_hi));
73 v128_t vt89AB = wasm_f32x4_add(vx89AB, wasm_f32x4_mul(vn89AB, vminus_ln2_hi));
74
75 vt0123 = wasm_f32x4_add(vt0123, wasm_f32x4_mul(vn0123, vminus_ln2_lo));
76 vt4567 = wasm_f32x4_add(vt4567, wasm_f32x4_mul(vn4567, vminus_ln2_lo));
77 vt89AB = wasm_f32x4_add(vt89AB, wasm_f32x4_mul(vn89AB, vminus_ln2_lo));
78
79 // Compute degree-5 polynomial approximation for exp(t) on [-log(2)/2, log(2)/2].
80 v128_t vp0123 = wasm_f32x4_add(vc4, wasm_f32x4_mul(vc5, vt0123));
81 v128_t vp4567 = wasm_f32x4_add(vc4, wasm_f32x4_mul(vc5, vt4567));
82 v128_t vp89AB = wasm_f32x4_add(vc4, wasm_f32x4_mul(vc5, vt89AB));
83
84 vp0123 = wasm_f32x4_add(vc3, wasm_f32x4_mul(vp0123, vt0123));
85 vp4567 = wasm_f32x4_add(vc3, wasm_f32x4_mul(vp4567, vt4567));
86 vp89AB = wasm_f32x4_add(vc3, wasm_f32x4_mul(vp89AB, vt89AB));
87
88 vp0123 = wasm_f32x4_add(vc2, wasm_f32x4_mul(vp0123, vt0123));
89 vp4567 = wasm_f32x4_add(vc2, wasm_f32x4_mul(vp4567, vt4567));
90 vp89AB = wasm_f32x4_add(vc2, wasm_f32x4_mul(vp89AB, vt89AB));
91
92 vp0123 = wasm_f32x4_add(vc1, wasm_f32x4_mul(vp0123, vt0123));
93 vp4567 = wasm_f32x4_add(vc1, wasm_f32x4_mul(vp4567, vt4567));
94 vp89AB = wasm_f32x4_add(vc1, wasm_f32x4_mul(vp89AB, vt89AB));
95
96 // Reconstruct the final f value:
97 // f = s * (1 + t * (c1 + t * (c2 + t * (c3 + t * (c4 + t * c5)))))
98 // = s + (t * s) * (c1 + t * (c2 + t * (c3 + t * (c4 + t * c5))))
99 // = s + (t * s) * p
100 vt0123 = wasm_f32x4_mul(vt0123, vs0123);
101 vt4567 = wasm_f32x4_mul(vt4567, vs4567);
102 vt89AB = wasm_f32x4_mul(vt89AB, vs89AB);
103
104 v128_t vf0123 = wasm_f32x4_add(vs0123, wasm_f32x4_mul(vt0123, vp0123));
105 v128_t vf4567 = wasm_f32x4_add(vs4567, wasm_f32x4_mul(vt4567, vp4567));
106 v128_t vf89AB = wasm_f32x4_add(vs89AB, wasm_f32x4_mul(vt89AB, vp89AB));
107
108 // For inputs below zero cutoff, replace output with +0.0f.
109 // Note that for NaN inputs, comparison result is false, and outputs are left unchanged.
110 vf0123 = wasm_v128_andnot(vf0123, wasm_f32x4_lt(vx0123, vdenorm_cutoff));
111 vf4567 = wasm_v128_andnot(vf4567, wasm_f32x4_lt(vx4567, vdenorm_cutoff));
112 vf89AB = wasm_v128_andnot(vf89AB, wasm_f32x4_lt(vx89AB, vdenorm_cutoff));
113
114 // Store 12 (3x4) outputs at a time.
115 wasm_v128_store(output, vf0123);
116 wasm_v128_store(output + 4, vf4567);
117 wasm_v128_store(output + 8, vf89AB);
118 output += 12;
119
120 // Accumulate computed exponents.
121 vacc0 = wasm_f32x4_add(vacc0, vf0123);
122 vacc0 = wasm_f32x4_add(vacc0, vf4567);
123 vacc0 = wasm_f32x4_add(vacc0, vf89AB);
124 }
125
126 v128_t vacc = vacc0;
127 for (; elements >= 4 * sizeof(float); elements -= 4 * sizeof(float)) {
128 // Load 4 inputs at a time.
129 const v128_t vi = wasm_v128_load(input);
130 input += 4;
131
132 // Subtract maximum input x := i - i_max. This implies x <= 0.
133 const v128_t vx = wasm_f32x4_sub(vi, vi_max);
134
135 // Compute reduced argument elements := round(x / log(2)).
136 v128_t vn = wasm_f32x4_add(vmagic_bias, wasm_f32x4_mul(vx, vlog2e));
137
138 // Create a floating-point number s (scale) such that s == 2**elements for inputs which don't cause underflow, i.e.
139 // -87.33642 <= x <= 0.0, and -126 <= elements <= 0 accordingly.
140 const v128_t vs = wasm_i32x4_shl(vn, 23);
141
142 // Subtract the large number back to get final elements := round(x / log(2)).
143 vn = wasm_f32x4_sub(vn, vmagic_bias);
144
145 // Compute reduced argument t := x - elements * log(2).
146 // Use Cody-Waite range reduction method (note two constants to represent log(2)) to improve accuracy.
147 v128_t vt = wasm_f32x4_add(vx, wasm_f32x4_mul(vn, vminus_ln2_hi));
148 vt = wasm_f32x4_add(vt, wasm_f32x4_mul(vn, vminus_ln2_lo));
149
150 // Compute degree-5 polynomial approximation for exp(t) on [-log(2)/2, log(2)/2].
151 v128_t vp = wasm_f32x4_add(vc4, wasm_f32x4_mul(vc5, vt));
152 vp = wasm_f32x4_add(vc3, wasm_f32x4_mul(vp, vt));
153 vp = wasm_f32x4_add(vc2, wasm_f32x4_mul(vp, vt));
154 vp = wasm_f32x4_add(vc1, wasm_f32x4_mul(vp, vt));
155
156 // Reconstruct the final f value:
157 // f = s * (1 + t * (c1 + t * (c2 + t * (c3 + t * (c4 + t * c5)))))
158 // = s + (t * s) * (c1 + t * (c2 + t * (c3 + t * (c4 + t * c5))))
159 // = s + (t * s) * p
160 vt = wasm_f32x4_mul(vt, vs);
161 v128_t vf = wasm_f32x4_add(vs, wasm_f32x4_mul(vt, vp));
162
163 // For inputs below zero cutoff, replace output with +0.0f.
164 // Note that for NaN inputs, comparison result is false, and outputs are left unchanged.
165 vf = wasm_v128_andnot(vf, wasm_f32x4_lt(vx, vdenorm_cutoff));
166
167 // Store 4 outputs at a time.
168 wasm_v128_store(output, vf);
169 output += 4;
170
171 // Accumulate computed exponents.
172 vacc = wasm_f32x4_add(vacc, vf);
173 }
174 vacc = wasm_f32x4_add(vacc, wasm_v32x4_shuffle(vacc, vacc, 2, 3, 2, 3));
175 float vsum = wasm_f32x4_extract_lane(vacc, 0) + wasm_f32x4_extract_lane(vacc, 1);
176 if (elements != 0) {
177 assert(elements >= 1 * sizeof(float));
178 assert(elements <= 3 * sizeof(float));
179 // Load 4 inputs at a time.
180 const v128_t vi = wasm_v128_load(input);
181
182 // Subtract maximum input x := i - i_max. This implies x <= 0.
183 const v128_t vx = wasm_f32x4_sub(vi, vi_max);
184
185 // Compute reduced argument elements := round(x / log(2)).
186 v128_t vn = wasm_f32x4_add(vmagic_bias, wasm_f32x4_mul(vx, vlog2e));
187
188 // Create a floating-point number s (scale) such that s == 2**elements for inputs which don't cause underflow, i.e.
189 // -87.33642 <= x <= 0.0, and -126 <= elements <= 0 accordingly.
190 const v128_t vs = wasm_i32x4_shl(vn, 23);
191
192 // Subtract the large number back to get final elements := round(x / log(2)).
193 vn = wasm_f32x4_sub(vn, vmagic_bias);
194
195 // Compute reduced argument t := x - elements * log(2).
196 // Use Cody-Waite range reduction method (note two constants to represent log(2)) to improve accuracy.
197 v128_t vt = wasm_f32x4_add(vx, wasm_f32x4_mul(vn, vminus_ln2_hi));
198 vt = wasm_f32x4_add(vt, wasm_f32x4_mul(vn, vminus_ln2_lo));
199
200 // Compute degree-5 polynomial approximation for exp(t) on [-log(2)/2, log(2)/2].
201 v128_t vp = wasm_f32x4_add(vc4, wasm_f32x4_mul(vc5, vt));
202 vp = wasm_f32x4_add(vc3, wasm_f32x4_mul(vp, vt));
203 vp = wasm_f32x4_add(vc2, wasm_f32x4_mul(vp, vt));
204 vp = wasm_f32x4_add(vc1, wasm_f32x4_mul(vp, vt));
205
206 // Reconstruct the final f value:
207 // f = s * (1 + t * (c1 + t * (c2 + t * (c3 + t * (c4 + t * c5)))))
208 // = s + (t * s) * (c1 + t * (c2 + t * (c3 + t * (c4 + t * c5))))
209 // = s + (t * s) * p
210 vt = wasm_f32x4_mul(vt, vs);
211 v128_t vf = wasm_f32x4_add(vs, wasm_f32x4_mul(vt, vp));
212
213 // For inputs below zero cutoff, replace output with +0.0f.
214 // Note that for NaN inputs, comparison result is false, and outputs are left unchanged.
215 vf = wasm_v128_andnot(vf, wasm_f32x4_lt(vx, vdenorm_cutoff));
216
217 if (elements & (2 * sizeof(float))) {
218 // Store and accumulate 2 outputs at a time.
219 const float vf0 = wasm_f32x4_extract_lane(vf, 0);
220 output[0] = vf0;
221 vsum += vf0;
222
223 const float vf1 = wasm_f32x4_extract_lane(vf, 1);
224 output[1] = vf1;
225 vsum += vf1;
226
227 vf = wasm_v32x4_shuffle(vf, vf, 2, 3, 2, 3);
228 output += 2;
229 }
230 if (elements & (1 * sizeof(float))) {
231 // Store 1 output at a time.
232 const float vf0 = wasm_f32x4_extract_lane(vf, 0);
233 *output = vf0;
234 vsum += vf0;
235 }
236 }
237 // Reduce 4 elements in the SIMD register
238 *sum = vsum;
239 }
240