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) XNN_DISABLE_TSAN 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 approximation 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 approximation 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 approximation 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