summaryrefslogtreecommitdiff
path: root/apps/fixedpoint.c
diff options
context:
space:
mode:
authorMaurus Cuelenaere <mcuelenaere@gmail.com>2009-07-07 13:36:34 +0000
committerMaurus Cuelenaere <mcuelenaere@gmail.com>2009-07-07 13:36:34 +0000
commit8d4d4610b6814c0c4a4abe9523b4b6b2b80ac9a4 (patch)
tree5823f6ea576834eb0b7fd9f1fc78a85997d0e224 /apps/fixedpoint.c
parent616905f965ed626c398624ec2ffa988f69899c32 (diff)
downloadrockbox-8d4d4610b6814c0c4a4abe9523b4b6b2b80ac9a4.tar.gz
rockbox-8d4d4610b6814c0c4a4abe9523b4b6b2b80ac9a4.zip
* FS#10411 - Fixed point math code is bloated by Jeffrey Goode
* Set svn:keywords properties git-svn-id: svn://svn.rockbox.org/rockbox/trunk@21701 a1c6a512-1295-4272-9138-f99709370657
Diffstat (limited to 'apps/fixedpoint.c')
-rw-r--r--apps/fixedpoint.c131
1 files changed, 41 insertions, 90 deletions
diff --git a/apps/fixedpoint.c b/apps/fixedpoint.c
index 917f624258..f9903f301f 100644
--- a/apps/fixedpoint.c
+++ b/apps/fixedpoint.c
@@ -5,7 +5,7 @@
5 * Jukebox | | ( <_> ) \___| < | \_\ ( <_> > < < 5 * Jukebox | | ( <_> ) \___| < | \_\ ( <_> > < <
6 * Firmware |____|_ /\____/ \___ >__|_ \|___ /\____/__/\_ \ 6 * Firmware |____|_ /\____/ \___ >__|_ \|___ /\____/__/\_ \
7 * \/ \/ \/ \/ \/ 7 * \/ \/ \/ \/ \/
8 * $Id: fixedpoint.c -1 $ 8 * $Id$
9 * 9 *
10 * Copyright (C) 2006 Jens Arnold 10 * Copyright (C) 2006 Jens Arnold
11 * 11 *
@@ -261,24 +261,18 @@ long fp16_log(int x) {
261 261
262#if (!defined(PLUGIN) && !defined(CODEC)) 262#if (!defined(PLUGIN) && !defined(CODEC))
263/** MODIFIED FROM replaygain.c */ 263/** MODIFIED FROM replaygain.c */
264/* These math routines have 64-bit internal precision to avoid overflows.
265 * Arguments and return values are 32-bit (long) precision.
266 */
267 264
268#define FP_MUL64(x, y) (((x) * (y)) >> (fracbits)) 265#define FP_MUL_FRAC(x, y) fp_mul(x, y, fracbits)
269#define FP_DIV64(x, y) (((x) << (fracbits)) / (y)) 266#define FP_DIV_FRAC(x, y) fp_div(x, y, fracbits)
270
271static long long fp_exp10(long long x, unsigned int fracbits);
272/* static long long fp_log10(long long n, unsigned int fracbits); */
273 267
274/* constants in fixed point format, 28 fractional bits */ 268/* constants in fixed point format, 28 fractional bits */
275#define FP28_LN2 (186065279LL) /* ln(2) */ 269#define FP28_LN2 (186065279L) /* ln(2) */
276#define FP28_LN2_INV (387270501LL) /* 1/ln(2) */ 270#define FP28_LN2_INV (387270501L) /* 1/ln(2) */
277#define FP28_EXP_ZERO (44739243LL) /* 1/6 */ 271#define FP28_EXP_ZERO (44739243L) /* 1/6 */
278#define FP28_EXP_ONE (-745654LL) /* -1/360 */ 272#define FP28_EXP_ONE (-745654L) /* -1/360 */
279#define FP28_EXP_TWO (12428LL) /* 1/21600 */ 273#define FP28_EXP_TWO (12428L) /* 1/21600 */
280#define FP28_LN10 (618095479LL) /* ln(10) */ 274#define FP28_LN10 (618095479L) /* ln(10) */
281#define FP28_LOG10OF2 (80807124LL) /* log10(2) */ 275#define FP28_LOG10OF2 (80807124L) /* log10(2) */
282 276
283#define TOL_BITS 2 /* log calculation tolerance */ 277#define TOL_BITS 2 /* log calculation tolerance */
284 278
@@ -290,24 +284,24 @@ static long long fp_exp10(long long x, unsigned int fracbits);
290/** FIXED POINT EXP10 284/** FIXED POINT EXP10
291 * Return 10^x as FP integer. Argument is FP integer. 285 * Return 10^x as FP integer. Argument is FP integer.
292 */ 286 */
293static long long fp_exp10(long long x, unsigned int fracbits) 287static long fp_exp10(long x, unsigned int fracbits)
294{ 288{
295 long long k; 289 long k;
296 long long z; 290 long z;
297 long long R; 291 long R;
298 long long xp; 292 long xp;
299 293
300 /* scale constants */ 294 /* scale constants */
301 const long long fp_one = (1 << fracbits); 295 const long fp_one = (1 << fracbits);
302 const long long fp_half = (1 << (fracbits - 1)); 296 const long fp_half = (1 << (fracbits - 1));
303 const long long fp_two = (2 << fracbits); 297 const long fp_two = (2 << fracbits);
304 const long long fp_mask = (fp_one - 1); 298 const long fp_mask = (fp_one - 1);
305 const long long fp_ln2_inv = (FP28_LN2_INV >> (28 - fracbits)); 299 const long fp_ln2_inv = (FP28_LN2_INV >> (28 - fracbits));
306 const long long fp_ln2 = (FP28_LN2 >> (28 - fracbits)); 300 const long fp_ln2 = (FP28_LN2 >> (28 - fracbits));
307 const long long fp_ln10 = (FP28_LN10 >> (28 - fracbits)); 301 const long fp_ln10 = (FP28_LN10 >> (28 - fracbits));
308 const long long fp_exp_zero = (FP28_EXP_ZERO >> (28 - fracbits)); 302 const long fp_exp_zero = (FP28_EXP_ZERO >> (28 - fracbits));
309 const long long fp_exp_one = (FP28_EXP_ONE >> (28 - fracbits)); 303 const long fp_exp_one = (FP28_EXP_ONE >> (28 - fracbits));
310 const long long fp_exp_two = (FP28_EXP_TWO >> (28 - fracbits)); 304 const long fp_exp_two = (FP28_EXP_TWO >> (28 - fracbits));
311 305
312 /* exp(0) = 1 */ 306 /* exp(0) = 1 */
313 if (x == 0) 307 if (x == 0)
@@ -316,21 +310,21 @@ static long long fp_exp10(long long x, unsigned int fracbits)
316 } 310 }
317 311
318 /* convert from base 10 to base e */ 312 /* convert from base 10 to base e */
319 x = FP_MUL64(x, fp_ln10); 313 x = FP_MUL_FRAC(x, fp_ln10);
320 314
321 /* calculate exp(x) */ 315 /* calculate exp(x) */
322 k = (FP_MUL64(abs(x), fp_ln2_inv) + fp_half) & ~fp_mask; 316 k = (FP_MUL_FRAC(abs(x), fp_ln2_inv) + fp_half) & ~fp_mask;
323 317
324 if (x < 0) 318 if (x < 0)
325 { 319 {
326 k = -k; 320 k = -k;
327 } 321 }
328 322
329 x -= FP_MUL64(k, fp_ln2); 323 x -= FP_MUL_FRAC(k, fp_ln2);
330 z = FP_MUL64(x, x); 324 z = FP_MUL_FRAC(x, x);
331 R = fp_two + FP_MUL64(z, fp_exp_zero + FP_MUL64(z, fp_exp_one 325 R = fp_two + FP_MUL_FRAC(z, fp_exp_zero + FP_MUL_FRAC(z, fp_exp_one
332 + FP_MUL64(z, fp_exp_two))); 326 + FP_MUL_FRAC(z, fp_exp_two)));
333 xp = fp_one + FP_DIV64(FP_MUL64(fp_two, x), R - x); 327 xp = fp_one + FP_DIV_FRAC(FP_MUL_FRAC(fp_two, x), R - x);
334 328
335 if (k < 0) 329 if (k < 0)
336 { 330 {
@@ -341,7 +335,7 @@ static long long fp_exp10(long long x, unsigned int fracbits)
341 k = fp_one << (k >> fracbits); 335 k = fp_one << (k >> fracbits);
342 } 336 }
343 337
344 return FP_MUL64(k, xp); 338 return FP_MUL_FRAC(k, xp);
345} 339}
346 340
347 341
@@ -349,13 +343,13 @@ static long long fp_exp10(long long x, unsigned int fracbits)
349/** FIXED POINT LOG10 343/** FIXED POINT LOG10
350 * Return log10(x) as FP integer. Argument is FP integer. 344 * Return log10(x) as FP integer. Argument is FP integer.
351 */ 345 */
352static long long fp_log10(long long n, unsigned int fracbits) 346static long fp_log10(long n, unsigned int fracbits)
353{ 347{
354 /* Calculate log2 of argument */ 348 /* Calculate log2 of argument */
355 349
356 long long log2, frac; 350 long log2, frac;
357 const long long fp_one = (1 << fracbits); 351 const long fp_one = (1 << fracbits);
358 const long long fp_two = (2 << fracbits); 352 const long fp_two = (2 << fracbits);
359 const long tolerance = (1 << ((fracbits / 2) + 2)); 353 const long tolerance = (1 << ((fracbits / 2) + 2));
360 354
361 if (n <=0) return FP_NEGINF; 355 if (n <=0) return FP_NEGINF;
@@ -378,7 +372,7 @@ static long long fp_log10(long long n, unsigned int fracbits)
378 while (frac > tolerance) 372 while (frac > tolerance)
379 { 373 {
380 frac >>= 1; 374 frac >>= 1;
381 n = FP_MUL64(n, n); 375 n = FP_MUL_FRAC(n, n);
382 if (n >= fp_two) 376 if (n >= fp_two)
383 { 377 {
384 n >>= 1; 378 n >>= 1;
@@ -387,31 +381,15 @@ static long long fp_log10(long long n, unsigned int fracbits)
387 } 381 }
388 382
389 /* convert log2 to log10 */ 383 /* convert log2 to log10 */
390 return FP_MUL64(log2, (FP28_LOG10OF2 >> (28 - fracbits))); 384 return FP_MUL_FRAC(log2, (FP28_LOG10OF2 >> (28 - fracbits)));
391} 385}
392 386
393 387
394/** CONVERT FACTOR TO DECIBELS */ 388/** CONVERT FACTOR TO DECIBELS */
395long fp_decibels(unsigned long factor, unsigned int fracbits) 389long fp_decibels(unsigned long factor, unsigned int fracbits)
396{ 390{
397 long long decibels;
398 long long f = (long long)factor;
399 bool neg;
400
401 /* keep factor in signed long range */
402 if (f >= (1LL << 31))
403 f = (1LL << 31) - 1;
404
405 /* decibels = 20 * log10(factor) */ 391 /* decibels = 20 * log10(factor) */
406 decibels = FP_MUL64((20LL << fracbits), fp_log10(f, fracbits)); 392 return FP_MUL_FRAC((20L << fracbits), fp_log10(factor, fracbits));
407
408 /* keep result in signed long range */
409 if ((neg = (decibels < 0)))
410 decibels = -decibels;
411 if (decibels >= (1LL << 31))
412 return neg ? FP_NEGINF : FP_INF;
413
414 return neg ? (long)-decibels : (long)decibels;
415} 393}
416#endif /* unused code */ 394#endif /* unused code */
417 395
@@ -419,34 +397,7 @@ long fp_decibels(unsigned long factor, unsigned int fracbits)
419/** CONVERT DECIBELS TO FACTOR */ 397/** CONVERT DECIBELS TO FACTOR */
420long fp_factor(long decibels, unsigned int fracbits) 398long fp_factor(long decibels, unsigned int fracbits)
421{ 399{
422 bool neg;
423 long long factor;
424 long long db = (long long)decibels;
425
426 /* if decibels is 0, factor is 1 */
427 if (db == 0)
428 return (1L << fracbits);
429
430 /* calculate for positive decibels only */
431 if ((neg = (db < 0)))
432 db = -db;
433
434 /* factor = 10 ^ (decibels / 20) */ 400 /* factor = 10 ^ (decibels / 20) */
435 factor = fp_exp10(FP_DIV64(db, (20LL << fracbits)), fracbits); 401 return fp_exp10(FP_DIV_FRAC(decibels, (20L << fracbits)), fracbits);
436
437 /* keep result in signed long range, return 0 if very small */
438 if (factor >= (1LL << 31))
439 {
440 if (neg)
441 return 0;
442 else
443 return FP_INF;
444 }
445
446 /* if negative argument, factor is 1 / result */
447 if (neg)
448 factor = FP_DIV64((1LL << fracbits), factor);
449
450 return (long)factor;
451} 402}
452#endif /* !PLUGIN and !CODEC */ 403#endif /* !PLUGIN and !CODEC */