blob: 784f90907af100428d438e402293c73ae388b955
| 1 | /* ***** BEGIN LICENSE BLOCK ***** |
| 2 | * Source last modified: $Id: sbrmath.c,v 1.1 2005/02/26 01:47:35 jrecker Exp $ |
| 3 | * |
| 4 | * Portions Copyright (c) 1995-2005 RealNetworks, Inc. All Rights Reserved. |
| 5 | * |
| 6 | * The contents of this file, and the files included with this file, |
| 7 | * are subject to the current version of the RealNetworks Public |
| 8 | * Source License (the "RPSL") available at |
| 9 | * http://www.helixcommunity.org/content/rpsl unless you have licensed |
| 10 | * the file under the current version of the RealNetworks Community |
| 11 | * Source License (the "RCSL") available at |
| 12 | * http://www.helixcommunity.org/content/rcsl, in which case the RCSL |
| 13 | * will apply. You may also obtain the license terms directly from |
| 14 | * RealNetworks. You may not use this file except in compliance with |
| 15 | * the RPSL or, if you have a valid RCSL with RealNetworks applicable |
| 16 | * to this file, the RCSL. Please see the applicable RPSL or RCSL for |
| 17 | * the rights, obligations and limitations governing use of the |
| 18 | * contents of the file. |
| 19 | * |
| 20 | * This file is part of the Helix DNA Technology. RealNetworks is the |
| 21 | * developer of the Original Code and owns the copyrights in the |
| 22 | * portions it created. |
| 23 | * |
| 24 | * This file, and the files included with this file, is distributed |
| 25 | * and made available on an 'AS IS' basis, WITHOUT WARRANTY OF ANY |
| 26 | * KIND, EITHER EXPRESS OR IMPLIED, AND REALNETWORKS HEREBY DISCLAIMS |
| 27 | * ALL SUCH WARRANTIES, INCLUDING WITHOUT LIMITATION, ANY WARRANTIES |
| 28 | * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, QUIET |
| 29 | * ENJOYMENT OR NON-INFRINGEMENT. |
| 30 | * |
| 31 | * Technology Compatibility Kit Test Suite(s) Location: |
| 32 | * http://www.helixcommunity.org/content/tck |
| 33 | * |
| 34 | * Contributor(s): |
| 35 | * |
| 36 | * ***** END LICENSE BLOCK ***** */ |
| 37 | |
| 38 | /************************************************************************************** |
| 39 | * Fixed-point HE-AAC decoder |
| 40 | * Jon Recker (jrecker@real.com) |
| 41 | * February 2005 |
| 42 | * |
| 43 | * sbrmath.c - fixed-point math functions for SBR |
| 44 | **************************************************************************************/ |
| 45 | |
| 46 | #include "sbr.h" |
| 47 | #include "assembly.h" |
| 48 | |
| 49 | #define Q28_2 0x20000000 /* Q28: 2.0 */ |
| 50 | #define Q28_15 0x30000000 /* Q28: 1.5 */ |
| 51 | |
| 52 | #define NUM_ITER_IRN 5 |
| 53 | |
| 54 | /************************************************************************************** |
| 55 | * Function: InvRNormalized |
| 56 | * |
| 57 | * Description: use Newton's method to solve for x = 1/r |
| 58 | * |
| 59 | * Inputs: r = Q31, range = [0.5, 1) (normalize your inputs to this range) |
| 60 | * |
| 61 | * Outputs: none |
| 62 | * |
| 63 | * Return: x = Q29, range ~= [1.0, 2.0] |
| 64 | * |
| 65 | * Notes: guaranteed to converge and not overflow for any r in [0.5, 1) |
| 66 | * |
| 67 | * xn+1 = xn - f(xn)/f'(xn) |
| 68 | * f(x) = 1/r - x = 0 (find root) |
| 69 | * = 1/x - r |
| 70 | * f'(x) = -1/x^2 |
| 71 | * |
| 72 | * so xn+1 = xn - (1/xn - r) / (-1/xn^2) |
| 73 | * = xn * (2 - r*xn) |
| 74 | * |
| 75 | * NUM_ITER_IRN = 2, maxDiff = 6.2500e-02 (precision of about 4 bits) |
| 76 | * NUM_ITER_IRN = 3, maxDiff = 3.9063e-03 (precision of about 8 bits) |
| 77 | * NUM_ITER_IRN = 4, maxDiff = 1.5288e-05 (precision of about 16 bits) |
| 78 | * NUM_ITER_IRN = 5, maxDiff = 3.0034e-08 (precision of about 24 bits) |
| 79 | **************************************************************************************/ |
| 80 | int InvRNormalized(int r) |
| 81 | { |
| 82 | int i, xn, t; |
| 83 | |
| 84 | /* r = [0.5, 1.0) |
| 85 | * 1/r = (1.0, 2.0] |
| 86 | * so use 1.5 as initial guess |
| 87 | */ |
| 88 | xn = Q28_15; |
| 89 | |
| 90 | /* xn = xn*(2.0 - r*xn) */ |
| 91 | for (i = NUM_ITER_IRN; i != 0; i--) { |
| 92 | t = MULSHIFT32(r, xn); /* Q31*Q29 = Q28 */ |
| 93 | t = Q28_2 - t; /* Q28 */ |
| 94 | xn = MULSHIFT32(xn, t) << 4; /* Q29*Q28 << 4 = Q29 */ |
| 95 | } |
| 96 | |
| 97 | return xn; |
| 98 | } |
| 99 | |
| 100 | #define NUM_TERMS_RPI 5 |
| 101 | #define LOG2_EXP_INV 0x58b90bfc /* 1/log2(e), Q31 */ |
| 102 | |
| 103 | /* invTab[x] = 1/(x+1), format = Q30 */ |
| 104 | static const int invTab[NUM_TERMS_RPI] = {0x40000000, 0x20000000, 0x15555555, 0x10000000, 0x0ccccccd}; |
| 105 | |
| 106 | /************************************************************************************** |
| 107 | * Function: RatioPowInv |
| 108 | * |
| 109 | * Description: use Taylor (MacLaurin) series expansion to calculate (a/b) ^ (1/c) |
| 110 | * |
| 111 | * Inputs: a = [1, 64], b = [1, 64], c = [1, 64], a >= b |
| 112 | * |
| 113 | * Outputs: none |
| 114 | * |
| 115 | * Return: y = Q24, range ~= [0.015625, 64] |
| 116 | **************************************************************************************/ |
| 117 | int RatioPowInv(int a, int b, int c) |
| 118 | { |
| 119 | int lna, lnb, i, p, t, y; |
| 120 | |
| 121 | if (a < 1 || b < 1 || c < 1 || a > 64 || b > 64 || c > 64 || a < b) { |
| 122 | return 0; |
| 123 | } |
| 124 | |
| 125 | lna = MULSHIFT32(log2Tab[a], LOG2_EXP_INV) << 1; /* ln(a), Q28 */ |
| 126 | lnb = MULSHIFT32(log2Tab[b], LOG2_EXP_INV) << 1; /* ln(b), Q28 */ |
| 127 | p = (lna - lnb) / c; /* Q28 */ |
| 128 | |
| 129 | /* sum in Q24 */ |
| 130 | y = (1 << 24); |
| 131 | t = p >> 4; /* t = p^1 * 1/1! (Q24)*/ |
| 132 | y += t; |
| 133 | |
| 134 | for (i = 2; i <= NUM_TERMS_RPI; i++) { |
| 135 | t = MULSHIFT32(invTab[i - 1], t) << 2; |
| 136 | t = MULSHIFT32(p, t) << 4; /* t = p^i * 1/i! (Q24) */ |
| 137 | y += t; |
| 138 | } |
| 139 | |
| 140 | return y; |
| 141 | } |
| 142 | |
| 143 | /************************************************************************************** |
| 144 | * Function: SqrtFix |
| 145 | * |
| 146 | * Description: use binary search to calculate sqrt(q) |
| 147 | * |
| 148 | * Inputs: q = Q30 |
| 149 | * number of fraction bits in input |
| 150 | * |
| 151 | * Outputs: number of fraction bits in output |
| 152 | * |
| 153 | * Return: lo = Q(fBitsOut) |
| 154 | * |
| 155 | * Notes: absolute precision varies depending on fBitsIn |
| 156 | * normalizes input to range [0x200000000, 0x7fffffff] and takes |
| 157 | * floor(sqrt(input)), and sets fBitsOut appropriately |
| 158 | **************************************************************************************/ |
| 159 | int SqrtFix(int q, int fBitsIn, int *fBitsOut) |
| 160 | { |
| 161 | int z, lo, hi, mid; |
| 162 | |
| 163 | if (q <= 0) { |
| 164 | *fBitsOut = fBitsIn; |
| 165 | return 0; |
| 166 | } |
| 167 | |
| 168 | /* force even fBitsIn */ |
| 169 | z = fBitsIn & 0x01; |
| 170 | q >>= z; |
| 171 | fBitsIn -= z; |
| 172 | |
| 173 | /* for max precision, normalize to [0x20000000, 0x7fffffff] */ |
| 174 | z = (CLZ(q) - 1); |
| 175 | z >>= 1; |
| 176 | q <<= (2 * z); |
| 177 | |
| 178 | /* choose initial bounds */ |
| 179 | lo = 1; |
| 180 | if (q >= 0x10000000) { |
| 181 | lo = 16384; /* (int)sqrt(0x10000000) */ |
| 182 | } |
| 183 | hi = 46340; /* (int)sqrt(0x7fffffff) */ |
| 184 | |
| 185 | /* do binary search with 32x32->32 multiply test */ |
| 186 | do { |
| 187 | mid = (lo + hi) >> 1; |
| 188 | if (mid * mid > q) { |
| 189 | hi = mid - 1; |
| 190 | } else { |
| 191 | lo = mid + 1; |
| 192 | } |
| 193 | } while (hi >= lo); |
| 194 | lo--; |
| 195 | |
| 196 | *fBitsOut = ((fBitsIn + 2 * z) >> 1); |
| 197 | return lo; |
| 198 | } |
| 199 |