• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
1// Copyright 2019 Google LLC
2//
3// This source code is licensed under the BSD-style license found in the
4// LICENSE file in the root directory of this source tree.
5
6$assert ELEMENTS_TILE % 4 == 0
7$assert ELEMENTS_TILE >= 4
8$SIMD_TILE = ELEMENTS_TILE // 4
9$ABC = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ"
10#include <assert.h>
11
12#include <emmintrin.h>
13
14#include <xnnpack/common.h>
15#include <xnnpack/raddstoreexpminusmax.h>
16
17
18void xnn_f32_raddstoreexpminusmax_ukernel__sse2_p5_x${ELEMENTS_TILE}${"" if ACCUMULATORS == 1 else "_acc%d" % ACCUMULATORS}(
19    size_t elements,
20    const float* input,
21    float* output,
22    float* sum,
23    float max)
24{
25  assert(elements % sizeof(float) == 0);
26
27  const __m128 vmagic_bias = _mm_set1_ps(0x1.8000FEp23f);
28  // The smallest x for which expf(x) is normalized.
29  const __m128 vdenorm_cutoff = _mm_set1_ps(-0x1.5D589Ep6f);
30  const __m128 vlog2e = _mm_set1_ps(0x1.715476p+0f);
31  // Last 7 bits are zeroes
32  const __m128 vminus_ln2_hi = _mm_set1_ps(-0x1.62E400p-1f);
33  const __m128 vminus_ln2_lo = _mm_set1_ps(-0x1.7F7D1Cp-20f);
34
35  const __m128 vc1 = _mm_set1_ps(0x1.FFFFF6p-1f);
36  const __m128 vc2 = _mm_set1_ps(0x1.FFFDC6p-2f);
37  const __m128 vc3 = _mm_set1_ps(0x1.555A80p-3f);
38  const __m128 vc4 = _mm_set1_ps(0x1.573A1Ap-5f);
39  const __m128 vc5 = _mm_set1_ps(0x1.0F9F9Cp-7f);
40
41  const __m128 vi_max = _mm_set1_ps(max);
42
43  $for K in range(ACCUMULATORS):
44    __m128 vacc${K} = _mm_setzero_ps();
45  for (; elements >= ${ELEMENTS_TILE} * sizeof(float); elements -= ${ELEMENTS_TILE} * sizeof(float)) {
46    // Load ${ELEMENTS_TILE} (${SIMD_TILE}x4) inputs at a time.
47    const __m128 vi${ABC[0:4]} = _mm_loadu_ps(input);
48    $for N in range(4, ELEMENTS_TILE, 4):
49      const __m128 vi${ABC[N:N+4]} = _mm_loadu_ps(input + ${N});
50    input += ${ELEMENTS_TILE};
51
52    // Subtract maximum input x := i - i_max. This implies x <= 0.
53    $for N in range(0, ELEMENTS_TILE, 4):
54      const __m128 vx${ABC[N:N+4]} = _mm_sub_ps(vi${ABC[N:N+4]}, vi_max);
55
56    // Compute reduced argument elements := round(x / log(2)).
57    $for N in range(0, ELEMENTS_TILE, 4):
58      __m128 vn${ABC[N:N+4]} = _mm_add_ps(_mm_mul_ps(vx${ABC[N:N+4]}, 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    $for N in range(0, ELEMENTS_TILE, 4):
63      const __m128 vs${ABC[N:N+4]} = _mm_castsi128_ps(_mm_slli_epi32(_mm_castps_si128(vn${ABC[N:N+4]}), 23));
64
65    // Subtract the large number back to get final elements := round(x / log(2)).
66    $for N in range(0, ELEMENTS_TILE, 4):
67      vn${ABC[N:N+4]} = _mm_sub_ps(vn${ABC[N:N+4]}, 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    $for N in range(0, ELEMENTS_TILE, 4):
72      __m128 vt${ABC[N:N+4]} = _mm_add_ps(_mm_mul_ps(vn${ABC[N:N+4]}, vminus_ln2_hi), vx${ABC[N:N+4]});
73
74    $for N in range(0, ELEMENTS_TILE, 4):
75      vt${ABC[N:N+4]} = _mm_add_ps(_mm_mul_ps(vn${ABC[N:N+4]}, vminus_ln2_lo), vt${ABC[N:N+4]});
76
77    // Compute degree-5 polynomial approxiatmion for exp(t) on [-log(2)/2, log(2)/2].
78    $for N in range(0, ELEMENTS_TILE, 4):
79      __m128 vp${ABC[N:N+4]} = _mm_add_ps(_mm_mul_ps(vc5, vt${ABC[N:N+4]}), vc4);
80
81    $for N in range(0, ELEMENTS_TILE, 4):
82      vp${ABC[N:N+4]} = _mm_add_ps(_mm_mul_ps(vp${ABC[N:N+4]}, vt${ABC[N:N+4]}), vc3);
83
84    $for N in range(0, ELEMENTS_TILE, 4):
85      vp${ABC[N:N+4]} = _mm_add_ps(_mm_mul_ps(vp${ABC[N:N+4]}, vt${ABC[N:N+4]}), vc2);
86
87    $for N in range(0, ELEMENTS_TILE, 4):
88      vp${ABC[N:N+4]} = _mm_add_ps(_mm_mul_ps(vp${ABC[N:N+4]}, vt${ABC[N:N+4]}), vc1);
89
90    // Reconstruct the final f value:
91    //   f = s * (1 + t * (c1 + t * (c2 + t * (c3 + t * (c4 + t * c5)))))
92    //     = s + (t * s) * (c1 + t * (c2 + t * (c3 + t * (c4 + t * c5))))
93    //     = s + (t * s) * p
94    $for N in range(0, ELEMENTS_TILE, 4):
95      vt${ABC[N:N+4]} = _mm_mul_ps(vt${ABC[N:N+4]}, vs${ABC[N:N+4]});
96
97    $for N in range(0, ELEMENTS_TILE, 4):
98      __m128 vf${ABC[N:N+4]} = _mm_add_ps(_mm_mul_ps(vt${ABC[N:N+4]}, vp${ABC[N:N+4]}), vs${ABC[N:N+4]});
99
100    // For inputs below zero cutoff, replace output with +0.0f.
101    // Note that for NaN inputs, comparison result is false, and outputs are left unchanged.
102    $for N in range(0, ELEMENTS_TILE, 4):
103      vf${ABC[N:N+4]} = _mm_andnot_ps(_mm_cmplt_ps(vx${ABC[N:N+4]}, vdenorm_cutoff), vf${ABC[N:N+4]});
104
105    // Store ${ELEMENTS_TILE} (${SIMD_TILE}x4) outputs at a time.
106    _mm_storeu_ps(output, vf${ABC[0:4]});
107    $for N in range(4, ELEMENTS_TILE, 4):
108      _mm_storeu_ps(output + ${N}, vf${ABC[N:N+4]});
109    output += ${ELEMENTS_TILE};
110
111    // Accumulate computed exponents.
112    $for N in range(0, ELEMENTS_TILE, 4):
113      vacc${N % ACCUMULATORS} = _mm_add_ps(vacc${N % ACCUMULATORS}, vf${ABC[N:N+4]});
114  }
115  $if ACCUMULATORS > 1:
116    // Add up all accumulators to vacc0
117    $ACC_SLICE = 1
118    $while ACC_SLICE < ACCUMULATORS:
119      $for A in range(0, ACCUMULATORS, ACC_SLICE * 2):
120        $if A + ACC_SLICE < ACCUMULATORS:
121          vacc${A} = _mm_add_ps(vacc${A}, vacc${A + ACC_SLICE});
122      $ACC_SLICE *= 2
123
124  __m128 vacc = vacc0;
125  for (; elements >= 4 * sizeof(float); elements -= 4 * sizeof(float)) {
126    // Load 4 inputs at a time.
127    const __m128 vi = _mm_loadu_ps(input);
128    input += 4;
129
130    // Subtract maximum input x := i - i_max. This implies x <= 0.
131    const __m128 vx = _mm_sub_ps(vi, vi_max);
132
133    // Compute reduced argument elements := round(x / log(2)).
134    __m128 vn = _mm_add_ps(_mm_mul_ps(vx, vlog2e), vmagic_bias);
135
136    // Create a floating-point number s (scale) such that s == 2**elements for inputs which don't cause underflow, i.e.
137    // -87.33642 <= x <= 0.0, and -126 <= elements <= 0 accordingly.
138    const __m128 vs = _mm_castsi128_ps(_mm_slli_epi32(_mm_castps_si128(vn), 23));
139
140    // Subtract the large number back to get final elements := round(x / log(2)).
141    vn = _mm_sub_ps(vn, vmagic_bias);
142
143    // Compute reduced argument t := x - elements * log(2).
144    // Use Cody-Waite range reduction method (note two constants to represent log(2)) to improve accuracy.
145    __m128 vt = _mm_add_ps(_mm_mul_ps(vn, vminus_ln2_hi), vx);
146    vt = _mm_add_ps(_mm_mul_ps(vn, vminus_ln2_lo), vt);
147
148    // Compute degree-5 polynomial approxiatmion for exp(t) on [-log(2)/2, log(2)/2].
149    __m128 vp = _mm_add_ps(_mm_mul_ps(vc5, vt), vc4);
150    vp = _mm_add_ps(_mm_mul_ps(vp, vt), vc3);
151    vp = _mm_add_ps(_mm_mul_ps(vp, vt), vc2);
152    vp = _mm_add_ps(_mm_mul_ps(vp, vt), vc1);
153
154    // Reconstruct the final f value:
155    //   f = s * (1 + t * (c1 + t * (c2 + t * (c3 + t * (c4 + t * c5)))))
156    //     = s + (t * s) * (c1 + t * (c2 + t * (c3 + t * (c4 + t * c5))))
157    //     = s + (t * s) * p
158    vt = _mm_mul_ps(vt, vs);
159    __m128 vf = _mm_add_ps(_mm_mul_ps(vt, vp), vs);
160
161    // For inputs below zero cutoff, replace output with +0.0f.
162    // Note that for NaN inputs, comparison result is false, and outputs are left unchanged.
163    vf = _mm_andnot_ps(_mm_cmplt_ps(vx, vdenorm_cutoff), vf);
164
165    // Store 4 outputs at a time.
166    _mm_storeu_ps(output, vf);
167    output += 4;
168
169    // Accumulate computed exponents.
170    vacc = _mm_add_ps(vacc, vf);
171  }
172  if (elements != 0) {
173    assert(elements >= 1 * sizeof(float));
174    assert(elements <= 3 * sizeof(float));
175    // Load 4 inputs at a time.
176    const __m128 vi = _mm_loadu_ps(input);
177
178    // Subtract maximum input x := i - i_max. This implies x <= 0.
179    const __m128 vx = _mm_sub_ps(vi, vi_max);
180
181    // Compute reduced argument elements := round(x / log(2)).
182    __m128 vn = _mm_add_ps(_mm_mul_ps(vx, vlog2e), vmagic_bias);
183
184    // Create a floating-point number s (scale) such that s == 2**elements for inputs which don't cause underflow, i.e.
185    // -87.33642 <= x <= 0.0, and -126 <= elements <= 0 accordingly.
186    const __m128 vs = _mm_castsi128_ps(_mm_slli_epi32(_mm_castps_si128(vn), 23));
187
188    // Subtract the large number back to get final elements := round(x / log(2)).
189    vn = _mm_sub_ps(vn, vmagic_bias);
190
191    // Compute reduced argument t := x - elements * log(2).
192    // Use Cody-Waite range reduction method (note two constants to represent log(2)) to improve accuracy.
193    __m128 vt = _mm_add_ps(_mm_mul_ps(vn, vminus_ln2_hi), vx);
194    vt = _mm_add_ps(_mm_mul_ps(vn, vminus_ln2_lo), vt);
195
196    // Compute degree-5 polynomial approxiatmion for exp(t) on [-log(2)/2, log(2)/2].
197    __m128 vp = _mm_add_ps(_mm_mul_ps(vc5, vt), vc4);
198    vp = _mm_add_ps(_mm_mul_ps(vp, vt), vc3);
199    vp = _mm_add_ps(_mm_mul_ps(vp, vt), vc2);
200    vp = _mm_add_ps(_mm_mul_ps(vp, vt), vc1);
201
202    // Reconstruct the final f value:
203    //   f = s * (1 + t * (c1 + t * (c2 + t * (c3 + t * (c4 + t * c5)))))
204    //     = s + (t * s) * (c1 + t * (c2 + t * (c3 + t * (c4 + t * c5))))
205    //     = s + (t * s) * p
206    vt = _mm_mul_ps(vt, vs);
207    __m128 vf = _mm_add_ps(_mm_mul_ps(vt, vp), vs);
208
209    // For inputs below zero cutoff, replace output with +0.0f.
210    // Note that for NaN inputs, comparison result is false, and outputs are left unchanged.
211    vf = _mm_andnot_ps(_mm_cmplt_ps(vx, vdenorm_cutoff), vf);
212
213    if (elements & (2 * sizeof(float))) {
214      // Store 2 outputs at a time.
215      _mm_storel_pi((__m64*) output, vf);
216      output += 2;
217
218      // Accumulate 2 computed exponents.
219      vacc = _mm_add_ps(vacc, _mm_movelh_ps(vf, _mm_setzero_ps()));
220
221      vf = _mm_movehl_ps(vf, vf);
222    }
223    if (elements & (1 * sizeof(float))) {
224      // Store 1 output at a time.
225      _mm_store_ss(output, vf);
226
227      // Accumulate 1 computed exponent.
228      vacc = _mm_add_ss(vacc, vf);
229    }
230  }
231  // Reduce 4 elements in the SIMD register
232  vacc = _mm_add_ps(vacc, _mm_movehl_ps(vacc, vacc));
233  vacc = _mm_add_ss(vacc, _mm_shuffle_ps(vacc, vacc, _MM_SHUFFLE(2, 3, 0, 1)));
234  _mm_store_ss(sum, vacc);
235}
236