src/common/RTMath.h

/***************************************************************************
 *                                                                         *
 *   LinuxSampler - modular, streaming capable sampler                     *
 *                                                                         *
 *   Copyright (C) 2003, 2004 by Benno Senoner and Christian Schoenebeck   *
 *   Copyright (C) 2005 - 2007 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:
        /**
         * Highly accurate time stamp.
         */
        typedef uint32_t time_stamp_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.
         */
        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]);
        }

        /**
         * 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);
        }

    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	schoenebeck	53	/***************************************************************************
2			* *
3			* LinuxSampler - modular, streaming capable sampler *
4			* *
5	schoenebeck	56	* Copyright (C) 2003, 2004 by Benno Senoner and Christian Schoenebeck *
6	schoenebeck	1212	* Copyright (C) 2005 - 2007 Christian Schoenebeck *
7	schoenebeck	53	* *
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	schoenebeck	328	#include <stdint.h>
29	schoenebeck	1424	#include "global_private.h"
30	schoenebeck	53
31			/// Needed for calculating frequency ratio used to pitch a sample
32			#define TWELVEHUNDREDTH_ROOT_OF_TWO 1.000577789506555
33
34	schoenebeck	319	enum implementation_t {
35	schoenebeck	361	CPP,
36			ASM_X86_MMX_SSE
37	schoenebeck	319	};
38
39	schoenebeck	563	/** @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	schoenebeck	319	class RTMathBase {
45	schoenebeck	53	public:
46			/**
47	schoenebeck	328	* Highly accurate time stamp.
48			*/
49			typedef uint32_t time_stamp_t;
50
51			/**
52			* We read the processor's cycle count register as a reference
53			* for the real time. These are of course only abstract values
54			* with arbitrary time entity, but that's not a problem as long
55			* as we calculate relatively.
56			*/
57			static time_stamp_t CreateTimeStamp();
58
59			/**
60	schoenebeck	53	* Calculates the frequency ratio for a pitch value given in cents
61			* (assuming equal tempered scale of course, divided into 12
62			* semitones per octave and 100 cents per semitone).
63			*
64	schoenebeck	554	* Note: CONFIG_MAX_PITCH (defined in config.h) has to be defined to an
65	schoenebeck	53	* appropriate value, otherwise the behavior of this function is
66	schoenebeck	554	* undefined, but most probably if CONFIG_MAX_PITCH is too small, the
67	schoenebeck	53	* application will crash due to segmentation fault here.
68			*
69			* @param cents - pitch value in cents (+1200 cents means +1 octave)
70			* @returns frequency ratio (e.g. +2.0 for +1 octave)
71			*/
72			inline static double CentsToFreqRatio(double Cents) {
73	schoenebeck	319	int index_int = (int) (Cents); // integer index
74	schoenebeck	53	float index_fract = Cents - index_int; // fractional part of index
75			return pCentsToFreqTable[index_int] + index_fract * (pCentsToFreqTable[index_int+1] - pCentsToFreqTable[index_int]);
76			}
77
78	schoenebeck	829	/**
79			* Inverse function to CentsToFreqRatio(). This function is a bit
80			* slow, so it should not be called too frequently.
81			*/
82			static double FreqRatioToCents(double FreqRatio) {
83			return log(FreqRatio) / log(TWELVEHUNDREDTH_ROOT_OF_TWO);
84			}
85
86	schoenebeck	319	private:
87			static float* pCentsToFreqTable;
88
89			static float* InitCentsToFreqTable();
90			};
91
92	schoenebeck	563	/** @brief Real Time Math
93	schoenebeck	319	*
94	schoenebeck	563	* This is a template which provides customized methods for the desired low
95			* level implementation. The ASM_X86_MMX_SSE implementation of each method
96			* for example doesn't use 387 FPU instruction. This is needed for MMX
97			* algorithms which do not allow mixed MMX and 387 instructions.
98	schoenebeck	319	*/
99			template<implementation_t IMPL = CPP>
100			class __RTMath : public RTMathBase {
101			public:
102			// conversion using truncate
103			inline static int Int(const float a) {
104			switch (IMPL) {
105	schoenebeck	617	#if CONFIG_ASM && ARCH_X86
106	schoenebeck	319	case ASM_X86_MMX_SSE: {
107			int ret;
108			asm (
109			"cvttss2si %1, %0 # convert to int\n\t"
110			: "=r" (ret)
111			: "m" (a)
112			);
113			return ret;
114			}
115	schoenebeck	617	#endif // CONFIG_ASM && ARCH_X86
116	persson	685	default: {
117			return (int) a;
118			}
119	schoenebeck	319	}
120			}
121
122			//for doubles and everything else except floats
123			template<class T_a> inline static int Int(const T_a a) {
124			return (int) a;
125			}
126
127			inline static float Float(const int a) {
128			switch (IMPL) {
129	schoenebeck	617	#if CONFIG_ASM && ARCH_X86
130	schoenebeck	319	case ASM_X86_MMX_SSE: {
131			float ret;
132			asm (
133			"cvtsi2ss %1, %%xmm0 # convert to float\n\t"
134			"movss %%xmm0,%0 # output\n\t"
135			: "=m" (ret)
136			: "r" (a)
137			);
138			return ret;
139			}
140	schoenebeck	617	#endif // CONFIG_ASM && ARCH_X86
141	persson	685	default: {
142			return (float) a;
143			}
144	schoenebeck	319	}
145			}
146
147			#if 0
148			//for everything except ints
149			template<class T_a> inline static float Float(T_a a) {
150			return (float) a;
151			}
152			#endif
153
154			inline static float Sum(const float& a, const float& b) {
155			switch (IMPL) {
156	schoenebeck	617	#if CONFIG_ASM && ARCH_X86
157	schoenebeck	319	case ASM_X86_MMX_SSE: {
158			float ret;
159			asm (
160			"movss %1, %%xmm0 # load a\n\t"
161			"addss %2, %%xmm0 # a + b\n\t"
162			"movss %%xmm0, %0 # output\n\t"
163			: "=m" (ret)
164			: "m" (a), "m" (b)
165			);
166			return ret;
167			}
168	schoenebeck	617	#endif // CONFIG_ASM && ARCH_X86
169	persson	685	default: {
170			return (a + b);
171			}
172	schoenebeck	319	}
173			}
174
175			template<class T_a, class T_b> inline static T_a Sum(const T_a a, const T_b b) {
176			return (a + b);
177			}
178
179			inline static float Sub(const float& a, const float& b) {
180			switch (IMPL) {
181	schoenebeck	617	#if CONFIG_ASM && ARCH_X86
182	schoenebeck	319	case ASM_X86_MMX_SSE: {
183			float ret;
184			asm (
185			"movss %1, %%xmm0 # load a\n\t"
186			"subss %2, %%xmm0 # a - b\n\t"
187			"movss %%xmm0, %0 # output\n\t"
188			: "=m" (ret)
189			: "m" (a), "m" (b)
190			);
191			return ret;
192			}
193	schoenebeck	617	#endif // CONFIG_ASM && ARCH_X86
194	persson	685	default: {
195			return (a - b);
196			}
197	schoenebeck	319	}
198			}
199
200			template<class T_a, class T_b> inline static T_a Sub(const T_a a, const T_b b) {
201			return (a - b);
202			}
203
204			inline static float Mul(const float a, const float b) {
205			switch (IMPL) {
206	schoenebeck	617	#if CONFIG_ASM && ARCH_X86
207	schoenebeck	319	case ASM_X86_MMX_SSE: {
208			float ret;
209			asm (
210			"movss %1, %%xmm0 # load a\n\t"
211			"mulss %2, %%xmm0 # a * b\n\t"
212			"movss %%xmm0, %0 # output\n\t"
213			: "=m" (ret)
214			: "m" (a), "m" (b)
215			);
216			return ret;
217			}
218	schoenebeck	617	#endif // CONFIG_ASM && ARCH_X86
219	persson	685	default: {
220			return (a * b);
221			}
222	schoenebeck	319	}
223			}
224
225			template<class T_a, class T_b> inline static T_a Mul(const T_a a, const T_b b) {
226			return (a * b);
227			}
228
229			inline static float Div(const float a, const float b) {
230			switch (IMPL) {
231	schoenebeck	617	#if CONFIG_ASM && ARCH_X86
232	schoenebeck	319	case ASM_X86_MMX_SSE: {
233			float ret;
234			asm (
235			"movss %1, %%xmm0 # load a\n\t"
236			"divss %2, %%xmm0 # a / b\n\t"
237			"movss %%xmm0, %0 # output\n\t"
238			: "=m" (ret)
239			: "m" (a), "m" (b)
240			);
241			return ret;
242			}
243	schoenebeck	617	#endif // CONFIG_ASM && ARCH_X86
244	persson	685	default: {
245			return (a / b);
246			}
247	schoenebeck	319	}
248			}
249
250			template<class T_a, class T_b> inline static T_a Div(const T_a a, const T_b b) {
251			return (a / b);
252			}
253
254			inline static float Min(const float a, const float b) {
255			switch (IMPL) {
256	schoenebeck	617	#if CONFIG_ASM && ARCH_X86
257	schoenebeck	319	case ASM_X86_MMX_SSE: {
258			float ret;
259			asm (
260			"movss %1, %%xmm0 # load a\n\t"
261			"minss %2, %%xmm0 # Minimum(a, b)\n\t"
262			"movss %%xmm0, %0 # output\n\t"
263			: "=m" (ret)
264			: "m" (a), "m" (b)
265			);
266			return ret;
267			}
268	schoenebeck	617	#endif // CONFIG_ASM && ARCH_X86
269	persson	685	default: {
270	wylder	818	return std::min(a, b);
271	persson	685	}
272	schoenebeck	319	}
273			}
274
275			template<class T_a, class T_b> inline static T_a Min(const T_a a, const T_b b) {
276	schoenebeck	53	return (b < a) ? b : a;
277			}
278
279	schoenebeck	319	inline static float Max(const float a, const float b) {
280			switch (IMPL) {
281	schoenebeck	617	#if CONFIG_ASM && ARCH_X86
282	schoenebeck	319	case ASM_X86_MMX_SSE: {
283			float ret;
284			asm (
285			"movss %1, %%xmm0 # load a\n\t"
286			"maxss %2, %%xmm0 # Maximum(a, b)\n\t"
287			"movss %%xmm0, %0 # output\n\t"
288			: "=m" (ret)
289			: "m" (a), "m" (b)
290			);
291			return ret;
292			}
293	schoenebeck	617	#endif // CONFIG_ASM && ARCH_X86
294	persson	685	default: {
295	wylder	818	return std::max(a, b);
296	persson	685	}
297	schoenebeck	319	}
298			}
299
300			template<class T_a, class T_b> inline static T_a Max(const T_a a, const T_b b) {
301	schoenebeck	53	return (b > a) ? b : a;
302			}
303
304	schoenebeck	319	inline static float Fmodf(const float &a, const float &b) {
305			switch (IMPL) {
306	schoenebeck	617	#if CONFIG_ASM && ARCH_X86
307	schoenebeck	319	case ASM_X86_MMX_SSE: {
308			float ret;
309			asm (
310			"movss %1, %%xmm0 # load a\n\t"
311			"movss %2, %%xmm1 # load b\n\t"
312			"movss %%xmm0,%%xmm2\n\t"
313			"divss %%xmm1, %%xmm2 # xmm2 = a / b\n\t"
314			"cvttss2si %%xmm2, %%ecx #convert to int\n\t"
315			"cvtsi2ss %%ecx, %%xmm2 #convert back to float\n\t"
316			"mulss %%xmm1, %%xmm2 # xmm2 = b * int(a/b)\n\t"
317			"subss %%xmm2, %%xmm0 #sub a\n\t"
318			"movss %%xmm0, %0 # output\n\t"
319			: "=m" (ret)
320			: "m" (a), "m" (b)
321			: "%ecx"
322			);
323			return ret;
324			}
325	schoenebeck	617	#endif // CONFIG_ASM && ARCH_X86
326	persson	685	default: {
327			return fmodf(a, b);
328			}
329	schoenebeck	319	}
330			}
331	schoenebeck	53	};
332
333	schoenebeck	319	/// convenience typedef for using the default implementation (which is CPP)
334			typedef __RTMath<> RTMath;
335
336	schoenebeck	53	#endif // __RT_MATH_H__