1 /* Copyright (C) 2005 Jean-Marc Valin */
2 /**
3 @file pseudofloat.h
4 @brief Pseudo-floating point
5 * This header file provides a lightweight floating point type for
6 * use on fixed-point platforms when a large dynamic range is
7 * required. The new type is not compatible with the 32-bit IEEE format,
8 * it is not even remotely as accurate as 32-bit floats, and is not
9 * even guaranteed to produce even remotely correct results for code
10 * other than Speex. It makes all kinds of shortcuts that are acceptable
11 * for Speex, but may not be acceptable for your application. You're
12 * quite welcome to reuse this code and improve it, but don't assume
13 * it works out of the box. Most likely, it doesn't.
14 */
15 /*
16 Redistribution and use in source and binary forms, with or without
17 modification, are permitted provided that the following conditions
18 are met:
19
20 - Redistributions of source code must retain the above copyright
21 notice, this list of conditions and the following disclaimer.
22
23 - Redistributions in binary form must reproduce the above copyright
24 notice, this list of conditions and the following disclaimer in the
25 documentation and/or other materials provided with the distribution.
26
27 - Neither the name of the Xiph.org Foundation nor the names of its
28 contributors may be used to endorse or promote products derived from
29 this software without specific prior written permission.
30
31 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
32 ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
33 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
34 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR
35 CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
36 EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
37 PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
38 PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
39 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
40 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
41 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
42 */
43
44 #ifndef PSEUDOFLOAT_H
45 #define PSEUDOFLOAT_H
46
47 #include "arch.h"
48 #include "os_support.h"
49 #include "math_approx.h"
50 #include <math.h>
51
52 #ifdef FIXED_POINT
53
54 typedef struct {
55 spx_int16_t m;
56 spx_int16_t e;
57 } spx_float_t;
58
59 static const spx_float_t FLOAT_ZERO = {0,0};
60 static const spx_float_t FLOAT_ONE = {16384,-14};
61 static const spx_float_t FLOAT_HALF = {16384,-15};
62
63 #define MIN(a,b) ((a)<(b)?(a):(b))
PSEUDOFLOAT(spx_int32_t x)64 static inline spx_float_t PSEUDOFLOAT(spx_int32_t x)
65 {
66 int e=0;
67 int sign=0;
68 if (x<0)
69 {
70 sign = 1;
71 x = -x;
72 }
73 if (x==0)
74 {
75 spx_float_t r = {0,0};
76 return r;
77 }
78 e = spx_ilog2(ABS32(x))-14;
79 x = VSHR32(x, e);
80 if (sign)
81 {
82 spx_float_t r;
83 r.m = -x;
84 r.e = e;
85 return r;
86 }
87 else
88 {
89 spx_float_t r;
90 r.m = x;
91 r.e = e;
92 return r;
93 }
94 }
95
96
FLOAT_ADD(spx_float_t a,spx_float_t b)97 static inline spx_float_t FLOAT_ADD(spx_float_t a, spx_float_t b)
98 {
99 spx_float_t r;
100 if (a.m==0)
101 return b;
102 else if (b.m==0)
103 return a;
104 if ((a).e > (b).e)
105 {
106 r.m = ((a).m>>1) + ((b).m>>MIN(15,(a).e-(b).e+1));
107 r.e = (a).e+1;
108 }
109 else
110 {
111 r.m = ((b).m>>1) + ((a).m>>MIN(15,(b).e-(a).e+1));
112 r.e = (b).e+1;
113 }
114 if (r.m>0)
115 {
116 if (r.m<16384)
117 {
118 r.m<<=1;
119 r.e-=1;
120 }
121 } else {
122 if (r.m>-16384)
123 {
124 r.m<<=1;
125 r.e-=1;
126 }
127 }
128 /*printf ("%f + %f = %f\n", REALFLOAT(a), REALFLOAT(b), REALFLOAT(r));*/
129 return r;
130 }
131
FLOAT_SUB(spx_float_t a,spx_float_t b)132 static inline spx_float_t FLOAT_SUB(spx_float_t a, spx_float_t b)
133 {
134 spx_float_t r;
135 if (a.m==0)
136 return b;
137 else if (b.m==0)
138 return a;
139 if ((a).e > (b).e)
140 {
141 r.m = ((a).m>>1) - ((b).m>>MIN(15,(a).e-(b).e+1));
142 r.e = (a).e+1;
143 }
144 else
145 {
146 r.m = ((a).m>>MIN(15,(b).e-(a).e+1)) - ((b).m>>1);
147 r.e = (b).e+1;
148 }
149 if (r.m>0)
150 {
151 if (r.m<16384)
152 {
153 r.m<<=1;
154 r.e-=1;
155 }
156 } else {
157 if (r.m>-16384)
158 {
159 r.m<<=1;
160 r.e-=1;
161 }
162 }
163 /*printf ("%f + %f = %f\n", REALFLOAT(a), REALFLOAT(b), REALFLOAT(r));*/
164 return r;
165 }
166
FLOAT_LT(spx_float_t a,spx_float_t b)167 static inline int FLOAT_LT(spx_float_t a, spx_float_t b)
168 {
169 if (a.m==0)
170 return b.m>0;
171 else if (b.m==0)
172 return a.m<0;
173 if ((a).e > (b).e)
174 return ((a).m>>1) < ((b).m>>MIN(15,(a).e-(b).e+1));
175 else
176 return ((b).m>>1) > ((a).m>>MIN(15,(b).e-(a).e+1));
177
178 }
179
FLOAT_GT(spx_float_t a,spx_float_t b)180 static inline int FLOAT_GT(spx_float_t a, spx_float_t b)
181 {
182 return FLOAT_LT(b,a);
183 }
184
FLOAT_MULT(spx_float_t a,spx_float_t b)185 static inline spx_float_t FLOAT_MULT(spx_float_t a, spx_float_t b)
186 {
187 spx_float_t r;
188 r.m = (spx_int16_t)((spx_int32_t)(a).m*(b).m>>15);
189 r.e = (a).e+(b).e+15;
190 if (r.m>0)
191 {
192 if (r.m<16384)
193 {
194 r.m<<=1;
195 r.e-=1;
196 }
197 } else {
198 if (r.m>-16384)
199 {
200 r.m<<=1;
201 r.e-=1;
202 }
203 }
204 /*printf ("%f * %f = %f\n", REALFLOAT(a), REALFLOAT(b), REALFLOAT(r));*/
205 return r;
206 }
207
FLOAT_AMULT(spx_float_t a,spx_float_t b)208 static inline spx_float_t FLOAT_AMULT(spx_float_t a, spx_float_t b)
209 {
210 spx_float_t r;
211 r.m = (spx_int16_t)((spx_int32_t)(a).m*(b).m>>15);
212 r.e = (a).e+(b).e+15;
213 return r;
214 }
215
216
FLOAT_SHL(spx_float_t a,int b)217 static inline spx_float_t FLOAT_SHL(spx_float_t a, int b)
218 {
219 spx_float_t r;
220 r.m = a.m;
221 r.e = a.e+b;
222 return r;
223 }
224
FLOAT_EXTRACT16(spx_float_t a)225 static inline spx_int16_t FLOAT_EXTRACT16(spx_float_t a)
226 {
227 if (a.e<0)
228 return EXTRACT16((EXTEND32(a.m)+(EXTEND32(1)<<(-a.e-1)))>>-a.e);
229 else
230 return a.m<<a.e;
231 }
232
FLOAT_EXTRACT32(spx_float_t a)233 static inline spx_int32_t FLOAT_EXTRACT32(spx_float_t a)
234 {
235 if (a.e<0)
236 return (EXTEND32(a.m)+(EXTEND32(1)<<(-a.e-1)))>>-a.e;
237 else
238 return EXTEND32(a.m)<<a.e;
239 }
240
FLOAT_MUL32(spx_float_t a,spx_word32_t b)241 static inline spx_int32_t FLOAT_MUL32(spx_float_t a, spx_word32_t b)
242 {
243 return VSHR32(MULT16_32_Q15(a.m, b),-a.e-15);
244 }
245
FLOAT_MUL32U(spx_word32_t a,spx_word32_t b)246 static inline spx_float_t FLOAT_MUL32U(spx_word32_t a, spx_word32_t b)
247 {
248 int e1, e2;
249 spx_float_t r;
250 if (a==0 || b==0)
251 {
252 return FLOAT_ZERO;
253 }
254 e1 = spx_ilog2(ABS32(a));
255 a = VSHR32(a, e1-14);
256 e2 = spx_ilog2(ABS32(b));
257 b = VSHR32(b, e2-14);
258 r.m = MULT16_16_Q15(a,b);
259 r.e = e1+e2-13;
260 return r;
261 }
262
263 /* Do NOT attempt to divide by a negative number */
FLOAT_DIV32_FLOAT(spx_word32_t a,spx_float_t b)264 static inline spx_float_t FLOAT_DIV32_FLOAT(spx_word32_t a, spx_float_t b)
265 {
266 int e=0;
267 spx_float_t r;
268 if (a==0)
269 {
270 return FLOAT_ZERO;
271 }
272 e = spx_ilog2(ABS32(a))-spx_ilog2(b.m-1)-15;
273 a = VSHR32(a, e);
274 if (ABS32(a)>=SHL32(EXTEND32(b.m-1),15))
275 {
276 a >>= 1;
277 e++;
278 }
279 r.m = DIV32_16(a,b.m);
280 r.e = e-b.e;
281 return r;
282 }
283
284
285 /* Do NOT attempt to divide by a negative number */
FLOAT_DIV32(spx_word32_t a,spx_word32_t b)286 static inline spx_float_t FLOAT_DIV32(spx_word32_t a, spx_word32_t b)
287 {
288 int e0=0,e=0;
289 spx_float_t r;
290 if (a==0)
291 {
292 return FLOAT_ZERO;
293 }
294 if (b>32767)
295 {
296 e0 = spx_ilog2(b)-14;
297 b = VSHR32(b, e0);
298 e0 = -e0;
299 }
300 e = spx_ilog2(ABS32(a))-spx_ilog2(b-1)-15;
301 a = VSHR32(a, e);
302 if (ABS32(a)>=SHL32(EXTEND32(b-1),15))
303 {
304 a >>= 1;
305 e++;
306 }
307 e += e0;
308 r.m = DIV32_16(a,b);
309 r.e = e;
310 return r;
311 }
312
313 /* Do NOT attempt to divide by a negative number */
FLOAT_DIVU(spx_float_t a,spx_float_t b)314 static inline spx_float_t FLOAT_DIVU(spx_float_t a, spx_float_t b)
315 {
316 int e=0;
317 spx_int32_t num;
318 spx_float_t r;
319 if (b.m<=0)
320 {
321 speex_warning_int("Attempted to divide by", b.m);
322 return FLOAT_ONE;
323 }
324 num = a.m;
325 a.m = ABS16(a.m);
326 while (a.m >= b.m)
327 {
328 e++;
329 a.m >>= 1;
330 }
331 num = num << (15-e);
332 r.m = DIV32_16(num,b.m);
333 r.e = a.e-b.e-15+e;
334 return r;
335 }
336
FLOAT_SQRT(spx_float_t a)337 static inline spx_float_t FLOAT_SQRT(spx_float_t a)
338 {
339 spx_float_t r;
340 spx_int32_t m;
341 m = SHL32(EXTEND32(a.m), 14);
342 r.e = a.e - 14;
343 if (r.e & 1)
344 {
345 r.e -= 1;
346 m <<= 1;
347 }
348 r.e >>= 1;
349 r.m = spx_sqrt(m);
350 return r;
351 }
352
353 #else
354
355 #define spx_float_t float
356 #define FLOAT_ZERO 0.f
357 #define FLOAT_ONE 1.f
358 #define FLOAT_HALF 0.5f
359 #define PSEUDOFLOAT(x) (x)
360 #define FLOAT_MULT(a,b) ((a)*(b))
361 #define FLOAT_AMULT(a,b) ((a)*(b))
362 #define FLOAT_MUL32(a,b) ((a)*(b))
363 #define FLOAT_DIV32(a,b) ((a)/(b))
364 #define FLOAT_EXTRACT16(a) (a)
365 #define FLOAT_EXTRACT32(a) (a)
366 #define FLOAT_ADD(a,b) ((a)+(b))
367 #define FLOAT_SUB(a,b) ((a)-(b))
368 #define REALFLOAT(x) (x)
369 #define FLOAT_DIV32_FLOAT(a,b) ((a)/(b))
370 #define FLOAT_MUL32U(a,b) ((a)*(b))
371 #define FLOAT_SHL(a,b) (a)
372 #define FLOAT_LT(a,b) ((a)<(b))
373 #define FLOAT_GT(a,b) ((a)>(b))
374 #define FLOAT_DIVU(a,b) ((a)/(b))
375 #define FLOAT_SQRT(a) (spx_sqrt(a))
376
377 #endif
378
379 #endif
380