1// Copyright 2020 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 <wasm_simd128.h> 13 14#include <xnnpack/common.h> 15#include <xnnpack/raddstoreexpminusmax.h> 16 17 18void xnn_f32_raddstoreexpminusmax_ukernel__wasmsimd_rr2_p5_x${ELEMENTS_TILE}${"" if ACCUMULATORS == 1 else "_acc%d" % ACCUMULATORS}( 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 v128_t vi_max = wasm_v128_load32_splat(max); 29 const v128_t vlog2e = wasm_v128_load64_splat(params->wasmsimd_rr2_p5.log2e); 30 const v128_t vmagic_bias = wasm_v128_load64_splat(params->wasmsimd_rr2_p5.magic_bias); 31 const v128_t vminus_ln2_hi = wasm_v128_load64_splat(params->wasmsimd_rr2_p5.minus_ln2_hi); 32 const v128_t vminus_ln2_lo = wasm_v128_load64_splat(params->wasmsimd_rr2_p5.minus_ln2_lo); 33 const v128_t vc5 = wasm_v128_load64_splat(params->wasmsimd_rr2_p5.c5); 34 const v128_t vc4 = wasm_v128_load64_splat(params->wasmsimd_rr2_p5.c4); 35 const v128_t vc3 = wasm_v128_load64_splat(params->wasmsimd_rr2_p5.c3); 36 const v128_t vc2 = wasm_v128_load64_splat(params->wasmsimd_rr2_p5.c2); 37 const v128_t vc1 = wasm_v128_load64_splat(params->wasmsimd_rr2_p5.c1); 38 const v128_t vdenorm_cutoff = wasm_v128_load64_splat(params->wasmsimd_rr2_p5.denorm_cutoff); 39 40 v128_t vacc0 = wasm_f32x4_const_splat(0.0f); 41 $for K in range(1, ACCUMULATORS): 42 v128_t vacc${K} = vacc0; 43 for (; elements >= ${ELEMENTS_TILE} * sizeof(float); elements -= ${ELEMENTS_TILE} * sizeof(float)) { 44 // Load ${ELEMENTS_TILE} (${SIMD_TILE}x4) inputs at a time. 45 const v128_t vi${ABC[0:4]} = wasm_v128_load(input); 46 $for N in range(4, ELEMENTS_TILE, 4): 47 const v128_t vi${ABC[N:N+4]} = wasm_v128_load(input + ${N}); 48 input += ${ELEMENTS_TILE}; 49 50 // Subtract maximum input x := i - i_max. This implies x <= 0. 51 $for N in range(0, ELEMENTS_TILE, 4): 52 const v128_t vx${ABC[N:N+4]} = wasm_f32x4_sub(vi${ABC[N:N+4]}, vi_max); 53 54 // Compute reduced argument elements := round(x / log(2)). 55 $for N in range(0, ELEMENTS_TILE, 4): 56 v128_t vn${ABC[N:N+4]} = wasm_f32x4_add(vmagic_bias, wasm_f32x4_mul(vx${ABC[N:N+4]}, vlog2e)); 57 58 // Create a floating-point number s (scale) such that s == 2**elements for inputs which don't cause underflow, i.e. 59 // -87.33642 <= x <= 0.0, and -126 <= elements <= 0 accordingly. 60 $for N in range(0, ELEMENTS_TILE, 4): 61 const v128_t vs${ABC[N:N+4]} = wasm_i32x4_shl(vn${ABC[N:N+4]}, 23); 62 63 // Subtract the large number back to get final elements := round(x / log(2)). 64 $for N in range(0, ELEMENTS_TILE, 4): 65 vn${ABC[N:N+4]} = wasm_f32x4_sub(vn${ABC[N:N+4]}, vmagic_bias); 66 67 // Compute reduced argument t := x - elements * log(2). 68 // Use Cody-Waite range reduction method (note two constants to represent log(2)) to improve accuracy. 69 $for N in range(0, ELEMENTS_TILE, 4): 70 v128_t vt${ABC[N:N+4]} = wasm_f32x4_add(vx${ABC[N:N+4]}, wasm_f32x4_mul(vn${ABC[N:N+4]}, vminus_ln2_hi)); 71 72 $for N in range(0, ELEMENTS_TILE, 4): 73 vt${ABC[N:N+4]} = wasm_f32x4_add(vt${ABC[N:N+4]}, wasm_f32x4_mul(vn${ABC[N:N+4]}, vminus_ln2_lo)); 74 75 // Compute degree-5 polynomial approximation for exp(t) on [-log(2)/2, log(2)/2]. 76 $for N in range(0, ELEMENTS_TILE, 4): 77 v128_t vp${ABC[N:N+4]} = wasm_f32x4_add(vc4, wasm_f32x4_mul(vc5, vt${ABC[N:N+4]})); 78 79 $for N in range(0, ELEMENTS_TILE, 4): 80 vp${ABC[N:N+4]} = wasm_f32x4_add(vc3, wasm_f32x4_mul(vp${ABC[N:N+4]}, vt${ABC[N:N+4]})); 81 82 $for N in range(0, ELEMENTS_TILE, 4): 83 vp${ABC[N:N+4]} = wasm_f32x4_add(vc2, wasm_f32x4_mul(vp${ABC[N:N+4]}, vt${ABC[N:N+4]})); 84 85 $for N in range(0, ELEMENTS_TILE, 4): 86 vp${ABC[N:N+4]} = wasm_f32x4_add(vc1, wasm_f32x4_mul(vp${ABC[N:N+4]}, vt${ABC[N:N+4]})); 87 88 // Reconstruct the final f value: 89 // f = s * (1 + t * (c1 + t * (c2 + t * (c3 + t * (c4 + t * c5))))) 90 // = s + (t * s) * (c1 + t * (c2 + t * (c3 + t * (c4 + t * c5)))) 91 // = s + (t * s) * p 92 $for N in range(0, ELEMENTS_TILE, 4): 93 vt${ABC[N:N+4]} = wasm_f32x4_mul(vt${ABC[N:N+4]}, vs${ABC[N:N+4]}); 94 95 $for N in range(0, ELEMENTS_TILE, 4): 96 v128_t vf${ABC[N:N+4]} = wasm_f32x4_add(vs${ABC[N:N+4]}, wasm_f32x4_mul(vt${ABC[N:N+4]}, vp${ABC[N:N+4]})); 97 98 // For inputs below zero cutoff, replace output with +0.0f. 99 // Note that for NaN inputs, comparison result is false, and outputs are left unchanged. 100 $for N in range(0, ELEMENTS_TILE, 4): 101 vf${ABC[N:N+4]} = wasm_v128_andnot(vf${ABC[N:N+4]}, wasm_f32x4_lt(vx${ABC[N:N+4]}, vdenorm_cutoff)); 102 103 // Store ${ELEMENTS_TILE} (${SIMD_TILE}x4) outputs at a time. 104 wasm_v128_store(output, vf${ABC[0:4]}); 105 $for N in range(4, ELEMENTS_TILE, 4): 106 wasm_v128_store(output + ${N}, vf${ABC[N:N+4]}); 107 output += ${ELEMENTS_TILE}; 108 109 // Accumulate computed exponents. 110 $for N in range(0, ELEMENTS_TILE, 4): 111 vacc${N % ACCUMULATORS} = wasm_f32x4_add(vacc${N % ACCUMULATORS}, vf${ABC[N:N+4]}); 112 } 113 $if ACCUMULATORS > 1: 114 // Add up all accumulators to vacc0 115 $ACC_SLICE = 1 116 $while ACC_SLICE < ACCUMULATORS: 117 $for A in range(0, ACCUMULATORS, ACC_SLICE * 2): 118 $if A + ACC_SLICE < ACCUMULATORS: 119 vacc${A} = wasm_f32x4_add(vacc${A}, vacc${A + ACC_SLICE}); 120 $ACC_SLICE *= 2 121 122 v128_t vacc = vacc0; 123 for (; elements >= 4 * sizeof(float); elements -= 4 * sizeof(float)) { 124 // Load 4 inputs at a time. 125 const v128_t vi = wasm_v128_load(input); 126 input += 4; 127 128 // Subtract maximum input x := i - i_max. This implies x <= 0. 129 const v128_t vx = wasm_f32x4_sub(vi, vi_max); 130 131 // Compute reduced argument elements := round(x / log(2)). 132 v128_t vn = wasm_f32x4_add(vmagic_bias, wasm_f32x4_mul(vx, vlog2e)); 133 134 // Create a floating-point number s (scale) such that s == 2**elements for inputs which don't cause underflow, i.e. 135 // -87.33642 <= x <= 0.0, and -126 <= elements <= 0 accordingly. 136 const v128_t vs = wasm_i32x4_shl(vn, 23); 137 138 // Subtract the large number back to get final elements := round(x / log(2)). 139 vn = wasm_f32x4_sub(vn, vmagic_bias); 140 141 // Compute reduced argument t := x - elements * log(2). 142 // Use Cody-Waite range reduction method (note two constants to represent log(2)) to improve accuracy. 143 v128_t vt = wasm_f32x4_add(vx, wasm_f32x4_mul(vn, vminus_ln2_hi)); 144 vt = wasm_f32x4_add(vt, wasm_f32x4_mul(vn, vminus_ln2_lo)); 145 146 // Compute degree-5 polynomial approximation for exp(t) on [-log(2)/2, log(2)/2]. 147 v128_t vp = wasm_f32x4_add(vc4, wasm_f32x4_mul(vc5, vt)); 148 vp = wasm_f32x4_add(vc3, wasm_f32x4_mul(vp, vt)); 149 vp = wasm_f32x4_add(vc2, wasm_f32x4_mul(vp, vt)); 150 vp = wasm_f32x4_add(vc1, wasm_f32x4_mul(vp, vt)); 151 152 // Reconstruct the final f value: 153 // f = s * (1 + t * (c1 + t * (c2 + t * (c3 + t * (c4 + t * c5))))) 154 // = s + (t * s) * (c1 + t * (c2 + t * (c3 + t * (c4 + t * c5)))) 155 // = s + (t * s) * p 156 vt = wasm_f32x4_mul(vt, vs); 157 v128_t vf = wasm_f32x4_add(vs, wasm_f32x4_mul(vt, vp)); 158 159 // For inputs below zero cutoff, replace output with +0.0f. 160 // Note that for NaN inputs, comparison result is false, and outputs are left unchanged. 161 vf = wasm_v128_andnot(vf, wasm_f32x4_lt(vx, vdenorm_cutoff)); 162 163 // Store 4 outputs at a time. 164 wasm_v128_store(output, vf); 165 output += 4; 166 167 // Accumulate computed exponents. 168 vacc = wasm_f32x4_add(vacc, vf); 169 } 170 vacc = wasm_f32x4_add(vacc, wasm_v32x4_shuffle(vacc, vacc, 2, 3, 2, 3)); 171 float vsum = wasm_f32x4_extract_lane(vacc, 0) + wasm_f32x4_extract_lane(vacc, 1); 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 v128_t vi = wasm_v128_load(input); 177 178 // Subtract maximum input x := i - i_max. This implies x <= 0. 179 const v128_t vx = wasm_f32x4_sub(vi, vi_max); 180 181 // Compute reduced argument elements := round(x / log(2)). 182 v128_t vn = wasm_f32x4_add(vmagic_bias, wasm_f32x4_mul(vx, vlog2e)); 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 v128_t vs = wasm_i32x4_shl(vn, 23); 187 188 // Subtract the large number back to get final elements := round(x / log(2)). 189 vn = wasm_f32x4_sub(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 v128_t vt = wasm_f32x4_add(vx, wasm_f32x4_mul(vn, vminus_ln2_hi)); 194 vt = wasm_f32x4_add(vt, wasm_f32x4_mul(vn, vminus_ln2_lo)); 195 196 // Compute degree-5 polynomial approximation for exp(t) on [-log(2)/2, log(2)/2]. 197 v128_t vp = wasm_f32x4_add(vc4, wasm_f32x4_mul(vc5, vt)); 198 vp = wasm_f32x4_add(vc3, wasm_f32x4_mul(vp, vt)); 199 vp = wasm_f32x4_add(vc2, wasm_f32x4_mul(vp, vt)); 200 vp = wasm_f32x4_add(vc1, wasm_f32x4_mul(vp, vt)); 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 = wasm_f32x4_mul(vt, vs); 207 v128_t vf = wasm_f32x4_add(vs, wasm_f32x4_mul(vt, vp)); 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 = wasm_v128_andnot(vf, wasm_f32x4_lt(vx, vdenorm_cutoff)); 212 213 if (elements & (2 * sizeof(float))) { 214 // Store and accumulate 2 outputs at a time. 215 const float vf0 = wasm_f32x4_extract_lane(vf, 0); 216 output[0] = vf0; 217 vsum += vf0; 218 219 const float vf1 = wasm_f32x4_extract_lane(vf, 1); 220 output[1] = vf1; 221 vsum += vf1; 222 223 vf = wasm_v32x4_shuffle(vf, vf, 2, 3, 2, 3); 224 output += 2; 225 } 226 if (elements & (1 * sizeof(float))) { 227 // Store 1 output at a time. 228 const float vf0 = wasm_f32x4_extract_lane(vf, 0); 229 *output = vf0; 230 vsum += vf0; 231 } 232 } 233 // Reduce 4 elements in the SIMD register 234 *sum = vsum; 235} 236