src/common/RTMath.h

/***************************************************************************
 *                                                                         *
 *   LinuxSampler - modular, streaming capable sampler                     *
 *                                                                         *
 *   Copyright (C) 2003, 2004 by Benno Senoner and Christian Schoenebeck   *
 *   Copyright (C) 2005 - 2019 Christian Schoenebeck                       *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the Free Software           *
 *   Foundation, Inc., 59 Temple Place, Suite 330, Boston,                 *
 *   MA  02111-1307  USA                                                   *
 ***************************************************************************/

#ifndef __RT_MATH_H__
#define __RT_MATH_H__

#include <math.h>
#include <stdint.h>
#include "global_private.h"

/// Needed for calculating frequency ratio used to pitch a sample
#define TWELVEHUNDREDTH_ROOT_OF_TWO     1.000577789506555

enum implementation_t {
    CPP,
    ASM_X86_MMX_SSE
};

/** @brief Real Time Math Base Class
 *
 * Math functions for real time operation. This base class contains all
 * non-template methods.
 */
class RTMathBase {
    public:
        /**
         * High resolution time stamp.
         */
        typedef uint32_t time_stamp_t;

        typedef uint64_t usecs_t;

        /**
         * We read the processor's cycle count register as a reference
         * for the real time. These are of course only abstract values
         * with arbitrary time entity, but that's not a problem as long
         * as we calculate relatively.
         *
         * @see unsafeMicroSeconds()
         */
        static time_stamp_t CreateTimeStamp();

        /**
         * Calculates the frequency ratio for a pitch value given in cents
         * (assuming equal tempered scale of course, divided into 12
         * semitones per octave and 100 cents per semitone).
         *
         * Note: CONFIG_MAX_PITCH (defined in config.h) has to be defined to an
         * appropriate value, otherwise the behavior of this function is
         * undefined, but most probably if CONFIG_MAX_PITCH is too small, the
         * application will crash due to segmentation fault here.
         *
         * @param cents - pitch value in cents (+1200 cents means +1 octave)
         * @returns  frequency ratio (e.g. +2.0 for +1 octave)
         */
        inline static double CentsToFreqRatio(double Cents) {
            int   index_int   = (int) (Cents);      // integer index
            float index_fract = Cents - index_int;  // fractional part of index
            return pCentsToFreqTable[index_int] + index_fract * (pCentsToFreqTable[index_int+1] - pCentsToFreqTable[index_int]);
        }

        /**
         * Slower version of CentsToFreqRatio, for big values.
         *
         * @param cents - pitch value in cents (+1200 cents means +1 octave)
         * @returns  frequency ratio (e.g. +2.0 for +1 octave)
         */
        static double CentsToFreqRatioUnlimited(double Cents) {
            int octaves = int(Cents / 1200);
            double x = CentsToFreqRatio(Cents - octaves * 1200);
            return  octaves < 0 ? x / (1 << -octaves) : x * (1 << octaves);
        }

        /**
         * Inverse function to CentsToFreqRatio(). This function is a bit
         * slow, so it should not be called too frequently.
         */
        static double FreqRatioToCents(double FreqRatio) {
            return log(FreqRatio) / log(TWELVEHUNDREDTH_ROOT_OF_TWO);
        }

        /**
         * Calculates the linear ratio value representation (linear scale)
         * of the @a decibel value provided (exponential scale).
         *
         * The context of audio acoustic sound pressure levels is assumed, and
         * hence the field version of the dB unit is used here (which uses a
         * linear factor of 20). This function is a bit slow, so it should
         * not be called too frequently.
         *
         * @param decibel - sound pressure level in dB
         * @returns linear ratio of the supplied dB value
         * @see LinRatioToDecibel() as inverse function
         */
        static float DecibelToLinRatio(float decibel) {
            return powf(10.f, decibel / 20.f);
        }

        /**
         * Calculates the decibel value (exponential scale) of the @a linear
         * ratio value representation (linear scale) provided.
         *
         * The context of audio acoustic sound pressure levels is assumed, and
         * hence the field version of the dB unit is used here (which uses a
         * linear factor of 20). This function is a bit slow, so it should
         * not be called too frequently.
         *
         * @param linear - sound pressure level as linear ratio value (linear scale)
         * @returns dB value representation
         * @see DecibelToLinRatio() as inverse function
         */
        static float LinRatioToDecibel(float linear) {
            return 20.f * log10f(linear);
        }

        /**
         * Calculates the relatively summed average of a set of values.
         *
         * @param current - the current avaerage value of all previously summed values
         * @param sample - new value to be applied as summed average to the existing values
         * @param n - amount of sample values applied so far
         * @returns new average value of all summed values (including the new @a sample)
         */
        template<typename T_int>
        inline static float RelativeSummedAvg(float current, float sample, T_int n) {
            return current + (sample - current) / float(n);
        }

        /**
         * Clock source to use for getting the current time.
         */
        enum clock_source_t {
            real_clock,    ///< Use this to measure time that passed in reality (no matter if process got suspended).
            process_clock, ///< Use this to measure only the CPU execution time of the current process (if the process got suspended, the clock is paused as well).
            thread_clock,  ///< Use this to measure only the CPU execution time of the current thread (if the process got suspended or another thread is executed, the clock is paused as well).
        };

        /**
         * Returns a time stamp of the current time in microseconds (in
         * probably real-time @b unsafe way). There is no guarantee about
         * what the returned amount of microseconds relates to (i.e.
         * microseconds since epoch, microseconds since system uptime, ...).
         * So you should only use it to calculate time differences between
         * values taken with this method.
         *
         * @b CAUTION: This method may not @b NOT be real-time safe! On some
         * systems it could be RT safe, but there is no guarantee whatsoever!
         * So this method should only be used for debugging, benchmarking and
         * other developing purposes !
         *
         * For creating time stamps in real-time context, use
         * CreateTimeStamp() instead.
         *
         * @param source - the actual clock to use for getting the current
         *                 time, note that the various clock sources may not
         *                 be implemented on all systems
         * @returns time stamp in microseconds
         *
         * @see CreateTimeStamp()
         */
        static usecs_t unsafeMicroSeconds(clock_source_t source);

        /**
         * 'Equalness' comparison between two 32 bit floating point numbers
         * (that is of single precision data type).
         *
         *  This methods deals with the expected tolerance of the floating point
         *  representation.
         *
         * @return true if the two numbers can be expected to be equal
         */
        static bool fEqual32(float a, float b);

        /**
         * 'Equalness' comparison between two 64 bit floating point numbers
         * (that is of double precision data type).
         *
         *  This methods deals with the expected tolerance of the floating point
         *  representation.
         *
         * @return true if the two numbers can be expected to be equal
         */
        static bool fEqual64(double a, double b);

    private:
        static float* pCentsToFreqTable;

        static float* InitCentsToFreqTable();
};

/** @brief Real Time Math
 *
 * This is a template which provides customized methods for the desired low
 * level implementation. The ASM_X86_MMX_SSE implementation of each method
 * for example doesn't use 387 FPU instruction. This is needed for MMX
 * algorithms which do not allow mixed MMX and 387 instructions.
 */
template<implementation_t IMPL = CPP>
class __RTMath : public RTMathBase {
    public:
        // conversion using truncate
        inline static int Int(const float a) {
            switch (IMPL) {
                #if CONFIG_ASM && ARCH_X86
                case ASM_X86_MMX_SSE: {
                    int ret;
                    asm (
                        "cvttss2si %1, %0  # convert to int\n\t"
                        : "=r" (ret)
                        : "m" (a)
                    );
                    return ret;
                }
                #endif // CONFIG_ASM && ARCH_X86
                default: {
                    return (int) a;
                }
            }
        }

        //for doubles and everything else except floats
        template<class T_a> inline static int Int(const T_a a) {
            return (int) a;
        }

        inline static float Float(const int a) {
            switch (IMPL) {
                #if CONFIG_ASM && ARCH_X86
                case ASM_X86_MMX_SSE: {
                    float ret;
                    asm (
                        "cvtsi2ss %1, %%xmm0  # convert to float\n\t"
                        "movss    %%xmm0,%0   # output\n\t"
                        : "=m" (ret)
                        : "r" (a)
                    );
                    return ret;
                }
                #endif // CONFIG_ASM && ARCH_X86
                default: {
                    return (float) a;
                }
            }
        }

#if 0
        //for everything except ints
        template<class T_a> inline static float Float(T_a a) {
            return (float) a;
        }
#endif

        inline static float Sum(const float& a, const float& b) {
            switch (IMPL) {
                #if CONFIG_ASM && ARCH_X86
                case ASM_X86_MMX_SSE: {
                    float ret;
                    asm (
                        "movss    %1, %%xmm0  # load a\n\t"
                        "addss    %2, %%xmm0  # a + b\n\t"
                        "movss    %%xmm0, %0  # output\n\t"
                        : "=m" (ret)
                        : "m" (a), "m" (b)
                    );
                    return ret;
                }
                #endif // CONFIG_ASM && ARCH_X86
                default: {
                    return (a + b);
                }
            }
        }

        template<class T_a, class T_b> inline static T_a Sum(const T_a a, const T_b b) {
            return (a + b);
        }

        inline static float Sub(const float& a, const float& b) {
            switch (IMPL) {
                #if CONFIG_ASM && ARCH_X86
                case ASM_X86_MMX_SSE: {
                    float ret;
                    asm (
                        "movss    %1, %%xmm0  # load a\n\t"
                        "subss    %2, %%xmm0  # a - b\n\t"
                        "movss    %%xmm0, %0  # output\n\t"
                        : "=m" (ret)
                        : "m" (a), "m" (b)
                    );
                    return ret;
                }
                #endif // CONFIG_ASM && ARCH_X86
                default: {
                    return (a - b);
                }
            }
        }

        template<class T_a, class T_b> inline static T_a Sub(const T_a a, const T_b b) {
            return (a - b);
        }

        inline static float Mul(const float a, const float b) {
            switch (IMPL) {
                #if CONFIG_ASM && ARCH_X86
                case ASM_X86_MMX_SSE: {
                    float ret;
                    asm (
                        "movss    %1, %%xmm0  # load a\n\t"
                        "mulss    %2, %%xmm0  # a * b\n\t"
                        "movss    %%xmm0, %0  # output\n\t"
                        : "=m" (ret)
                        : "m" (a), "m" (b)
                    );
                    return ret;
                }
                #endif // CONFIG_ASM && ARCH_X86
                default: {
                    return (a * b);
                }
            }
        }

        template<class T_a, class T_b> inline static T_a Mul(const T_a a, const T_b b) {
            return (a * b);
        }

        inline static float Div(const float a, const float b) {
            switch (IMPL) {
                #if CONFIG_ASM && ARCH_X86
                case ASM_X86_MMX_SSE: {
                    float ret;
                    asm (
                        "movss    %1, %%xmm0  # load a\n\t"
                        "divss    %2, %%xmm0  # a / b\n\t"
                        "movss    %%xmm0, %0  # output\n\t"
                        : "=m" (ret)
                        : "m" (a), "m" (b)
                    );
                    return ret;
                }
                #endif // CONFIG_ASM && ARCH_X86
                default: {
                    return (a / b);
                }
            }
        }

        template<class T_a, class T_b> inline static T_a Div(const T_a a, const T_b b) {
            return (a / b);
        }

        inline static float Min(const float a, const float b) {
            switch (IMPL) {
                #if CONFIG_ASM && ARCH_X86
                case ASM_X86_MMX_SSE: {
                    float ret;
                    asm (
                        "movss    %1, %%xmm0  # load a\n\t"
                        "minss    %2, %%xmm0  # Minimum(a, b)\n\t"
                        "movss    %%xmm0, %0  # output\n\t"
                        : "=m" (ret)
                        : "m" (a), "m" (b)
                    );
                    return ret;
                }
                #endif // CONFIG_ASM && ARCH_X86
                default: {
                    return std::min(a, b);
                }
            }
        }

        template<class T_a, class T_b> inline static T_a Min(const T_a a, const T_b b) {
            return (b < a) ? b : a;
        }

        inline static float Max(const float a, const float b) {
            switch (IMPL) {
                #if CONFIG_ASM && ARCH_X86
                case ASM_X86_MMX_SSE: {
                    float ret;
                    asm (
                        "movss    %1, %%xmm0  # load a\n\t"
                        "maxss    %2, %%xmm0  # Maximum(a, b)\n\t"
                        "movss    %%xmm0, %0  # output\n\t"
                        : "=m" (ret)
                        : "m" (a), "m" (b)
                    );
                    return ret;
                }
                #endif // CONFIG_ASM && ARCH_X86
                default: {
                    return std::max(a, b);
                }
            }
        }

        template<class T_a, class T_b> inline static T_a Max(const T_a a, const T_b b) {
            return (b > a) ? b : a;
        }

        inline static float Fmodf(const float &a, const float &b) {
            switch (IMPL) {
                #if CONFIG_ASM && ARCH_X86
                case ASM_X86_MMX_SSE: {
                    float ret;
                    asm (
                        "movss    %1, %%xmm0  # load a\n\t"
                        "movss    %2, %%xmm1  # load b\n\t"
                        "movss    %%xmm0,%%xmm2\n\t"
                        "divss    %%xmm1, %%xmm2  # xmm2 = a / b\n\t"
                        "cvttss2si %%xmm2, %%ecx  #convert to int\n\t"
                        "cvtsi2ss %%ecx, %%xmm2  #convert back to float\n\t"
                        "mulss    %%xmm1, %%xmm2  # xmm2 = b * int(a/b)\n\t"
                        "subss    %%xmm2, %%xmm0  #sub a\n\t"
                        "movss    %%xmm0, %0  # output\n\t"
                        : "=m" (ret)
                        : "m" (a), "m" (b)
                        : "%ecx"
                    );
                    return ret;
                }
                #endif // CONFIG_ASM && ARCH_X86
                default: {
                    return fmodf(a, b);
                }
            }
        }
};

/// convenience typedef for using the default implementation (which is CPP)
typedef __RTMath<> RTMath;

#endif // __RT_MATH_H__
1	/***************************************************************************
2	* *
3	* LinuxSampler - modular, streaming capable sampler *
4	* *
5	* Copyright (C) 2003, 2004 by Benno Senoner and Christian Schoenebeck *
6	* Copyright (C) 2005 - 2019 Christian Schoenebeck *
7	* *
8	* This program is free software; you can redistribute it and/or modify *
9	* it under the terms of the GNU General Public License as published by *
10	* the Free Software Foundation; either version 2 of the License, or *
11	* (at your option) any later version. *
12	* *
13	* This program is distributed in the hope that it will be useful, *
14	* but WITHOUT ANY WARRANTY; without even the implied warranty of *
15	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
16	* GNU General Public License for more details. *
17	* *
18	* You should have received a copy of the GNU General Public License *
19	* along with this program; if not, write to the Free Software *
20	* Foundation, Inc., 59 Temple Place, Suite 330, Boston, *
21	* MA 02111-1307 USA *
22	***************************************************************************/
23
24	#ifndef __RT_MATH_H__
25	#define __RT_MATH_H__
26
27	#include <math.h>
28	#include <stdint.h>
29	#include "global_private.h"
30
31	/// Needed for calculating frequency ratio used to pitch a sample
32	#define TWELVEHUNDREDTH_ROOT_OF_TWO 1.000577789506555
33
34	enum implementation_t {
35	CPP,
36	ASM_X86_MMX_SSE
37	};
38
39	/** @brief Real Time Math Base Class
40	*
41	* Math functions for real time operation. This base class contains all
42	* non-template methods.
43	*/
44	class RTMathBase {
45	public:
46	/**
47	* High resolution time stamp.
48	*/
49	typedef uint32_t time_stamp_t;
50
51	typedef uint64_t usecs_t;
52
53	/**
54	* We read the processor's cycle count register as a reference
55	* for the real time. These are of course only abstract values
56	* with arbitrary time entity, but that's not a problem as long
57	* as we calculate relatively.
58	*
59	* @see unsafeMicroSeconds()
60	*/
61	static time_stamp_t CreateTimeStamp();
62
63	/**
64	* Calculates the frequency ratio for a pitch value given in cents
65	* (assuming equal tempered scale of course, divided into 12
66	* semitones per octave and 100 cents per semitone).
67	*
68	* Note: CONFIG_MAX_PITCH (defined in config.h) has to be defined to an
69	* appropriate value, otherwise the behavior of this function is
70	* undefined, but most probably if CONFIG_MAX_PITCH is too small, the
71	* application will crash due to segmentation fault here.
72	*
73	* @param cents - pitch value in cents (+1200 cents means +1 octave)
74	* @returns frequency ratio (e.g. +2.0 for +1 octave)
75	*/
76	inline static double CentsToFreqRatio(double Cents) {
77	int index_int = (int) (Cents); // integer index
78	float index_fract = Cents - index_int; // fractional part of index
79	return pCentsToFreqTable[index_int] + index_fract * (pCentsToFreqTable[index_int+1] - pCentsToFreqTable[index_int]);
80	}
81
82	/**
83	* Slower version of CentsToFreqRatio, for big values.
84	*
85	* @param cents - pitch value in cents (+1200 cents means +1 octave)
86	* @returns frequency ratio (e.g. +2.0 for +1 octave)
87	*/
88	static double CentsToFreqRatioUnlimited(double Cents) {
89	int octaves = int(Cents / 1200);
90	double x = CentsToFreqRatio(Cents - octaves * 1200);
91	return octaves < 0 ? x / (1 << -octaves) : x * (1 << octaves);
92	}
93
94	/**
95	* Inverse function to CentsToFreqRatio(). This function is a bit
96	* slow, so it should not be called too frequently.
97	*/
98	static double FreqRatioToCents(double FreqRatio) {
99	return log(FreqRatio) / log(TWELVEHUNDREDTH_ROOT_OF_TWO);
100	}
101
102	/**
103	* Calculates the linear ratio value representation (linear scale)
104	* of the @a decibel value provided (exponential scale).
105	*
106	* The context of audio acoustic sound pressure levels is assumed, and
107	* hence the field version of the dB unit is used here (which uses a
108	* linear factor of 20). This function is a bit slow, so it should
109	* not be called too frequently.
110	*
111	* @param decibel - sound pressure level in dB
112	* @returns linear ratio of the supplied dB value
113	* @see LinRatioToDecibel() as inverse function
114	*/
115	static float DecibelToLinRatio(float decibel) {
116	return powf(10.f, decibel / 20.f);
117	}
118
119	/**
120	* Calculates the decibel value (exponential scale) of the @a linear
121	* ratio value representation (linear scale) provided.
122	*
123	* The context of audio acoustic sound pressure levels is assumed, and
124	* hence the field version of the dB unit is used here (which uses a
125	* linear factor of 20). This function is a bit slow, so it should
126	* not be called too frequently.
127	*
128	* @param linear - sound pressure level as linear ratio value (linear scale)
129	* @returns dB value representation
130	* @see DecibelToLinRatio() as inverse function
131	*/
132	static float LinRatioToDecibel(float linear) {
133	return 20.f * log10f(linear);
134	}
135
136	/**
137	* Calculates the relatively summed average of a set of values.
138	*
139	* @param current - the current avaerage value of all previously summed values
140	* @param sample - new value to be applied as summed average to the existing values
141	* @param n - amount of sample values applied so far
142	* @returns new average value of all summed values (including the new @a sample)
143	*/
144	template<typename T_int>
145	inline static float RelativeSummedAvg(float current, float sample, T_int n) {
146	return current + (sample - current) / float(n);
147	}
148
149	/**
150	* Clock source to use for getting the current time.
151	*/
152	enum clock_source_t {
153	real_clock, ///< Use this to measure time that passed in reality (no matter if process got suspended).
154	process_clock, ///< Use this to measure only the CPU execution time of the current process (if the process got suspended, the clock is paused as well).
155	thread_clock, ///< Use this to measure only the CPU execution time of the current thread (if the process got suspended or another thread is executed, the clock is paused as well).
156	};
157
158	/**
159	* Returns a time stamp of the current time in microseconds (in
160	* probably real-time @b unsafe way). There is no guarantee about
161	* what the returned amount of microseconds relates to (i.e.
162	* microseconds since epoch, microseconds since system uptime, ...).
163	* So you should only use it to calculate time differences between
164	* values taken with this method.
165	*
166	* @b CAUTION: This method may not @b NOT be real-time safe! On some
167	* systems it could be RT safe, but there is no guarantee whatsoever!
168	* So this method should only be used for debugging, benchmarking and
169	* other developing purposes !
170	*
171	* For creating time stamps in real-time context, use
172	* CreateTimeStamp() instead.
173	*
174	* @param source - the actual clock to use for getting the current
175	* time, note that the various clock sources may not
176	* be implemented on all systems
177	* @returns time stamp in microseconds
178	*
179	* @see CreateTimeStamp()
180	*/
181	static usecs_t unsafeMicroSeconds(clock_source_t source);
182
183	/**
184	* 'Equalness' comparison between two 32 bit floating point numbers
185	* (that is of single precision data type).
186	*
187	* This methods deals with the expected tolerance of the floating point
188	* representation.
189	*
190	* @return true if the two numbers can be expected to be equal
191	*/
192	static bool fEqual32(float a, float b);
193
194	/**
195	* 'Equalness' comparison between two 64 bit floating point numbers
196	* (that is of double precision data type).
197	*
198	* This methods deals with the expected tolerance of the floating point
199	* representation.
200	*
201	* @return true if the two numbers can be expected to be equal
202	*/
203	static bool fEqual64(double a, double b);
204
205	private:
206	static float* pCentsToFreqTable;
207
208	static float* InitCentsToFreqTable();
209	};
210
211	/** @brief Real Time Math
212	*
213	* This is a template which provides customized methods for the desired low
214	* level implementation. The ASM_X86_MMX_SSE implementation of each method
215	* for example doesn't use 387 FPU instruction. This is needed for MMX
216	* algorithms which do not allow mixed MMX and 387 instructions.
217	*/
218	template<implementation_t IMPL = CPP>
219	class __RTMath : public RTMathBase {
220	public:
221	// conversion using truncate
222	inline static int Int(const float a) {
223	switch (IMPL) {
224	#if CONFIG_ASM && ARCH_X86
225	case ASM_X86_MMX_SSE: {
226	int ret;
227	asm (
228	"cvttss2si %1, %0 # convert to int\n\t"
229	: "=r" (ret)
230	: "m" (a)
231	);
232	return ret;
233	}
234	#endif // CONFIG_ASM && ARCH_X86
235	default: {
236	return (int) a;
237	}
238	}
239	}
240
241	//for doubles and everything else except floats
242	template<class T_a> inline static int Int(const T_a a) {
243	return (int) a;
244	}
245
246	inline static float Float(const int a) {
247	switch (IMPL) {
248	#if CONFIG_ASM && ARCH_X86
249	case ASM_X86_MMX_SSE: {
250	float ret;
251	asm (
252	"cvtsi2ss %1, %%xmm0 # convert to float\n\t"
253	"movss %%xmm0,%0 # output\n\t"
254	: "=m" (ret)
255	: "r" (a)
256	);
257	return ret;
258	}
259	#endif // CONFIG_ASM && ARCH_X86
260	default: {
261	return (float) a;
262	}
263	}
264	}
265
266	#if 0
267	//for everything except ints
268	template<class T_a> inline static float Float(T_a a) {
269	return (float) a;
270	}
271	#endif
272
273	inline static float Sum(const float& a, const float& b) {
274	switch (IMPL) {
275	#if CONFIG_ASM && ARCH_X86
276	case ASM_X86_MMX_SSE: {
277	float ret;
278	asm (
279	"movss %1, %%xmm0 # load a\n\t"
280	"addss %2, %%xmm0 # a + b\n\t"
281	"movss %%xmm0, %0 # output\n\t"
282	: "=m" (ret)
283	: "m" (a), "m" (b)
284	);
285	return ret;
286	}
287	#endif // CONFIG_ASM && ARCH_X86
288	default: {
289	return (a + b);
290	}
291	}
292	}
293
294	template<class T_a, class T_b> inline static T_a Sum(const T_a a, const T_b b) {
295	return (a + b);
296	}
297
298	inline static float Sub(const float& a, const float& b) {
299	switch (IMPL) {
300	#if CONFIG_ASM && ARCH_X86
301	case ASM_X86_MMX_SSE: {
302	float ret;
303	asm (
304	"movss %1, %%xmm0 # load a\n\t"
305	"subss %2, %%xmm0 # a - b\n\t"
306	"movss %%xmm0, %0 # output\n\t"
307	: "=m" (ret)
308	: "m" (a), "m" (b)
309	);
310	return ret;
311	}
312	#endif // CONFIG_ASM && ARCH_X86
313	default: {
314	return (a - b);
315	}
316	}
317	}
318
319	template<class T_a, class T_b> inline static T_a Sub(const T_a a, const T_b b) {
320	return (a - b);
321	}
322
323	inline static float Mul(const float a, const float b) {
324	switch (IMPL) {
325	#if CONFIG_ASM && ARCH_X86
326	case ASM_X86_MMX_SSE: {
327	float ret;
328	asm (
329	"movss %1, %%xmm0 # load a\n\t"
330	"mulss %2, %%xmm0 # a * b\n\t"
331	"movss %%xmm0, %0 # output\n\t"
332	: "=m" (ret)
333	: "m" (a), "m" (b)
334	);
335	return ret;
336	}
337	#endif // CONFIG_ASM && ARCH_X86
338	default: {
339	return (a * b);
340	}
341	}
342	}
343
344	template<class T_a, class T_b> inline static T_a Mul(const T_a a, const T_b b) {
345	return (a * b);
346	}
347
348	inline static float Div(const float a, const float b) {
349	switch (IMPL) {
350	#if CONFIG_ASM && ARCH_X86
351	case ASM_X86_MMX_SSE: {
352	float ret;
353	asm (
354	"movss %1, %%xmm0 # load a\n\t"
355	"divss %2, %%xmm0 # a / b\n\t"
356	"movss %%xmm0, %0 # output\n\t"
357	: "=m" (ret)
358	: "m" (a), "m" (b)
359	);
360	return ret;
361	}
362	#endif // CONFIG_ASM && ARCH_X86
363	default: {
364	return (a / b);
365	}
366	}
367	}
368
369	template<class T_a, class T_b> inline static T_a Div(const T_a a, const T_b b) {
370	return (a / b);
371	}
372
373	inline static float Min(const float a, const float b) {
374	switch (IMPL) {
375	#if CONFIG_ASM && ARCH_X86
376	case ASM_X86_MMX_SSE: {
377	float ret;
378	asm (
379	"movss %1, %%xmm0 # load a\n\t"
380	"minss %2, %%xmm0 # Minimum(a, b)\n\t"
381	"movss %%xmm0, %0 # output\n\t"
382	: "=m" (ret)
383	: "m" (a), "m" (b)
384	);
385	return ret;
386	}
387	#endif // CONFIG_ASM && ARCH_X86
388	default: {
389	return std::min(a, b);
390	}
391	}
392	}
393
394	template<class T_a, class T_b> inline static T_a Min(const T_a a, const T_b b) {
395	return (b < a) ? b : a;
396	}
397
398	inline static float Max(const float a, const float b) {
399	switch (IMPL) {
400	#if CONFIG_ASM && ARCH_X86
401	case ASM_X86_MMX_SSE: {
402	float ret;
403	asm (
404	"movss %1, %%xmm0 # load a\n\t"
405	"maxss %2, %%xmm0 # Maximum(a, b)\n\t"
406	"movss %%xmm0, %0 # output\n\t"
407	: "=m" (ret)
408	: "m" (a), "m" (b)
409	);
410	return ret;
411	}
412	#endif // CONFIG_ASM && ARCH_X86
413	default: {
414	return std::max(a, b);
415	}
416	}
417	}
418
419	template<class T_a, class T_b> inline static T_a Max(const T_a a, const T_b b) {
420	return (b > a) ? b : a;
421	}
422
423	inline static float Fmodf(const float &a, const float &b) {
424	switch (IMPL) {
425	#if CONFIG_ASM && ARCH_X86
426	case ASM_X86_MMX_SSE: {
427	float ret;
428	asm (
429	"movss %1, %%xmm0 # load a\n\t"
430	"movss %2, %%xmm1 # load b\n\t"
431	"movss %%xmm0,%%xmm2\n\t"
432	"divss %%xmm1, %%xmm2 # xmm2 = a / b\n\t"
433	"cvttss2si %%xmm2, %%ecx #convert to int\n\t"
434	"cvtsi2ss %%ecx, %%xmm2 #convert back to float\n\t"
435	"mulss %%xmm1, %%xmm2 # xmm2 = b * int(a/b)\n\t"
436	"subss %%xmm2, %%xmm0 #sub a\n\t"
437	"movss %%xmm0, %0 # output\n\t"
438	: "=m" (ret)
439	: "m" (a), "m" (b)
440	: "%ecx"
441	);
442	return ret;
443	}
444	#endif // CONFIG_ASM && ARCH_X86
445	default: {
446	return fmodf(a, b);
447	}
448	}
449	}
450	};
451
452	/// convenience typedef for using the default implementation (which is CPP)
453	typedef __RTMath<> RTMath;
454
455	#endif // __RT_MATH_H__