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