2026-08-31

Printing a double without printf

tpl2c's runtime prints numbers without printf. A typed tag such as <?=%.2f price ?> goes through the runtime's own formatter, which handles every printf conversion and must print the same text as glibc.

I rebuilt every version for this article and timed it on one x86-64 machine under Debian 13, with gcc 14.2 at -O2 unless noted, against glibc 2.41 on the same machine. Times are in nanoseconds per call, the best of five runs of at least 0.2 seconds on one core. Two runs of the same build can differ by up to 20%.

Why not printf

For one value, printf is overkill. It is an interpreter for format strings: on every call it parses the format, walks the flags, the width, the precision and the length, picks the conversion, and fetches the argument through varargs. In a template, the type is known when the page is compiled. tpl2c checks the conversion at build time. <?=%d n ?> becomes a direct call to hc_put_int_esc, with no format string at all. <?=%.2f x ?> still goes through the runtime's own formatter, which reads its one conversion on each call.

Dave Plummer tells a similar story in Why I Deleted printf() from Windows COM in 1994!: as a Windows developer, he found a routine that COM called constantly and that was written very inefficiently, and took printf out of it.

Integers are easy. A double is a binary fraction, and %.15g of it means the decimal digits of its exact value, rounded half to even. The hard cases are the values that sit almost exactly halfway between two outputs.

What a double holds

A float has 32 bits and a double 64. Both have the same three fields, from the highest bit down: a sign, an exponent and a fraction.

        sign exponent fraction
float:  1    8        23 bits
double: 1    11       52 bits

For a normal number, the value is 1.fraction × 2^{exponent − bias}, where the bias is 127 for a float and 1023 for a double. With the leading 1 put in front of its 52 bits, the fraction is a 53-bit integer m, and the value is m * 2^e, with e the exponent field minus 1075. Here is 0.1, the double 0x3FB999999999999A:

0 | 01111111011 | 1001100110011001100110011001100110011001100110011010

The exponent field is 1019, so e is -56, and m is 7205759403792794. Halving m and adding 1 to e gives 3602879701896397 × 2^-55, the form used below. When the exponent field is 0, there is no leading 1: those are the subnormals, down to 5e-324. Over all doubles, e goes from -1074 to 971. A float goes through the same code, converted to a double first, as printf does.

Written in decimal, that value can be long: 1e-300 has more than 700 significant digits. They are kept in a big number: an array of unsigned long whose elements, called limbs, each hold a few decimal digits. With 8 digits per limb, each limb is one digit in base 10⁸. 123456789012 is stored as two limbs, the lowest first:

unsigned long a[2];

a[0] = 56789012; /* lowest */
a[1] = 1234;

Multiplying the big number by a small factor is done limb by limb with a carry, from a[0] up, as on paper.

For a longer introduction, Fabien Sanglard's Floating Point Visually Explained starts from the beginning, with pictures. Bartosz Ciechanowski's Exposing Floating Point goes further, and his float.exposed shows the bits of any value.

First: exact, and slow

The first version followed the method of Steele and White: compute the exact decimal value with a big number, then round. I kept m in limbs of 8 decimal digits and multiplied it by 25 or by 32, one pass at a time, until the exponent was gone. It passed every test from the start, but it was slow.

The idea fits in one example. 0.1 is 3602879701896397 × 2^-55, and dividing by 2⁵⁵ is the same as multiplying by 5⁵⁵ and moving the decimal point 55 places. So the exact value of 0.1 is 3602879701896397 × 5⁵⁵, with the point 55 places from the right:

0.1000000000000000055511151231257827021181583404541015625

%.15g keeps 15 of those digits and rounds. The big number does the multiplication one pass at a time. A pass multiplies every limb by a small factor m and carries the rest into the next limb:

#define HC_BASE 100000000UL /* 10^8 */

static size_t hc_mul(unsigned long *a,
                     size_t n, unsigned long m) {
    unsigned long c = 0;
    size_t i;

    for (i = 0; i < n; i++) {
        unsigned long t = a[i] * m + c;

        a[i] = t % HC_BASE;
        c = t / HC_BASE;
    }
    if (c != 0) {
        a[n++] = c;
    }
    return n;
}

A limb is below 10⁸, so a limb times 32, plus the carry, stays below 2³² and fits in a 32-bit unsigned long. The passes multiply by 32 until the power of 2 is done, or by 25 until the power of 5 is done:

for (j = e; j >= 5; j -= 5) {
    n = hc_mul(a, n, 32);
}
if (j > 0) {
    n = hc_mul(a, n, 1UL << j);
}
for (j = -e; j >= 2; j -= 2) {
    n = hc_mul(a, n, 25);
}
if (j > 0) {
    n = hc_mul(a, n, 5);
}

For the value 0.1, that is 27 passes by 25 and one by 5. For the value 1e-300, it is 525 passes, over a number that grows to 94 limbs.

The first version was two to four times slower than glibc on values such as 0.1, 3.14159 or 2.5e-7, about as fast on 1234.5678 or a price, and took about 70 microseconds for 1e-300, where the big number has hundreds of digits. It won on %f of huge values, where glibc itself takes 26 microseconds, and on whole numbers: the runtime prints a double such as 123456789.0 through its integer code, in 46 ns, where glibc takes 751.

Second: bigger limbs, then a window

Two changes made a huge improvement on the extreme values: from 70 microseconds to only 3.

The first optimization is bigger limbs. Where unsigned long has 64 bits, a limb can hold 9 digits instead of 8. The same 123456789012 becomes:

a[0] = 456789012; /* lowest */
a[1] = 123;

A limb below 10⁹ times 5¹⁴ or 2³³, plus the carry, still fits in 64 bits, so one pass can multiply by those instead of 25 or 32:

#define HC_BASE 1000000000UL /* 10^9 */
#define HC_M5   6103515625UL /* 5^14 */
#define HC_K2   33

for (j = e; j >= HC_K2; j -= HC_K2) {
    n = hc_mul(a, n, 1UL << HC_K2);
}
for (j = -e; j >= 14; j -= 14) {
    n = hc_mul(a, n, HC_M5);
}

The last pass multiplies by what is left, such as 5¹³. 5¹⁴ is bigger than a limb, so the carry at the end of a pass can now fill two limbs, and the if at the end of hc_mul became a loop:

for (; c != 0; c /= HC_BASE) {
    a[n++] = c % HC_BASE;
}

For the value 0.1, that is 4 passes instead of 28. For the value 1e-300, it is 75 passes over 84 limbs instead of 525 over 94, and it went from 69 to 12 microseconds. Where unsigned long has 32 bits, the code keeps 8 digits and 25 or 32.

The second optimization is to stop computing digits when they are not needed anymore. %.15g needs 15 digits, not the 700-odd of 1e-300, so the big number keeps only a window of its top limbs, plus a few spare ones, and drops the rest as it goes. Each drop makes the value a little smaller, but the error stays below one unit of the third lowest limb. If the digits after the rounding place are neither all zeros nor all nines, rounding the window gives the right answer; if they are, the exact computation runs instead. 1e-300 went from 12 microseconds to 3.

Dropping limbs only makes the value smaller, never larger, so the question is whether the missing part could change the rounding. d holds the window's digits, and d[keep] is the first digit dropped, the one that decides the rounding. The missing part can change that decision in two ways. It can carry into d[keep], which needs every digit after it to be 9. Or it can turn an exact half, a 5 followed by zeros, into slightly more than half, which needs them all to be 0. So the check looks between d[keep] and the last digit for one digit that is not 0 and one that is not 9:

static int hc_digits_sure(const char *d, size_t n, long keep) {
    size_t i = keep > 0 ? (size_t)keep : 0;
    int z = 0, nine = 0;

    for (i++; i + 1 < n; i++) {
        z |= d[i] != '0';
        nine |= d[i] != '9';
    }
    return z && nine;
}

Third: one multiplication

Most double numbers a program prints are prices, coordinates or ratios: a few digits and a moderate exponent. For those, the big number is overkill. The idea of Loitsch's Grisu is to compute the digits with machine integers, check whether the result could be wrong, and fall back to the exact method only when it could.

With a 64-bit long, when the digits to print fit in one long, they are m * 5^k shifted by a power of 2. One 64 by 64 bit product, kept in 128 bits, gives the digits, and the bits below the rounding place tell whether rounding them is safe. For scales beyond 10²⁷, a table of 26 powers of 5 rounded to 64 bits covers every exponent, and the rounding is accepted only when the lowest and the highest value the table's error allows round the same way.

Take %.15g of 0.1 again. Its 15 digits are 0.1 × 10¹⁵, rounded, which is 3602879701896397 × 5¹⁵ × 2^-40. The product with 5¹⁵ has 87 bits; shifted right by 40 it gives 100000000000000, and the 40 bits shifted out decide the rounding. In the code, slightly shortened, with p = 5^k and s the shift:

hc_mul128(m, p, &hi, &lo);
...
*q = (lo >> s) | (s == 0 ? 0 : hi << (64 - s));
...
if (hc_bit128(hi, lo, s - 1) && (hc_below128(hi, lo, s - 1) || (*q & 1))) {
    (*q)++;
}

The last test is round half to even: round up when the first bit shifted out is set, and either a lower one is set or the result is odd.

The result is an integer, so it becomes text through the routine that prints integers, the one in An alternative to utoa in C.

The single product made ordinary values faster than glibc, and the table did the same for extreme exponents. Over 1000 doubles with random bit patterns, which are mostly very large or very small, %.2f takes 619 ns on average, against 7,028 in glibc. On typical values the gain is smaller than on extreme ones: %.2f of 0.1 takes 113 ns, against 200 in glibc.

It is not faster everywhere. %.17g of 1e-300 still takes 2.4 microseconds, against 273 ns for %.15g of the same value. My guess is that at 17 digits the check cannot decide for that value, and it goes to the big number. I have not investigated it.

No 128-bit type

C89 has no 128-bit type. gcc and clang have unsigned __int128 as an extension, but gcc rejects it under -pedantic, so the product is built by hand with every compiler. Each 64-bit number is split into its high and low 32 bits. The four products of these parts fit in 64 bits, and adding them at their place, with the carries, gives the 128-bit result:

unsigned long a0 = a & 0xFFFFFFFFUL, a1 = a >> 32, b0 = b & 0xFFFFFFFFUL, b1 = b >> 32;
unsigned long p00 = a0 * b0, p01 = a0 * b1, p10 = a1 * b0;
unsigned long mid = (p00 >> 32) + (p01 & 0xFFFFFFFFUL) + (p10 & 0xFFFFFFFFUL);
*lo = (p00 & 0xFFFFFFFFUL) | (mid << 32);
*hi = a1 * b1 + (p01 >> 32) + (p10 >> 32) + (mid >> 32);

I hoped the compilers would see through it. They do not. gcc 14 turns it into 33 instructions with four imul, clang 19 into 27 with four imul. The same product written with unsigned __int128 compiles to one mul. The portable version still wins, because it replaces the big number for nearly every double, but it pays four multiplications where the hardware needs one.

32 bits

Where unsigned long has 32 bits, there is no 64-bit product to build on. The runtime falls back to the big number for every double, with 8-digit limbs. Against 32-bit glibc:

Slower than glibc on ordinary values, up to 14 times slower on the smallest ones, and faster only on whole numbers and huge %f. Building the same narrow algorithm on x86-64 gives about the same times, so it is the algorithm, not the instruction set.

Testing

I built a small unit test that compares it with the system's printf for 4000 random doubles at every precision from %.0e to %.17e, and at every binary exponent a double can have.

I also built every version with gcc, clang and tcc in 32 and 64 bits, and compared it with glibc on three sets of doubles. The first is the values I timed: nine chosen ones, 1000 random bit patterns and 1000 prices, in six formats. The second is the 44 values that Vern Paxson lists in A Program for Testing IEEE Decimal–Binary Conversion (1991) as the hardest to print with 1 to 22 digits: their exact value lies so close to halfway between two outputs that rounding them right takes more than 60 bits beyond the 53 of the double. The third is the 155,114 doubles of the testbase trace that the Go project publishes, each in %.0e to %.18e, as Go's own tests do. No build printed a different result.

I found one case where it differs from glibc, and I keep the difference on purpose. In C89, %#g keeps its trailing zeros, so 999999.5 prints as 1.00000e+06. glibc prints 1.e+06.

Next

The 32-bit path is the one I would improve next, with a fast path for a 32-bit long.