NEMath.h
1 /* 2 * Copyright (c) 2016, 2019 ARM Limited. 3 * 4 * SPDX-License-Identifier: MIT 5 * 6 * Permission is hereby granted, free of charge, to any person obtaining a copy 7 * of this software and associated documentation files (the "Software"), to 8 * deal in the Software without restriction, including without limitation the 9 * rights to use, copy, modify, merge, publish, distribute, sublicense, and/or 10 * sell copies of the Software, and to permit persons to whom the Software is 11 * furnished to do so, subject to the following conditions: 12 * 13 * The above copyright notice and this permission notice shall be included in all 14 * copies or substantial portions of the Software. 15 * 16 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR 17 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 18 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE 19 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 20 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, 21 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE 22 * SOFTWARE. 23 */ 24 #ifndef __ARM_COMPUTE_NEMATH_H__ 25 #define __ARM_COMPUTE_NEMATH_H__ 26 27 28 #if defined(ARM_MATH_NEON) 29 /** Calculate floor of a vector. 30 * 31 * @param[in] val Input vector value in F32 format. 32 * 33 * @return The calculated floor vector. 34 */ 35 static inline float32x4_t vfloorq_f32(float32x4_t val); 36 37 /** Calculate inverse square root. 38 * 39 * @param[in] x Input value. 40 * 41 * @return The calculated inverse square root. 42 */ 43 static inline float32x2_t vinvsqrt_f32(float32x2_t x); 44 45 /** Calculate inverse square root. 46 * 47 * @param[in] x Input value. 48 * 49 * @return The calculated inverse square root. 50 */ 51 static inline float32x4_t vinvsqrtq_f32(float32x4_t x); 52 53 /** Calculate reciprocal. 54 * 55 * @param[in] x Input value. 56 * 57 * @return The calculated reciprocal. 58 */ 59 static inline float32x2_t vinv_f32(float32x2_t x); 60 61 /** Calculate reciprocal. 62 * 63 * @param[in] x Input value. 64 * 65 * @return The calculated reciprocal. 66 */ 67 static inline float32x4_t vinvq_f32(float32x4_t x); 68 69 /** Perform a 7th degree polynomial approximation using Estrin's method. 70 * 71 * @param[in] x Input vector value in F32 format. 72 * @param[in] coeffs Polynomial coefficients table. (array of flattened float32x4_t vectors) 73 * 74 * @return The calculated approximation. 75 */ 76 static inline float32x4_t vtaylor_polyq_f32(float32x4_t x, const float32_t *coeffs); 77 78 /** Calculate exponential 79 * 80 * @param[in] x Input vector value in F32 format. 81 * 82 * @return The calculated exponent. 83 */ 84 static inline float32x4_t vexpq_f32(float32x4_t x); 85 86 /** Calculate logarithm 87 * 88 * @param[in] x Input vector value in F32 format. 89 * 90 * @return The calculated logarithm. 91 */ 92 static inline float32x4_t vlogq_f32(float32x4_t x); 93 94 /** Calculate hyperbolic tangent. 95 * 96 * tanh(x) = (e^2x - 1)/(e^2x + 1) 97 * 98 * @note We clamp x to [-5,5] to avoid overflowing issues. 99 * 100 * @param[in] val Input vector value in F32 format. 101 * 102 * @return The calculated Hyperbolic Tangent. 103 */ 104 static inline float32x4_t vtanhq_f32(float32x4_t val); 105 106 /** Calculate n power of a number. 107 * 108 * pow(x,n) = e^(n*log(x)) 109 * 110 * @param[in] val Input vector value in F32 format. 111 * @param[in] n Powers to raise the input to. 112 * 113 * @return The calculated power. 114 */ 115 static inline float32x4_t vpowq_f32(float32x4_t val, float32x4_t n); 116 117 #ifdef __ARM_FEATURE_FP16_VECTOR_ARITHMETIC 118 /** Calculate hyperbolic tangent. 119 * 120 * tanh(x) = (e^2x - 1)/(e^2x + 1) 121 * 122 * @note We clamp x to [-5,5] to avoid overflowing issues. 123 * 124 * @param[in] val Input vector value in F32 format. 125 * 126 * @return The calculated Hyperbolic Tangent. 127 */ 128 static inline float16x8_t vtanhq_f16(float16x8_t val); 129 130 /** Calculate reciprocal. 131 * 132 * @param[in] x Input value. 133 * 134 * @return The calculated reciprocal. 135 */ 136 static inline float16x4_t vinv_f16(float16x4_t x); 137 138 /** Calculate reciprocal. 139 * 140 * @param[in] x Input value. 141 * 142 * @return The calculated reciprocal. 143 */ 144 static inline float16x8_t vinvq_f16(float16x8_t x); 145 146 /** Calculate inverse square root. 147 * 148 * @param[in] x Input value. 149 * 150 * @return The calculated inverse square root. 151 */ 152 static inline float16x4_t vinvsqrt_f16(float16x4_t x); 153 154 /** Calculate inverse square root. 155 * 156 * @param[in] x Input value. 157 * 158 * @return The calculated inverse square root. 159 */ 160 static inline float16x8_t vinvsqrtq_f16(float16x8_t x); 161 162 /** Calculate exponential 163 * 164 * @param[in] x Input vector value in F16 format. 165 * 166 * @return The calculated exponent. 167 */ 168 static inline float16x8_t vexpq_f16(float16x8_t x); 169 170 /** Calculate n power of a number. 171 * 172 * pow(x,n) = e^(n*log(x)) 173 * 174 * @param[in] val Input vector value in F16 format. 175 * @param[in] n Powers to raise the input to. 176 * 177 * @return The calculated power. 178 */ 179 static inline float16x8_t vpowq_f16(float16x8_t val, float16x8_t n); 180 #endif /* __ARM_FEATURE_FP16_VECTOR_ARITHMETIC */ 181 182 /** Exponent polynomial coefficients */ 183 extern const float32_t exp_tab[4*8]; 184 185 186 /** Logarithm polynomial coefficients */ 187 extern const float32_t log_tab[4*8]; 188 189 #ifndef DOXYGEN_SKIP_THIS 190 inline float32x4_t vfloorq_f32(float32x4_t val) 191 { 192 static const float32_t CONST_1[4] = {1.f,1.f,1.f,1.f}; 193 194 const int32x4_t z = vcvtq_s32_f32(val); 195 const float32x4_t r = vcvtq_f32_s32(z); 196 197 return vbslq_f32(vcgtq_f32(r, val), vsubq_f32(r, vld1q_f32(CONST_1)), r); 198 } 199 200 inline float32x2_t vinvsqrt_f32(float32x2_t x) 201 { 202 float32x2_t sqrt_reciprocal = vrsqrte_f32(x); 203 sqrt_reciprocal = vmul_f32(vrsqrts_f32(vmul_f32(x, sqrt_reciprocal), sqrt_reciprocal), sqrt_reciprocal); 204 sqrt_reciprocal = vmul_f32(vrsqrts_f32(vmul_f32(x, sqrt_reciprocal), sqrt_reciprocal), sqrt_reciprocal); 205 206 return sqrt_reciprocal; 207 } 208 209 inline float32x4_t vinvsqrtq_f32(float32x4_t x) 210 { 211 float32x4_t sqrt_reciprocal = vrsqrteq_f32(x); 212 sqrt_reciprocal = vmulq_f32(vrsqrtsq_f32(vmulq_f32(x, sqrt_reciprocal), sqrt_reciprocal), sqrt_reciprocal); 213 sqrt_reciprocal = vmulq_f32(vrsqrtsq_f32(vmulq_f32(x, sqrt_reciprocal), sqrt_reciprocal), sqrt_reciprocal); 214 215 return sqrt_reciprocal; 216 } 217 218 inline float32x2_t vinv_f32(float32x2_t x) 219 { 220 float32x2_t recip = vrecpe_f32(x); 221 recip = vmul_f32(vrecps_f32(x, recip), recip); 222 recip = vmul_f32(vrecps_f32(x, recip), recip); 223 return recip; 224 } 225 226 inline float32x4_t vinvq_f32(float32x4_t x) 227 { 228 float32x4_t recip = vrecpeq_f32(x); 229 recip = vmulq_f32(vrecpsq_f32(x, recip), recip); 230 recip = vmulq_f32(vrecpsq_f32(x, recip), recip); 231 return recip; 232 } 233 234 inline float32x4_t vtaylor_polyq_f32(float32x4_t x, const float32_t *coeffs) 235 { 236 float32x4_t A = vmlaq_f32(vld1q_f32(&coeffs[4*0]), vld1q_f32(&coeffs[4*4]), x); 237 float32x4_t B = vmlaq_f32(vld1q_f32(&coeffs[4*2]), vld1q_f32(&coeffs[4*6]), x); 238 float32x4_t C = vmlaq_f32(vld1q_f32(&coeffs[4*1]), vld1q_f32(&coeffs[4*5]), x); 239 float32x4_t D = vmlaq_f32(vld1q_f32(&coeffs[4*3]), vld1q_f32(&coeffs[4*7]), x); 240 float32x4_t x2 = vmulq_f32(x, x); 241 float32x4_t x4 = vmulq_f32(x2, x2); 242 float32x4_t res = vmlaq_f32(vmlaq_f32(A, B, x2), vmlaq_f32(C, D, x2), x4); 243 return res; 244 } 245 246 inline float32x4_t vexpq_f32(float32x4_t x) 247 { 248 static const float32_t CONST_LN2[4] = {0.6931471805f,0.6931471805f,0.6931471805f,0.6931471805f}; // ln(2) 249 static const float32_t CONST_INV_LN2[4] = {1.4426950408f,1.4426950408f,1.4426950408f,1.4426950408f}; // 1/ln(2) 250 static const float32_t CONST_0[4] = {0.f,0.f,0.f,0.f}; 251 static const int32_t CONST_NEGATIVE_126[4] = {-126,-126,-126,-126}; 252 253 // Perform range reduction [-log(2),log(2)] 254 int32x4_t m = vcvtq_s32_f32(vmulq_f32(x, vld1q_f32(CONST_INV_LN2))); 255 float32x4_t val = vmlsq_f32(x, vcvtq_f32_s32(m), vld1q_f32(CONST_LN2)); 256 257 // Polynomial Approximation 258 float32x4_t poly = vtaylor_polyq_f32(val, exp_tab); 259 260 // Reconstruct 261 poly = vreinterpretq_f32_s32(vqaddq_s32(vreinterpretq_s32_f32(poly), vqshlq_n_s32(m, 23))); 262 poly = vbslq_f32(vcltq_s32(m, vld1q_s32(CONST_NEGATIVE_126)), vld1q_f32(CONST_0), poly); 263 264 return poly; 265 } 266 267 inline float32x4_t vlogq_f32(float32x4_t x) 268 { 269 static const int32_t CONST_127[4] = {127,127,127,127}; // 127 270 static const float32_t CONST_LN2[4] = {0.6931471805f,0.6931471805f,0.6931471805f,0.6931471805f}; // ln(2) 271 272 // Extract exponent 273 int32x4_t m = vsubq_s32(vreinterpretq_s32_u32(vshrq_n_u32(vreinterpretq_u32_f32(x), 23)), vld1q_s32(CONST_127)); 274 float32x4_t val = vreinterpretq_f32_s32(vsubq_s32(vreinterpretq_s32_f32(x), vshlq_n_s32(m, 23))); 275 276 // Polynomial Approximation 277 float32x4_t poly = vtaylor_polyq_f32(val, log_tab); 278 279 // Reconstruct 280 poly = vmlaq_f32(poly, vcvtq_f32_s32(m), vld1q_f32(CONST_LN2)); 281 282 return poly; 283 } 284 285 inline float32x4_t vtanhq_f32(float32x4_t val) 286 { 287 static const float32_t CONST_1[4] = {1.f,1.f,1.f,1.f}; 288 static const float32_t CONST_2[4] = {2.f,2.f,2.f,2.f}; 289 static const float32_t CONST_MIN_TANH[4] = {-10.f,-10.f,-10.f,-10.f}; 290 static const float32_t CONST_MAX_TANH[4] = {10.f,10.f,10.f,10.f}; 291 292 float32x4_t x = vminq_f32(vmaxq_f32(val, vld1q_f32(CONST_MIN_TANH)), vld1q_f32(CONST_MAX_TANH)); 293 float32x4_t exp2x = vexpq_f32(vmulq_f32(vld1q_f32(CONST_2), x)); 294 float32x4_t num = vsubq_f32(exp2x, vld1q_f32(CONST_1)); 295 float32x4_t den = vaddq_f32(exp2x, vld1q_f32(CONST_1)); 296 float32x4_t tanh = vmulq_f32(num, vinvq_f32(den)); 297 return tanh; 298 } 299 300 inline float32x4_t vpowq_f32(float32x4_t val, float32x4_t n) 301 { 302 return vexpq_f32(vmulq_f32(n, vlogq_f32(val))); 303 } 304 #endif /* DOXYGEN_SKIP_THIS */ 305 306 #ifdef __ARM_FEATURE_FP16_VECTOR_ARITHMETIC 307 /** Exponent polynomial coefficients */ 308 /** Logarithm polynomial coefficients */ 309 #ifndef DOXYGEN_SKIP_THIS 310 inline float16x8_t vfloorq_f16(float16x8_t val) 311 { 312 static const float16_t CONST_1[8] = {1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f}; 313 314 const int16x8_t z = vcvtq_s16_f16(val); 315 const float16x8_t r = vcvtq_f16_s16(z); 316 317 return vbslq_f16(vcgtq_f16(r, val), vsubq_f16(r, vld1q_f16(CONST_1)), r); 318 } 319 inline float16x4_t vinvsqrt_f16(float16x4_t x) 320 { 321 float16x4_t sqrt_reciprocal = vrsqrte_f16(x); 322 sqrt_reciprocal = vmul_f16(vrsqrts_f16(vmul_f16(x, sqrt_reciprocal), sqrt_reciprocal), sqrt_reciprocal); 323 sqrt_reciprocal = vmul_f16(vrsqrts_f16(vmul_f16(x, sqrt_reciprocal), sqrt_reciprocal), sqrt_reciprocal); 324 return sqrt_reciprocal; 325 } 326 327 inline float16x8_t vinvsqrtq_f16(float16x8_t x) 328 { 329 float16x8_t sqrt_reciprocal = vrsqrteq_f16(x); 330 sqrt_reciprocal = vmulq_f16(vrsqrtsq_f16(vmulq_f16(x, sqrt_reciprocal), sqrt_reciprocal), sqrt_reciprocal); 331 sqrt_reciprocal = vmulq_f16(vrsqrtsq_f16(vmulq_f16(x, sqrt_reciprocal), sqrt_reciprocal), sqrt_reciprocal); 332 return sqrt_reciprocal; 333 } 334 335 inline float16x4_t vinv_f16(float16x4_t x) 336 { 337 float16x4_t recip = vrecpe_f16(x); 338 recip = vmul_f16(vrecps_f16(x, recip), recip); 339 recip = vmul_f16(vrecps_f16(x, recip), recip); 340 return recip; 341 } 342 343 inline float16x8_t vinvq_f16(float16x8_t x) 344 { 345 float16x8_t recip = vrecpeq_f16(x); 346 recip = vmulq_f16(vrecpsq_f16(x, recip), recip); 347 recip = vmulq_f16(vrecpsq_f16(x, recip), recip); 348 return recip; 349 } 350 351 inline float16x8_t vtanhq_f16(float16x8_t val) 352 { 353 const float16_t CONST_1[8] = {1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f}; 354 const float16_t CONST_2[8] = {2.f,2.f,2.f,2.f,2.f,2.f,2.f,2.f}; 355 const float16_t CONST_MIN_TANH[8] = {-10.f,-10.f,-10.f,-10.f,-10.f,-10.f,-10.f,-10.f}; 356 const float16_t CONST_MAX_TANH[8] = {10.f,10.f,10.f,10.f,10.f,10.f,10.f,10.f}; 357 358 const float16x8_t x = vminq_f16(vmaxq_f16(val, vld1q_f16(CONST_MIN_TANH)), vld1q_f16(CONST_MAX_TANH)); 359 const float16x8_t exp2x = vexpq_f16(vmulq_f16(vld1q_f16(CONST_2), x)); 360 const float16x8_t num = vsubq_f16(exp2x, vld1q_f16(CONST_1)); 361 const float16x8_t den = vaddq_f16(exp2x, vld1q_f16(CONST_1)); 362 const float16x8_t tanh = vmulq_f16(num, vinvq_f16(den)); 363 return tanh; 364 } 365 366 inline float16x8_t vtaylor_polyq_f16(float16x8_t x, const float16_t *coeffs) 367 { 368 const float16x8_t A = vaddq_f16(vld1q_f16(&coeffs[8*0]), vmulq_f16(vld1q_f16(&coeffs[8*4]), x)); 369 const float16x8_t B = vaddq_f16(vld1q_f16(&coeffs[8*2]), vmulq_f16(vld1q_f16(&coeffs[8*6]), x)); 370 const float16x8_t C = vaddq_f16(vld1q_f16(&coeffs[8*1]), vmulq_f16(vld1q_f16(&coeffs[8*5]), x)); 371 const float16x8_t D = vaddq_f16(vld1q_f16(&coeffs[8*3]), vmulq_f16(vld1q_f16(&coeffs[8*7]), x)); 372 const float16x8_t x2 = vmulq_f16(x, x); 373 const float16x8_t x4 = vmulq_f16(x2, x2); 374 const float16x8_t res = vaddq_f16(vaddq_f16(A, vmulq_f16(B, x2)), vmulq_f16(vaddq_f16(C, vmulq_f16(D, x2)), x4)); 375 return res; 376 } 377 378 inline float16x8_t vexpq_f16(float16x8_t x) 379 { 380 // TODO (COMPMID-1535) : Revisit FP16 approximations 381 const float32x4_t x_high = vcvt_f32_f16(vget_high_f16(x)); 382 const float32x4_t x_low = vcvt_f32_f16(vget_low_f16(x)); 383 384 const float16x8_t res = vcvt_high_f16_f32(vcvt_f16_f32(vexpq_f32(x_low)), vexpq_f32(x_high)); 385 return res; 386 } 387 388 inline float16x8_t vlogq_f16(float16x8_t x) 389 { 390 // TODO (COMPMID-1535) : Revisit FP16 approximations 391 const float32x4_t x_high = vcvt_f32_f16(vget_high_f16(x)); 392 const float32x4_t x_low = vcvt_f32_f16(vget_low_f16(x)); 393 394 const float16x8_t res = vcvt_high_f16_f32(vcvt_f16_f32(vlogq_f32(x_low)), vlogq_f32(x_high)); 395 return res; 396 } 397 398 inline float16x8_t vpowq_f16(float16x8_t val, float16x8_t n) 399 { 400 // TODO (giaiod01) - COMPMID-1535 401 float32x4_t n0_f32 = vcvt_f32_f16(vget_low_f16(n)); 402 float32x4_t n1_f32 = vcvt_f32_f16(vget_high_f16(n)); 403 float32x4_t val0_f32 = vcvt_f32_f16(vget_low_f16(val)); 404 float32x4_t val1_f32 = vcvt_f32_f16(vget_high_f16(val)); 405 406 float32x4_t res0_f32 = vexpq_f32(vmulq_f32(n0_f32, vlogq_f32(val0_f32))); 407 float32x4_t res1_f32 = vexpq_f32(vmulq_f32(n1_f32, vlogq_f32(val1_f32))); 408 409 return vcombine_f16(vcvt_f16_f32(res0_f32), vcvt_f16_f32(res1_f32)); 410 } 411 #endif /* DOXYGEN_SKIP_THIS */ 412 #endif /* __ARM_FEATURE_FP16_VECTOR_ARITHMETIC */ 413 #endif 414 #endif /* __ARM_COMPUTE_NEMATH_H__ */