Cleanup lmms_math.h (#7382)

* simplified fraction and absfraction functions * removed unused fastSqrt() and fastPow() functions * unused absMin() and absMax() * move roundAt to math header * Code review from saker Co-authored-by: saker <sakertooth@gmail.com> * use std::trunc() * fixup after fixing merge conflicts * remove unused fastFma and fastFmal functions. * remove lmms_basics include, not needed * use signedPowf from lmms_math in NES * removed fastRand function, unused * remove unused sinc function * cleanup signedPowf * code review * further simplify random number math * removed static from lmms_math file --------- Co-authored-by: saker <sakertooth@gmail.com>
2024-08-28 12:48:56 +05:30
parent ff8c47062f
commit a992019626
9 changed files with 75 additions and 233 deletions
--- a/include/interpolation.h
+++ b/include/interpolation.h
@@ -71,11 +71,11 @@ inline float cubicInterpolate( float v0, float v1, float v2, float v3, float x )
 {
 	float frsq = x*x;
 	float frcu = frsq*v0;
-	float t1 = v3 + 3*v1;
+	float t1 = std::fma(v1, 3, v3);

-	return( v1 + fastFmaf( 0.5f, frcu, x ) * ( v2 - frcu * ( 1.0f/6.0f ) -
-		fastFmaf( t1, ( 1.0f/6.0f ), -v0 ) * ( 1.0f/3.0f ) ) + frsq * x * ( t1 *
-		( 1.0f/6.0f ) - 0.5f * v2 ) + frsq * fastFmaf( 0.5f, v2, -v1 ) );
+	return (v1 + std::fma(0.5f, frcu, x) * (v2 - frcu * (1.0f / 6.0f) -
+		std::fma(t1, (1.0f / 6.0f), -v0) * (1.0f / 3.0f)) + frsq * x * (t1 *
+		(1.0f / 6.0f) - 0.5f * v2) + frsq * std::fma(0.5f, v2, -v1));
 }


@@ -83,13 +83,13 @@ inline float cubicInterpolate( float v0, float v1, float v2, float v3, float x )
 inline float cosinusInterpolate( float v0, float v1, float x )
 {
 	const float f = ( 1.0f - cosf( x * F_PI ) ) * 0.5f;
-	return fastFmaf( f, v1-v0, v0 );
+	return std::fma(f, v1 - v0, v0);
 }


 inline float linearInterpolate( float v0, float v1, float x )
 {
-	return fastFmaf( x, v1-v0, v0 );
+	return std::fma(x, v1 - v0, v0);
 }


@@ -104,7 +104,7 @@ inline float optimalInterpolate( float v0, float v1, float x )
 	const float c2 = even * -0.004541102062639801;
 	const float c3 = odd * -1.57015627178718420;
 	
-	return fastFmaf( fastFmaf( fastFmaf( c3, z, c2 ), z, c1 ), z, c0 );
+	return std::fma(std::fma(std::fma(c3, z, c2), z, c1), z, c0);
 }


@@ -121,7 +121,7 @@ inline float optimal4pInterpolate( float v0, float v1, float v2, float v3, float
 	const float c2 = even1 * -0.246185007019907091 + even2 * 0.24614027139700284;
 	const float c3 = odd1 * -0.36030925263849456 + odd2 * 0.10174985775982505;

-	return fastFmaf( fastFmaf( fastFmaf( c3, z, c2 ), z, c1 ), z, c0 );
+	return std::fma(std::fma(std::fma(c3, z, c2), z, c1), z, c0);
 }


@@ -132,7 +132,7 @@ inline float lagrangeInterpolate( float v0, float v1, float v2, float v3, float
 	const float c1 = v2 - v0 * ( 1.0f / 3.0f ) - v1 * 0.5f - v3 * ( 1.0f / 6.0f );
 	const float c2 = 0.5f * (v0 + v2) - v1;
 	const float c3 = ( 1.0f/6.0f ) * ( v3 - v0 ) + 0.5f * ( v1 - v2 );
-	return fastFmaf( fastFmaf( fastFmaf( c3, x, c2 ), x, c1 ), x, c0 );
+	return std::fma(std::fma(std::fma(c3, x, c2), x, c1), x, c0);
 }


--- a/include/lmms_math.h
+++ b/include/lmms_math.h
@@ -37,40 +37,11 @@
 namespace lmms
 {

-static inline bool approximatelyEqual(float x, float y)
+inline bool approximatelyEqual(float x, float y)
 {
 	return x == y ? true : std::abs(x - y) < F_EPSILON;
 }

-#ifdef __INTEL_COMPILER
-
-static inline float absFraction( const float _x )
-{
-	return( _x - floorf( _x ) );
-}
-
-static inline float fraction( const float _x )
-{
-	return( _x - floorf( _x ) - ( _x >= 0.0f ? 0.0 : 1.0 ) );
-}
-
-#else
-
-/*!
- * @brief Returns the wrapped fractional part of a float, a value between 0.0f and 1.0f.
- *
- * absFraction( 2.3) =>  0.3
- * absFraction(-2.3) =>  0.7
- *
- * Note that this not the same as the absolute value of the fraction (as the function name suggests).
- * If the result is interpreted as a phase of an oscillator, it makes that negative phases are
- * converted to positive phases.
- */
-static inline float absFraction(const float x)
-{
-	return x - std::floor(x);
-}
-
 /*!
 * @brief Returns the fractional part of a float, a value between -1.0f and 1.0f.
 *
@@ -80,143 +51,75 @@ static inline float absFraction(const float x)
 * Note that if the return value is used as a phase of an oscillator, that the oscillator must support
 * negative phases.
 */
-static inline float fraction( const float _x )
+inline float fraction(const float x)
 {
-	return( _x - static_cast<int>( _x ) );
+	return x - std::trunc(x);
+}
+
+/*!
+ * @brief Returns the wrapped fractional part of a float, a value between 0.0f and 1.0f.
+ *
+ * absFraction( 2.3) =>  0.3
+ * absFraction(-2.3) =>  0.7
+ *
+ * Note that this not the same as the absolute value of the fraction (as the function name suggests).
+ * If the result is interpreted as a phase of an oscillator, it makes that negative phases are
+ * converted to positive phases.
+ */
+inline float absFraction(const float x)
+{
+	return x - std::floor(x);
 }


-#if 0
-// SSE3-version
-static inline float absFraction( float _x )
-{
-	unsigned int tmp;
-	asm(
-		"fld %%st\n\t"
-		"fisttp %1\n\t"
-		"fild %1\n\t"
-		"ftst\n\t"
-		"sahf\n\t"
-		"jae 1f\n\t"
-		"fld1\n\t"
-		"fsubrp %%st, %%st(1)\n\t"
-	"1:\n\t"
-		"fsubrp %%st, %%st(1)"
-		: "+t"( _x ), "=m"( tmp )
-		:
-		: "st(1)", "cc" );
-	return( _x );
-}
-
-static inline float absFraction( float _x )
-{
-	unsigned int tmp;
-	asm(
-		"fld %%st\n\t"
-		"fisttp %1\n\t"
-		"fild %1\n\t"
-		"fsubrp %%st, %%st(1)"
-		: "+t"( _x ), "=m"( tmp )
-		:
-		: "st(1)" );
-	return( _x );
-}
-#endif
-
-#endif // __INTEL_COMPILER
-
-
-
-constexpr int FAST_RAND_MAX = 32767;
-static inline int fast_rand()
+constexpr float FAST_RAND_RATIO = 1.0f / 32767;
+inline int fast_rand()
 {
 	static unsigned long next = 1;
 	next = next * 1103515245 + 12345;
 	return( (unsigned)( next / 65536 ) % 32768 );
 }

-static inline double fastRand( double range )
+inline float fastRandf(float range)
 {
-	static const double fast_rand_ratio = 1.0 / FAST_RAND_MAX;
-	return fast_rand() * range * fast_rand_ratio;
+	return fast_rand() * range * FAST_RAND_RATIO;
 }

-static inline float fastRandf( float range )
-{
-	static const float fast_rand_ratio = 1.0f / FAST_RAND_MAX;
-	return fast_rand() * range * fast_rand_ratio;
-}

-//! @brief Takes advantage of fmal() function if present in hardware
-static inline long double fastFmal( long double a, long double b, long double c ) 
+//! Round `value` to `where` depending on step size
+template<class T>
+static void roundAt(T& value, const T& where, const T& stepSize)
 {
-#ifdef FP_FAST_FMAL
-	#ifdef __clang__
-		return fma( a, b, c );
-	#else
-		return fmal( a, b, c );
-	#endif
-#else
-	return a * b + c;
-#endif // FP_FAST_FMAL
-}
-
-//! @brief Takes advantage of fmaf() function if present in hardware
-static inline float fastFmaf( float a, float b, float c ) 
-{
-#ifdef FP_FAST_FMAF
-	#ifdef __clang__
-		return fma( a, b, c );
-	#else
-		return fmaf( a, b, c );
-	#endif
-#else
-	return a * b + c;
-#endif // FP_FAST_FMAF
-}
-
-//! @brief Takes advantage of fma() function if present in hardware
-static inline double fastFma( double a, double b, double c ) 
-{
-#ifdef FP_FAST_FMA
-	return fma( a, b, c );
-#else
-	return a * b + c;
-#endif
-}
-
-// source: http://martin.ankerl.com/2007/10/04/optimized-pow-approximation-for-java-and-c-c/
-static inline double fastPow( double a, double b )
-{
-	union
+	if (std::abs(value - where) < F_EPSILON * std::abs(stepSize))
 	{
-		double d;
-		int32_t x[2];
-	} u = { a };
-	u.x[1] = static_cast<int32_t>( b * ( u.x[1] - 1072632447 ) + 1072632447 );
-	u.x[0] = 0;
-	return u.d;
+		value = where;
+	}
 }

-// sinc function
-static inline double sinc( double _x )
+
+//! returns 1.0f if val >= 0.0f, -1.0 else
+inline float sign(float val) 
+{ 
+	return val >= 0.0f ? 1.0f : -1.0f; 
+}
+
+
+//! if val >= 0.0f, returns sqrtf(val), else: -sqrtf(-val)
+inline float sqrt_neg(float val) 
 {
-	return _x == 0.0 ? 1.0 : sin( F_PI * _x ) / ( F_PI * _x );
+	return std::sqrt(std::abs(val)) * sign(val);
 }

-
 //! @brief Exponential function that deals with negative bases
-static inline float signedPowf( float v, float e )
+inline float signedPowf(float v, float e)
 {
-	return v < 0 
-		? powf( -v, e ) * -1.0f
-		: powf( v, e );
+	return std::pow(std::abs(v), e) * sign(v);
 }


 //! @brief Scales @value from linear to logarithmic.
 //! Value should be within [0,1]
-static inline float logToLinearScale( float min, float max, float value )
+inline float logToLinearScale(float min, float max, float value)
 {
 	if( min < 0 )
 	{
@@ -231,7 +134,7 @@ static inline float logToLinearScale( float min, float max, float value )


 //! @brief Scales value from logarithmic to linear. Value should be in min-max range.
-static inline float linearToLogScale( float min, float max, float value )
+inline float linearToLogScale(float min, float max, float value)
 {
 	static const float EXP = 1.0f / F_E;
 	const float valueLimited = std::clamp(value, min, max);
@@ -252,7 +155,7 @@ static inline float linearToLogScale( float min, float max, float value )
 //! @brief Converts linear amplitude (0-1.0) to dBFS scale. Handles zeroes as -inf.
 //! @param amp Linear amplitude, where 1.0 = 0dBFS. 
 //! @return Amplitude in dBFS. -inf for 0 amplitude.
-static inline float safeAmpToDbfs( float amp )
+inline float safeAmpToDbfs(float amp)
 {
 	return amp == 0.0f
 		? -INFINITY
@@ -263,7 +166,7 @@ static inline float safeAmpToDbfs( float amp )
 //! @brief Converts dBFS-scale to linear amplitude with 0dBFS = 1.0. Handles infinity as zero.
 //! @param dbfs The dBFS value to convert: all infinites are treated as -inf and result in 0
 //! @return Linear amplitude
-static inline float safeDbfsToAmp( float dbfs )
+inline float safeDbfsToAmp(float dbfs)
 {
 	return std::isinf( dbfs )
 		? 0.0f
@@ -274,7 +177,7 @@ static inline float safeDbfsToAmp( float dbfs )
 //! @brief Converts linear amplitude (>0-1.0) to dBFS scale. 
 //! @param amp Linear amplitude, where 1.0 = 0dBFS. ** Must be larger than zero! **
 //! @return Amplitude in dBFS. 
-static inline float ampToDbfs(float amp)
+inline float ampToDbfs(float amp)
 {
 	return log10f(amp) * 20.0f;
 }
@@ -283,54 +186,12 @@ static inline float ampToDbfs(float amp)
 //! @brief Converts dBFS-scale to linear amplitude with 0dBFS = 1.0
 //! @param dbfs The dBFS value to convert. ** Must be a real number - not inf/nan! **
 //! @return Linear amplitude
-static inline float dbfsToAmp(float dbfs)
+inline float dbfsToAmp(float dbfs)
 {
 	return std::pow(10.f, dbfs * 0.05f);
 }


-
-//! returns 1.0f if val >= 0.0f, -1.0 else
-static inline float sign( float val ) 
-{ 
-	return val >= 0.0f ? 1.0f : -1.0f; 
-}
-
-
-//! if val >= 0.0f, returns sqrtf(val), else: -sqrtf(-val)
-static inline float sqrt_neg( float val ) 
-{
-	return sqrtf( fabs( val ) ) * sign( val );
-}
-
-
-// fast approximation of square root
-static inline float fastSqrt( float n )
-{
-	union 
-	{
-		int32_t i;
-		float f;
-	} u;
-	u.f = n;
-	u.i = ( u.i + ( 127 << 23 ) ) >> 1;
-	return u.f;
-}
-
-//! returns value furthest from zero
-template<class T>
-static inline T absMax( T a, T b )
-{
-	return std::abs(a) > std::abs(b) ? a : b;
-}
-
-//! returns value nearest to zero
-template<class T>
-static inline T absMin( T a, T b )
-{
-	return std::abs(a) < std::abs(b) ? a : b;
-}
-
 //! Returns the linear interpolation of the two values
 template<class T, class F>
 constexpr T lerp(T a, T b, F t)
@@ -340,13 +201,12 @@ constexpr T lerp(T a, T b, F t)

 // @brief Calculate number of digits which LcdSpinBox would show for a given number
 // @note Once we upgrade to C++20, we could probably use std::formatted_size
-static inline int numDigitsAsInt(float f)
+inline int numDigitsAsInt(float f)
 {
 	// use rounding:
-	// LcdSpinBox sometimes uses roundf(), sometimes cast rounding
+	// LcdSpinBox sometimes uses std::round(), sometimes cast rounding
 	// we use rounding to be on the "safe side"
-	const float rounded = roundf(f);
-	int asInt = static_cast<int>(rounded);
+	int asInt = static_cast<int>(std::round(f));
 	int digits = 1; // always at least 1
 	if(asInt < 0)
 	{
@@ -354,11 +214,11 @@ static inline int numDigitsAsInt(float f)
 		asInt = -asInt;
 	}
 	// "asInt" is positive from now
-	int32_t power = 1;
-	for(int32_t i = 1; i<10; ++i)
+	int power = 1;
+	for (int i = 1; i < 10; ++i)
 	{
 		power *= 10;
-		if(static_cast<int32_t>(asInt) >= power) { ++digits; } // 2 digits for >=10, 3 for >=100
+		if (asInt >= power) { ++digits; } // 2 digits for >=10, 3 for >=100
 		else { break; }
 	}
 	return digits;
--- a/plugins/Flanger/Noise.cpp
+++ b/plugins/Flanger/Noise.cpp
@@ -31,7 +31,7 @@ namespace lmms

 Noise::Noise()
 {
-	inv_randmax = 1.0/FAST_RAND_MAX; /* for range of 0 - 1.0 */
+	inv_randmax = FAST_RAND_RATIO; /* for range of 0 - 1.0 */
 }


@@ -39,7 +39,7 @@ Noise::Noise()

 float Noise::tick()
 {
-	return (float) ((2.0 * fast_rand() * inv_randmax) - 1.0);
+	return fastRandf(2.0f) - 1.0f;
 }


--- a/plugins/GranularPitchShifter/GranularPitchShifterEffect.cpp
+++ b/plugins/GranularPitchShifter/GranularPitchShifterEffect.cpp
@@ -157,17 +157,17 @@ bool GranularPitchShifterEffect::processAudioBuffer(SampleFrame* buf, const fpp_
 		if (++m_timeSinceLastGrain >= m_nextWaitRandomization * waitMult)
 		{
 			m_timeSinceLastGrain = 0;
-			double randThing = (fast_rand()/static_cast<double>(FAST_RAND_MAX) * 2. - 1.);
+			double randThing = fast_rand() * static_cast<double>(FAST_RAND_RATIO) * 2. - 1.;
 			m_nextWaitRandomization = std::exp2(randThing * twitch);
 			double grainSpeed = 1. / std::exp2(randThing * jitter);

 			std::array<float, 2> sprayResult = {0, 0};
 			if (spray > 0)
 			{
-				sprayResult[0] = (fast_rand() / static_cast<float>(FAST_RAND_MAX)) * spray * m_sampleRate;
+				sprayResult[0] = fast_rand() * FAST_RAND_RATIO * spray * m_sampleRate;
 				sprayResult[1] = linearInterpolate(
 					sprayResult[0],
-					(fast_rand() / static_cast<float>(FAST_RAND_MAX)) * spray * m_sampleRate,
+					fast_rand() * FAST_RAND_RATIO * spray * m_sampleRate,
 					spraySpread);
 			}
 			
--- a/plugins/Kicker/KickerOsc.h
+++ b/plugins/Kicker/KickerOsc.h
@@ -64,7 +64,7 @@ public:
 	{
 		for( fpp_t frame = 0; frame < frames; ++frame )
 		{
-			const double gain = ( 1 - fastPow( ( m_counter < m_length ) ? m_counter / m_length : 1, m_env ) );
+			const double gain = 1 - std::pow((m_counter < m_length) ? m_counter / m_length : 1, m_env);
 			const sample_t s = ( Oscillator::sinSample( m_phase ) * ( 1 - m_noise ) ) + ( Oscillator::noiseSample( 0 ) * gain * gain * m_noise );
 			buf[frame][0] = s * gain;
 			buf[frame][1] = s * gain;
@@ -80,7 +80,7 @@ public:
 			m_FX.nextSample( buf[frame][0], buf[frame][1] );
 			m_phase += m_freq / sampleRate;

-			const double change = ( m_counter < m_length ) ? ( ( m_startFreq - m_endFreq ) * ( 1 - fastPow( m_counter / m_length, m_slope ) ) ) : 0;
+			const double change = (m_counter < m_length) ? ((m_startFreq - m_endFreq) * (1 - std::pow(m_counter / m_length, m_slope))) : 0;
 			m_freq = m_endFreq + change;
 			++m_counter;
 		}
--- a/plugins/Monstro/Monstro.cpp
+++ b/plugins/Monstro/Monstro.cpp
@@ -858,7 +858,7 @@ inline sample_t MonstroSynth::calcSlope( int slope, sample_t s )
 {
 	if( m_parent->m_slope[slope] == 1.0f ) return s;
 	if( s == 0.0f ) return s;
-	return fastPow( s, m_parent->m_slope[slope] );
+	return std::pow(s, m_parent->m_slope[slope]);
 }


--- a/plugins/Nes/Nes.cpp
+++ b/plugins/Nes/Nes.cpp
@@ -34,6 +34,7 @@
 #include "Oscillator.h"

 #include "embed.h"
+#include "lmms_math.h"
 #include "plugin_export.h"

 namespace lmms
@@ -392,7 +393,7 @@ void NesObject::renderOutput( SampleFrame* buf, fpp_t frames )
 		pin1 *= 1.0 + ( Oscillator::noiseSample( 0.0f ) * DITHER_AMP );		
 		pin1 = pin1 / 30.0f;
 		
-		pin1 = signedPow( pin1, NES_DIST );
+		pin1 = signedPowf(pin1, NES_DIST);
 		
 		pin1 = pin1 * 2.0f - 1.0f;
 		
@@ -401,7 +402,7 @@ void NesObject::renderOutput( SampleFrame* buf, fpp_t frames )
 		m_12Last = pin1;

 		// compensate DC offset
-		pin1 += 1.0f - signedPow( static_cast<float>( ch1Level + ch2Level ) / 30.0f, NES_DIST );
+		pin1 += 1.0f - signedPowf(static_cast<float>(ch1Level + ch2Level) / 30.0f, NES_DIST);
 		
 		pin1 *= NES_MIXING_12;

@@ -410,7 +411,7 @@ void NesObject::renderOutput( SampleFrame* buf, fpp_t frames )
 		pin2 *= 1.0 + ( Oscillator::noiseSample( 0.0f ) * DITHER_AMP );		
 		pin2 = pin2 / 30.0f;
 		
-		pin2 = signedPow( pin2, NES_DIST );
+		pin2 = signedPowf(pin2, NES_DIST);
 		
 		pin2 = pin2 * 2.0f - 1.0f;

@@ -419,7 +420,7 @@ void NesObject::renderOutput( SampleFrame* buf, fpp_t frames )
 		m_34Last = pin2;
 		
 		// compensate DC offset
-		pin2 += 1.0f - signedPow( static_cast<float>( ch3Level + ch4Level ) / 30.0f, NES_DIST );
+		pin2 += 1.0f - signedPowf(static_cast<float>(ch3Level + ch4Level) / 30.0f, NES_DIST);
 		
 		pin2 *= NES_MIXING_34;
 		
--- a/plugins/Nes/Nes.h
+++ b/plugins/Nes/Nes.h
@@ -136,13 +136,6 @@ public:
 		return static_cast<int>( m_samplerate / freq );
 	}
 	
-	inline float signedPow( float f, float e )
-	{
-		return f < 0 
-			? powf( qAbs( f ), e ) * -1.0f
-			: powf( f, e );
-	}
-	
 	inline int nearestNoiseFreq( float f )
 	{
 		int n = 15;
--- a/src/core/AutomatableModel.cpp
+++ b/src/core/AutomatableModel.cpp
@@ -350,18 +350,6 @@ float AutomatableModel::inverseScaledValue( float value ) const



-//! @todo: this should be moved into a maths header
-template<class T>
-void roundAt( T& value, const T& where, const T& step_size )
-{
-	if (std::abs(value - where) < F_EPSILON * std::abs(step_size))
-	{
-		value = where;
-	}
-}
-
-
-

 template<class T>
 void AutomatableModel::roundAt( T& value, const T& where ) const