1 /*
2 * Musepack audio compression
3 *
4 * This library is free software; you can redistribute it and/or
5 * modify it under the terms of the GNU Lesser General Public
6 * License as published by the Free Software Foundation; either
7 * version 2.1 of the License, or (at your option) any later version.
8 *
9 * This library is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 * Lesser General Public License for more details.
13 *
14 * You should have received a copy of the GNU Lesser General Public
15 * License along with this library; If not, see <http://www.gnu.org/licenses/>
16 */
17
18 #include <math.h>
19
20 #include "mpc_types.h"
21
22 typedef union mpc_floatint
23 {
24 float f;
25 mpc_int32_t n;
26 } mpc_floatint;
27
28 typedef union mpc_doubleint
29 {
30 double d;
31 mpc_int32_t n[2];
32 } mpc_doubleint;
33
mpc_lrintf(float fVal)34 static mpc_inline mpc_int32_t mpc_lrintf(float fVal)
35 {
36 mpc_floatint tmp;
37 tmp.f = fVal + 0x00FF8000;
38 return tmp.n - 0x4B7F8000;
39 }
40
41 #define mpc_round32 mpc_lrintf
42 #define mpc_nearbyintf mpc_lrintf
43
44
45 #ifndef M_PI
46 # define M_PI 3.1415926535897932384626433832795029 // 4*atan(1)
47 # define M_PIl 3.1415926535897932384626433832795029L
48 # define M_LN2 0.6931471805599453094172321214581766 // ln(2)
49 # define M_LN2l 0.6931471805599453094172321214581766L
50 # define M_LN10 2.3025850929940456840179914546843642 // ln 10 */
51 # define M_LN10l 2.3025850929940456840179914546843642L
52 #endif
53
54 // fast but maybe more inaccurate, use if you need speed
55 #if defined(__GNUC__) && !defined(__APPLE__)
56 # define SIN(x) sinf ((float)(x))
57 # define COS(x) cosf ((float)(x))
58 # define ATAN2(x,y) atan2f ((float)(x), (float)(y))
59 # define SQRT(x) sqrtf ((float)(x))
60 # define LOG(x) logf ((float)(x))
61 # define LOG10(x) log10f ((float)(x))
62 # define POW(x,y) expf (logf(x) * (y))
63 # define POW10(x) expf (M_LN10 * (x))
64 # define FLOOR(x) floorf ((float)(x))
65 # define IFLOOR(x) (int) floorf ((float)(x))
66 # define FABS(x) fabsf ((float)(x))
67 #else
68 # define SIN(x) (float) sin (x)
69 # define COS(x) (float) cos (x)
70 # define ATAN2(x,y) (float) atan2 (x, y)
71 # define SQRT(x) (float) sqrt (x)
72 # define LOG(x) (float) log (x)
73 # define LOG10(x) (float) log10 (x)
74 # define POW(x,y) (float) pow (x,y)
75 # define POW10(x) (float) pow (10., (x))
76 # define FLOOR(x) (float) floor (x)
77 # define IFLOOR(x) (int) floor (x)
78 # define FABS(x) (float) fabs (x)
79 #endif
80
81 #define SQRTF(x) SQRT (x)
82 #ifdef FAST_MATH
83 # define TABSTEP 64
84 # define COSF(x) my_cos ((float)(x))
85 # define ATAN2F(x,y) my_atan2 ((float)(x), (float)(y))
86 # define IFLOORF(x) my_ifloor ((float)(x))
87
88 void Init_FastMath ( void );
89 extern const float tabatan2 [] [2];
90 extern const float tabcos [] [2];
91 extern const float tabsqrt_ex [];
92 extern const float tabsqrt_m [] [2];
93
my_atan2(float x,float y)94 static mpc_inline float my_atan2 ( float x, float y )
95 {
96 float t, ret; int i; mpc_floatint mx, my;
97
98 mx.f = x;
99 my.f = y;
100 if ( (mx.n & 0x7FFFFFFF) < (my.n & 0x7FFFFFFF) ) {
101 i = mpc_round32 (t = TABSTEP * (mx.f / my.f));
102 ret = tabatan2 [1*TABSTEP+i][0] + tabatan2 [1*TABSTEP+i][1] * (t-i);
103 if ( my.n < 0 )
104 ret = (float)(ret - M_PI);
105 }
106 else if ( mx.n < 0 ) {
107 i = mpc_round32 (t = TABSTEP * (my.f / mx.f));
108 ret = - M_PI/2 - tabatan2 [1*TABSTEP+i][0] + tabatan2 [1*TABSTEP+i][1] * (i-t);
109 }
110 else if ( mx.n > 0 ) {
111 i = mpc_round32 (t = TABSTEP * (my.f / mx.f));
112 ret = + M_PI/2 - tabatan2 [1*TABSTEP+i][0] + tabatan2 [1*TABSTEP+i][1] * (i-t);
113 }
114 else {
115 ret = 0.;
116 }
117 return ret;
118 }
119
120
my_cos(float x)121 static mpc_inline float my_cos ( float x )
122 {
123 float t, ret; int i;
124 i = mpc_round32 (t = TABSTEP * x);
125 ret = tabcos [13*TABSTEP+i][0] + tabcos [13*TABSTEP+i][1] * (t-i);
126 return ret;
127 }
128
129
my_ifloor(float x)130 static mpc_inline int my_ifloor ( float x )
131 {
132 mpc_floatint mx;
133 mx.f = (float) (x + (0x0C00000L + 0.500000001));
134 return mx.n - 1262485505;
135 }
136
137
my_sqrt(float x)138 static mpc_inline float my_sqrt ( float x )
139 {
140 float ret; int i, ex; mpc_floatint mx;
141 mx.f = x;
142 ex = mx.n >> 23; // get the exponent
143 mx.n = (mx.n & 0x7FFFFF) | 0x42800000; // delete the exponent
144 i = mpc_round32 (mx.f); // Integer-part of the mantissa (round ????????????)
145 ret = tabsqrt_m [i-TABSTEP][0] + tabsqrt_m [i-TABSTEP][1] * (mx.f-i); // calculate value
146 ret *= tabsqrt_ex [ex];
147 return ret;
148 }
149 #else
150 # define COSF(x) COS (x)
151 # define ATAN2F(x,y) ATAN2 (x,y)
152 # define IFLOORF(x) IFLOOR (x)
153 #endif
154
155