log-space, quasipolynomial-time floating point base converter

ftoa.c

#define __STDC_WANT_IEC_60559_BFP_EXT__
#define __STDC_WANT_IEC_60559_TYPES_EXT__
#include <stdlib.h>
#include <stdio.h>
#include <string.h>

int memo_max_exp;
#ifdef CHEAT
static char **memo[6];
#endif

/* Calculate floor(base**exponent / 10**position) mod 10
   base should be 2 or 5, exponent should be nonnegative */
static int gen_power_digit(int base, int exponent, int position) {
    int d = 0;
    int rv = 0;
    int cur_exp;
    /* largest input exponent is a bit over 2^15, so if d=15 we'll never get
       past the base cases; pos[15] needs to exist for pos[d+1] */
    /* these types are slightly larger than they need to be, e.g. level[8]
       and higher would fit in a char */
    int accumulator[15];
    short level[15], pos[16];
    char flag[15];

    pos[0] = position;
nest:
    cur_exp = exponent >> d;

    if (pos[d] < 0) {
        rv = 0;
        goto out;
    }
    if (cur_exp == 0) {
        rv = pos[d] == 0;
        goto out;
    }
    if (cur_exp == 1) {
        rv = (pos[d] == 0) * base;
        goto out;
    }
    /* optimization opportunity: multiply exponent by an upper bound of
       log10(base) */
    if (pos[d] > cur_exp) {
        rv = 0;
        goto out;
    }
#ifdef CHEAT
    if (cur_exp < memo_max_exp) {
        if (!memo[base]) {
            memo[base] = calloc(memo_max_exp, sizeof(char*));
        }
        if (!memo[base][cur_exp]) {
            memo[base][cur_exp] = calloc(exponent + 1, 1);
        }
        if (memo[base][cur_exp][pos[d]]) {
            rv = memo[base][cur_exp][pos[d]] - '0';
            goto out;
        }
    }
#endif

    accumulator[d] = 0;
    for (level[d] = 0; level[d] <= cur_exp; level[d]++) {
        for (pos[d+1] = 0; pos[d+1] <= level[d]; pos[d+1]++) {
            flag[d] = 0;
            for (;;) {
                d++;
                goto nest;
unnest:
                d--;
                cur_exp = exponent >> d;
                pos[d+1] = level[d] - pos[d+1];
                if (flag[d]) break;
                flag[d] = 1 + rv;
            }
            accumulator[d] += rv * (flag[d] - 1) * ((cur_exp & 1) ? base : 1);
        }

#ifdef CHEAT
        if (cur_exp < memo_max_exp) {
            memo[base][cur_exp][level[d]] = '0' + (accumulator[d] % 10);
        } else {
#endif
            if (level[d] == pos[d]) {
                rv = accumulator[d] % 10;
                goto out;
            }
#ifdef CHEAT
        }
#endif

        accumulator[d] /= 10;
    }

#ifdef CHEAT
    /* reachable only if memoizing */
    rv = memo[base][cur_exp][pos[d]] - '0';
#endif
out:
    if (d) goto unnest;
    return rv;
}

/* exp=112 gives an integer */
static inline int get_float_exp(unsigned char *fp) {
    int top = fp[14] + (fp[15] << 8);
    return (top & 0x7fff) - 16383;
}

static inline int min_pos(unsigned char *fp) {
    int exp = get_float_exp(fp);
    return exp < 112 ? exp - 112 : 0;
}

static void get_float_predig(int *acc, int level, unsigned char *fp) {
    int pos2;
    int exp = get_float_exp(fp);
    for (pos2 = 0; pos2 <= 112; pos2++) {
        /* considering bit pos2 in mantissa, 2**(exp+pos2-112) */
        /* skip if 0 */
        if (pos2 == 112 && exp == -16383) continue; /* implied bit */
        if (pos2 < 112 && !((fp[pos2 >> 3] >> (pos2 & 7)) & 1))
            continue; /* explicit bit */
        *acc += gen_power_digit(
            exp + pos2 >= 112 ? 2 : 5,
            exp + pos2 >= 112 ? exp + pos2 - 112 : 112 - exp - pos2,
            exp + pos2 >= 112 ? level : level + (112 - exp - pos2));
    }
}

static inline int max_pos(unsigned char *fp) {
    int exp = get_float_exp(fp);
    return exp > 0 ? exp : 0;
}

int my_strfromf128(char *out, size_t n, char *format, unsigned char *fp) {
    int field_width = 0;
    int precision = -1;
    size_t written = 0;
    enum { f_zeros=64, f_left=128, f_sign=8, f_space=16,
        f_radix=32, f_e=1, f_f=2, f_g=4 };
    enum { f_suppress=1, f_upcase=2, f_thou=4 };
    int flags = 0, flags2 = 0;

    for (; *format; format++) {
        switch (*format) {
            case '%': continue;
            case '#': flags |= f_radix; continue;
            case '+': flags |= f_sign; continue;
            case '-': flags |= f_left; continue;
            case ' ': flags |= f_space; continue;
            case '\'': flags2 |= f_thou; continue;
            case '.': precision = 0; continue;
            case 'e': case 'f': case 'g':
                flags |= 1 << (*format - 'e'); format = "!";
                continue;
            case 'E': case 'F': case 'G':
                flags |= 1 << (*format - 'E'); flags2 |= f_upcase;
                format = "!";
                continue;
            case '0':
                if (!field_width && precision < 0) {
                    flags |= f_zeros; continue;
                }
                /*fallthrough*/
            default:
                if (*format >= '0' && *format <= '9') {
                    if (precision >= 0) {
                        precision = precision * 10 + *format - '0';
                    } else {
                        field_width = field_width * 10 + *format - '0';
                    }
                }
                continue;
        }
    }

#define OUT(ch) do { if (written < n) out[written] = ch; written++; } while(0)
    if (get_float_exp(fp) == 16384) {
        /* inf/nan handling */
        int is_nan = 0;
        for (int i = 0; i < 14; i++) is_nan = fp[i] ? 1 : is_nan;
        int pad = field_width - 3 -
            ((fp[15] & 128) || flags & f_sign || (flags & f_space));
        if (!(flags & f_left)) while (pad-- > 0) OUT(' ');
        if (fp[15] & 128) OUT('-');
        else if (flags & f_sign) OUT('+');
        else if (flags & f_space) OUT(' ');
        for (int i = 0; i < 3; i++)
            OUT("infnanINFNAN"[i + 3*is_nan + 6*!!(flags2 & f_upcase)]);
        if (flags & f_left) while (pad-- > 0) OUT(' ');
        OUT(0);
        return (int)written;
    }

    if (precision == -1) precision = 6;
    if ((flags & f_g) && precision == 0) precision = 1;

    /* find true top */
    short top = -20000, level;
    int acc = 0;
    /*int rstate;*/
    int pass;

    /* pass 0: finds true exponent, rounding disabled, skip for f */
    /* pass 1: find rounded exponent for e/f/g, changes top */
    /* after pass 1: resolve g to e/f */
    /* pass 2: reduce precision to remove trailing zeros if orig g */
    /* after pass 2: print length is known, emit pad */
    /* pass 3+: print digits; the pass number counts digits */

    for (pass = 0; ; pass++) {
        if (pass == 0 && (flags & f_f)) continue;
        if (pass == 2 && !(flags2 & f_suppress)) continue;

        if (pass == 3) {
            /* one-time convert field_width to the number of padding
               characters needed */

            if (flags & f_sign || (flags & f_space) || fp[15] & 128) field_width--;
            if ((flags & f_f) && top > 0) {
                field_width -= top + 1;
                if (flags2 & f_thou) field_width -= top / 3;
            } else {
                field_width--;
            }
            if ((flags & f_radix) || precision > 0) field_width--;
            field_width -= precision;
            if (flags & f_e) {
                field_width -= 2;
                field_width -= (top < 100 && top > -100) ? 2 :
                    (top < 1000 && top > -1000) ? 3 :
                    (top < 10000 && top > -10000) ? 4 : 5;
            }
        }

        /* output left-padding (if not leftjust) and sign immediately before
           first digit */
        if (pass == 3) {
            if (!(flags & f_left) && !(flags & f_zeros)) {
                while (field_width-- > 0) OUT(' ');
            }
            if (fp[15] & 128) OUT('-');
            else if (flags & f_sign) OUT('+');
            else if (flags & f_space) OUT(' ');
            if (!(flags & f_left) && (flags & f_zeros)) {
                while (field_width-- > 0) OUT('0');
            }
        }

        /* we might be done? */
        if (pass == 3 + ((flags & f_f) && top > 0 ? top + 1 : 1) +
                precision) {
            if (flags & f_e) {
                OUT((flags2 & f_upcase) ? 'E' : 'e');
                OUT(top < 0 ? '-' : '+');
                if (top < 0) top = -top;
                if (top >= 10000) OUT('0' + top / 10000);
                if (top >= 1000) OUT('0' + (top / 1000) % 10);
                if (top >= 100) OUT('0' + (top / 100) % 10);
                OUT('0' + (top / 10) % 10);
                OUT('0' + top % 10);
            }
            if (flags & f_left) {
                while (field_width-- > 0) OUT(' ');
            }
            OUT(0);
            return (int)written;
        }

        /* rstate = 0; */
        level = 0;
        if (pass >= 3) {
            level = ((flags & f_f) && top < 0 ? 0 : top) - (pass - 3);
        }
        if (level > min_pos(fp)) level = min_pos(fp);
        for (; level <= max_pos(fp) || acc; level++) {
            get_float_predig(&acc, level, fp);
            /* actual rounding happens here, while looking at the first
               hidden digit */
            if (pass && acc % 10 >= 5) {
                if ((flags & f_g) && top - level == precision) acc += 10;
                if ((flags & f_f) && level == -1 - precision) acc += 10;
                if ((flags & f_e) && top - level == precision + 1)
                    acc += 10;
            }
            if (pass == 2 && precision > 0 && (acc % 10) == 0) {
                if (level == ((flags & f_e) ? top : 0) - precision) {
                    /* remove a trailing 0. this shouldn't change rounding */
                    /* may happen more than once in pass 2 */
                    precision--;
                }
            }
            /* sticky bit - needed for configurable rounding */
            /* if (acc % 10) rstate = 1; */
            if (acc && level > top) {
                /* carry out makes the top higher */
                /* note, this moves the rounding point in subsequent passes!
                   this should be safe because we've just established the (new)
                   first hidden digit is a 9 */
                top = level;
            }
            if (pass >= 3 && level == ((flags & f_f) && top < 0 ? 0 : top)
                    - (pass - 3)) {
                OUT('0' + acc % 10);
                int rlevel = level - ((flags & f_e) ? top : 0);
                if (rlevel == 0 && ((flags & f_radix) || precision > 0)) {
                    OUT('.');
                }
                if ((flags2 & f_thou) && rlevel > 0 && (rlevel % 3) == 0) {
                    OUT(',');
                }
            }
            acc /= 10;
        }
        /* fix top for 0 */
        if (top == -20000) top = 0;
        /* top is the e-mode exponent if pass >= 1 */

        if (pass == 1 && (flags & f_g)) {
            if (!(flags & f_radix)) flags2 |= f_suppress;
            /* both of these have the same maximum number of digits,
               so we can reuse the rounding calculation */
            if (precision > top && top >= -4) {
                flags = (flags & ~f_g) | f_f;
                precision -= top + 1;
            } else {
                flags = (flags & ~f_g) | f_e;
                precision--;
            }
        }
    }
}

#ifdef TEST
int main(int argc, char **argv) {
    if (argc != 4) {
        printf("usage: ./a.out (0-99999) (format string) (float)\n");
        return 1;
    }
    memo_max_exp = atoi(argv[1]);
    _Float128 fpval = strtof128(argv[3], 0);
    int len = my_strfromf128(0, 0, argv[2], (unsigned char*)&fpval);
    char tmp[len + 1];
    my_strfromf128(tmp, len+1, argv[2], (unsigned char*)&fpval);

    printf("DUT <%s>\n", tmp);

    if (strlen(argv[2]) > 2) {
        int len2 = snprintf(0, 0, argv[2], (double)fpval);
        char tmp2[len2 + 1];
        snprintf(tmp2, len2 + 1, argv[2], (double)fpval);

        printf("REF <%s>\n", tmp2);
    } else {
        int len2 = strfromf128(0, 0, argv[2], fpval);
        char tmp2[len2 + 1];
        strfromf128(tmp2, len2 + 1, argv[2], fpval);

        printf("REF <%s>\n", tmp2);
    }
    return 0;
}
#endif