How do I convert double to string using only math.h?

Question

I am trying to convert a double to a string in a native NT application, i.e. an application that only depends on ntdll.dll. Unfortunately, ntdll's version of vsnprintf does not support %f et al., forcing me to implement the conversion on my own.

The aforementioned ntdll.dll exports only a few of the math.h functions (floor, ceil, log, pow, ...). However, I am reasonably sure that I can implement any of the unavailable math.h functions if necessary.

There is an implementation of floating point conversion in GNU's libc, but the code is extremely dense and difficult to comprehent (the GNU indentation style does not help here).

I've already implemented the conversion by normalizing the number (i.e. multiplying/dividing the number by 10 until it's in the interval [1, 10)) and then generating each digit by cutting the integral part off with modf and multiplying the fractional part by 10. This works, but there is a loss of precision (only the first 15 digits are correct). The loss of precision is, of course, inherent to the algorithm.

I'd settle with 17 digits, but an algorithm that would be able to generate an arbitrary number of digits correctly would be preferred.

Could you please suggest an algorithm or point me to a good resource?

Answer 1

Double-precision numbers do not have more than 15 significant (decimal) figures of precision. There is absolutely no way you can get "an arbitrary number of digits correctly"; doubles are not bignums.

Since you say you're happy with 17 significant figures, use long double; on Windows, I think, that will give you 19 significant figures.

Answer 2

I've thought about this a bit more. You lose precision because you normalize by multiplying by some power of 10 (you chose [1,10) rather than [0,1), but that's a minor detail). If you did so with a power of 2, you'd lose no precision, but then you'd get "decimal digits"*2^e; you could implement bcd arithmetic and compute the product yourself, but that doesn't sound like fun.

I'm pretty confident that you could split the double g=m*2^e into two parts: h=floor(g*10^k) and i=modf(g*10^k) for some k, and then separately convert to decimal digits and then stitch them together, but how about a simpler approach: use "long double" (80 bits, but I've heard that Visual C++ may not support it?) with your current approach and stop after 17 digits.

_gcvt should do it (edit - it's not in ntdll.dll, it's in some msvcrt*.dll?)

As for decimal digits of precision, IEEE binary64 has 52 binary digits. 52*log10(2)=15.65... (edit: as you pointed out, to round trip, you need more than 16 digits)

Answer 3

Does vsnprintf supports I64?

double x = SOME_VAL; // allowed to be from -1.e18 to 1.e18
bool sign = (SOME_VAL < 0);
if ( sign ) x = -x;
__int64 i = static_cast<__int64>( x );
double xm = x - static_cast<double>( i );
__int64 w = static_cast<__int64>( xm*pow(10.0, DIGITS_VAL) ); // DIGITS_VAL indicates how many digits after the decimal point you want to get

char out[100];
vsnprintf( out, sizeof out, "%s%I64.%I64", (sign?"-":""), i, w );

Another option is to try to find implementation of gcvt.

Answer 4

// --------------------------------------------------------------------------
// return number of decimal-digits of a given unsigned-integer
// N is BYTE/WORD/UINT/ULONGLONG
template <class N> inline BYTE GetUnsignedDecDigits(const N n)
{
    if ((N)-1 < (N)1)              // if type of N is floating-point / signed-integer
        if (::IsDebuggerPresent())
        {
            ::OutputDebugString(_T("GetUnsignedDecDigits: Incorrect type passed\n"));
            ::DebugBreak();
        }

    const BYTE anMaxDigits[]= {3, 5, 8, 10, 13, 15, 17, 20};
    const BYTE nMaxDigits   = anMaxDigits[sizeof(N)-1];

    BYTE nDigits=  1;
    N    nRoof  = 10;

    while ((n >= nRoof) && (nDigits<nMaxDigits))
    {
        nDigits++;
        nRoof*= 10;
    }

    return nDigits;
}
// --------------------------------------------------------------------------
// convert floating-point value to NULL-terminated string represention
// T is char or wchar_t
template <class T> T* DoubleToStr(double f       ,  // [i  ]
                                  T*     pczStr  ,  // [i/o]
                                  int    nDigitsI,  // [i  ] digits of integer    part including sign / <1: auto
                                  int    nDigitsF ) // [i  ] digits of fractional part                / <0: auto
{
    switch (_fpclass(f))
    {
    case _FPCLASS_SNAN:
    case _FPCLASS_QNAN: str_cpy(pczStr, 5, _T("NaN" )); return pczStr;
    case _FPCLASS_NINF: str_cpy(pczStr, 5, _T("-INF")); return pczStr;
    case _FPCLASS_PINF: str_cpy(pczStr, 5, _T("+INF")); return pczStr;
    }

    if (nDigitsI> 18) nDigitsI= 18;  if (nDigitsI< 1) nDigitsI= -1;
    if (nDigitsF> 18) nDigitsF= 18;  if (nDigitsF< 0) nDigitsF= -1;

    bool bNeg= (f<0);
    if (f<0)
        f= -f;

    int nE= 0;                                  // exponent (displayed if != 0)

    if ( ((-1 == nDigitsI) && (f >= 1e18              )) ||   // large value: switch to scientific representation
         ((-1 != nDigitsI) && (f >= pow(10., nDigitsI)))    )
    {
       nE= (int)log10(f);
       f/= (double)pow(10., nE);

       if (-1 != nDigitsF)
           nDigitsF= max(nDigitsF, nDigitsI+nDigitsF-(bNeg?2:1)-4);

       nDigitsI= (bNeg?2:1);
    }
    else if (f>0)
    if ((-1 == nDigitsF) && (f <= 1e-10))       // small value: switch to scientific representation
    {
        nE= (int)log10(f)-1;
        f/= (double)pow(10., nE);

       if (-1 != nDigitsF)
           nDigitsF= max(nDigitsF, nDigitsI+nDigitsF-(bNeg?2:1)-4);

        nDigitsI= (bNeg?2:1);
    }

    double fI;
    double fF= modf(f, &fI);                    // fI: integer part, fF: fractional part

    if (-1 == nDigitsF)                         // figure out number of meaningfull digits in fF
    {
        double fG, fGI, fGF;
        do
        {
            nDigitsF++;
            fG = fF*pow(10., nDigitsF);
            fGF= modf(fG, &fGI);
        }
        while (fGF > 1e-10);
    }

    ULONGLONG uI= Round<ULONGLONG>(fI                    );
    ULONGLONG uF= Round<ULONGLONG>(fF*afPower10[nDigitsF]);

    if (uF)
        if (GetUnsignedDecDigits(uF) > nDigitsF)    // X.99999 was rounded to X+1
        {
            uF= 0;
            uI++;

            if (nE)
            {
                uI/= 10;
                nE++;
            }
        }

    BYTE nRealDigitsI= GetUnsignedDecDigits(uI);
    if (bNeg)
        nRealDigitsI++;

    int nPads= 0;

    if (-1 != nDigitsI)
    {
        nPads= nDigitsI-nRealDigitsI;

        for (int i= nPads-1; i>=0; i--)         // leading spaces
            pczStr[i]= _T(' ');
    }

    if (bNeg)                                   // minus sign
    {
        pczStr[nPads]= _T('-');
        nRealDigitsI--;
        nPads++;
    }

    for (int j= nRealDigitsI-1; j>=0; j--)      // digits of integer    part
    {
        pczStr[nPads+j]= (BYTE)(uI%10) + _T('0');
        uI /= 10;
    }

    nPads+= nRealDigitsI;

    if (nDigitsF)
    {
        pczStr[nPads++]= _T('.');               // decimal point
        for (int k= nDigitsF-1; k>=0; k--)      // digits of fractional part
        {
            pczStr[nPads+k]= (BYTE)(uF%10)+ _T('0');
            uF /= 10;
        }
    }

    nPads+= nDigitsF;

    if (nE)
    {
        pczStr[nPads++]= _T('e');               // exponent sign

        if (nE<0)
        {
            pczStr[nPads++]= _T('-');
            nE= -nE;
        }
        else
            pczStr[nPads++]= _T('+');

        for (int l= 2; l>=0; l--)               // digits of exponent
        {
            pczStr[nPads+l]= (BYTE)(nE%10) + _T('0');
            nE /= 10;
        }

        pczStr[nPads+3]= 0;
    }
    else
        pczStr[nPads]= 0;

    return pczStr;
}

Answer 5

Have you looked at the uClibc implementation of printf?

Answer 6

After a lot of research, I found a paper titled Printing Floating-Point Numbers Quickly and Accurately. It uses exact rational arithmetic to avoid precision loss. It cites a little older paper: How to Print Floating-Point Numbers Accurately, which however seems to require ACM subscription to access.

Since the former paper was reprinted in 2006, I am inclined to believe that it is still current. The exact rational arithmetic (which requires dynamic allocation) seems to be a necessary evil.