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 % 16 == 0 7$assert ELEMENTS_TILE >= 16 8$SIMD_TILE = ELEMENTS_TILE // 16 9$ABC = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ" 10#include <assert.h> 11#include <math.h> 12 13#include <immintrin.h> 14 15#include <xnnpack/common.h> 16#include <xnnpack/intrinsics-polyfill.h> 17#include <xnnpack/raddextexp.h> 18 19 20void xnn_f32_raddextexp_ukernel__avx512f_p5_scalef_x${ELEMENTS_TILE}${"" if ACCUMULATORS == 1 else "_acc%d" % ACCUMULATORS}( 21 size_t elements, 22 const float* x, 23 float* sum) 24{ 25 assert(elements % sizeof(float) == 0); 26 27 const __m512 vlog2e = _mm512_set1_ps(0x1.715476p+0f); 28 const __m512 vminus_ln2_hi = _mm512_set1_ps(-0x1.62E43p-1f); 29 const __m512 vminus_ln2_lo = _mm512_set1_ps(0x1.05C61p-29f); 30 31 const __m512 vc0 = _mm512_set1_ps(1.0f); 32 const __m512 vc1 = _mm512_set1_ps(0x1.FFFFF6p-1f); 33 const __m512 vc2 = _mm512_set1_ps(0x1.FFFDC6p-2f); 34 const __m512 vc3 = _mm512_set1_ps(0x1.555A80p-3f); 35 const __m512 vc4 = _mm512_set1_ps(0x1.573A1Ap-5f); 36 const __m512 vc5 = _mm512_set1_ps(0x1.0F9F9Cp-7f); 37 38 const __m512 vminus_inf = _mm512_set1_ps(-INFINITY); 39 40 $for K in range(ACCUMULATORS): 41 __m512 vaccv${K} = _mm512_setzero_ps(); 42 $for K in range(ACCUMULATORS): 43 __m512 vacce${K} = vminus_inf; 44 for (; elements >= ${ELEMENTS_TILE} * sizeof(float); elements -= ${ELEMENTS_TILE} * sizeof(float)) { 45 // Load ${ELEMENTS_TILE} (${SIMD_TILE}x16) inputs at a time. 46 const __m512 vx0 = _mm512_loadu_ps(x); 47 $for N in range(1, SIMD_TILE): 48 const __m512 vx${N} = _mm512_loadu_ps(x + ${N * 16}); 49 x += ${ELEMENTS_TILE}; 50 51 // Compute reduced argument elements := round(x / log(2)). 52 $for N in range(SIMD_TILE): 53 const __m512 vn${N} = _mm512_roundscale_ps(_mm512_mul_ps(vx${N}, vlog2e), 0); 54 55 // Compute reduced argument t := x - elements * log(2). 56 // Use Cody-Waite range reduction method (note two constants to represent log(2)) to improve accuracy. 57 $for N in range(SIMD_TILE): 58 __m512 vt${N} = _mm512_fmadd_ps(vn${N}, vminus_ln2_hi, vx${N}); 59 60 $for N in range(SIMD_TILE): 61 vt${N} = _mm512_fmadd_ps(vn${N}, vminus_ln2_lo, vt${N}); 62 63 // Compute degree-5 polynomial approxiatmion for exp(t) on [-log(2)/2, log(2)/2]. 64 $for N in range(SIMD_TILE): 65 __m512 vp${N} = _mm512_fmadd_ps(vc5, vt${N}, vc4); 66 67 $for N in range(SIMD_TILE): 68 vp${N} = _mm512_fmadd_ps(vp${N}, vt${N}, vc3); 69 70 $for N in range(SIMD_TILE): 71 vp${N} = _mm512_fmadd_ps(vp${N}, vt${N}, vc2); 72 73 $for N in range(SIMD_TILE): 74 vp${N} = _mm512_fmadd_ps(vp${N}, vt${N}, vc1); 75 76 $for N in range(SIMD_TILE): 77 vp${N} = _mm512_fmadd_ps(vp${N}, vt${N}, vc0); 78 79 // Accumulate "extended" floating-point numbers in ("mantissa", "exponent") representation where 80 // - vnX is "exponent" 81 // - vpX is "mantissa" 82 // 83 // exp2(ae) * av + exp2(be) * bv = 84 // = exp2(max(ae, be)) * exp2(ae - max(ae, be)) * av + exp2(max(ae, be)) * exp2(be - max(ae, be)) * bv 85 // = exp2(max_e) * (exp2(ae - max_e) * av + exp2(be - max_e) * bv) 86 // = exp2(max_e) * (exp2(delta_ae) * av + exp2(delta_be) * bv) 87 // 88 // For computational efficiency we add three "extended" floating-point numbers at a time. 89 $for N in range(SIMD_TILE): 90 $if N < ACCUMULATORS: 91 __m512 vmax_e${N} = _mm512_max_ps(vacce${N}, vn${N}); 92 $else: 93 vmax_e${N % ACCUMULATORS} = _mm512_max_ps(vmax_e${N % ACCUMULATORS}, vn${N}); 94 95 $for K in range(ACCUMULATORS): 96 const __m512 vdelta_acce${K} = _mm512_sub_ps(vacce${K}, vmax_e${K}); 97 $for N in range(SIMD_TILE): 98 const __m512 vdelta_e${N} = _mm512_sub_ps(vn${N}, vmax_e${N % ACCUMULATORS}); 99 100 // Update accumulated "mantissa" and "exponent" values 101 $for K in range(ACCUMULATORS): 102 vaccv${K} = _mm512_scalef_ps(vaccv${K}, vdelta_acce${K}); 103 $for N in range(SIMD_TILE): 104 vaccv${N % ACCUMULATORS} = _mm512_add_ps(vaccv${N % ACCUMULATORS}, _mm512_scalef_ps(vp${N}, vdelta_e${N})); 105 106 $for K in range(ACCUMULATORS): 107 vacce${K} = vmax_e${K}; 108 } 109 110 // Reduce partial sums of "extended" floating-point numbers into a single "extended" SIMD vector of sums. 111 $if ACCUMULATORS > 1: 112 $for A in range(0, ACCUMULATORS, 2): 113 $if A + 1 < ACCUMULATORS: 114 const __m512 vmax_acce${ABC[A:A+2]} = _mm512_max_ps(vacce${A}, vacce${A+1}); 115 $else: 116 const __m512 vmax_acce${ABC[A]} = vacce${A}; 117 $ACC_SLICE = 2 118 $while ACC_SLICE < ACCUMULATORS: 119 $for A in range(0, ACCUMULATORS, ACC_SLICE * 2): 120 $if A + ACC_SLICE < ACCUMULATORS: 121 const __m512 vmax_acce${ABC[A:min(A+ACC_SLICE*2, ACCUMULATORS)]} = _mm512_max_ps(vmax_acce${ABC[A:A+ACC_SLICE]}, vmax_acce${ABC[A+ACC_SLICE:min(ACCUMULATORS,A+ACC_SLICE*2)]}); 122 $ACC_SLICE *= 2 123 124 $for K in range(ACCUMULATORS): 125 const __m512 vdelta_acce${K} = _mm512_sub_ps(vacce${K}, vmax_acce${ABC[0:ACCUMULATORS]}); 126 127 __m512 vaccv = _mm512_scalef_ps(vaccv0, vdelta_acce0); 128 $for K in range(1, ACCUMULATORS): 129 vaccv = _mm512_add_ps(vaccv, _mm512_scalef_ps(vaccv${K}, vdelta_acce${K})); 130 __m512 vacce = vmax_acce${ABC[0:ACCUMULATORS]}; 131 $else: 132 __m512 vaccv = vaccv0; 133 __m512 vacce = vacce0; 134 135 for (; elements >= 16 * sizeof(float); elements -= 16 * sizeof(float)) { 136 // Load 16 inputs at a time. 137 const __m512 vx = _mm512_loadu_ps(x); 138 x += 16; 139 140 // Compute reduced argument elements := round(x / log(2)). 141 const __m512 vn = _mm512_roundscale_ps(_mm512_mul_ps(vx, vlog2e), 0); 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 __m512 vt = _mm512_fmadd_ps(vn, vminus_ln2_hi, vx); 146 vt = _mm512_fmadd_ps(vn, vminus_ln2_lo, vt); 147 148 // Compute degree-5 polynomial approxiatmion for exp(t) on [-log(2)/2, log(2)/2]. 149 __m512 vp = _mm512_fmadd_ps(vc5, vt, vc4); 150 vp = _mm512_fmadd_ps(vp, vt, vc3); 151 vp = _mm512_fmadd_ps(vp, vt, vc2); 152 vp = _mm512_fmadd_ps(vp, vt, vc1); 153 vp = _mm512_fmadd_ps(vp, vt, vc0); 154 155 // Accumulate "extended" floating-point numbers in ("mantissa", "exponent") representation. 156 const __m512 vmax_e = _mm512_max_ps(vacce, vn); 157 const __m512 vdelta_acce = _mm512_sub_ps(vacce, vmax_e); 158 const __m512 vdelta_e = _mm512_sub_ps(vn, vmax_e); 159 vaccv = _mm512_scalef_ps(vaccv, vdelta_acce); 160 vaccv = _mm512_add_ps(vaccv, _mm512_scalef_ps(vp, vdelta_e)); 161 162 vacce = vmax_e; 163 } 164 if XNN_UNLIKELY(elements != 0) { 165 // Prepare mask for valid 32-bit elements (depends on elements). 166 elements >>= 2 /* log2(sizeof(float)) */; 167 const __mmask16 vmask = _cvtu32_mask16((uint16_t) ((uint32_t) (UINT32_C(1) << elements) - UINT32_C(1))); 168 169 // Load up to 15 inputs at a time. 170 const __m512 vx = _mm512_maskz_loadu_ps(vmask, x); 171 172 // Compute reduced argument elements := round(x / log(2)). 173 const __m512 vn = _mm512_roundscale_ps(_mm512_mul_ps(vx, vlog2e), 0); 174 175 // Compute reduced argument t := x - elements * log(2). 176 // Use Cody-Waite range reduction method (note two constants to represent log(2)) to improve accuracy. 177 __m512 vt = _mm512_fmadd_ps(vn, vminus_ln2_hi, vx); 178 vt = _mm512_fmadd_ps(vn, vminus_ln2_lo, vt); 179 180 // Compute degree-5 polynomial approxiatmion for exp(t) on [-log(2)/2, log(2)/2]. 181 __m512 vp = _mm512_fmadd_ps(vc5, vt, vc4); 182 vp = _mm512_fmadd_ps(vp, vt, vc3); 183 vp = _mm512_fmadd_ps(vp, vt, vc2); 184 vp = _mm512_fmadd_ps(vp, vt, vc1); 185 vp = _mm512_fmadd_ps(vp, vt, vc0); 186 187 // Accumulate "extended" floating-point numbers in ("mantissa", "exponent") representation. 188 const __m512 vmax_e = _mm512_mask_max_ps(vacce, vmask, vacce, vn); 189 const __m512 vdelta_acce = _mm512_sub_ps(vacce, vmax_e); 190 const __m512 vdelta_e = _mm512_sub_ps(vn, vmax_e); 191 vaccv = _mm512_mask_scalef_ps(vaccv, vmask, vaccv, vdelta_acce); 192 vaccv = _mm512_mask_add_ps(vaccv, vmask, vaccv, _mm512_maskz_scalef_ps(vmask, vp, vdelta_e)); 193 vacce = vmax_e; 194 } 195 196 // Reduce partial sums of "extended" floating-point numbers into a single "extended" floating-point sum. 197 const float vmax_acce = _mm512_reduce_max_ps(vacce); 198 const __m512 vdelta_acce = _mm512_sub_ps(vacce, _mm512_set1_ps(vmax_acce)); 199 200 sum[0] = _mm512_reduce_add_ps(_mm512_scalef_ps(vaccv, vdelta_acce)); 201 sum[1] = vmax_acce; 202 203 _mm256_zeroupper(); 204} 205