diff options
Diffstat (limited to 'src/lib/std_dtoa.c')
| -rw-r--r-- | src/lib/std_dtoa.c | 496 |
1 files changed, 496 insertions, 0 deletions
diff --git a/src/lib/std_dtoa.c b/src/lib/std_dtoa.c new file mode 100644 index 0000000..880f3c6 --- /dev/null +++ b/src/lib/std_dtoa.c @@ -0,0 +1,496 @@ +/* +* Copyright (c) 2017, The Linux Foundation. All rights reserved. +* All rights reserved. +* +* Redistribution and use in source and binary forms, with or without +* modification, are permitted provided that the following conditions are met: +* +* 1. Redistributions of source code must retain the above copyright notice, +* this list of conditions and the following disclaimer. +* +* 2. Redistributions in binary form must reproduce the above copyright notice, +* this list of conditions and the following disclaimer in the documentation +* and/or other materials provided with the distribution. +* +* 3. Neither the name of the copyright holder nor the names of its contributors +* may be used to endorse or promote products derived from this software without +* specific prior written permission. +* +* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +* ARE DISCLAIMED. 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 OR OTHERWISE) +* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +* POSSIBILITY OF SUCH DAMAGE. +*/ +#include "AEEStdDef.h" +#include "AEEstd.h" +#include "AEEStdErr.h" +#include "std_dtoa.h" +#include "math.h" + +#define FAILED(b) ( (b) != AEE_SUCCESS ? TRUE : FALSE ) +#define CLEANUP_ON_ERROR(b,l) if( FAILED( b ) ) { goto l; } +#define FP_POW_10(n) fp_pow_10(n) + +/* Returns the number of leading zeroes in a uint32. + This is a naive implementation that uses binary search. This could be + replaced by an optimized inline assembly code. */ +static __inline uint32 std_dtoa_clz32(uint32 ulVal) +{ + if ((int)ulVal <= 0) + { + return (ulVal == 0) ? 32 : 0; + } + else + { + uint32 uRet = 28; + uint32 uTmp = 0; + uTmp = (ulVal > 0xFFFF) * 16; + ulVal >>= uTmp, uRet -= uTmp; + uTmp = (ulVal > 0xFF) * 8; + ulVal >>= uTmp, uRet -= uTmp; + uTmp = (ulVal > 0xF) * 4; + ulVal >>= uTmp, uRet -= uTmp; + return uRet + ((0x55AF >> (ulVal * 2)) & 3); + } +} + +/* Returns the number of leading zeroes in a uint64. */ +static __inline uint32 std_dtoa_clz64(uint64 ulVal) + +{ + uint32 ulCount = 0; + + if (!(ulVal >> 32)) + { + ulCount += 32; + } + else + { + ulVal >>= 32; + } + + return ulCount + std_dtoa_clz32((uint32)ulVal); +} + +double fp_pow_10(int nPow) +{ + double dRet = 1.0; + int nI = 0; + boolean bNegative = FALSE; + double aTablePos[] = { 0, 1e1, 1e2, 1e4, 1e8, 1e16, 1e32, 1e64, 1e128, + 1e256 + }; + double aTableNeg[] = { 0, 1e-1, 1e-2, 1e-4, 1e-8, 1e-16, 1e-32, 1e-64, 1e-128, + 1e-256 + }; + double *pTable = aTablePos; + int nTableSize = STD_ARRAY_SIZE(aTablePos); + + if (0 == nPow) + { + return 1.0; + } + + if (nPow < 0) + { + bNegative = TRUE; + nPow = -nPow; + pTable = aTableNeg; + nTableSize = STD_ARRAY_SIZE(aTableNeg); + } + + for (nI = 1; nPow && (nI < nTableSize); nI++) + { + if (nPow & 1) + { + dRet *= pTable[nI]; + } + + nPow >>= 1; + } + + if (nPow) + { + /* Overflow. Trying to compute a large power value. */ + uint64 ulInf = STD_DTOA_FP_POSITIVE_INF; + dRet = bNegative ? 0 : UINT64_TO_DOUBLE(ulInf); + } + + return dRet; +} + +/* Rounds dNumber to the specified precision nPrecision. + For example: + fp_round(2.34553, 3) = 2.346 + fp_round(2.34553, 4) = 2.3455 */ +double fp_round(double dNumber, int nPrecision) +{ + double dResult = dNumber; + double dRoundingFactor = FP_POW_10(-nPrecision) * 0.5; + + if (dNumber < 0) + { + dResult = dNumber - dRoundingFactor; + } + else + { + dResult = dNumber + dRoundingFactor; + } + + return dResult; +} + +/* Finds the integer part of the log_10( dNumber ). + The function assumes that dNumber != 0. */ +int fp_log_10(double dNumber) +{ + /* Absorb the negative sign */ + if (dNumber < 0) + { + dNumber = -dNumber; + } + + return (int)(floor(log10(dNumber))); +} + +/* Evaluates the input floating-point number dNumber to check for + following special cases: NaN, +/-Infinity. + The evaluation is based on the IEEE Standard 754 for Floating Point Numbers. */ +int fp_check_special_cases(double dNumber, FloatingPointType *pNumberType) +{ + int nError = AEE_SUCCESS; + FloatingPointType NumberType = FP_TYPE_UNKOWN; + uint64 ullValue = 0; + uint64 ullSign = 0; + int64 n64Exponent = 0; + uint64 ullMantissa = 0; + + ullValue = DOUBLE_TO_UINT64(dNumber); + + /* Extract the sign, exponent and mantissa */ + ullSign = FP_SIGN(ullValue); + n64Exponent = FP_EXPONENT_BIASED(ullValue); + ullMantissa = FP_MANTISSA_DENORM(ullValue); + + /* Rules for special cases are listed below: + For Infinity, the following needs to be true: + 1. Exponent should have all bits set to 1. + 2. Mantissa should have all bits set to 0. + + For NaN, the following needs to be true: + 1. Exponent should have all bits set to 1. + 2. Mantissa should be non-zero. + Note that we do not differentiate between QNaNs and SNaNs. */ + + if (STD_DTOA_DP_INFINITY_EXPONENT_ID == n64Exponent) + { + if (0 == ullMantissa) + { + /* Inifinity. */ + if (ullSign) + { + NumberType = FP_TYPE_NEGATIVE_INF; + } + else + { + NumberType = FP_TYPE_POSITIVE_INF; + } + } + else + { + /* NaN */ + NumberType = FP_TYPE_NAN; + } + } + else + { + /* A normal number */ + NumberType = FP_TYPE_GENERAL; + } + + /* Set the output value */ + *pNumberType = NumberType; + + return nError; +} + +int std_dtoa_decimal(double dNumber, int nPrecision, + char acIntegerPart[ STD_DTOA_FORMAT_INTEGER_SIZE ], + char acFractionPart[ STD_DTOA_FORMAT_FRACTION_SIZE ]) +{ + int nError = AEE_SUCCESS; + boolean bNegativeNumber = FALSE; + double dIntegerPart = 0.0; + double dFractionPart = 0.0; + double dTempIp = 0.0; + double dTempFp = 0.0; + int nMaxIntDigs = STD_DTOA_FORMAT_INTEGER_SIZE; + uint32 ulI = 0; + int nIntStartPos = 0; + + /* Optimization: Special case an input of 0. */ + if (0.0 == dNumber) + { + acIntegerPart[0] = '0'; + acIntegerPart[1] = '\0'; + + for (ulI = 0; (ulI < STD_DTOA_FORMAT_FRACTION_SIZE - 1) && (nPrecision > 0); + ulI++, nPrecision--) + { + acFractionPart[ulI] = '0'; + } + acFractionPart[ ulI ] = '\0'; + + goto bail; + } + + /* Absorb the negative sign */ + if (dNumber < 0) + { + acIntegerPart[0] = '-'; + nIntStartPos = 1; + dNumber = -dNumber; + bNegativeNumber = TRUE; + } + + /* Split the input number into it's integer and fraction parts. */ + dFractionPart = modf(dNumber, &dIntegerPart); + + /* First up, convert the integer part */ + if (0.0 == dIntegerPart) + { + acIntegerPart[ nIntStartPos ] = '0'; + } + else + { + double dRoundingConst = FP_POW_10(-STD_DTOA_PRECISION_ROUNDING_VALUE); + int nIntDigs = 0; + int nI = 0; + + /* Compute the number of digits in the integer part of the number*/ + nIntDigs = fp_log_10(dIntegerPart) + 1; + + /* For negative numbers, a '-' sign has already been written. */ + if (TRUE == bNegativeNumber) + { + nIntDigs++; + } + + /* Check for overflow */ + if (nIntDigs >= nMaxIntDigs) + { + /* Overflow! + Note that currently, we return a simple AEE_EFAILED for all + errors. */ + nError = AEE_EFAILED; + goto bail; + } + + /* Null Terminate the string */ + acIntegerPart[ nIntDigs ] = '\0'; + + for (nI = nIntDigs - 1; nI >= nIntStartPos; nI--) + { + dIntegerPart = dIntegerPart / 10.0; + dTempFp = modf(dIntegerPart, &dTempIp); + + /* Round it to the a specific precision */ + dTempFp = dTempFp + dRoundingConst; + + /* Convert the digit to a character */ + acIntegerPart[ nI ] = (int)(dTempFp * 10) + '0'; + if (!MY_ISDIGIT(acIntegerPart[ nI ])) + { + /* Overflow! + Note that currently, we return a simple AEE_EFAILED for all + errors. */ + nError = AEE_EFAILED; + goto bail; + } + dIntegerPart = dTempIp; + } + } + + /* Just a double check for integrity sake. This should ideally never happen. + Out of bounds scenario. That is, the integer part of the input number is + too large. */ + if (dIntegerPart != 0.0) + { + /* Note that currently, we return a simple AEE_EFAILED for all + errors. */ + nError = AEE_EFAILED; + goto bail; + } + + /* Now, convert the fraction part */ + for (ulI = 0; (nPrecision > 0) && (ulI < STD_DTOA_FORMAT_FRACTION_SIZE - 1); + nPrecision--, ulI++) + { + if (0.0 == dFractionPart) + { + acFractionPart[ ulI ] = '0'; + } + else + { + double dRoundingValue = FP_POW_10(-(nPrecision + + STD_DTOA_PRECISION_ROUNDING_VALUE)); + acFractionPart[ ulI ] = (int)((dFractionPart + dRoundingValue) * 10.0) + '0'; + if (!MY_ISDIGIT(acFractionPart[ ulI ])) + { + /* Overflow! + Note that currently, we return a simple AEE_EFAILED for all + errors. */ + nError = AEE_EFAILED; + goto bail; + } + + dFractionPart = (dFractionPart * 10.0) - + (int)((dFractionPart + FP_POW_10(-nPrecision - 6)) * 10.0); + } + } + + +bail: + + return nError; +} + +int std_dtoa_hex(double dNumber, int nPrecision, char cFormat, + char acIntegerPart[ STD_DTOA_FORMAT_INTEGER_SIZE ], + char acFractionPart[ STD_DTOA_FORMAT_FRACTION_SIZE ], + int *pnExponent) +{ + int nError = AEE_SUCCESS; + uint64 ullMantissa = 0; + uint64 ullSign = 0; + int64 n64Exponent = 0; + static const char HexDigitsU[] = "0123456789ABCDEF"; + static const char HexDigitsL[] = "0123456789abcde"; + boolean bFirstDigit = TRUE; + int nI = 0; + int nF = 0; + uint64 ullValue = DOUBLE_TO_UINT64(dNumber); + int nManShift = 0; + const char *pcDigitArray = (cFormat == 'A') ? HexDigitsU : HexDigitsL; + boolean bPrecisionSpecified = TRUE; + + /* If no precision is specified, then set the precision to be fairly + large. */ + if (nPrecision < 0) + { + nPrecision = STD_DTOA_FORMAT_FRACTION_SIZE; + bPrecisionSpecified = FALSE; + } + else + { + bPrecisionSpecified = TRUE; + } + + /* Extract the sign, exponent and mantissa */ + ullSign = FP_SIGN(ullValue); + n64Exponent = FP_EXPONENT(ullValue); + ullMantissa = FP_MANTISSA(ullValue); + + /* Write out the sign */ + if (ullSign) + { + acIntegerPart[ nI++ ] = '-'; + } + + /* Optimization: Special case an input of 0 */ + if (0.0 == dNumber) + { + acIntegerPart[0] = '0'; + acIntegerPart[1] = '\0'; + + for (nF = 0; (nF < STD_DTOA_FORMAT_FRACTION_SIZE - 1) && (nPrecision > 0); + nF++, nPrecision--) + { + acFractionPart[nF] = '0'; + } + acFractionPart[nF] = '\0'; + + goto bail; + } + + /* The mantissa is in lower 53 bits (52 bits + an implicit 1). + If we are dealing with a denormalized number, then the implicit 1 + is absent. The above macros would have then set that bit to 0. + Shift the mantisaa on to the highest bits. */ + + if (0 == (n64Exponent + STD_DTOA_DP_EXPONENT_BIAS)) + { + /* DENORMALIZED NUMBER. + A denormalized number is of the form: + 0.bbb...bbb x 2^Exponent + Shift the mantissa to the higher bits while discarding the leading 0 */ + ullMantissa <<= 12; + + /* Lets update the exponent so as to make sure that the first hex value + in the mantissa is non-zero, i.e., at least one of the first 4 bits is + non-zero. */ + nManShift = std_dtoa_clz64(ullMantissa) - 3; + if (nManShift > 0) + { + ullMantissa <<= nManShift; + n64Exponent -= nManShift; + } + } + else + { + /* NORMALIZED NUMBER. + A normalized number has the following form: + 1.bbb...bbb x 2^Exponent + Shift the mantissa to the higher bits while retaining the leading 1 */ + ullMantissa <<= 11; + } + + /* Now, lets get the decimal point out of the picture by shifting the + exponent by 1. */ + n64Exponent++; + + /* Read the mantissa four bits at a time to form the hex output */ + for (nI = 0, nF = 0, bFirstDigit = TRUE; ullMantissa != 0; + ullMantissa <<= 4) + { + uint64 ulHexVal = ullMantissa & 0xF000000000000000uLL; + ulHexVal >>= 60; + if (bFirstDigit) + { + // Write to the integral part of the number + acIntegerPart[ nI++ ] = pcDigitArray[ulHexVal]; + bFirstDigit = FALSE; + } + else if (nF < nPrecision) + { + // Write to the fractional part of the number + acFractionPart[ nF++ ] = pcDigitArray[ulHexVal]; + } + } + + /* Pad the fraction with trailing zeroes upto the specified precision */ + for (; bPrecisionSpecified && (nF < nPrecision); nF++) + { + acFractionPart[ nF ] = '0'; + } + + /* Now the output is of the form; + h.hhh x 2^Exponent + where h is a non-zero hexadecimal number. + But we were dealing with a binary fraction 0.bbb...bbb x 2^Exponent. + Therefore, we need to subtract 4 from the exponent (since the shift + was to the base 16 and the exponent is to the base 2). */ + n64Exponent -= 4; + *pnExponent = (int)n64Exponent; + +bail: + return nError; +} + |