• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
1 // Auto-generated file. Do not edit!
2 //   Template: src/f32-sigmoid/avx2-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 <immintrin.h>
13 
14 #include <xnnpack/common.h>
15 #include <xnnpack/vunary.h>
16 
17 
18 static const int32_t mask_table[14] = {-1, -1, -1, -1, -1, -1, -1, 0, 0, 0, 0, 0, 0, 0};
19 
xnn_f32_sigmoid_ukernel__avx2_rr1_p5_div_x16(size_t n,const float * x,float * y,const void * params)20 void xnn_f32_sigmoid_ukernel__avx2_rr1_p5_div_x16(
21     size_t n,
22     const float* x,
23     float* y,
24     const void* params)
25 {
26   assert(n % sizeof(float) == 0);
27 
28   const __m256 vmagic_bias = _mm256_set1_ps(0x1.8000FEp23f);
29   // The smallest x for which sigmoidf(x) is normalized.
30   // This number is also the smallest x for which expf(x) is normalized.
31   const __m256 vdenorm_cutoff = _mm256_set1_ps(-0x1.5D589Ep+6f);
32   const __m256 vlog2e = _mm256_set1_ps(0x1.715476p+0f);
33   const __m256 vminus_ln2 = _mm256_set1_ps(-0x1.62E43p-1f);
34   const __m256 vone = _mm256_set1_ps(1.0f);
35   const __m256 vsign_mask = _mm256_set1_ps(-0.0f);
36 
37   const __m256 vc1 = _mm256_set1_ps(0x1.FFFFF6p-1f);
38   const __m256 vc2 = _mm256_set1_ps(0x1.FFFDC6p-2f);
39   const __m256 vc3 = _mm256_set1_ps(0x1.555A80p-3f);
40   const __m256 vc4 = _mm256_set1_ps(0x1.573A1Ap-5f);
41   const __m256 vc5 = _mm256_set1_ps(0x1.0F9F9Cp-7f);
42 
43   for (; n >= 16 * sizeof(float); n -= 16 * sizeof(float)) {
44     const __m256 vx0 = _mm256_loadu_ps(x);
45     const __m256 vx1 = _mm256_loadu_ps(x + 8);
46     x += 16;
47 
48     // General structure of the algorithm:
49     //           / exp(x) / (1 + exp(x)) if x <= 0
50     //   f[x] :=
51     //           \ 1 - f[-x] if x >= 0
52     //
53     // First we compute f[z] := exp(z) / (1 + exp(z)) where z = -abs(x),
54     // then replace result with 1 - f[z] if x >= 0.
55     const __m256 vz0 = _mm256_or_ps(vx0, vsign_mask);
56     const __m256 vz1 = _mm256_or_ps(vx1, vsign_mask);
57 
58     // Compute reduced argument n := round(z / log(2)).
59     // We do it by adding a large number (magic bias) to the product z * (1/log(2)), which cause rounding of the result
60     // to an integer, then subtracing the large number back. The trick with adding large number is valid only within
61     // certain bounds (|x| <= 2**22), but thats ok, because inputs x outside of [-87.336544, 17.328678] (i.e. z outsize
62     // [0, 87.336544]) underflow or saturate sigmoidf(x) anyway. We fixup the result for such inputs at the very end of
63     // the algorithm.
64     __m256 vn0 = _mm256_fmadd_ps(vz0, vlog2e, vmagic_bias);
65     __m256 vn1 = _mm256_fmadd_ps(vz1, 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 <= z <= 0.0, and -126 <= n <= 0 accordingly.
69     const __m256 vs0 = _mm256_castsi256_ps(_mm256_slli_epi32(_mm256_castps_si256(vn0), 23));
70     const __m256 vs1 = _mm256_castsi256_ps(_mm256_slli_epi32(_mm256_castps_si256(vn1), 23));
71 
72     // Subtract the large number back to get final n := round(z / log(2)).
73     vn0 = _mm256_sub_ps(vn0, vmagic_bias);
74     vn1 = _mm256_sub_ps(vn1, vmagic_bias);
75 
76     // Compute reduced argument t := z - n * log(2).
77     __m256 vt0 = _mm256_fmadd_ps(vn0, vminus_ln2, vz0);
78     __m256 vt1 = _mm256_fmadd_ps(vn1, vminus_ln2, vz1);
79 
80     // Compute degree-5 polynomial approxiatmion for exp(t) on [-log(2)/2, log(2)/2].
81     __m256 vp0 = _mm256_fmadd_ps(vc5, vt0, vc4);
82     __m256 vp1 = _mm256_fmadd_ps(vc5, vt1, vc4);
83 
84     vp0 = _mm256_fmadd_ps(vp0, vt0, vc3);
85     vp1 = _mm256_fmadd_ps(vp1, vt1, vc3);
86 
87     vp0 = _mm256_fmadd_ps(vp0, vt0, vc2);
88     vp1 = _mm256_fmadd_ps(vp1, vt1, vc2);
89 
90     vp0 = _mm256_fmadd_ps(vp0, vt0, vc1);
91     vp1 = _mm256_fmadd_ps(vp1, vt1, vc1);
92 
93     // Reconstruct the exp(z) value:
94     //   e = s * (1 + t * (c1 + t * (c2 + t * (c3 + t * (c4 + t * c5)))))
95     //     = s + (t * s) * (c1 + t * (c2 + t * (c3 + t * (c4 + t * c5))))
96     //     = s + (t * s) * p
97     vt0 = _mm256_mul_ps(vt0, vs0);
98     vt1 = _mm256_mul_ps(vt1, vs1);
99 
100     const __m256 ve0 = _mm256_fmadd_ps(vt0, vp0, vs0);
101     const __m256 ve1 = _mm256_fmadd_ps(vt1, vp1, vs1);
102 
103     // Denominator of the sigmoid fraction: 1.0 + exp(z)
104     const __m256 vd0 = _mm256_add_ps(ve0, vone);
105     const __m256 vd1 = _mm256_add_ps(ve1, vone);
106 
107     // Reconstruct sigmoid(z) = exp(z) / (1.0 + exp(z))
108     __m256 vf0 = _mm256_div_ps(ve0, vd0);
109     __m256 vf1 = _mm256_div_ps(ve1, vd1);
110 
111     // For inputs below denormal cutoff, replace output with +0.0f.
112     // Note that for NaN inputs, comparison result is false, and outputs are left unchanged.
113     vf0 = _mm256_andnot_ps(_mm256_cmp_ps(vz0, vdenorm_cutoff, _CMP_LT_OS), vf0);
114     vf1 = _mm256_andnot_ps(_mm256_cmp_ps(vz1, vdenorm_cutoff, _CMP_LT_OS), vf1);
115 
116     // Reconstruct sigmoid(x) = x < 0 ? sigmoid(z) : 1.0 - sigmoid(z)
117     vf0 = _mm256_blendv_ps(_mm256_sub_ps(vone, vf0), vf0, vx0);
118     vf1 = _mm256_blendv_ps(_mm256_sub_ps(vone, vf1), vf1, vx1);
119 
120     _mm256_storeu_ps(y, vf0);
121     _mm256_storeu_ps(y + 8, vf1);
122     y += 16;
123   }
124   for (; n >= 8 * sizeof(float); n -= 8 * sizeof(float)) {
125     const __m256 vx = _mm256_loadu_ps(x);
126     x += 8;
127 
128     // General structure of the algorithm:
129     //           / exp(x) / (1 + exp(x)) if x <= 0
130     //   f[x] :=
131     //           \ 1 - f[-x] if x >= 0
132     //
133     // First we compute f[z] := exp(z) / (1 + exp(z)) where z = -abs(x),
134     // then replace result with 1 - f[z] if x >= 0.
135     const __m256 vz = _mm256_or_ps(vx, vsign_mask);
136 
137     // Compute reduced argument n := round(z / log(2)).
138     // We do it by adding a large number (magic bias) to the product z * (1/log(2)), which cause rounding of the result
139     // to an integer, then subtracing the large number back. The trick with adding large number is valid only within
140     // certain bounds (|x| <= 2**22), but thats ok, because inputs x outside of [-87.336544, 17.328678] (i.e. z outsize
141     // [0, 87.336544]) underflow or saturate sigmoidf(x) anyway. We fixup the result for such inputs at the very end of
142     // the algorithm.
143     __m256 vn = _mm256_fmadd_ps(vz, vlog2e, vmagic_bias);
144 
145     // Create a floating-point number s (scale) such that s == 2**n for inputs which don't cause underflow, i.e.
146     // -87.33642 <= z <= 0.0, and -126 <= n <= 0 accordingly.
147     const __m256 vs = _mm256_castsi256_ps(_mm256_slli_epi32(_mm256_castps_si256(vn), 23));
148 
149     // Subtract the large number back to get final n := round(z / log(2)).
150     vn = _mm256_sub_ps(vn, vmagic_bias);
151 
152     // Compute reduced argument t := z - n * log(2).
153     __m256 vt = _mm256_fmadd_ps(vn, vminus_ln2, vz);
154 
155     // Compute degree-5 polynomial approxiatmion for exp(t) on [-log(2)/2, log(2)/2].
156     __m256 vp = _mm256_fmadd_ps(vc5, vt, vc4);
157     vp = _mm256_fmadd_ps(vp, vt, vc3);
158     vp = _mm256_fmadd_ps(vp, vt, vc2);
159     vp = _mm256_fmadd_ps(vp, vt, vc1);
160 
161     // Reconstruct the exp(z) value:
162     //   e = 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 = _mm256_mul_ps(vt, vs);
166     const __m256 ve = _mm256_fmadd_ps(vt, vp, vs);
167 
168     // Denominator of the sigmoid fraction: 1.0 + exp(z)
169     const __m256 vd = _mm256_add_ps(ve, vone);
170 
171     // Reconstruct sigmoid(z) = exp(z) / (1.0 + exp(z))
172     __m256 vf = _mm256_div_ps(ve, vd);
173 
174     // For inputs below denormal cutoff, replace output with +0.0f.
175     // Note that for NaN inputs, comparison result is false, and outputs are left unchanged.
176     vf = _mm256_andnot_ps(_mm256_cmp_ps(vz, vdenorm_cutoff, _CMP_LT_OS), vf);
177 
178     // Reconstruct sigmoid(x) = x < 0 ? sigmoid(z) : 1.0 - sigmoid(z)
179     vf = _mm256_blendv_ps(_mm256_sub_ps(vone, vf), vf, vx);
180 
181     _mm256_storeu_ps(y, vf);
182     y += 8;
183   }
184   if XNN_UNLIKELY(n != 0) {
185     assert(n >= 1 * sizeof(float));
186     assert(n <= 7 * sizeof(float));
187     __m256i vmask = _mm256_loadu_si256((const __m256i*) ((uintptr_t) &mask_table[7] - n));
188 
189     const __m256 vx = _mm256_maskload_ps(x, vmask);
190 
191     // General structure of the algorithm:
192     //           / exp(x) / (1 + exp(x)) if x <= 0
193     //   f[x] :=
194     //           \ 1 - f[-x] if x >= 0
195     //
196     // First we compute f[z] := exp(z) / (1 + exp(z)) where z = -abs(x),
197     // then replace result with 1 - f[z] if x >= 0.
198     const __m256 vz = _mm256_or_ps(vx, vsign_mask);
199 
200     // Compute reduced argument n := round(z / log(2)).
201     // We do it by adding a large number (magic bias) to the product z * (1/log(2)), which cause rounding of the result
202     // to an integer, then subtracing the large number back. The trick with adding large number is valid only within
203     // certain bounds (|x| <= 2**22), but thats ok, because inputs x outside of [-87.336544, 17.328678] (i.e. z outsize
204     // [0, 87.336544]) underflow or saturate sigmoidf(x) anyway. We fixup the result for such inputs at the very end of
205     // the algorithm.
206     __m256 vn = _mm256_fmadd_ps(vz, vlog2e, vmagic_bias);
207 
208     // Create a floating-point number s (scale) such that s == 2**n for inputs which don't cause underflow, i.e.
209     // -87.33642 <= z <= 0.0, and -126 <= n <= 0 accordingly.
210     const __m256 vs = _mm256_castsi256_ps(_mm256_slli_epi32(_mm256_castps_si256(vn), 23));
211 
212     // Subtract the large number back to get final n := round(z / log(2)).
213     vn = _mm256_sub_ps(vn, vmagic_bias);
214 
215     // Compute reduced argument t := z - n * log(2).
216     __m256 vt = _mm256_fmadd_ps(vn, vminus_ln2, vz);
217 
218     // Compute degree-5 polynomial approxiatmion for exp(t) on [-log(2)/2, log(2)/2].
219     __m256 vp = _mm256_fmadd_ps(vc5, vt, vc4);
220     vp = _mm256_fmadd_ps(vp, vt, vc3);
221     vp = _mm256_fmadd_ps(vp, vt, vc2);
222     vp = _mm256_fmadd_ps(vp, vt, vc1);
223 
224     // Reconstruct the exp(z) value:
225     //   e = s * (1 + t * (c1 + t * (c2 + t * (c3 + t * (c4 + t * c5)))))
226     //     = s + (t * s) * (c1 + t * (c2 + t * (c3 + t * (c4 + t * c5))))
227     //     = s + (t * s) * p
228     vt = _mm256_mul_ps(vt, vs);
229     const __m256 ve = _mm256_fmadd_ps(vt, vp, vs);
230 
231     // Denominator of the sigmoid fraction: 1.0 + exp(z)
232     const __m256 vd = _mm256_add_ps(ve, vone);
233 
234     // Reconstruct sigmoid(z) = exp(z) / (1.0 + exp(z))
235     __m256 vf = _mm256_div_ps(ve, vd);
236 
237     // For inputs below denormal cutoff, replace output with +0.0f.
238     // Note that for NaN inputs, comparison result is false, and outputs are left unchanged.
239     vf = _mm256_andnot_ps(_mm256_cmp_ps(vz, vdenorm_cutoff, _CMP_LT_OS), vf);
240 
241     // Reconstruct sigmoid(x) = x < 0 ? sigmoid(z) : 1.0 - sigmoid(z)
242     vf = _mm256_blendv_ps(_mm256_sub_ps(vone, vf), vf, vx);
243 
244     // _mm256_maskstore_ps(y, vmask, vf) could be used here, but triggers msan failures (probably an msan bug).
245     __m128 vf_lo = _mm256_castps256_ps128(vf);
246     if (n & (4 * sizeof(float))) {
247       _mm_storeu_ps(y, vf_lo);
248       vf_lo = _mm256_extractf128_ps(vf, 1);
249       y += 4;
250     }
251     if (n & (2 * sizeof(float))) {
252       _mm_storel_pi((__m64*) y, vf_lo);
253       vf_lo = _mm_movehl_ps(vf_lo, vf_lo);
254       y += 2;
255     }
256     if (n & (1 * sizeof(float))) {
257       _mm_store_ss(y, vf_lo);
258     }
259   }
260 }
261