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 BATCH_TILE >= 1 7$ABC = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ" 8#include <assert.h> 9#include <math.h> 10 11#include <xnnpack/common.h> 12#include <xnnpack/vunary.h> 13 14#include <fp16/bitcasts.h> 15 16 17// Note redefine as uint32[] to avoid redundant bitcasts. 18extern XNN_INTERNAL const uint32_t xnn_table_exp2_k_over_2048[2048]; 19 20void xnn_f32_sigmoid_ukernel__scalar_lut2048_p1_div_x${BATCH_TILE}( 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 float vmagic_bias = 0x1.800000p23f; 29 // The largest z for which sigmoidf(-z) is normalized. 30 // This number is also the largest z for which expf(-z) is normalized. 31 const float vdenorm_cutoff = 0x1.5D589Ep+6f; 32 const float vminus_log2e_x2048 = -0x1.715476p11f; 33 // Last 18 bits are zeroes 34 const float vln2_o2048_hi = 0x1.600000p-12f; 35 const float vln2_o2048_lo = 0x1.7217F8p-19f; 36 const float vone = 1.0f; 37 38 const float vc1 = -0x1.FFFFFEp-1f; 39 40 const uint32_t vindex_mask = UINT32_C(0x7FF); 41 42 $if BATCH_TILE > 1: 43 for (; n >= ${BATCH_TILE} * sizeof(float); n -= ${BATCH_TILE} * sizeof(float)) { 44 $for N in range(BATCH_TILE): 45 const float vx${N} = x[${N}]; 46 x += ${BATCH_TILE}; 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 $for N in range(BATCH_TILE): 56 const float vz${N} = fabsf(vx${N}); 57 58 // Compute reduced argument n := round(-z * 2048 / log(2)). 59 // We do it by adding a large number (magic bias), which cause rounding of the result to an integer, then subtracing 60 // the large number back. The first addition is combined with multiplication by log2e into a single FMA instruction. 61 // The trick with adding large number is valid only within certain bounds (|z * 2048 / log(2)| <= 2**22, i.e. 62 // |z| <= 0x1.62E43p+10 = 1419.5654296875), but that is acceptable, because inputs x outside of 63 // [-87.336544, 17.328678] (i.e. z outsize [0, 87.336544]) underflow or saturate sigmoidf(x). We fixup the result 64 // for such inputs at the very end of the algorithm. 65 $for N in range(BATCH_TILE): 66 float vn${N} = vz${N} * vminus_log2e_x2048 + vmagic_bias; 67 68 // Create a floating-point number s (scale) such that s := 2**(n / 2048) for such inputs that sigmoidf(-z) is 69 // normalized, i.e. 0 <= z <= 87.33642. As n has 11 fractional bits, we split s == 2**(n / 2048) = 70 // = 2**e * 2**(n / 2048 - e), where e := int(n / 2048). We create s in two steps: 71 // 1. Fetch 2**(n / 2048 - e) = 2**(n % 2048) from table using the 6 low bits of n, as integer. 72 // Note that the fetched values are in the [1.0, 2.0) range, i.e. their floating-point exponent is 0. 73 // 2. Adjust fecthed value by addition of e to its floating-point exponent. The result is always a normalized 74 // number, because for 0 <= z <= 87.33642 (inputs for which sigmoidf(-z) is normalized) we have -126 <= e <= 0, 75 // and thus the adjusted exponent is not lower than -126. 76 // 77 // Extract e from bits 11:19 of n and shift it into bits 23:31 (position of floating-point exponent). 78 $for N in range(BATCH_TILE): 79 const uint32_t ve${N} = (fp32_to_bits(vn${N}) & ~vindex_mask) << 12; 80 81 // Use bits 0:11 bits of n, as integer, as an index for table lookup of l := 2**(n % 2048). 82 $for N in range(BATCH_TILE): 83 const uint32_t vidx${N} = fp32_to_bits(vn${N}) & vindex_mask; 84 // Adjust exponent of the value l fetched from the table to get the final s value. 85 $for N in range(BATCH_TILE): 86 const float vs${N} = fp32_from_bits(xnn_table_exp2_k_over_2048[vidx${N}] + ve${N}); 87 88 // Subtract the large number back to get the final n := round(-z * 2048 / log(2)) as a floating-point number. 89 $for N in range(BATCH_TILE): 90 vn${N} -= vmagic_bias; 91 92 // Compute reduced argument t := (z + n * log(2) / 2048). Note that -t = -z - n * log(2) / 2048. 93 // Use Cody-Waite range reduction method (note two constants to represent log(2) / 2048) to improve accuracy. 94 $for N in range(BATCH_TILE): 95 float vt${N} = vn${N} * vln2_o2048_hi + vz${N}; 96 97 $for N in range(BATCH_TILE): 98 vt${N} = vn${N} * vln2_o2048_lo + vt${N}; 99 100 // Compute degree-1 polynomial approximation for exp(-t) on [-log(2)/4096, log(2)/4096]: 101 // P1(t) = 1 + t * c1 102 $for N in range(BATCH_TILE): 103 const float vp${N} = vt${N} * vc1; 104 105 // Reconstruct the exp(-z) value: 106 // y = s * (1 + t * c1) 107 // = s + s * (t * c1)) 108 // = s + s * p 109 $for N in range(BATCH_TILE): 110 const float vy${N} = vp${N} * vs${N} + vs${N}; 111 112 // Reconstruct sigmoid(-z) = exp(-z) / (1.0 + exp(-z)) 113 $for N in range(BATCH_TILE): 114 float vf${N} = vy${N} / (vy${N} + vone); 115 116 // For inputs above denormal cutoff, replace output with +0.0f. 117 // Note that for NaN inputs, comparison result is false, and outputs are left unchanged. 118 $for N in range(BATCH_TILE): 119 if XNN_UNPREDICTABLE(vz${N} > vdenorm_cutoff) { 120 vf${N} = 0.0f; 121 } 122 123 // Reconstruct sigmoid(x) = x < 0 ? sigmoid(-z) : 1.0 - sigmoid(-z) 124 $for N in range(BATCH_TILE): 125 if XNN_UNPREDICTABLE(vx${N} > 0.0f) { 126 vf${N} = vone - vf${N}; 127 } 128 129 $for N in range(BATCH_TILE): 130 y[${N}] = vf${N}; 131 y += ${BATCH_TILE}; 132 } 133 $if BATCH_TILE == 1: 134 do { 135 const float vx = *x++; 136 137 // General structure of the algorithm: 138 // / exp(x) / (1 + exp(x)) if x <= 0 139 // f[x] := 140 // \ 1 - f[-x] if x >= 0 141 // 142 // First we compute f[-z] := exp(-z) / (1 + exp(-z)) where z = abs(x), 143 // then replace result with 1 - f[-z] if x >= 0. 144 const float vz = fabsf(vx); 145 146 // Compute reduced argument n := round(-z * 2048 / log(2)). 147 // We do it by adding a large number (magic bias), which cause rounding of the result to an integer, then subtracing 148 // the large number back. The first addition is combined with multiplication by log2e into a single FMA instruction. 149 // The trick with adding large number is valid only within certain bounds (|z * 2048 / log(2)| <= 2**22, i.e. 150 // |z| <= 0x1.62E43p+10 = 1419.5654296875), but that is acceptable, because inputs x outside of 151 // [-87.336544, 17.328678] (i.e. z outsize [0, 87.336544]) underflow or saturate sigmoidf(x). We fixup the result 152 // for such inputs at the very end of the algorithm. 153 float vn = vz * vminus_log2e_x2048 + vmagic_bias; 154 155 // Create a floating-point number s (scale) such that s := 2**(n / 2048) for such inputs that sigmoidf(-z) is 156 // normalized, i.e. 0 <= z <= 87.33642. As n has 11 fractional bits, we split s == 2**(n / 2048) = 157 // = 2**e * 2**(n / 2048 - e), where e := int(n / 2048). We create s in two steps: 158 // 1. Fetch 2**(n / 2048 - e) = 2**(n % 2048) from table using the 6 low bits of n, as integer. 159 // Note that the fetched values are in the [1.0, 2.0) range, i.e. their floating-point exponent is 0. 160 // 2. Adjust fecthed value by addition of e to its floating-point exponent. The result is always a normalized 161 // number, because for 0 <= z <= 87.33642 (inputs for which sigmoidf(-z) is normalized) we have -126 <= e <= 0, 162 // and thus the adjusted exponent is not lower than -126. 163 // 164 // Extract e from bits 11:19 of n and shift it into bits 23:31 (position of floating-point exponent). 165 const uint32_t ve = (fp32_to_bits(vn) & ~vindex_mask) << 12; 166 167 // Use bits 0:11 bits of n, as integer, as an index for table lookup of l := 2**(n % 2048). 168 const uint32_t vidx = fp32_to_bits(vn) & vindex_mask; 169 // Adjust exponent of the value l fetched from the table to get the final s value. 170 const float vs = fp32_from_bits(xnn_table_exp2_k_over_2048[vidx] + ve); 171 172 // Subtract the large number back to get the final n := round(-z * 2048 / log(2)) as a floating-point number. 173 vn -= vmagic_bias; 174 175 // Compute reduced argument t := (z + n * log(2) / 2048). Note that -t = -z - n * log(2) / 2048. 176 // Use Cody-Waite range reduction method (note two constants to represent log(2) / 2048) to improve accuracy. 177 float vt = vn * vln2_o2048_hi + vz; 178 vt = vn * vln2_o2048_lo + vt; 179 180 // Compute degree-1 polynomial approximation for exp(-t) on [-log(2)/4096, log(2)/4096]: 181 // P1(t) = 1 + t * c1 182 const float vp = vt * vc1; 183 184 // Reconstruct the exp(-z) value: 185 // y = s * (1 + t * c1) 186 // = s + s * (t * c1)) 187 // = s + s * p 188 const float vy = vp * vs + vs; 189 190 // Reconstruct sigmoid(-z) = exp(-z) / (1.0 + exp(-z)) 191 float vf = vy / (vy + vone); 192 193 // For inputs above denormal cutoff, replace output with +0.0f. 194 // Note that for NaN inputs, comparison result is false, and outputs are left unchanged. 195 if XNN_UNPREDICTABLE(vz > vdenorm_cutoff) { 196 vf = 0.0f; 197 } 198 199 // Reconstruct sigmoid(x) = x < 0 ? sigmoid(-z) : 1.0 - sigmoid(-z) 200 if XNN_UNPREDICTABLE(vx > 0.0f) { 201 vf = vone - vf; 202 } 203 204 *y++ = vf; 205 206 n -= sizeof(float); 207 } while (n != 0); 208 $elif BATCH_TILE == 2: 209 if XNN_UNLIKELY(n != 0) { 210 const float vx = *x; 211 212 // General structure of the algorithm: 213 // / exp(x) / (1 + exp(x)) if x <= 0 214 // f[x] := 215 // \ 1 - f[-x] if x >= 0 216 // 217 // First we compute f[-z] := exp(-z) / (1 + exp(-z)) where z = abs(x), 218 // then replace result with 1 - f[-z] if x >= 0. 219 const float vz = fabsf(vx); 220 221 // Compute reduced argument n := round(-z * 2048 / log(2)). 222 // We do it by adding a large number (magic bias), which cause rounding of the result to an integer, then subtracing 223 // the large number back. The first addition is combined with multiplication by log2e into a single FMA instruction. 224 // The trick with adding large number is valid only within certain bounds (|z * 2048 / log(2)| <= 2**22, i.e. 225 // |z| <= 0x1.62E43p+10 = 1419.5654296875), but that is acceptable, because inputs x outside of 226 // [-87.336544, 17.328678] (i.e. z outsize [0, 87.336544]) underflow or saturate sigmoidf(x). We fixup the result 227 // for such inputs at the very end of the algorithm. 228 float vn = vz * vminus_log2e_x2048 + vmagic_bias; 229 230 // Create a floating-point number s (scale) such that s := 2**(n / 2048) for such inputs that sigmoidf(-z) is 231 // normalized, i.e. 0 <= z <= 87.33642. As n has 11 fractional bits, we split s == 2**(n / 2048) = 232 // = 2**e * 2**(n / 2048 - e), where e := int(n / 2048). We create s in two steps: 233 // 1. Fetch 2**(n / 2048 - e) = 2**(n % 2048) from table using the 6 low bits of n, as integer. 234 // Note that the fetched values are in the [1.0, 2.0) range, i.e. their floating-point exponent is 0. 235 // 2. Adjust fecthed value by addition of e to its floating-point exponent. The result is always a normalized 236 // number, because for 0 <= z <= 87.33642 (inputs for which sigmoidf(-z) is normalized) we have -126 <= e <= 0, 237 // and thus the adjusted exponent is not lower than -126. 238 // 239 // Extract e from bits 11:19 of n and shift it into bits 23:31 (position of floating-point exponent). 240 const uint32_t ve = (fp32_to_bits(vn) & ~vindex_mask) << 12; 241 242 // Use bits 0:11 bits of n, as integer, as an index for table lookup of l := 2**(n % 2048). 243 const uint32_t vidx = fp32_to_bits(vn) & vindex_mask; 244 // Adjust exponent of the value l fetched from the table to get the final s value. 245 const float vs = fp32_from_bits(xnn_table_exp2_k_over_2048[vidx] + ve); 246 247 // Subtract the large number back to get the final n := round(-z * 2048 / log(2)) as a floating-point number. 248 vn -= vmagic_bias; 249 250 // Compute reduced argument t := (z + n * log(2) / 2048). Note that -t = -z - n * log(2) / 2048. 251 // Use Cody-Waite range reduction method (note two constants to represent log(2) / 2048) to improve accuracy. 252 float vt = vn * vln2_o2048_hi + vz; 253 vt = vn * vln2_o2048_lo + vt; 254 255 // Compute degree-1 polynomial approximation for exp(-t) on [-log(2)/4096, log(2)/4096]: 256 // P1(t) = 1 + t * c1 257 const float vp = vt * vc1; 258 259 // Reconstruct the exp(-z) value: 260 // y = s * (1 + t * c1) 261 // = s + s * (t * c1)) 262 // = s + s * p 263 const float vy = vp * vs + vs; 264 265 // Reconstruct sigmoid(-z) = exp(-z) / (1.0 + exp(-z)) 266 float vf = vy / (vy + vone); 267 268 // For inputs above denormal cutoff, replace output with +0.0f. 269 // Note that for NaN inputs, comparison result is false, and outputs are left unchanged. 270 if XNN_UNPREDICTABLE(vz > vdenorm_cutoff) { 271 vf = 0.0f; 272 } 273 274 // Reconstruct sigmoid(x) = x < 0 ? sigmoid(-z) : 1.0 - sigmoid(-z) 275 if XNN_UNPREDICTABLE(vx > 0.0f) { 276 vf = vone - vf; 277 } 278 279 *y = vf; 280 } 281 $else: 282 if XNN_UNLIKELY(n != 0) { 283 do { 284 const float vx = *x++; 285 286 // General structure of the algorithm: 287 // / exp(x) / (1 + exp(x)) if x <= 0 288 // f[x] := 289 // \ 1 - f[-x] if x >= 0 290 // 291 // First we compute f[-z] := exp(-z) / (1 + exp(-z)) where z = abs(x), 292 // then replace result with 1 - f[-z] if x >= 0. 293 const float vz = fabsf(vx); 294 295 // Compute reduced argument n := round(-z * 2048 / log(2)). 296 // We do it by adding a large number (magic bias), which cause rounding of the result to an integer, then subtracing 297 // the large number back. The first addition is combined with multiplication by log2e into a single FMA instruction. 298 // The trick with adding large number is valid only within certain bounds (|z * 2048 / log(2)| <= 2**22, i.e. 299 // |z| <= 0x1.62E43p+10 = 1419.5654296875), but that is acceptable, because inputs x outside of 300 // [-87.336544, 17.328678] (i.e. z outsize [0, 87.336544]) underflow or saturate sigmoidf(x). We fixup the result 301 // for such inputs at the very end of the algorithm. 302 float vn = vz * vminus_log2e_x2048 + vmagic_bias; 303 304 // Create a floating-point number s (scale) such that s := 2**(n / 2048) for such inputs that sigmoidf(-z) is 305 // normalized, i.e. 0 <= z <= 87.33642. As n has 11 fractional bits, we split s == 2**(n / 2048) = 306 // = 2**e * 2**(n / 2048 - e), where e := int(n / 2048). We create s in two steps: 307 // 1. Fetch 2**(n / 2048 - e) = 2**(n % 2048) from table using the 6 low bits of n, as integer. 308 // Note that the fetched values are in the [1.0, 2.0) range, i.e. their floating-point exponent is 0. 309 // 2. Adjust fecthed value by addition of e to its floating-point exponent. The result is always a normalized 310 // number, because for 0 <= z <= 87.33642 (inputs for which sigmoidf(-z) is normalized) we have -126 <= e <= 0, 311 // and thus the adjusted exponent is not lower than -126. 312 // 313 // Extract e from bits 11:19 of n and shift it into bits 23:31 (position of floating-point exponent). 314 const uint32_t ve = (fp32_to_bits(vn) & ~vindex_mask) << 12; 315 316 // Use bits 0:11 bits of n, as integer, as an index for table lookup of l := 2**(n % 2048). 317 const uint32_t vidx = fp32_to_bits(vn) & vindex_mask; 318 // Adjust exponent of the value l fetched from the table to get the final s value. 319 const float vs = fp32_from_bits(xnn_table_exp2_k_over_2048[vidx] + ve); 320 321 // Subtract the large number back to get the final n := round(-z * 2048 / log(2)) as a floating-point number. 322 vn -= vmagic_bias; 323 324 // Compute reduced argument t := (z + n * log(2) / 2048). Note that -t = -z - n * log(2) / 2048. 325 // Use Cody-Waite range reduction method (note two constants to represent log(2) / 2048) to improve accuracy. 326 float vt = vn * vln2_o2048_hi + vz; 327 vt = vn * vln2_o2048_lo + vt; 328 329 // Compute degree-1 polynomial approximation for exp(-t) on [-log(2)/4096, log(2)/4096]: 330 // P1(t) = 1 + t * c1 331 const float vp = vt * vc1; 332 333 // Reconstruct the exp(-z) value: 334 // y = s * (1 + t * c1) 335 // = s + s * (t * c1)) 336 // = s + s * p 337 const float vy = vp * vs + vs; 338 339 // Reconstruct sigmoid(-z) = exp(-z) / (1.0 + exp(-z)) 340 float vf = vy / (vy + vone); 341 342 // For inputs above denormal cutoff, replace output with +0.0f. 343 // Note that for NaN inputs, comparison result is false, and outputs are left unchanged. 344 if XNN_UNPREDICTABLE(vz > vdenorm_cutoff) { 345 vf = 0.0f; 346 } 347 348 // Reconstruct sigmoid(x) = x < 0 ? sigmoid(-z) : 1.0 - sigmoid(-z) 349 if XNN_UNPREDICTABLE(vx > 0.0f) { 350 vf = vone - vf; 351 } 352 353 *y++ = vf; 354 355 n -= sizeof(float); 356 } while (n != 0); 357 } 358} 359