`atof` revisited

Question

In an answer to this question I mentioned best effort. Here I try to explain what I meant. Please keep in mind that the implementation is intentionally incomplete (missing features such as NAN are trivial to add) and focuses specifically on the numerical stability.

It was surprisingly hard to get it right, but it paid off: the results match stdlib implementation even for the most frivolous inputs I could come up with.

However, normalize is the messiest piece of code I have ever written. It computes three results - that alone is enough to raise some brows. And memmove looks particularly out of place.

All suggestions are welcome.

#include <stdlib.h>
#include <stdbool.h>
#include <string.h>
#include <math.h>

static char * drop_leading_zeroes(char * start)
{
    return start + strspn(start, "0");
}

static char * digit_span(char * start)
{
    return start + strspn(start, "0123456789");
}

static int normalize(char ** start, char ** end)
{
    int shift = 0;
    *start = drop_leading_zeroes(*start);
    if (isdigit(**start)) {
        *end = digit_span(*start);
        shift = *end - *start;
        if (**end == '.') {
            *end = digit_span(*end + 1);
            memmove(*start + 1, *start, shift);
            *start += 1;
        }
    } else if (**start == '.') {
        *start += 1;
        *end = drop_leading_zeroes(*start);
        shift = *start - *end;
        *start = *end;
        *end = digit_span(*end);
    } else {
        *end = *start;
    }
    return shift;
}

static double compute_mantissa(char * start, char * end)
{
    double result = 0.0;
    while (end != start) {
        result = (result + (*--end - '0')) / 10;
    }
    return result;
}

double my_atof(char * s)
{
    bool minus = false;
    switch(*s) {
        case '-': minus = true;
        case '+': ++s; break;
    }

    char * end;
    int exponent = normalize(&s, &end);

    double mantissa = compute_mantissa(s, end);
    if (minus) {
        mantissa = -mantissa;
    }

    if (tolower(*end) == 'e') {
        exponent += atoi(end + 1);
    }

    return mantissa * pow(10, exponent);
}

Answer 1

Modifying the input string

The real atof() does not modify the input string and it seems wrong for yours to do so. It could lead to all sorts of unexpected behavior. For example, your program could crash if you pass in a string literal that is in a read-only section:

// Segmentation violation!
double val = my_atof("55.5");

Or you might do this and get a surprising result:

double val1 = my_atof(str);
double val2 = my_atof(str);

// val1 != val2 because the string mutated

You can make a few small changes to make your atof() work without modifying the input. First remove the lines that modify the string:

        memmove(*start + 1, *start, shift);
        *start += 1;

The only purpose of those lines is to erase the '.' separating the whole part from the fractional part. Then modify compute_mantissa() to skip any '.' characters:

static double compute_mantissa(const char * start, const char * end)
{
    double result = 0.0;
    while (end != start) {
        char c = *--end;
        if (c == '.')
            continue;
        result = (result + (c - '0')) / 10;
    }
    return result;
}

The last step is to change every char * to a const char *, now that you aren't going to modify any strings. You can see I already did that to compute_mantissa() above.

Reduce levels of indirection

In normalize(), you operate on pointers to pointers throughout the function. I find that if you use temporary variables, you can reduce the levels of indirection and the code will be easier to read. So for example:

static int normalize(const char ** pStart, const char ** pEnd)
{
    int shift = 0;
    const char * start = *pStart;
    const char * end;

    start = drop_leading_zeroes(start);
    if (isdigit(*start)) {
        end = digit_span(start);
        shift = end - start;
        if (*end == '.') {
            end = digit_span(end + 1);
        }
    } else if (*start == '.') {
        start += 1;
        end = drop_leading_zeroes(start);
        shift = start - end;
        start = end;
        end = digit_span(end);
    } else {
        end = start;
    }
    *pStart = start;
    *pEnd   = end;
    return shift;
}

Answer 2

1. Any math related function deserves an assessment of its correctness. A typical metric is in ULP. An OK atof() should have a worst case 1.0 ULP or less and an average certainly less than 0.5 ULP. Yet this code simple asserts "results match stdlib implementation even for the most frivolous inputs I could come up with" without providing sample inputs nor statistics. It was good a comparison was done, yet a quantitative assessment is deserved. To rectify this takes work, possible more than writing the code.

2. The modification of the input string breaks the contract of the standard double atof(const char *nptr). That alone greatly degrades the code's value. Similar noted by @JS1. To correct this, form "matissa" as in #3, skipping the optional decimal point. No need to change the source string.

3. compute_mantissa() incurs round-off error with each iteration as dividing by 10 is rarely exact. This repeated error accumulation is easily avoided. Form result by multiplying from start to end. No error is expected with less than DBL_DIG (e. g. 15) digits. The power-of-ten of end-start can be accounted for in the exponent calculation of compute_mantissa(s, end)

4. mantissa * pow(10, exponent) suffers when pow(10, exponent) is outside double normal range even when the mathematical product is within range. Simple use (mantissa * pow(5, exponent)) * pow(2, exponent) as suggested in this late review of the earlier code.

5. Good use of static for helper functions and interesting use of strspn().

6. Minor: isdigit() is not defined for all negative inputs like those that may come from char ** start ... isdigit(**start). Use isdigit((unsigned char) (**start))

7. Minor: Sign detection is OK, but simpler code exists. Example:

   char *sign = s;
   if (*s == '-' || *s == '+') s++;
   ...
   if (*sign == '-') {
     mantissa = -mantissa;
   }
   

8. If you continue with case '-': minus = true;, recommend a comment after to indicate the intentional drop though: case '-': minus = true; // drop though as at first glance, it looks like a bug missing a break before the next case.

9. Very minor: Pedantically int shift is too narrow a type for *end - *start which returns type ptrdiff_t. Code could use shift = (int) (*end - *start); to quiet warnings.

10. Code is using a hard-coded decimal point '.'. The decimal point is locale sensitive. Could use char dp = localeconv()->decimal_point;

Overall: Descent implementation but breaking the const function signature is a non-no.