src/fastmath.h - third_party/github.com/google/liblc3 - Git at Google

 /******************************************************************************
  *
  *  Copyright 2022 Google LLC
  *
  *  Licensed under the Apache License, Version 2.0 (the "License");
  *  you may not use this file except in compliance with the License.
  *  You may obtain a copy of the License at:
  *
  *  http://www.apache.org/licenses/LICENSE-2.0
  *
  *  Unless required by applicable law or agreed to in writing, software
  *  distributed under the License is distributed on an "AS IS" BASIS,
  *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  *  See the License for the specific language governing permissions and
  *  limitations under the License.
  *
  ******************************************************************************/

 #ifndef __LC3_FASTMATH_H
 #define __LC3_FASTMATH_H

 #include <stdint.h>
 #include <float.h>
 #include <math.h>


 /**
  * IEEE 754 Floating point representation
  */

 #define LC3_IEEE754_SIGN_SHL   (31)
 #define LC3_IEEE754_SIGN_MASK  (1 << LC3_IEEE754_SIGN_SHL)

 #define LC3_IEEE754_EXP_SHL    (23)
 #define LC3_IEEE754_EXP_MASK   (0xff << LC3_IEEE754_EXP_SHL)
 #define LC3_IEEE754_EXP_BIAS   (127)


 /**
  * Fast multiply floating-point number by integral power of 2
  * x               Operand, finite number
  * exp             Exponent such that 2^x is finite
  * return          2^exp
  */
 static inline float lc3_ldexpf(float _x, int exp) {
     union { float f; int32_t s; } x = { .f = _x };

     if (x.s & LC3_IEEE754_EXP_MASK)
         x.s += exp << LC3_IEEE754_EXP_SHL;

     return x.f;
 }

 /**
  * Fast convert floating-point number to fractional and integral components
  * x               Operand, finite number
  * exp             Return the exponent part
  * return          The normalized fraction in [0.5:1[
  */
 static inline float lc3_frexpf(float _x, int *exp) {
     union { float f; uint32_t u; } x = { .f = _x };

     int e = (x.u & LC3_IEEE754_EXP_MASK) >> LC3_IEEE754_EXP_SHL;
     *exp = e - (LC3_IEEE754_EXP_BIAS - 1);

     x.u = (x.u & ~LC3_IEEE754_EXP_MASK) |
         ((LC3_IEEE754_EXP_BIAS - 1) << LC3_IEEE754_EXP_SHL);

     return x.f;
 }

 /**
  * Fast 2^n approximation
  * x               Operand, range -100 to 100
  * return          2^x approximation (max relative error ~ 4e-7)
  */
 static inline float lc3_exp2f(float x)
 {
     /* --- 2^(i/8) for i from 0 to 7 --- */

     static const float e[] = {
         1.00000000e+00, 1.09050773e+00, 1.18920712e+00, 1.29683955e+00,
         1.41421356e+00, 1.54221083e+00, 1.68179283e+00, 1.83400809e+00 };

     /* --- Polynomial approx in range 0 to 1/8 --- */

     static const float p[] = {
         1.00448128e-02, 5.54563260e-02, 2.40228756e-01, 6.93147140e-01 };

     /* --- Split the operand ---
      *
      * Such as x = k/8 + y, with k an integer, and |y| < 0.5/8
      *
      * Note that `fast-math` compiler option leads to rounding error,
      * disable optimisation with `volatile`. */

     volatile union { float f; int32_t s; } v;

     v.f = x + 0x1.8p20f;
     int k = v.s;
     x -= v.f - 0x1.8p20f;

     /* --- Compute 2^x, with |x| < 1 ---
      * Perform polynomial approximation in range -0.5/8 to 0.5/8,
      * and muplity by precomputed value of 2^(i/8), i in [0:7] */

     union { float f; int32_t s; } y;

     y.f = (      p[0]) * x;
     y.f = (y.f + p[1]) * x;
     y.f = (y.f + p[2]) * x;
     y.f = (y.f + p[3]) * x;
     y.f = (y.f +  1.f) * e[k & 7];

     /* --- Add the exponent --- */

     y.s += (k >> 3) << LC3_IEEE754_EXP_SHL;

     return y.f;
 }

 /**
  * Fast log2(x) approximation
  * x               Operand, greater than 0
  * return          log2(x) approximation (max absolute error ~ 1e-4)
  */
 static inline float lc3_log2f(float x)
 {
     float y;
     int e;

     /* --- Polynomial approx in range 0.5 to 1 --- */

     static const float c[] = {
         -1.29479677, 5.11769018, -8.42295281, 8.10557963, -3.50567360 };

     x = lc3_frexpf(x, &e);

     y = (    c[0]) * x;
     y = (y + c[1]) * x;
     y = (y + c[2]) * x;
     y = (y + c[3]) * x;
     y = (y + c[4]);

     /* --- Add log2f(2^e) and return --- */

     return e + y;
 }

 /**
  * Fast log10(x) approximation
  * x               Operand, greater than 0
  * return          log10(x) approximation (max absolute error ~ 1e-4)
  */
 static inline float lc3_log10f(float x)
 {
     return log10f(2) * lc3_log2f(x);
 }

 /**
  * Fast `10 * log10(x)` (or dB) approximation in fixed Q16
  * x               Operand, in range 2^-63 to 2^63 (1e-19 to 1e19)
  * return          10 * log10(x) in fixed Q16 (-190 to 192 dB)
  *
  * - The 0 value is accepted and return the minimum value ~ -191dB
  * - This function assumed that float 32 bits is coded IEEE 754
  */
 static inline int32_t lc3_db_q16(float x)
 {
     /* --- Table in Q15 --- */

     static const uint16_t t[][2] = {

         /* [n][0] = 10 * log10(2) * log2(1 + n/32), with n = [0..15]     */
         /* [n][1] = [n+1][0] - [n][0] (while defining [16][0])           */

         {     0, 4379 }, {  4379, 4248 }, {  8627, 4125 }, { 12753, 4009 },
         { 16762, 3899 }, { 20661, 3795 }, { 24456, 3697 }, { 28153, 3603 },
         { 31755, 3514 }, { 35269, 3429 }, { 38699, 3349 }, { 42047, 3272 },
         { 45319, 3198 }, { 48517, 3128 }, { 51645, 3061 }, { 54705, 2996 },

         /* [n][0] = 10 * log10(2) * log2(1 + n/32) - 10 * log10(2) / 2,  */
         /*     with n = [16..31]                                         */
         /* [n][1] = [n+1][0] - [n][0] (while defining [32][0])           */

         {  8381, 2934 }, { 11315, 2875 }, { 14190, 2818 }, { 17008, 2763 },
         { 19772, 2711 }, { 22482, 2660 }, { 25142, 2611 }, { 27754, 2564 },
         { 30318, 2519 }, { 32837, 2475 }, { 35312, 2433 }, { 37744, 2392 },
         { 40136, 2352 }, { 42489, 2314 }, { 44803, 2277 }, { 47080, 2241 },

     };

     /* --- Approximation ---
      *
      *   10 * log10(x^2) = 10 * log10(2) * log2(x^2)
      *
      *   And log2(x^2) = 2 * log2( (1 + m) * 2^e )
      *                 = 2 * (e + log2(1 + m)) , with m in range [0..1]
      *
      * Split the float values in :
      *   e2  Double value of the exponent (2 * e + k)
      *   hi  High 5 bits of mantissa, for precalculated result `t[hi][0]`
      *   lo  Low 16 bits of mantissa, for linear interpolation `t[hi][1]`
      *
      * Two cases, from the range of the mantissa :
      *   0 to 0.5   `k = 0`, use 1st part of the table
      *   0.5 to 1   `k = 1`, use 2nd part of the table  */

     union { float f; uint32_t u; } x2 = { .f = x*x };

     int e2 = (int)(x2.u >> 22) - 2*127;
     int hi = (x2.u >> 18) & 0x1f;
     int lo = (x2.u >>  2) & 0xffff;

     return e2 * 49321 + t[hi][0] + ((t[hi][1] * lo) >> 16);
 }


 #endif /* __LC3_FASTMATH_H */
	/******************************************************************************
	*
	* Copyright 2022 Google LLC
	*
	* Licensed under the Apache License, Version 2.0 (the "License");
	* you may not use this file except in compliance with the License.
	* You may obtain a copy of the License at:
	*
	* http://www.apache.org/licenses/LICENSE-2.0
	*
	* Unless required by applicable law or agreed to in writing, software
	* distributed under the License is distributed on an "AS IS" BASIS,
	* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
	* See the License for the specific language governing permissions and
	* limitations under the License.
	*
	******************************************************************************/

	#ifndef __LC3_FASTMATH_H
	#define __LC3_FASTMATH_H

	#include <stdint.h>
	#include <float.h>
	#include <math.h>


	/**
	* IEEE 754 Floating point representation
	*/

	#define LC3_IEEE754_SIGN_SHL (31)
	#define LC3_IEEE754_SIGN_MASK (1 << LC3_IEEE754_SIGN_SHL)

	#define LC3_IEEE754_EXP_SHL (23)
	#define LC3_IEEE754_EXP_MASK (0xff << LC3_IEEE754_EXP_SHL)
	#define LC3_IEEE754_EXP_BIAS (127)


	/**
	* Fast multiply floating-point number by integral power of 2
	* x Operand, finite number
	* exp Exponent such that 2^x is finite
	* return 2^exp
	*/
	static inline float lc3_ldexpf(float _x, int exp) {
	union { float f; int32_t s; } x = { .f = _x };

	if (x.s & LC3_IEEE754_EXP_MASK)
	x.s += exp << LC3_IEEE754_EXP_SHL;

	return x.f;
	}

	/**
	* Fast convert floating-point number to fractional and integral components
	* x Operand, finite number
	* exp Return the exponent part
	* return The normalized fraction in [0.5:1[
	*/
	static inline float lc3_frexpf(float _x, int *exp) {
	union { float f; uint32_t u; } x = { .f = _x };

	int e = (x.u & LC3_IEEE754_EXP_MASK) >> LC3_IEEE754_EXP_SHL;
	*exp = e - (LC3_IEEE754_EXP_BIAS - 1);

	x.u = (x.u & ~LC3_IEEE754_EXP_MASK) \|
	((LC3_IEEE754_EXP_BIAS - 1) << LC3_IEEE754_EXP_SHL);

	return x.f;
	}

	/**
	* Fast 2^n approximation
	* x Operand, range -100 to 100
	* return 2^x approximation (max relative error ~ 4e-7)
	*/
	static inline float lc3_exp2f(float x)
	{
	/* --- 2^(i/8) for i from 0 to 7 --- */

	static const float e[] = {
	1.00000000e+00, 1.09050773e+00, 1.18920712e+00, 1.29683955e+00,
	1.41421356e+00, 1.54221083e+00, 1.68179283e+00, 1.83400809e+00 };

	/* --- Polynomial approx in range 0 to 1/8 --- */

	static const float p[] = {
	1.00448128e-02, 5.54563260e-02, 2.40228756e-01, 6.93147140e-01 };

	/* --- Split the operand ---
	*
	* Such as x = k/8 + y, with k an integer, and \|y\| < 0.5/8
	*
	* Note that `fast-math` compiler option leads to rounding error,
	* disable optimisation with `volatile`. */

	volatile union { float f; int32_t s; } v;

	v.f = x + 0x1.8p20f;
	int k = v.s;
	x -= v.f - 0x1.8p20f;

	/* --- Compute 2^x, with \|x\| < 1 ---
	* Perform polynomial approximation in range -0.5/8 to 0.5/8,
	* and muplity by precomputed value of 2^(i/8), i in [0:7] */

	union { float f; int32_t s; } y;

	y.f = ( p[0]) * x;
	y.f = (y.f + p[1]) * x;
	y.f = (y.f + p[2]) * x;
	y.f = (y.f + p[3]) * x;
	y.f = (y.f + 1.f) * e[k & 7];

	/* --- Add the exponent --- */

	y.s += (k >> 3) << LC3_IEEE754_EXP_SHL;

	return y.f;
	}

	/**
	* Fast log2(x) approximation
	* x Operand, greater than 0
	* return log2(x) approximation (max absolute error ~ 1e-4)
	*/
	static inline float lc3_log2f(float x)
	{
	float y;
	int e;

	/* --- Polynomial approx in range 0.5 to 1 --- */

	static const float c[] = {
	-1.29479677, 5.11769018, -8.42295281, 8.10557963, -3.50567360 };

	x = lc3_frexpf(x, &e);

	y = ( c[0]) * x;
	y = (y + c[1]) * x;
	y = (y + c[2]) * x;
	y = (y + c[3]) * x;
	y = (y + c[4]);

	/* --- Add log2f(2^e) and return --- */

	return e + y;
	}

	/**
	* Fast log10(x) approximation
	* x Operand, greater than 0
	* return log10(x) approximation (max absolute error ~ 1e-4)
	*/
	static inline float lc3_log10f(float x)
	{
	return log10f(2) * lc3_log2f(x);
	}

	/**
	* Fast `10 * log10(x)` (or dB) approximation in fixed Q16
	* x Operand, in range 2^-63 to 2^63 (1e-19 to 1e19)
	* return 10 * log10(x) in fixed Q16 (-190 to 192 dB)
	*
	* - The 0 value is accepted and return the minimum value ~ -191dB
	* - This function assumed that float 32 bits is coded IEEE 754
	*/
	static inline int32_t lc3_db_q16(float x)
	{
	/* --- Table in Q15 --- */

	static const uint16_t t[][2] = {

	/* [n][0] = 10 * log10(2) * log2(1 + n/32), with n = [0..15] */
	/* [n][1] = [n+1][0] - [n][0] (while defining [16][0]) */

	{ 0, 4379 }, { 4379, 4248 }, { 8627, 4125 }, { 12753, 4009 },
	{ 16762, 3899 }, { 20661, 3795 }, { 24456, 3697 }, { 28153, 3603 },
	{ 31755, 3514 }, { 35269, 3429 }, { 38699, 3349 }, { 42047, 3272 },
	{ 45319, 3198 }, { 48517, 3128 }, { 51645, 3061 }, { 54705, 2996 },

	/* [n][0] = 10 * log10(2) * log2(1 + n/32) - 10 * log10(2) / 2, */
	/* with n = [16..31] */
	/* [n][1] = [n+1][0] - [n][0] (while defining [32][0]) */

	{ 8381, 2934 }, { 11315, 2875 }, { 14190, 2818 }, { 17008, 2763 },
	{ 19772, 2711 }, { 22482, 2660 }, { 25142, 2611 }, { 27754, 2564 },
	{ 30318, 2519 }, { 32837, 2475 }, { 35312, 2433 }, { 37744, 2392 },
	{ 40136, 2352 }, { 42489, 2314 }, { 44803, 2277 }, { 47080, 2241 },

	};

	/* --- Approximation ---
	*
	* 10 * log10(x^2) = 10 * log10(2) * log2(x^2)
	*
	* And log2(x^2) = 2 * log2( (1 + m) * 2^e )
	* = 2 * (e + log2(1 + m)) , with m in range [0..1]
	*
	* Split the float values in :
	* e2 Double value of the exponent (2 * e + k)
	* hi High 5 bits of mantissa, for precalculated result `t[hi][0]`
	* lo Low 16 bits of mantissa, for linear interpolation `t[hi][1]`
	*
	* Two cases, from the range of the mantissa :
	* 0 to 0.5 `k = 0`, use 1st part of the table
	* 0.5 to 1 `k = 1`, use 2nd part of the table */

	union { float f; uint32_t u; } x2 = { .f = x*x };

	int e2 = (int)(x2.u >> 22) - 2*127;
	int hi = (x2.u >> 18) & 0x1f;
	int lo = (x2.u >> 2) & 0xffff;

	return e2 * 49321 + t[hi][0] + ((t[hi][1] * lo) >> 16);
	}


	#endif /* __LC3_FASTMATH_H */