src/common/RTMath.h

/***************************************************************************
 *                                                                         *
 *   LinuxSampler - modular, streaming capable sampler                     *
 *                                                                         *
 *   Copyright (C) 2003, 2004 by Benno Senoner and Christian Schoenebeck   *
 *   Copyright (C) 2005, 2006 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.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  CentsToFreqTable[CONFIG_MAX_PITCH * 1200 * 2 + 1];
        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, 2006 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.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	* 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	* 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	* Note: CONFIG_MAX_PITCH (defined in config.h) has to be defined to an
65	* appropriate value, otherwise the behavior of this function is
66	* undefined, but most probably if CONFIG_MAX_PITCH is too small, the
67	* 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	int index_int = (int) (Cents); // integer index
74	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	/**
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	private:
87	static float CentsToFreqTable[CONFIG_MAX_PITCH * 1200 * 2 + 1];
88	static float* pCentsToFreqTable;
89
90	static float* InitCentsToFreqTable();
91	};
92
93	/** @brief Real Time Math
94	*
95	* This is a template which provides customized methods for the desired low
96	* level implementation. The ASM_X86_MMX_SSE implementation of each method
97	* for example doesn't use 387 FPU instruction. This is needed for MMX
98	* algorithms which do not allow mixed MMX and 387 instructions.
99	*/
100	template<implementation_t IMPL = CPP>
101	class __RTMath : public RTMathBase {
102	public:
103	// conversion using truncate
104	inline static int Int(const float a) {
105	switch (IMPL) {
106	#if CONFIG_ASM && ARCH_X86
107	case ASM_X86_MMX_SSE: {
108	int ret;
109	asm (
110	"cvttss2si %1, %0 # convert to int\n\t"
111	: "=r" (ret)
112	: "m" (a)
113	);
114	return ret;
115	}
116	#endif // CONFIG_ASM && ARCH_X86
117	default: {
118	return (int) a;
119	}
120	}
121	}
122
123	//for doubles and everything else except floats
124	template<class T_a> inline static int Int(const T_a a) {
125	return (int) a;
126	}
127
128	inline static float Float(const int a) {
129	switch (IMPL) {
130	#if CONFIG_ASM && ARCH_X86
131	case ASM_X86_MMX_SSE: {
132	float ret;
133	asm (
134	"cvtsi2ss %1, %%xmm0 # convert to float\n\t"
135	"movss %%xmm0,%0 # output\n\t"
136	: "=m" (ret)
137	: "r" (a)
138	);
139	return ret;
140	}
141	#endif // CONFIG_ASM && ARCH_X86
142	default: {
143	return (float) a;
144	}
145	}
146	}
147
148	#if 0
149	//for everything except ints
150	template<class T_a> inline static float Float(T_a a) {
151	return (float) a;
152	}
153	#endif
154
155	inline static float Sum(const float& a, const float& b) {
156	switch (IMPL) {
157	#if CONFIG_ASM && ARCH_X86
158	case ASM_X86_MMX_SSE: {
159	float ret;
160	asm (
161	"movss %1, %%xmm0 # load a\n\t"
162	"addss %2, %%xmm0 # a + b\n\t"
163	"movss %%xmm0, %0 # output\n\t"
164	: "=m" (ret)
165	: "m" (a), "m" (b)
166	);
167	return ret;
168	}
169	#endif // CONFIG_ASM && ARCH_X86
170	default: {
171	return (a + b);
172	}
173	}
174	}
175
176	template<class T_a, class T_b> inline static T_a Sum(const T_a a, const T_b b) {
177	return (a + b);
178	}
179
180	inline static float Sub(const float& a, const float& b) {
181	switch (IMPL) {
182	#if CONFIG_ASM && ARCH_X86
183	case ASM_X86_MMX_SSE: {
184	float ret;
185	asm (
186	"movss %1, %%xmm0 # load a\n\t"
187	"subss %2, %%xmm0 # a - b\n\t"
188	"movss %%xmm0, %0 # output\n\t"
189	: "=m" (ret)
190	: "m" (a), "m" (b)
191	);
192	return ret;
193	}
194	#endif // CONFIG_ASM && ARCH_X86
195	default: {
196	return (a - b);
197	}
198	}
199	}
200
201	template<class T_a, class T_b> inline static T_a Sub(const T_a a, const T_b b) {
202	return (a - b);
203	}
204
205	inline static float Mul(const float a, const float b) {
206	switch (IMPL) {
207	#if CONFIG_ASM && ARCH_X86
208	case ASM_X86_MMX_SSE: {
209	float ret;
210	asm (
211	"movss %1, %%xmm0 # load a\n\t"
212	"mulss %2, %%xmm0 # a * b\n\t"
213	"movss %%xmm0, %0 # output\n\t"
214	: "=m" (ret)
215	: "m" (a), "m" (b)
216	);
217	return ret;
218	}
219	#endif // CONFIG_ASM && ARCH_X86
220	default: {
221	return (a * b);
222	}
223	}
224	}
225
226	template<class T_a, class T_b> inline static T_a Mul(const T_a a, const T_b b) {
227	return (a * b);
228	}
229
230	inline static float Div(const float a, const float b) {
231	switch (IMPL) {
232	#if CONFIG_ASM && ARCH_X86
233	case ASM_X86_MMX_SSE: {
234	float ret;
235	asm (
236	"movss %1, %%xmm0 # load a\n\t"
237	"divss %2, %%xmm0 # a / b\n\t"
238	"movss %%xmm0, %0 # output\n\t"
239	: "=m" (ret)
240	: "m" (a), "m" (b)
241	);
242	return ret;
243	}
244	#endif // CONFIG_ASM && ARCH_X86
245	default: {
246	return (a / b);
247	}
248	}
249	}
250
251	template<class T_a, class T_b> inline static T_a Div(const T_a a, const T_b b) {
252	return (a / b);
253	}
254
255	inline static float Min(const float a, const float b) {
256	switch (IMPL) {
257	#if CONFIG_ASM && ARCH_X86
258	case ASM_X86_MMX_SSE: {
259	float ret;
260	asm (
261	"movss %1, %%xmm0 # load a\n\t"
262	"minss %2, %%xmm0 # Minimum(a, b)\n\t"
263	"movss %%xmm0, %0 # output\n\t"
264	: "=m" (ret)
265	: "m" (a), "m" (b)
266	);
267	return ret;
268	}
269	#endif // CONFIG_ASM && ARCH_X86
270	default: {
271	return std::min(a, b);
272	}
273	}
274	}
275
276	template<class T_a, class T_b> inline static T_a Min(const T_a a, const T_b b) {
277	return (b < a) ? b : a;
278	}
279
280	inline static float Max(const float a, const float b) {
281	switch (IMPL) {
282	#if CONFIG_ASM && ARCH_X86
283	case ASM_X86_MMX_SSE: {
284	float ret;
285	asm (
286	"movss %1, %%xmm0 # load a\n\t"
287	"maxss %2, %%xmm0 # Maximum(a, b)\n\t"
288	"movss %%xmm0, %0 # output\n\t"
289	: "=m" (ret)
290	: "m" (a), "m" (b)
291	);
292	return ret;
293	}
294	#endif // CONFIG_ASM && ARCH_X86
295	default: {
296	return std::max(a, b);
297	}
298	}
299	}
300
301	template<class T_a, class T_b> inline static T_a Max(const T_a a, const T_b b) {
302	return (b > a) ? b : a;
303	}
304
305	inline static float Fmodf(const float &a, const float &b) {
306	switch (IMPL) {
307	#if CONFIG_ASM && ARCH_X86
308	case ASM_X86_MMX_SSE: {
309	float ret;
310	asm (
311	"movss %1, %%xmm0 # load a\n\t"
312	"movss %2, %%xmm1 # load b\n\t"
313	"movss %%xmm0,%%xmm2\n\t"
314	"divss %%xmm1, %%xmm2 # xmm2 = a / b\n\t"
315	"cvttss2si %%xmm2, %%ecx #convert to int\n\t"
316	"cvtsi2ss %%ecx, %%xmm2 #convert back to float\n\t"
317	"mulss %%xmm1, %%xmm2 # xmm2 = b * int(a/b)\n\t"
318	"subss %%xmm2, %%xmm0 #sub a\n\t"
319	"movss %%xmm0, %0 # output\n\t"
320	: "=m" (ret)
321	: "m" (a), "m" (b)
322	: "%ecx"
323	);
324	return ret;
325	}
326	#endif // CONFIG_ASM && ARCH_X86
327	default: {
328	return fmodf(a, b);
329	}
330	}
331	}
332	};
333
334	/// convenience typedef for using the default implementation (which is CPP)
335	typedef __RTMath<> RTMath;
336
337	#endif // __RT_MATH_H__