sbrmath.c

/* ***** BEGIN LICENSE BLOCK *****  
 * Source last modified: $Id: sbrmath.c,v 1.1 2005/02/26 01:47:35 jrecker Exp $ 
 *   
 * Portions Copyright (c) 1995-2005 RealNetworks, Inc. All Rights Reserved.  
 *       
 * The contents of this file, and the files included with this file, 
 * are subject to the current version of the RealNetworks Public 
 * Source License (the "RPSL") available at 
 * http://www.helixcommunity.org/content/rpsl unless you have licensed 
 * the file under the current version of the RealNetworks Community 
 * Source License (the "RCSL") available at 
 * http://www.helixcommunity.org/content/rcsl, in which case the RCSL 
 * will apply. You may also obtain the license terms directly from 
 * RealNetworks.  You may not use this file except in compliance with 
 * the RPSL or, if you have a valid RCSL with RealNetworks applicable 
 * to this file, the RCSL.  Please see the applicable RPSL or RCSL for 
 * the rights, obligations and limitations governing use of the 
 * contents of the file. 
 *   
 * This file is part of the Helix DNA Technology. RealNetworks is the 
 * developer of the Original Code and owns the copyrights in the 
 * portions it created. 
 *   
 * This file, and the files included with this file, is distributed 
 * and made available on an 'AS IS' basis, WITHOUT WARRANTY OF ANY 
 * KIND, EITHER EXPRESS OR IMPLIED, AND REALNETWORKS HEREBY DISCLAIMS 
 * ALL SUCH WARRANTIES, INCLUDING WITHOUT LIMITATION, ANY WARRANTIES 
 * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, QUIET 
 * ENJOYMENT OR NON-INFRINGEMENT. 
 *  
 * Technology Compatibility Kit Test Suite(s) Location:  
 *    http://www.helixcommunity.org/content/tck  
 *  
 * Contributor(s):  
 *   
 * ***** END LICENSE BLOCK ***** */  

/**************************************************************************************
 * Fixed-point HE-AAC decoder
 * Jon Recker (jrecker@real.com)
 * February 2005
 *
 * sbrmath.c - fixed-point math functions for SBR
 **************************************************************************************/

#include "sbr.h"
#include "assembly.h"

#define Q28_2	0x20000000	/* Q28: 2.0 */
#define Q28_15	0x30000000	/* Q28: 1.5 */

#define NUM_ITER_IRN		5

/**************************************************************************************
 * Function:    InvRNormalized
 *
 * Description: use Newton's method to solve for x = 1/r
 *
 * Inputs:      r = Q31, range = [0.5, 1) (normalize your inputs to this range)
 *
 * Outputs:     none
 *
 * Return:      x = Q29, range ~= [1.0, 2.0]
 *
 * Notes:       guaranteed to converge and not overflow for any r in [0.5, 1)
 * 
 *              xn+1  = xn - f(xn)/f'(xn)
 *              f(x)  = 1/r - x = 0 (find root)
 *                    = 1/x - r
 *              f'(x) = -1/x^2
 *
 *              so xn+1 = xn - (1/xn - r) / (-1/xn^2)
 *                      = xn * (2 - r*xn)
 *
 *              NUM_ITER_IRN = 2, maxDiff = 6.2500e-02 (precision of about 4 bits)
 *              NUM_ITER_IRN = 3, maxDiff = 3.9063e-03 (precision of about 8 bits)
 *              NUM_ITER_IRN = 4, maxDiff = 1.5288e-05 (precision of about 16 bits)
 *              NUM_ITER_IRN = 5, maxDiff = 3.0034e-08 (precision of about 24 bits)
 **************************************************************************************/
int InvRNormalized(int r)
{
	int i, xn, t;

	/* r =   [0.5, 1.0) 
	 * 1/r = (1.0, 2.0] 
	 *   so use 1.5 as initial guess 
	 */
	xn = Q28_15;

	/* xn = xn*(2.0 - r*xn) */
	for (i = NUM_ITER_IRN; i != 0; i--) {
		t = MULSHIFT32(r, xn);			/* Q31*Q29 = Q28 */
		t = Q28_2 - t;					/* Q28 */
		xn = MULSHIFT32(xn, t) << 4;	/* Q29*Q28 << 4 = Q29 */
	}

	return xn;
}

#define NUM_TERMS_RPI	5
#define LOG2_EXP_INV	0x58b90bfc	/* 1/log2(e), Q31 */

/* invTab[x] = 1/(x+1), format = Q30 */
static const int invTab[NUM_TERMS_RPI] PROGMEM = {0x40000000, 0x20000000, 0x15555555, 0x10000000, 0x0ccccccd};

/**************************************************************************************
 * Function:    RatioPowInv
 *
 * Description: use Taylor (MacLaurin) series expansion to calculate (a/b) ^ (1/c)
 *
 * Inputs:      a = [1, 64], b = [1, 64], c = [1, 64], a >= b
 *
 * Outputs:     none
 *
 * Return:      y = Q24, range ~= [0.015625, 64]
 **************************************************************************************/
int RatioPowInv(int a, int b, int c)
{
	int lna, lnb, i, p, t, y;

	if (a < 1 || b < 1 || c < 1 || a > 64 || b > 64 || c > 64 || a < b)
		return 0;

	lna = MULSHIFT32(log2Tab[a], LOG2_EXP_INV) << 1;	/* ln(a), Q28 */
	lnb = MULSHIFT32(log2Tab[b], LOG2_EXP_INV) << 1;	/* ln(b), Q28 */
	p = (lna - lnb) / c;	/* Q28 */

	/* sum in Q24 */
	y = (1 << 24);
	t = p >> 4;		/* t = p^1 * 1/1! (Q24)*/
	y += t;

	for (i = 2; i <= NUM_TERMS_RPI; i++) {
		t = MULSHIFT32(invTab[i-1], t) << 2;
		t = MULSHIFT32(p, t) << 4;	/* t = p^i * 1/i! (Q24) */
		y += t;
	}

	return y;
}

/**************************************************************************************
 * Function:    SqrtFix
 *
 * Description: use binary search to calculate sqrt(q)
 *
 * Inputs:      q = Q30
 *              number of fraction bits in input
 *
 * Outputs:     number of fraction bits in output
 *
 * Return:      lo = Q(fBitsOut)
 *
 * Notes:       absolute precision varies depending on fBitsIn
 *              normalizes input to range [0x200000000, 0x7fffffff] and takes 
 *                floor(sqrt(input)), and sets fBitsOut appropriately
 **************************************************************************************/
int SqrtFix(int q, int fBitsIn, int *fBitsOut)
{
	int z, lo, hi, mid;

	if (q <= 0) {
		*fBitsOut = fBitsIn;
		return 0;
	}

	/* force even fBitsIn */
	z = fBitsIn & 0x01;
	q >>= z;
	fBitsIn -= z;

	/* for max precision, normalize to [0x20000000, 0x7fffffff] */
	z = (CLZ(q) - 1);
	z >>= 1;
	q <<= (2*z);

	/* choose initial bounds */
	lo = 1;
	if (q >= 0x10000000)
		lo = 16384;	/* (int)sqrt(0x10000000) */
	hi = 46340;		/* (int)sqrt(0x7fffffff) */

	/* do binary search with 32x32->32 multiply test */
	do {
		mid = (lo + hi) >> 1;
		if (mid*mid > q)
			hi = mid - 1;
		else
			lo = mid + 1;
	} while (hi >= lo);
	lo--;

	*fBitsOut = ((fBitsIn + 2*z) >> 1);
	return lo;
}