How GLIBC’s strtod() Works

The string to double function, strtod(), converts decimal numbers represented as strings into binary numbers represented in IEEE double-precision floating-point. Many programming environments implement their string to double conversions with David Gay’s strtod(); glibc, the GNU C Library, does not.

Like David Gay’s strtod(), glibc’s strtod() produces correctly rounded conversions. But it uses a simpler algorithm: it doesn’t have a floating-point only fast path for small inputs; it doesn’t compute a floating-point approximation to the correct result; it doesn’t check the approximation with big integers; it doesn’t adjust the approximation and recheck it; it doesn’t have an optimization for really long inputs. Instead, it handles all inputs uniformly, converting their integer and fractional parts separately, using only big integers. I will give an overview of how glibc’s strtod() works.

Overview

glibc’s strtod() is written in C. It lives in /stdlib/strtod_l.c, at entry point ____STRTOF_INTERNAL(). It is invoked from entry point STRTOF() in /stdlib/strtod.c, which acts as a wrapper. strtod() uses the low-level mpn_* GMP functions for its binary big integer operations.

strtod() starts by parsing its decimal input to determine its integer and fractional parts. For example:

• 1.789e6 has an integer part of 1789000 and no fractional part.
• 123e-8 has no integer part and a fractional part of 0.00000123.
• 456.833e2 has an integer part of 45683 and a fractional part 0.3.
• 9007199254740992.4567 has an integer part of 9007199254740992 and a fractional part of 0.4567.

It then converts these parts separately, first the integer part, and then the fractional part.

Since strtod() returns a double-precision floating-point number, it must produce 53-bits of precision (subnormal numbers require 52 bits or less, but I won’t consider that further); its algorithm for obtaining these bits is simple:

• Get the initial bits from the binary representation of the integer part, if there is one.
• Get the remaining bits, if more are needed, from the fractional part, if there is one.

Obtaining the bits from the integer part is easy: strtod() converts it to a big integer and takes its bits (just the first 53 if there are more than that). If there is a fractional part, and the integer part is less than 53 bits (any integer less than 2⁵³), then the remaining bits are obtained from the fractional part. strtod() calculates how many bits it needs and does big integer based division to obtain them.

strtod() rounds the 53-bit result according to the current IEEE 754 rounding mode (with round-to-nearest/ties-to-even being the default). It does this by inspecting the rounding bit (significant bit 54) and bits beyond. It creates the double by setting its sign, exponent, and significand fields directly, using a union of a double and a struct of integer fields.

Converting the Integer Part

strtod() creates a binary big integer from the digits of the integer part, using the straightforward conversion algorithm for integers: scan the next most significant decimal digit and convert it to binary based on its ASCII code; multiply the accumulating result by ten in binary and add the digit value; repeat until all digits are consumed. For example, 1789000 becomes the 21-bit binary integer 110110100110001001000, and 95514225908761452 becomes the 57-bit binary integer 101010011010101011010111010110010010010011100001101101100. In the case of the latter, there are more than 53 bits; strtod() rounds it to 10101001101010101101011101011001001001001110000110111 (assuming round-half-to-even rounding).

Powers Of Ten

If there is an exponent in the input string, strtod() multiplies the integer digits by the appropriate power of ten. For example, for input 1.789e6, strtod() creates integer 1789, and then multiplies it by 10³, getting 1789000.

Precomputed powers of ten — in big integer form — are stored in two tables. One table, called _tens_in_limb[], stores all powers of ten that fit in one GMP limb. (A limb is a 64-bit word on 64-bit systems and a 32-bit word on 32-bit systems.) For 64-bit limbs, _tens_in_limb[] stores powers of ten 10¹ through 10¹⁹. When the exponent is small enough, a power of ten from _tens_in_limb[] is used.

A second table, __tens[], stores the big integer representations of powers of ten of the form 10²ⁱ, for i = 0 to 14; that is, powers of ten 10²⁰ through 10²¹⁴. (Only powers of ten up to 10²¹⁰ are needed to convert doubles.) These are used to form the power of ten factor using binary exponentiation. The big integer representation of the integer digits is multiplied by the appropriate factors from this table, forming the complete big integer representation of the input.

Converting the Fractional Part

To represent the fractional part, strtod() creates a fraction with a big integer numerator and big integer denominator. For example, 0.00000123 is represented as 123/10⁸. (The power of ten is computed as described above.) Simply dividing the numerator by the denominator using multiple-precision floating-point binary division would give it the bits it needs. But like other conversion algorithms, strtod() uses integer division instead. To make integer division work, the numerator must be scaled up by a power of two that’s big enough so that the resulting integer quotient has the required number of bits.

strtod() could scale the numerator and then do a full-blown multiple-precision integer division. But that would be overkill, since only up to 53 bits are needed. Instead, it does the division manually, using only one or two native integer divisions (e.g., the div instruction on Intel processors).

When the numerator and denominator of the fraction are small enough, one native division is used to produce a full quotient. One division may also suffice when the numerator or denominator are too big; strtod() will truncate the numerator and possibly the denominator to compute a partial quotient. In some cases, two divisions are required; strtod() has to construct two numerator and denominator pairs that taken together produce the desired leading bits of the quotient. In all cases, the scaling is built into the process; the numerator is placed “high enough” to produce the desired number of bits.

This division process accounts for much of the complexity of strtod().

Deciding How Many Fractional Digits Are Necessary

For the fraction, strtod() takes only as many digits as it needs in order to fill out the required 53 bits. Unlike for the integer part, the approximately 3.32 bits per digit that you get should not be your guide (that was the assumption in the original code, but it led to incorrect conversions; it has since been fixed). In the worst case, you’ll need as many digits as you need bits.

For example, consider this decimal value:

0.000976562500000000325260651745651330202235840260982513427734375

In binary, it’s

0.000000000100000000000000000000000000000000000000000000000000011

That’s 2^-10 + 2^-62 + 2^-63, a value that is halfway between two floating-point numbers.

If you drop the last decimal digit (‘5’) from the input and convert the remaining 59 significant digits, you get a decimal value to about 195 bits of precision. But that’s still not enough to get bit 54 correct (it comes out as 0 instead of 1). You need all 60 input decimal significant digits to get it right — to round up in this case.

The Fix

For inputs that have only fractional parts, you have to find out where the significant bits start and gather your “one digit per bit” from there. In the example above, both representations have 63 digits, but they have different numbers of which are significant. The decimal value has 60 significant digits, offset by 3 leading zeros, and the binary value has 54 significant digits, offset by 9 leading zeros. So even though the number of digits are the same, the significant digits don’t correspond. From the decimal value, how do you know where the significant bits start, so that you can be sure to take enough decimal digits to include them?

The simple answer is to use a logarithm. If d is the decimal number, then abs(floor(log₂(d))) gives the starting position of the first significant bit. For our example, that’s 10. But computing a logarithm of an arbitrary number is expensive — and for our purposes, unnecessary — so an approximation is used instead.

That’s where “approximately 3.32 bits per decimal digit” fits properly into the picture. Think of it as the ratio of binary leading zeros, z₂ , to decimal leading zeros, z₁₀ . We know z₁₀ , so we can approximate z₂ as ceil(z₁₀ * 3.32). For our example, z₂ is 10. strtod() actually does a more conservative approximation: ((z₁₀ + 1) * 10) / 3 + 1. (Multiplying by 10 and dividing by 3 is the integer division way of multiplying by 10/3, or approximately 3.3.) In our example, that gives 14.

To create the numerator, strtod() takes the decimal digits from the first significant decimal digit up to the approximate location of the first significant bit and continues for as many digits as it needs bits (up to 53 — actually 54, for rounding). Any digits remaining after the ones taken are not made part of the numerator, but are considered in rounding (they are reflected in sort of a “sticky” bit, like in IEEE arithmetic).

Examples

Here are a few examples (remember, strtod() needs to generate 54 bits; the 54th bit is for rounding):

Example 1

For input 1234.56789012345678901234567890123456789, strtod() gets 11 bits from the integer part, and takes all 35 decimal digits from the fractional part. It creates the fraction 56789012345678901234567890123456789/10³⁵, and does one native integer division.

Example 2

For input 123456789012345.678901234567890123456789, strtod() gets 47 bits from the integer part, and takes the first 7 decimal digits from the fractional part. It creates the fraction 6789012/10⁷, and does one native integer division.

Example 3

For input 123456789012345.0034375, strtod() gets 47 bits from the integer part, and takes all 7 decimal digits from the fractional part. It creates the fraction 34375/10⁷, and does one native integer division.

This is like the prior example but with leading fractional zeros. These zeros are significant; the division gives a quotient with 7 leading 0s, 6 of which become part of the result. In other words, the resulting double has no fractional part; it’s just the integer 123456789012345.

Example 4

For input 123456789012345.0234375, strtod() gets 47 bits from the integer part, and takes all 7 decimal digits from the fractional part. It creates the fraction 234375/10⁷, and does one native integer division.

This is like the prior example except that the fractional part, 0.0234375, is exactly representable in binary: 0.0000011. All 7 of those bits figure into the answer; properly rounded, it’s binary 0.00001, or decimal 0.03125. The result is 123456789012345.03125 — it has 20 significant digits and a fractional part!

Example 5

For input 12345678901234567.8901234567890123456789, strtod() gets all 54 bits from the integer part; it ignores the fractional part.

Example 6

For input 9007199254740991.05, strtod() gets 53 bits from the integer part, and takes 1 decimal digit from the fractional part. (The single fractional digit is 0, but interestingly, strtod() creates a fraction of 0/10 and does the division anyway.)

Example 7

For this input

0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000022250738585072008890245868760858598876504231122409594654935248025624400092282356951787758888037591552642309780950434312085877387158357291821993020294379224223559819827501242041788969571311791082261043971979604000454897391938079198936081525613113376149842043271751033627391549782731594143828136275113838604094249464942286316695429105080201815926642134996606517803095075913058719846423906068637102005108723282784678843631944515866135041223479014792369585208321597621066375401613736583044193603714778355306682834535634005074073040135602968046375918583163124224521599262546494300836851861719422417646455137135420132217031370496583210154654068035397417906022589503023501937519773030945763173210852507299305089761582519159720757232455434770912461317493580281734466552734375

strtod() gets all 54 bits from the fractional part. This number is the largest subnormal double-precision floating-point number, 2^-1022 – 2^-1074. It has 1074 digits — 307 leading zeros and 767 significant digits. (I first saw the calculation of maximum significant digits done by Bruce Dawson, for single-precision floating-point; the largest float subnormal has 149 digits — 37 leading zeros and 112 significant digits. The number of digits in the largest double subnormal is shown on page 41 of the “Handbook of Floating-Point Arithmetic”.)

To get the 54 bits it needs, strtod() writes the decimal as a fraction with a big integer numerator consisting of the 767 significant digits and a big integer denominator of 10¹⁰⁷⁴. Fortunately though, it requires only two native integer divisions.

(The binary representation has 1022 leading zeros; strtod() conservatively puts it at 1027. Also, notice that the ratio of binary to decimal leading zeros is 1022/307, which is approximately 3.33.)