From 6cfabd35363c3ef5e3b209b867169a500b3ccc3c Mon Sep 17 00:00:00 2001 From: Fraunhofer IIS FDK Date: Mon, 26 Feb 2018 20:17:00 +0100 Subject: Upgrade to FDKv2 Bug: 71430241 Test: CTS DecoderTest and DecoderTestAacDrc original-Change-Id: Iaa20f749b8a04d553b20247cfe1a8930ebbabe30 Apply clang-format also on header files. original-Change-Id: I14de1ef16bbc79ec0283e745f98356a10efeb2e4 Fixes for MPEG-D DRC original-Change-Id: If1de2d74bbbac84b3f67de3b88b83f6a23b8a15c Catch unsupported tw_mdct at an early stage original-Change-Id: Ied9dd00d754162a0e3ca1ae3e6b854315d818afe Fixing PVC transition frames original-Change-Id: Ib75725abe39252806c32d71176308f2c03547a4e Move qmf bands sanity check original-Change-Id: Iab540c3013c174d9490d2ae100a4576f51d8dbc4 Initialize scaling variable original-Change-Id: I3c4087101b70e998c71c1689b122b0d7762e0f9e Add 16 qmf band configuration to getSlotNrgHQ() original-Change-Id: I49a5d30f703a1b126ff163df9656db2540df21f1 Always apply byte alignment at the end of the AudioMuxElement original-Change-Id: I42d560287506d65d4c3de8bfe3eb9a4ebeb4efc7 Setup SBR element only if no parse error exists original-Change-Id: I1915b73704bc80ab882b9173d6bec59cbd073676 Additional array index check in HCR original-Change-Id: I18cc6e501ea683b5009f1bbee26de8ddd04d8267 Fix fade-in index selection in concealment module original-Change-Id: Ibf802ed6ed8c05e9257e1f3b6d0ac1162e9b81c1 Enable explicit backward compatible parser for AAC_LD original-Change-Id: I27e9c678dcb5d40ed760a6d1e06609563d02482d Skip spatial specific config in explicit backward compatible ASC original-Change-Id: Iff7cc365561319e886090cedf30533f562ea4d6e Update flags description in decoder API original-Change-Id: I9a5b4f8da76bb652f5580cbd3ba9760425c43830 Add QMF domain reset function original-Change-Id: I4f89a8a2c0277d18103380134e4ed86996e9d8d6 DRC upgrade v2.1.0 original-Change-Id: I5731c0540139dab220094cd978ef42099fc45b74 Fix integer overflow in sqrtFixp_lookup() original-Change-Id: I429a6f0d19aa2cc957e0f181066f0ca73968c914 Fix integer overflow in invSqrtNorm2() original-Change-Id: I84de5cbf9fb3adeb611db203fe492fabf4eb6155 Fix integer overflow in GenerateRandomVector() original-Change-Id: I3118a641008bd9484d479e5b0b1ee2b5d7d44d74 Fix integer overflow in adjustTimeSlot_EldGrid() original-Change-Id: I29d503c247c5c8282349b79df940416a512fb9d5 Fix integer overflow in FDKsbrEnc_codeEnvelope() original-Change-Id: I6b34b61ebb9d525b0c651ed08de2befc1f801449 Follow-up on: Fix integer overflow in adjustTimeSlot_EldGrid() original-Change-Id: I6f8f578cc7089e5eb7c7b93e580b72ca35ad689a Fix integer overflow in get_pk_v2() original-Change-Id: I63375bed40d45867f6eeaa72b20b1f33e815938c Fix integer overflow in Syn_filt_zero() original-Change-Id: Ie0c02fdfbe03988f9d3b20d10cd9fe4c002d1279 Fix integer overflow in CFac_CalcFacSignal() original-Change-Id: Id2d767c40066c591b51768e978eb8af3b803f0c5 Fix integer overflow in FDKaacEnc_FDKaacEnc_calcPeNoAH() original-Change-Id: Idcbd0f4a51ae2550ed106aa6f3d678d1f9724841 Fix integer overflow in sbrDecoder_calculateGainVec() original-Change-Id: I7081bcbe29c5cede9821b38d93de07c7add2d507 Fix integer overflow in CLpc_SynthesisLattice() original-Change-Id: I4a95ddc18de150102352d4a1845f06094764c881 Fix integer overflow in Pred_Lt4() original-Change-Id: I4dbd012b2de7d07c3e70a47b92e3bfae8dbc750a Fix integer overflow in FDKsbrEnc_InitSbrFastTransientDetector() original-Change-Id: I788cbec1a4a00f44c2f3a72ad7a4afa219807d04 Fix unsigned integer overflow in FDKaacEnc_WriteBitstream() original-Change-Id: I68fc75166e7d2cd5cd45b18dbe3d8c2a92f1822a Fix unsigned integer overflow in FDK_MetadataEnc_Init() original-Change-Id: Ie8d025f9bcdb2442c704bd196e61065c03c10af4 Fix overflow in pseudo random number generators original-Change-Id: I3e2551ee01356297ca14e3788436ede80bd5513c Fix unsigned integer overflow in sbrDecoder_Parse() original-Change-Id: I3f231b2f437e9c37db4d5b964164686710eee971 Fix unsigned integer overflow in longsub() original-Change-Id: I73c2bc50415cac26f1f5a29e125bbe75f9180a6e Fix unsigned integer overflow in CAacDecoder_DecodeFrame() original-Change-Id: Ifce2db4b1454b46fa5f887e9d383f1cc43b291e4 Fix overflow at CLpdChannelStream_Read() original-Change-Id: Idb9d822ce3a4272e4794b643644f5434e2d4bf3f Fix unsigned integer overflow in Hcr_State_BODY_SIGN_ESC__ESC_WORD() original-Change-Id: I1ccf77c0015684b85534c5eb97162740a870b71c Fix unsigned integer overflow in UsacConfig_Parse() original-Change-Id: Ie6d27f84b6ae7eef092ecbff4447941c77864d9f Fix unsigned integer overflow in aacDecoder_drcParse() original-Change-Id: I713f28e883eea3d70b6fa56a7b8f8c22bcf66ca0 Fix unsigned integer overflow in aacDecoder_drcReadCompression() original-Change-Id: Ia34dfeb88c4705c558bce34314f584965cafcf7a Fix unsigned integer overflow in CDataStreamElement_Read() original-Change-Id: Iae896cc1d11f0a893d21be6aa90bd3e60a2c25f0 Fix unsigned integer overflow in transportDec_AdjustEndOfAccessUnit() original-Change-Id: I64cf29a153ee784bb4a16fdc088baabebc0007dc Fix unsigned integer overflow in transportDec_GetAuBitsRemaining() original-Change-Id: I975b3420faa9c16a041874ba0db82e92035962e4 Fix unsigned integer overflow in extractExtendedData() original-Change-Id: I2a59eb09e2053cfb58dfb75fcecfad6b85a80a8f Fix signed integer overflow in CAacDecoder_ExtPayloadParse() original-Change-Id: I4ad5ca4e3b83b5d964f1c2f8c5e7b17c477c7929 Fix unsigned integer overflow in CAacDecoder_DecodeFrame() original-Change-Id: I29a39df77d45c52a0c9c5c83c1ba81f8d0f25090 Follow-up on: Fix integer overflow in CLpc_SynthesisLattice() original-Change-Id: I8fb194ffc073a3432a380845be71036a272d388f Fix signed integer overflow in _interpolateDrcGain() original-Change-Id: I879ec9ab14005069a7c47faf80e8bc6e03d22e60 Fix unsigned integer overflow in FDKreadBits() original-Change-Id: I1f47a6a8037ff70375aa8844947d5681bb4287ad Fix unsigned integer overflow in FDKbyteAlign() original-Change-Id: Id5f3a11a0c9e50fc6f76ed6c572dbd4e9f2af766 Fix unsigned integer overflow in FDK_get32() original-Change-Id: I9d33b8e97e3d38cbb80629cb859266ca0acdce96 Fix unsigned integer overflow in FDK_pushBack() original-Change-Id: Ic87f899bc8c6acf7a377a8ca7f3ba74c3a1e1c19 Fix unsigned integer overflow in FDK_pushForward() original-Change-Id: I3b754382f6776a34be1602e66694ede8e0b8effc Fix unsigned integer overflow in ReadPsData() original-Change-Id: I25361664ba8139e32bbbef2ca8c106a606ce9c37 Fix signed integer overflow in E_UTIL_residu() original-Change-Id: I8c3abd1f437ee869caa8fb5903ce7d3d641b6aad REVERT: Follow-up on: Integer overflow in CLpc_SynthesisLattice(). original-Change-Id: I3d340099acb0414795c8dfbe6362bc0a8f045f9b Follow-up on: Fix integer overflow in CLpc_SynthesisLattice() original-Change-Id: I4aedb8b3a187064e9f4d985175aa55bb99cc7590 Follow-up on: Fix unsigned integer overflow in aacDecoder_drcParse() original-Change-Id: I2aa2e13916213bf52a67e8b0518e7bf7e57fb37d Fix integer overflow in acelp original-Change-Id: Ie6390c136d84055f8b728aefbe4ebef6e029dc77 Fix unsigned integer overflow in aacDecoder_UpdateBitStreamCounters() original-Change-Id: I391ffd97ddb0b2c184cba76139bfb356a3b4d2e2 Adjust concealment default settings original-Change-Id: I6a95db935a327c47df348030bcceafcb29f54b21 Saturate estimatedStartPos original-Change-Id: I27be2085e0ae83ec9501409f65e003f6bcba1ab6 Negative shift exponent in _interpolateDrcGain() original-Change-Id: I18edb26b26d002aafd5e633d4914960f7a359c29 Negative shift exponent in calculateICC() original-Change-Id: I3dcd2ae98d2eb70ee0d59750863cbb2a6f4f8aba Too large shift exponent in FDK_put() original-Change-Id: Ib7d9aaa434d2d8de4a13b720ca0464b31ca9b671 Too large shift exponent in CalcInvLdData() original-Change-Id: I43e6e78d4cd12daeb1dcd5d82d1798bdc2550262 Member access within null pointer of type SBR_CHANNEL original-Change-Id: Idc5e4ea8997810376d2f36bbdf628923b135b097 Member access within null pointer of type CpePersistentData original-Change-Id: Ib6c91cb0d37882768e5baf63324e429589de0d9d Member access within null pointer FDKaacEnc_psyMain() original-Change-Id: I7729b7f4479970531d9dc823abff63ca52e01997 Member access within null pointer FDKaacEnc_GetPnsParam() original-Change-Id: I9aa3b9f3456ae2e0f7483dbd5b3dde95fc62da39 Member access within null pointer FDKsbrEnc_EnvEncodeFrame() original-Change-Id: I67936f90ea714e90b3e81bc0dd1472cc713eb23a Add HCR sanity check original-Change-Id: I6c1d9732ebcf6af12f50b7641400752f74be39f7 Fix memory issue for HBE edge case with 8:3 SBR original-Change-Id: I11ea58a61e69fbe8bf75034b640baee3011e63e9 Additional SBR parametrization sanity check for ELD original-Change-Id: Ie26026fbfe174c2c7b3691f6218b5ce63e322140 Add MPEG-D DRC channel layout check original-Change-Id: Iea70a74f171b227cce636a9eac4ba662777a2f72 Additional out-of-bounds checks in MPEG-D DRC original-Change-Id: Ife4a8c3452c6fde8a0a09e941154a39a769777d4 Change-Id: Ic63cb2f628720f54fe9b572b0cb528e2599c624e --- libFDK/include/fixpoint_math.h | 902 ++++++++++++++++++++++++++++++----------- 1 file changed, 664 insertions(+), 238 deletions(-) (limited to 'libFDK/include/fixpoint_math.h') diff --git a/libFDK/include/fixpoint_math.h b/libFDK/include/fixpoint_math.h index 0d50f0a..3805892 100644 --- a/libFDK/include/fixpoint_math.h +++ b/libFDK/include/fixpoint_math.h @@ -1,74 +1,85 @@ - -/* ----------------------------------------------------------------------------------------------------------- +/* ----------------------------------------------------------------------------- Software License for The Fraunhofer FDK AAC Codec Library for Android -© Copyright 1995 - 2015 Fraunhofer-Gesellschaft zur Förderung der angewandten Forschung e.V. - All rights reserved. +© Copyright 1995 - 2018 Fraunhofer-Gesellschaft zur Förderung der angewandten +Forschung e.V. All rights reserved. 1. INTRODUCTION -The Fraunhofer FDK AAC Codec Library for Android ("FDK AAC Codec") is software that implements -the MPEG Advanced Audio Coding ("AAC") encoding and decoding scheme for digital audio. -This FDK AAC Codec software is intended to be used on a wide variety of Android devices. - -AAC's HE-AAC and HE-AAC v2 versions are regarded as today's most efficient general perceptual -audio codecs. AAC-ELD is considered the best-performing full-bandwidth communications codec by -independent studies and is widely deployed. AAC has been standardized by ISO and IEC as part -of the MPEG specifications. - -Patent licenses for necessary patent claims for the FDK AAC Codec (including those of Fraunhofer) -may be obtained through Via Licensing (www.vialicensing.com) or through the respective patent owners -individually for the purpose of encoding or decoding bit streams in products that are compliant with -the ISO/IEC MPEG audio standards. Please note that most manufacturers of Android devices already license -these patent claims through Via Licensing or directly from the patent owners, and therefore FDK AAC Codec -software may already be covered under those patent licenses when it is used for those licensed purposes only. - -Commercially-licensed AAC software libraries, including floating-point versions with enhanced sound quality, -are also available from Fraunhofer. Users are encouraged to check the Fraunhofer website for additional -applications information and documentation. +The Fraunhofer FDK AAC Codec Library for Android ("FDK AAC Codec") is software +that implements the MPEG Advanced Audio Coding ("AAC") encoding and decoding +scheme for digital audio. This FDK AAC Codec software is intended to be used on +a wide variety of Android devices. + +AAC's HE-AAC and HE-AAC v2 versions are regarded as today's most efficient +general perceptual audio codecs. AAC-ELD is considered the best-performing +full-bandwidth communications codec by independent studies and is widely +deployed. AAC has been standardized by ISO and IEC as part of the MPEG +specifications. + +Patent licenses for necessary patent claims for the FDK AAC Codec (including +those of Fraunhofer) may be obtained through Via Licensing +(www.vialicensing.com) or through the respective patent owners individually for +the purpose of encoding or decoding bit streams in products that are compliant +with the ISO/IEC MPEG audio standards. Please note that most manufacturers of +Android devices already license these patent claims through Via Licensing or +directly from the patent owners, and therefore FDK AAC Codec software may +already be covered under those patent licenses when it is used for those +licensed purposes only. + +Commercially-licensed AAC software libraries, including floating-point versions +with enhanced sound quality, are also available from Fraunhofer. Users are +encouraged to check the Fraunhofer website for additional applications +information and documentation. 2. COPYRIGHT LICENSE -Redistribution and use in source and binary forms, with or without modification, are permitted without -payment of copyright license fees provided that you satisfy the following conditions: +Redistribution and use in source and binary forms, with or without modification, +are permitted without payment of copyright license fees provided that you +satisfy the following conditions: -You must retain the complete text of this software license in redistributions of the FDK AAC Codec or -your modifications thereto in source code form. +You must retain the complete text of this software license in redistributions of +the FDK AAC Codec or your modifications thereto in source code form. -You must retain the complete text of this software license in the documentation and/or other materials -provided with redistributions of the FDK AAC Codec or your modifications thereto in binary form. -You must make available free of charge copies of the complete source code of the FDK AAC Codec and your +You must retain the complete text of this software license in the documentation +and/or other materials provided with redistributions of the FDK AAC Codec or +your modifications thereto in binary form. You must make available free of +charge copies of the complete source code of the FDK AAC Codec and your modifications thereto to recipients of copies in binary form. -The name of Fraunhofer may not be used to endorse or promote products derived from this library without -prior written permission. +The name of Fraunhofer may not be used to endorse or promote products derived +from this library without prior written permission. -You may not charge copyright license fees for anyone to use, copy or distribute the FDK AAC Codec -software or your modifications thereto. +You may not charge copyright license fees for anyone to use, copy or distribute +the FDK AAC Codec software or your modifications thereto. -Your modified versions of the FDK AAC Codec must carry prominent notices stating that you changed the software -and the date of any change. For modified versions of the FDK AAC Codec, the term -"Fraunhofer FDK AAC Codec Library for Android" must be replaced by the term -"Third-Party Modified Version of the Fraunhofer FDK AAC Codec Library for Android." +Your modified versions of the FDK AAC Codec must carry prominent notices stating +that you changed the software and the date of any change. For modified versions +of the FDK AAC Codec, the term "Fraunhofer FDK AAC Codec Library for Android" +must be replaced by the term "Third-Party Modified Version of the Fraunhofer FDK +AAC Codec Library for Android." 3. NO PATENT LICENSE -NO EXPRESS OR IMPLIED LICENSES TO ANY PATENT CLAIMS, including without limitation the patents of Fraunhofer, -ARE GRANTED BY THIS SOFTWARE LICENSE. Fraunhofer provides no warranty of patent non-infringement with -respect to this software. +NO EXPRESS OR IMPLIED LICENSES TO ANY PATENT CLAIMS, including without +limitation the patents of Fraunhofer, ARE GRANTED BY THIS SOFTWARE LICENSE. +Fraunhofer provides no warranty of patent non-infringement with respect to this +software. -You may use this FDK AAC Codec software or modifications thereto only for purposes that are authorized -by appropriate patent licenses. +You may use this FDK AAC Codec software or modifications thereto only for +purposes that are authorized by appropriate patent licenses. 4. DISCLAIMER -This FDK AAC Codec software is provided by Fraunhofer on behalf of the copyright holders and contributors -"AS IS" and WITHOUT ANY EXPRESS OR IMPLIED WARRANTIES, including but not limited to the implied warranties -of merchantability and fitness for a particular purpose. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR -CONTRIBUTORS BE LIABLE for any direct, indirect, incidental, special, exemplary, or consequential damages, -including but not limited to procurement of substitute goods or services; loss of use, data, or profits, -or business interruption, however caused and on any theory of liability, whether in contract, strict -liability, or tort (including negligence), arising in any way out of the use of this software, even if -advised of the possibility of such damage. +This FDK AAC Codec software is provided by Fraunhofer on behalf of the copyright +holders and contributors "AS IS" and WITHOUT ANY EXPRESS OR IMPLIED WARRANTIES, +including but not limited to the implied warranties of merchantability and +fitness for a particular purpose. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR +CONTRIBUTORS BE LIABLE for any direct, indirect, incidental, special, exemplary, +or consequential damages, including but not limited to procurement of substitute +goods or services; loss of use, data, or profits, or business interruption, +however caused and on any theory of liability, whether in contract, strict +liability, or tort (including negligence), arising in any way out of the use of +this software, even if advised of the possibility of such damage. 5. CONTACT INFORMATION @@ -79,21 +90,77 @@ Am Wolfsmantel 33 www.iis.fraunhofer.de/amm amm-info@iis.fraunhofer.de ------------------------------------------------------------------------------------------------------------ */ +----------------------------------------------------------------------------- */ -/*************************** Fraunhofer IIS FDK Tools ********************** +/******************* Library for basic calculation routines ******************** Author(s): M. Gayer - Description: Fixed point specific mathematical functions -******************************************************************************/ + Description: Fixed point specific mathematical functions -#ifndef __fixpoint_math_H -#define __fixpoint_math_H +*******************************************************************************/ +#ifndef FIXPOINT_MATH_H +#define FIXPOINT_MATH_H #include "common_fix.h" +#include "scale.h" + +/* + * Data definitions + */ + +#define LD_DATA_SCALING (64.0f) +#define LD_DATA_SHIFT 6 /* pow(2, LD_DATA_SHIFT) = LD_DATA_SCALING */ + +#define MAX_LD_PRECISION 10 +#define LD_PRECISION 10 + +/* Taylor series coefficients for ln(1-x), centered at 0 (MacLaurin polynomial). + */ +#ifndef LDCOEFF_16BIT +LNK_SECTION_CONSTDATA_L1 +static const FIXP_DBL ldCoeff[MAX_LD_PRECISION] = { + FL2FXCONST_DBL(-1.0), FL2FXCONST_DBL(-1.0 / 2.0), + FL2FXCONST_DBL(-1.0 / 3.0), FL2FXCONST_DBL(-1.0 / 4.0), + FL2FXCONST_DBL(-1.0 / 5.0), FL2FXCONST_DBL(-1.0 / 6.0), + FL2FXCONST_DBL(-1.0 / 7.0), FL2FXCONST_DBL(-1.0 / 8.0), + FL2FXCONST_DBL(-1.0 / 9.0), FL2FXCONST_DBL(-1.0 / 10.0)}; +#else /* LDCOEFF_16BIT */ +LNK_SECTION_CONSTDATA_L1 +static const FIXP_SGL ldCoeff[MAX_LD_PRECISION] = { + FL2FXCONST_SGL(-1.0), FL2FXCONST_SGL(-1.0 / 2.0), + FL2FXCONST_SGL(-1.0 / 3.0), FL2FXCONST_SGL(-1.0 / 4.0), + FL2FXCONST_SGL(-1.0 / 5.0), FL2FXCONST_SGL(-1.0 / 6.0), + FL2FXCONST_SGL(-1.0 / 7.0), FL2FXCONST_SGL(-1.0 / 8.0), + FL2FXCONST_SGL(-1.0 / 9.0), FL2FXCONST_SGL(-1.0 / 10.0)}; +#endif /* LDCOEFF_16BIT */ + +/***************************************************************************** + functionname: invSqrtNorm2 + description: delivers 1/sqrt(op) normalized to .5...1 and the shift value +of the OUTPUT + +*****************************************************************************/ +#define SQRT_BITS 7 +#define SQRT_VALUES (128 + 2) +#define SQRT_BITS_MASK 0x7f +#define SQRT_FRACT_BITS_MASK 0x007FFFFF + +extern const FIXP_DBL invSqrtTab[SQRT_VALUES]; + +/* + * Hardware specific implementations + */ + +#if defined(__x86__) +#include "x86/fixpoint_math_x86.h" +#endif /* target architecture selector */ + +/* + * Fallback implementations + */ #if !defined(FUNCTION_fIsLessThan) /** * \brief Compares two fixpoint values incl. scaling. @@ -103,61 +170,83 @@ amm-info@iis.fraunhofer.de * \param b_e exponent of the second input value. * \return non-zero if (a_m*2^a_e) < (b_m*2^b_e), 0 otherwise */ -FDK_INLINE INT fIsLessThan(FIXP_DBL a_m, INT a_e, FIXP_DBL b_m, INT b_e) -{ +FDK_INLINE INT fIsLessThan(FIXP_DBL a_m, INT a_e, FIXP_DBL b_m, INT b_e) { if (a_e > b_e) { - return (b_m >> fMin(a_e-b_e, DFRACT_BITS-1) > a_m); + return ((b_m >> fMin(a_e - b_e, DFRACT_BITS - 1)) > a_m); } else { - return (a_m >> fMin(b_e-a_e, DFRACT_BITS-1) < b_m); + return ((a_m >> fMin(b_e - a_e, DFRACT_BITS - 1)) < b_m); } } -FDK_INLINE INT fIsLessThan(FIXP_SGL a_m, INT a_e, FIXP_SGL b_m, INT b_e) -{ +FDK_INLINE INT fIsLessThan(FIXP_SGL a_m, INT a_e, FIXP_SGL b_m, INT b_e) { if (a_e > b_e) { - return (b_m >> fMin(a_e-b_e, FRACT_BITS-1) > a_m); + return ((b_m >> fMin(a_e - b_e, FRACT_BITS - 1)) > a_m); } else { - return (a_m >> fMin(b_e-a_e, FRACT_BITS-1) < b_m); + return ((a_m >> fMin(b_e - a_e, FRACT_BITS - 1)) < b_m); } } #endif - - -#define LD_DATA_SCALING (64.0f) -#define LD_DATA_SHIFT 6 /* pow(2, LD_DATA_SHIFT) = LD_DATA_SCALING */ - /** * \brief deprecated. Use fLog2() instead. */ -FIXP_DBL CalcLdData(FIXP_DBL op); +#define CalcLdData(op) fLog2(op, 0) void LdDataVector(FIXP_DBL *srcVector, FIXP_DBL *destVector, INT number); -FIXP_DBL CalcInvLdData(FIXP_DBL op); +extern const UINT exp2_tab_long[32]; +extern const UINT exp2w_tab_long[32]; +extern const UINT exp2x_tab_long[32]; + +LNK_SECTION_CODE_L1 +FDK_INLINE FIXP_DBL CalcInvLdData(const FIXP_DBL x) { + int set_zero = (x < FL2FXCONST_DBL(-31.0 / 64.0)) ? 0 : 1; + int set_max = (x >= FL2FXCONST_DBL(31.0 / 64.0)) | (x == FL2FXCONST_DBL(0.0)); + + FIXP_SGL frac = (FIXP_SGL)((LONG)x & 0x3FF); + UINT index3 = (UINT)(LONG)(x >> 10) & 0x1F; + UINT index2 = (UINT)(LONG)(x >> 15) & 0x1F; + UINT index1 = (UINT)(LONG)(x >> 20) & 0x1F; + int exp = fMin(31, ((x > FL2FXCONST_DBL(0.0f)) ? (31 - (int)(x >> 25)) + : (int)(-(x >> 25)))); + + UINT lookup1 = exp2_tab_long[index1] * set_zero; + UINT lookup2 = exp2w_tab_long[index2]; + UINT lookup3 = exp2x_tab_long[index3]; + UINT lookup3f = + lookup3 + (UINT)(LONG)fMultDiv2((FIXP_DBL)(0x0016302F), (FIXP_SGL)frac); + + UINT lookup12 = (UINT)(LONG)fMult((FIXP_DBL)lookup1, (FIXP_DBL)lookup2); + UINT lookup = (UINT)(LONG)fMult((FIXP_DBL)lookup12, (FIXP_DBL)lookup3f); + FIXP_DBL retVal = (lookup << 3) >> exp; -void InitLdInt(); + if (set_max) { + retVal = (FIXP_DBL)MAXVAL_DBL; + } + + return retVal; +} + +void InitLdInt(); FIXP_DBL CalcLdInt(INT i); extern const USHORT sqrt_tab[49]; -inline FIXP_DBL sqrtFixp_lookup(FIXP_DBL x) -{ +inline FIXP_DBL sqrtFixp_lookup(FIXP_DBL x) { UINT y = (INT)x; - UCHAR is_zero=(y==0); - INT zeros=fixnormz_D(y) & 0x1e; - y<<=zeros; - UINT idx=(y>>26)-16; - USHORT frac=(y>>10)&0xffff; - USHORT nfrac=0xffff^frac; - UINT t=nfrac*sqrt_tab[idx]+frac*sqrt_tab[idx+1]; - t=t>>(zeros>>1); - return(is_zero ? 0 : t); + UCHAR is_zero = (y == 0); + INT zeros = fixnormz_D(y) & 0x1e; + y <<= zeros; + UINT idx = (y >> 26) - 16; + USHORT frac = (y >> 10) & 0xffff; + USHORT nfrac = 0xffff ^ frac; + UINT t = (UINT)nfrac * sqrt_tab[idx] + (UINT)frac * sqrt_tab[idx + 1]; + t = t >> (zeros >> 1); + return (is_zero ? 0 : t); } -inline FIXP_DBL sqrtFixp_lookup(FIXP_DBL x, INT *x_e) -{ +inline FIXP_DBL sqrtFixp_lookup(FIXP_DBL x, INT *x_e) { UINT y = (INT)x; INT e; @@ -166,106 +255,153 @@ inline FIXP_DBL sqrtFixp_lookup(FIXP_DBL x, INT *x_e) } /* Normalize */ - e=fixnormz_D(y); - y<<=e; - e = *x_e - e + 2; + e = fixnormz_D(y); + y <<= e; + e = *x_e - e + 2; /* Correct odd exponent. */ if (e & 1) { y >>= 1; - e ++; + e++; } /* Get square root */ - UINT idx=(y>>26)-16; - USHORT frac=(y>>10)&0xffff; - USHORT nfrac=0xffff^frac; - UINT t=nfrac*sqrt_tab[idx]+frac*sqrt_tab[idx+1]; + UINT idx = (y >> 26) - 16; + USHORT frac = (y >> 10) & 0xffff; + USHORT nfrac = 0xffff ^ frac; + UINT t = (UINT)nfrac * sqrt_tab[idx] + (UINT)frac * sqrt_tab[idx + 1]; /* Write back exponent */ *x_e = e >> 1; - return (FIXP_DBL)(LONG)(t>>1); + return (FIXP_DBL)(LONG)(t >> 1); } - - -FIXP_DBL sqrtFixp(FIXP_DBL op); - void InitInvSqrtTab(); -FIXP_DBL invSqrtNorm2(FIXP_DBL op, INT *shift); - +#ifndef FUNCTION_invSqrtNorm2 +/** + * \brief calculate 1.0/sqrt(op) + * \param op_m mantissa of input value. + * \param result_e pointer to return the exponent of the result + * \return mantissa of the result + */ /***************************************************************************** - - functionname: invFixp - description: delivers 1/(op) - + delivers 1/sqrt(op) normalized to .5...1 and the shift value of the OUTPUT, + i.e. the denormalized result is 1/sqrt(op) = invSqrtNorm(op) * 2^(shift) + uses Newton-iteration for approximation + Q(n+1) = Q(n) + Q(n) * (0.5 - 2 * V * Q(n)^2) + with Q = 0.5* V ^-0.5; 0.5 <= V < 1.0 *****************************************************************************/ -inline FIXP_DBL invFixp(FIXP_DBL op) -{ - INT tmp_exp ; - FIXP_DBL tmp_inv = invSqrtNorm2(op, &tmp_exp) ; - FDK_ASSERT((31-(2*tmp_exp+1))>=0) ; - return ( fPow2Div2( (FIXP_DBL)tmp_inv ) >> (31-(2*tmp_exp+1)) ) ; -} +static FDK_FORCEINLINE FIXP_DBL invSqrtNorm2(FIXP_DBL op, INT *shift) { + FIXP_DBL val = op; + FIXP_DBL reg1, reg2; + if (val == FL2FXCONST_DBL(0.0)) { + *shift = 16; + return ((LONG)MAXVAL_DBL); /* maximum positive value */ + } +#define INVSQRTNORM2_LINEAR_INTERPOLATE +#define INVSQRTNORM2_LINEAR_INTERPOLATE_HQ + + /* normalize input, calculate shift value */ + FDK_ASSERT(val > FL2FXCONST_DBL(0.0)); + *shift = fNormz(val) - 1; /* CountLeadingBits() is not necessary here since + test value is always > 0 */ + val <<= *shift; /* normalized input V */ + *shift += 2; /* bias for exponent */ + +#if defined(INVSQRTNORM2_LINEAR_INTERPOLATE) + INT index = + (INT)(val >> (DFRACT_BITS - 1 - (SQRT_BITS + 1))) & SQRT_BITS_MASK; + FIXP_DBL Fract = + (FIXP_DBL)(((INT)val & SQRT_FRACT_BITS_MASK) << (SQRT_BITS + 1)); + FIXP_DBL diff = invSqrtTab[index + 1] - invSqrtTab[index]; + reg1 = invSqrtTab[index] + (fMultDiv2(diff, Fract) << 1); +#if defined(INVSQRTNORM2_LINEAR_INTERPOLATE_HQ) + /* reg1 = t[i] + (t[i+1]-t[i])*fract ... already computed ... + + (1-fract)fract*(t[i+2]-t[i+1])/2 */ + if (Fract != (FIXP_DBL)0) { + /* fract = fract * (1 - fract) */ + Fract = fMultDiv2(Fract, (FIXP_DBL)((ULONG)0x80000000 - (ULONG)Fract)) << 1; + diff = diff - (invSqrtTab[index + 2] - invSqrtTab[index + 1]); + reg1 = fMultAddDiv2(reg1, Fract, diff); + } +#endif /* INVSQRTNORM2_LINEAR_INTERPOLATE_HQ */ +#else +#error \ + "Either define INVSQRTNORM2_NEWTON_ITERATE or INVSQRTNORM2_LINEAR_INTERPOLATE" +#endif + /* calculate the output exponent = input exp/2 */ + if (*shift & 0x00000001) { /* odd shift values ? */ + /* Note: Do not use rounded value 0x5A82799A to avoid overflow with + * shift-by-2 */ + reg2 = (FIXP_DBL)0x5A827999; + /* FL2FXCONST_DBL(0.707106781186547524400844362104849f);*/ /* 1/sqrt(2); + */ + reg1 = fMultDiv2(reg1, reg2) << 2; + } -#if defined(__mips__) && (__GNUC__==2) + *shift = *shift >> 1; -#define FUNCTION_schur_div -inline FIXP_DBL schur_div(FIXP_DBL num,FIXP_DBL denum, INT count) -{ - INT result, tmp ; - __asm__ ("srl %1, %2, 15\n" - "div %3, %1\n" : "=lo" (result) - : "%d" (tmp), "d" (denum) , "d" (num) - : "hi" ) ; - return result<<16 ; + return (reg1); } +#endif /* FUNCTION_invSqrtNorm2 */ -/*###########################################################################################*/ -#elif defined(__mips__) && (__GNUC__==3) - -#define FUNCTION_schur_div -inline FIXP_DBL schur_div(FIXP_DBL num,FIXP_DBL denum, INT count) -{ - INT result, tmp; +#ifndef FUNCTION_sqrtFixp +static FDK_FORCEINLINE FIXP_DBL sqrtFixp(FIXP_DBL op) { + INT tmp_exp = 0; + FIXP_DBL tmp_inv = invSqrtNorm2(op, &tmp_exp); - __asm__ ("srl %[tmp], %[denum], 15\n" - "div %[result], %[num], %[tmp]\n" - : [tmp] "+r" (tmp), [result]"=r"(result) - : [denum]"r"(denum), [num]"r"(num) - : "hi", "lo"); - return result << (DFRACT_BITS-16); + FDK_ASSERT(tmp_exp > 0); + return ((FIXP_DBL)(fMultDiv2((op << (tmp_exp - 1)), tmp_inv) << 2)); } +#endif /* FUNCTION_sqrtFixp */ -/*###########################################################################################*/ -#elif defined(SIMULATE_MIPS_DIV) - -#define FUNCTION_schur_div -inline FIXP_DBL schur_div(FIXP_DBL num, FIXP_DBL denum, INT count) -{ - FDK_ASSERT (count<=DFRACT_BITS-1); - FDK_ASSERT (num>=(FIXP_DBL)0); - FDK_ASSERT (denum>(FIXP_DBL)0); - FDK_ASSERT (num <= denum); +#ifndef FUNCTION_invFixp +/** + * \brief calculate 1.0/op + * \param op mantissa of the input value. + * \return mantissa of the result with implicit exponent of 31 + * \exceptions are provided for op=0,1 setting max. positive value + */ +static inline FIXP_DBL invFixp(FIXP_DBL op) { + if ((op == (FIXP_DBL)0x00000000) || (op == (FIXP_DBL)0x00000001)) { + return ((LONG)MAXVAL_DBL); + } + INT tmp_exp; + FIXP_DBL tmp_inv = invSqrtNorm2(op, &tmp_exp); + FDK_ASSERT((31 - (2 * tmp_exp + 1)) >= 0); + int shift = 31 - (2 * tmp_exp + 1); + tmp_inv = fPow2Div2(tmp_inv); + if (shift) { + tmp_inv = ((tmp_inv >> (shift - 1)) + (FIXP_DBL)1) >> 1; + } + return tmp_inv; +} - INT tmp = denum >> (count-1); - INT result = 0; +/** + * \brief calculate 1.0/(op_m * 2^op_e) + * \param op_m mantissa of the input value. + * \param op_e pointer into were the exponent of the input value is stored, and + * the result will be stored into. + * \return mantissa of the result + */ +static inline FIXP_DBL invFixp(FIXP_DBL op_m, int *op_e) { + if ((op_m == (FIXP_DBL)0x00000000) || (op_m == (FIXP_DBL)0x00000001)) { + *op_e = 31 - *op_e; + return ((LONG)MAXVAL_DBL); + } - while (num > tmp) - { - num -= tmp; - result++; - } + INT tmp_exp; + FIXP_DBL tmp_inv = invSqrtNorm2(op_m, &tmp_exp); - return result << (DFRACT_BITS-count); + *op_e = (tmp_exp << 1) - *op_e + 1; + return fPow2Div2(tmp_inv); } +#endif /* FUNCTION_invFixp */ -/*###########################################################################################*/ -#endif /* target architecture selector */ +#ifndef FUNCTION_schur_div -#if !defined(FUNCTION_schur_div) /** * \brief Divide two FIXP_DBL values with given precision. * \param num dividend @@ -273,31 +409,34 @@ inline FIXP_DBL schur_div(FIXP_DBL num, FIXP_DBL denum, INT count) * \param count amount of significant bits of the result (starting to the MSB) * \return num/divisor */ -FIXP_DBL schur_div(FIXP_DBL num,FIXP_DBL denum, INT count); -#endif +FIXP_DBL schur_div(FIXP_DBL num, FIXP_DBL denum, INT count); +#endif /* FUNCTION_schur_div */ -FIXP_DBL mul_dbl_sgl_rnd (const FIXP_DBL op1, - const FIXP_SGL op2); +FIXP_DBL mul_dbl_sgl_rnd(const FIXP_DBL op1, const FIXP_SGL op2); +#ifndef FUNCTION_fMultNorm /** * \brief multiply two values with normalization, thus max precision. * Author: Robert Weidner * * \param f1 first factor - * \param f2 secod factor - * \param result_e pointer to an INT where the exponent of the result is stored into + * \param f2 second factor + * \param result_e pointer to an INT where the exponent of the result is stored + * into * \return mantissa of the product f1*f2 */ -FIXP_DBL fMultNorm( - FIXP_DBL f1, - FIXP_DBL f2, - INT *result_e - ); +FIXP_DBL fMultNorm(FIXP_DBL f1, FIXP_DBL f2, INT *result_e); -inline FIXP_DBL fMultNorm(FIXP_DBL f1, FIXP_DBL f2) -{ +/** + * \brief Multiply 2 values using maximum precision. The exponent of the result + * is 0. + * \param f1_m mantissa of factor 1 + * \param f2_m mantissa of factor 2 + * \return mantissa of the result with exponent equal to 0 + */ +inline FIXP_DBL fMultNorm(FIXP_DBL f1, FIXP_DBL f2) { FIXP_DBL m; INT e; @@ -308,52 +447,253 @@ inline FIXP_DBL fMultNorm(FIXP_DBL f1, FIXP_DBL f2) return m; } +/** + * \brief Multiply 2 values with exponent and use given exponent for the + * mantissa of the result. + * \param f1_m mantissa of factor 1 + * \param f1_e exponent of factor 1 + * \param f2_m mantissa of factor 2 + * \param f2_e exponent of factor 2 + * \param result_e exponent for the returned mantissa of the result + * \return mantissa of the result with exponent equal to result_e + */ +inline FIXP_DBL fMultNorm(FIXP_DBL f1_m, INT f1_e, FIXP_DBL f2_m, INT f2_e, + INT result_e) { + FIXP_DBL m; + INT e; + + m = fMultNorm(f1_m, f2_m, &e); + + m = scaleValueSaturate(m, e + f1_e + f2_e - result_e); + + return m; +} +#endif /* FUNCTION_fMultNorm */ + +#ifndef FUNCTION_fMultI +/** + * \brief Multiplies a fractional value and a integer value and performs + * rounding to nearest + * \param a fractional value + * \param b integer value + * \return integer value + */ +inline INT fMultI(FIXP_DBL a, INT b) { + FIXP_DBL m, mi; + INT m_e; + + m = fMultNorm(a, (FIXP_DBL)b, &m_e); + + if (m_e < (INT)0) { + if (m_e > (INT)-DFRACT_BITS) { + m = m >> ((-m_e) - 1); + mi = (m + (FIXP_DBL)1) >> 1; + } else { + mi = (FIXP_DBL)0; + } + } else { + mi = scaleValueSaturate(m, m_e); + } + + return ((INT)mi); +} +#endif /* FUNCTION_fMultI */ + +#ifndef FUNCTION_fMultIfloor +/** + * \brief Multiplies a fractional value and a integer value and performs floor + * rounding + * \param a fractional value + * \param b integer value + * \return integer value + */ +inline INT fMultIfloor(FIXP_DBL a, INT b) { + FIXP_DBL m, mi; + INT m_e; + + m = fMultNorm(a, (FIXP_DBL)b, &m_e); + + if (m_e < (INT)0) { + if (m_e > (INT)-DFRACT_BITS) { + mi = m >> (-m_e); + } else { + mi = (FIXP_DBL)0; + if (m < (FIXP_DBL)0) { + mi = (FIXP_DBL)-1; + } + } + } else { + mi = scaleValueSaturate(m, m_e); + } + + return ((INT)mi); +} +#endif /* FUNCTION_fMultIfloor */ + +#ifndef FUNCTION_fMultIceil +/** + * \brief Multiplies a fractional value and a integer value and performs ceil + * rounding + * \param a fractional value + * \param b integer value + * \return integer value + */ +inline INT fMultIceil(FIXP_DBL a, INT b) { + FIXP_DBL m, mi; + INT m_e; + + m = fMultNorm(a, (FIXP_DBL)b, &m_e); + + if (m_e < (INT)0) { + if (m_e > (INT)-DFRACT_BITS) { + mi = (m >> (-m_e)); + if ((LONG)m & ((1 << (-m_e)) - 1)) { + mi = mi + (FIXP_DBL)1; + } + } else { + mi = (FIXP_DBL)1; + if (m < (FIXP_DBL)0) { + mi = (FIXP_DBL)0; + } + } + } else { + mi = scaleValueSaturate(m, m_e); + } + + return ((INT)mi); +} +#endif /* FUNCTION_fMultIceil */ + +#ifndef FUNCTION_fDivNorm /** * \brief Divide 2 FIXP_DBL values with normalization of input values. * \param num numerator - * \param denum denomintator - * \return num/denum with exponent = 0 + * \param denum denominator + * \param result_e pointer to an INT where the exponent of the result is stored + * into + * \return num/denum with exponent = *result_e */ FIXP_DBL fDivNorm(FIXP_DBL num, FIXP_DBL denom, INT *result_e); /** - * \brief Divide 2 FIXP_DBL values with normalization of input values. + * \brief Divide 2 positive FIXP_DBL values with normalization of input values. * \param num numerator - * \param denum denomintator - * \param result_e pointer to an INT where the exponent of the result is stored into - * \return num/denum with exponent = *result_e + * \param denum denominator + * \return num/denum with exponent = 0 */ FIXP_DBL fDivNorm(FIXP_DBL num, FIXP_DBL denom); /** - * \brief Divide 2 FIXP_DBL values with normalization of input values. + * \brief Divide 2 signed FIXP_DBL values with normalization of input values. * \param num numerator - * \param denum denomintator + * \param denum denominator + * \param result_e pointer to an INT where the exponent of the result is stored + * into + * \return num/denum with exponent = *result_e + */ +FIXP_DBL fDivNormSigned(FIXP_DBL L_num, FIXP_DBL L_denum, INT *result_e); + +/** + * \brief Divide 2 signed FIXP_DBL values with normalization of input values. + * \param num numerator + * \param denum denominator * \return num/denum with exponent = 0 */ -FIXP_DBL fDivNormHighPrec(FIXP_DBL L_num, FIXP_DBL L_denum, INT *result_e); +FIXP_DBL fDivNormSigned(FIXP_DBL num, FIXP_DBL denom); +#endif /* FUNCTION_fDivNorm */ /** - * \brief Calculate log(argument)/log(2) (logarithm with base 2). deprecated. Use fLog2() instead. - * \param arg mantissa of the argument - * \param arg_e exponent of the argument - * \param result_e pointer to an INT to store the exponent of the result - * \return the mantissa of the result. - * \param + * \brief Adjust mantissa to exponent -1 + * \param a_m mantissa of value to be adjusted + * \param pA_e pointer to the exponen of a_m + * \return adjusted mantissa */ -FIXP_DBL CalcLog2(FIXP_DBL arg, INT arg_e, INT *result_e); +inline FIXP_DBL fAdjust(FIXP_DBL a_m, INT *pA_e) { + INT shift; + shift = fNorm(a_m) - 1; + *pA_e -= shift; + + return scaleValue(a_m, shift); +} + +#ifndef FUNCTION_fAddNorm /** - * \brief return 2 ^ (exp * 2^exp_e) + * \brief Add two values with normalization + * \param a_m mantissa of first summand + * \param a_e exponent of first summand + * \param a_m mantissa of second summand + * \param a_e exponent of second summand + * \param pResult_e pointer to where the exponent of the result will be stored + * to. + * \return mantissa of result + */ +inline FIXP_DBL fAddNorm(FIXP_DBL a_m, INT a_e, FIXP_DBL b_m, INT b_e, + INT *pResult_e) { + INT result_e; + FIXP_DBL result_m; + + /* If one of the summands is zero, return the other. + This is necessary for the summation of a very small number to zero */ + if (a_m == (FIXP_DBL)0) { + *pResult_e = b_e; + return b_m; + } + if (b_m == (FIXP_DBL)0) { + *pResult_e = a_e; + return a_m; + } + + a_m = fAdjust(a_m, &a_e); + b_m = fAdjust(b_m, &b_e); + + if (a_e > b_e) { + result_m = a_m + (b_m >> fMin(a_e - b_e, DFRACT_BITS - 1)); + result_e = a_e; + } else { + result_m = (a_m >> fMin(b_e - a_e, DFRACT_BITS - 1)) + b_m; + result_e = b_e; + } + + *pResult_e = result_e; + return result_m; +} + +inline FIXP_DBL fAddNorm(FIXP_DBL a_m, INT a_e, FIXP_DBL b_m, INT b_e, + INT result_e) { + FIXP_DBL result_m; + + a_m = scaleValue(a_m, a_e - result_e); + b_m = scaleValue(b_m, b_e - result_e); + + result_m = a_m + b_m; + + return result_m; +} +#endif /* FUNCTION_fAddNorm */ + +/** + * \brief Divide 2 FIXP_DBL values with normalization of input values. + * \param num numerator + * \param denum denomintator + * \return num/denum with exponent = 0 + */ +FIXP_DBL fDivNormHighPrec(FIXP_DBL L_num, FIXP_DBL L_denum, INT *result_e); + +#ifndef FUNCTION_fPow +/** + * \brief return 2 ^ (exp_m * 2^exp_e) * \param exp_m mantissa of the exponent to 2.0f * \param exp_e exponent of the exponent to 2.0f - * \param result_e pointer to a INT where the exponent of the result will be stored into + * \param result_e pointer to a INT where the exponent of the result will be + * stored into * \return mantissa of the result */ FIXP_DBL f2Pow(const FIXP_DBL exp_m, const INT exp_e, INT *result_e); /** - * \brief return 2 ^ (exp_m * 2^exp_e). This version returns only the mantissa with implicit exponent of zero. + * \brief return 2 ^ (exp_m * 2^exp_e). This version returns only the mantissa + * with implicit exponent of zero. * \param exp_m mantissa of the exponent to 2.0f * \param exp_e exponent of the exponent to 2.0f * \return mantissa of the result @@ -361,57 +701,70 @@ FIXP_DBL f2Pow(const FIXP_DBL exp_m, const INT exp_e, INT *result_e); FIXP_DBL f2Pow(const FIXP_DBL exp_m, const INT exp_e); /** - * \brief return x ^ (exp * 2^exp_e), where log2(x) = baseLd_m * 2^(baseLd_e). This saves - * the need to compute log2() of constant values (when x is a constant). - * \param ldx_m mantissa of log2() of x. - * \param ldx_e exponent of log2() of x. + * \brief return x ^ (exp_m * 2^exp_e), where log2(x) = baseLd_m * 2^(baseLd_e). + * This saves the need to compute log2() of constant values (when x is a + * constant). + * \param baseLd_m mantissa of log2() of x. + * \param baseLd_e exponent of log2() of x. * \param exp_m mantissa of the exponent to 2.0f * \param exp_e exponent of the exponent to 2.0f - * \param result_e pointer to a INT where the exponent of the result will be stored into + * \param result_e pointer to a INT where the exponent of the result will be + * stored into * \return mantissa of the result */ -FIXP_DBL fLdPow( - FIXP_DBL baseLd_m, - INT baseLd_e, - FIXP_DBL exp_m, INT exp_e, - INT *result_e - ); +FIXP_DBL fLdPow(FIXP_DBL baseLd_m, INT baseLd_e, FIXP_DBL exp_m, INT exp_e, + INT *result_e); /** - * \brief return x ^ (exp * 2^exp_e), where log2(x) = baseLd_m * 2^(baseLd_e). This saves - * the need to compute log2() of constant values (when x is a constant). This version - * does not return an exponent, which is implicitly 0. - * \param ldx_m mantissa of log2() of x. - * \param ldx_e exponent of log2() of x. + * \brief return x ^ (exp_m * 2^exp_e), where log2(x) = baseLd_m * 2^(baseLd_e). + * This saves the need to compute log2() of constant values (when x is a + * constant). This version does not return an exponent, which is + * implicitly 0. + * \param baseLd_m mantissa of log2() of x. + * \param baseLd_e exponent of log2() of x. * \param exp_m mantissa of the exponent to 2.0f * \param exp_e exponent of the exponent to 2.0f * \return mantissa of the result */ -FIXP_DBL fLdPow( - FIXP_DBL baseLd_m, INT baseLd_e, - FIXP_DBL exp_m, INT exp_e - ); +FIXP_DBL fLdPow(FIXP_DBL baseLd_m, INT baseLd_e, FIXP_DBL exp_m, INT exp_e); /** - * \brief return (base * 2^base_e) ^ (exp * 2^exp_e). Use fLdPow() instead whenever possible. + * \brief return (base_m * 2^base_e) ^ (exp * 2^exp_e). Use fLdPow() instead + * whenever possible. * \param base_m mantissa of the base. * \param base_e exponent of the base. * \param exp_m mantissa of power to be calculated of the base. * \param exp_e exponent of power to be calculated of the base. - * \param result_e pointer to a INT where the exponent of the result will be stored into. + * \param result_e pointer to a INT where the exponent of the result will be + * stored into. * \return mantissa of the result. */ -FIXP_DBL fPow(FIXP_DBL base_m, INT base_e, FIXP_DBL exp_m, INT exp_e, INT *result_e); +FIXP_DBL fPow(FIXP_DBL base_m, INT base_e, FIXP_DBL exp_m, INT exp_e, + INT *result_e); /** - * \brief return (base * 2^base_e) ^ N - * \param base mantissa of the base + * \brief return (base_m * 2^base_e) ^ N + * \param base_m mantissa of the base * \param base_e exponent of the base - * \param power to be calculated of the base - * \param result_e pointer to a INT where the exponent of the result will be stored into + * \param N power to be calculated of the base + * \param result_e pointer to a INT where the exponent of the result will be + * stored into * \return mantissa of the result */ FIXP_DBL fPowInt(FIXP_DBL base_m, INT base_e, INT N, INT *result_e); +#endif /* #ifndef FUNCTION_fPow */ + +#ifndef FUNCTION_fLog2 +/** + * \brief Calculate log(argument)/log(2) (logarithm with base 2). deprecated. + * Use fLog2() instead. + * \param arg mantissa of the argument + * \param arg_e exponent of the argument + * \param result_e pointer to an INT to store the exponent of the result + * \return the mantissa of the result. + * \param + */ +FIXP_DBL CalcLog2(FIXP_DBL arg, INT arg_e, INT *result_e); /** * \brief calculate logarithm of base 2 of x_m * 2^(x_e) @@ -420,7 +773,68 @@ FIXP_DBL fPowInt(FIXP_DBL base_m, INT base_e, INT N, INT *result_e); * \param pointer to an INT where the exponent of the result is returned into. * \return mantissa of the result. */ -FIXP_DBL fLog2(FIXP_DBL x_m, INT x_e, INT *result_e); +FDK_INLINE FIXP_DBL fLog2(FIXP_DBL x_m, INT x_e, INT *result_e) { + FIXP_DBL result_m; + + /* Short cut for zero and negative numbers. */ + if (x_m <= FL2FXCONST_DBL(0.0f)) { + *result_e = DFRACT_BITS - 1; + return FL2FXCONST_DBL(-1.0f); + } + + /* Calculate log2() */ + { + FIXP_DBL x2_m; + + /* Move input value x_m * 2^x_e toward 1.0, where the taylor approximation + of the function log(1-x) centered at 0 is most accurate. */ + { + INT b_norm; + + b_norm = fNormz(x_m) - 1; + x2_m = x_m << b_norm; + x_e = x_e - b_norm; + } + + /* map x from log(x) domain to log(1-x) domain. */ + x2_m = -(x2_m + FL2FXCONST_DBL(-1.0)); + + /* Taylor polynomial approximation of ln(1-x) */ + { + FIXP_DBL px2_m; + result_m = FL2FXCONST_DBL(0.0); + px2_m = x2_m; + for (int i = 0; i < LD_PRECISION; i++) { + result_m = fMultAddDiv2(result_m, ldCoeff[i], px2_m); + px2_m = fMult(px2_m, x2_m); + } + } + /* Multiply result with 1/ln(2) = 1.0 + 0.442695040888 (get log2(x) from + * ln(x) result). */ + result_m = + fMultAddDiv2(result_m, result_m, + FL2FXCONST_DBL(2.0 * 0.4426950408889634073599246810019)); + + /* Add exponent part. log2(x_m * 2^x_e) = log2(x_m) + x_e */ + if (x_e != 0) { + int enorm; + + enorm = DFRACT_BITS - fNorm((FIXP_DBL)x_e); + /* The -1 in the right shift of result_m compensates the fMultDiv2() above + * in the taylor polynomial evaluation loop.*/ + result_m = (result_m >> (enorm - 1)) + + ((FIXP_DBL)x_e << (DFRACT_BITS - 1 - enorm)); + + *result_e = enorm; + } else { + /* 1 compensates the fMultDiv2() above in the taylor polynomial evaluation + * loop.*/ + *result_e = 1; + } + } + + return result_m; +} /** * \brief calculate logarithm of base 2 of x_m * 2^(x_e) @@ -428,16 +842,27 @@ FIXP_DBL fLog2(FIXP_DBL x_m, INT x_e, INT *result_e); * \param x_e exponent of the input value. * \return mantissa of the result with implicit exponent of LD_DATA_SHIFT. */ -FIXP_DBL fLog2(FIXP_DBL x_m, INT x_e); +FDK_INLINE FIXP_DBL fLog2(FIXP_DBL x_m, INT x_e) { + if (x_m <= FL2FXCONST_DBL(0.0f)) { + x_m = FL2FXCONST_DBL(-1.0f); + } else { + INT result_e; + x_m = fLog2(x_m, x_e, &result_e); + x_m = scaleValue(x_m, result_e - LD_DATA_SHIFT); + } + return x_m; +} +#endif /* FUNCTION_fLog2 */ + +#ifndef FUNCTION_fAddSaturate /** * \brief Add with saturation of the result. * \param a first summand * \param b second summand * \return saturated sum of a and b. */ -inline FIXP_SGL fAddSaturate(const FIXP_SGL a, const FIXP_SGL b) -{ +inline FIXP_SGL fAddSaturate(const FIXP_SGL a, const FIXP_SGL b) { LONG sum; sum = (LONG)(SHORT)a + (LONG)(SHORT)b; @@ -451,19 +876,26 @@ inline FIXP_SGL fAddSaturate(const FIXP_SGL a, const FIXP_SGL b) * \param b second summand * \return saturated sum of a and b. */ -inline FIXP_DBL fAddSaturate(const FIXP_DBL a, const FIXP_DBL b) -{ +inline FIXP_DBL fAddSaturate(const FIXP_DBL a, const FIXP_DBL b) { LONG sum; - sum = (LONG)(a>>1) + (LONG)(b>>1); - sum = fMax(fMin((INT)sum, (INT)(MAXVAL_DBL>>1)), (INT)(MINVAL_DBL>>1)); - return (FIXP_DBL)(LONG)(sum<<1); + sum = (LONG)(a >> 1) + (LONG)(b >> 1); + sum = fMax(fMin((INT)sum, (INT)(MAXVAL_DBL >> 1)), (INT)(MINVAL_DBL >> 1)); + return (FIXP_DBL)(LONG)(sum << 1); } +#endif /* FUNCTION_fAddSaturate */ -//#define TEST_ROUNDING +INT fixp_floorToInt(FIXP_DBL f_inp, INT sf); +FIXP_DBL fixp_floor(FIXP_DBL f_inp, INT sf); +INT fixp_ceilToInt(FIXP_DBL f_inp, INT sf); +FIXP_DBL fixp_ceil(FIXP_DBL f_inp, INT sf); +INT fixp_truncateToInt(FIXP_DBL f_inp, INT sf); +FIXP_DBL fixp_truncate(FIXP_DBL f_inp, INT sf); +INT fixp_roundToInt(FIXP_DBL f_inp, INT sf); +FIXP_DBL fixp_round(FIXP_DBL f_inp, INT sf); /***************************************************************************** @@ -471,25 +903,19 @@ inline FIXP_DBL fAddSaturate(const FIXP_DBL a, const FIXP_DBL b) ****************************************************************************/ - extern const FIXP_DBL invCount[80]; - - LNK_SECTION_INITCODE - inline void InitInvInt(void) {} +extern const FIXP_DBL invCount[80]; +LNK_SECTION_INITCODE +inline void InitInvInt(void) {} /** * \brief Calculate the value of 1/i where i is a integer value. It supports - * input values from 1 upto 80. + * input values from 1 upto (80-1). * \param intValue Integer input value. * \param FIXP_DBL representation of 1/intValue */ -inline FIXP_DBL GetInvInt(int intValue) -{ - FDK_ASSERT((intValue > 0) && (intValue < 80)); - FDK_ASSERT(intValue<80); - return invCount[intValue]; +inline FIXP_DBL GetInvInt(int intValue) { + return invCount[fMin(fMax(intValue, 0), 80 - 1)]; } - -#endif - +#endif /* FIXPOINT_MATH_H */ -- cgit v1.2.3