# 实现sprintf--浮点数打印字符串

TAG：无

## The Basics基本

constintc_maxDigits = 256;

// input binary floating-point number to convert from
doublevalue = 122.5;   // the number to convert

// output decimal representation to convert to
char  decimalDigits[c_maxDigits];// buffer to put the decimal representation in
int   numDigits = 0;              // this will be set to the number of digits in the buffer
int   firstDigitExponent = 0;     // this will be set to the the base 10 exponent of the first
//  digit in the buffer

// Compute the first digit's exponent in the equation digit*10^exponent
// (e.g. 122.5 would compute a 2 because its first digit is in the hundreds place)
firstDigitExponent = (int)floor(log10(value) );

// Scale the input value such that the first digit is in the ones place
// (e.g. 122.5 would become 1.225).
value = value / pow(10, firstDigitExponent);

// while there is a non-zero value to print and we have room in the buffer
while(value > 0.0 && numDigits < c_maxDigits)
{
// Output the current digit.
doubledigit = floor(value);
decimalDigits[numDigits] = '0'+ (char)digit;// convert to an ASCII character
++numDigits;

// Compute the remainder by subtracting the current digit
// (e.g. 1.225 would becom 0.225)
value -= digit;

// Scale the next digit into the ones place.
// (e.g. 0.225 would becom 2.25)
value *= 10.0;
}

?

• 精度的要求（比如只是打印小数点后6位）

• 如果数字因为宽度或者buffer不够而被截断，最后的数字应该要正确的舍入。

• 我们应该可以正确转换这个浮点数（包括舍入），使得这个数字是唯一表示的，可以精确从text转换回来。Integers整数

type sign bits exponent bits mantissa bits 32-bit float 1 8 23 64-bit double 1 11 52

exponent mantissa usage 非规格化 0 任何值 非规格化数字是一些太小而不能fit in the normalized range. 他们真正的值是这样计算的： mantissa2numMantissaBits?21+1?2numExponentBits?1. 32位的浮点数(代入32)就变成这样： mantissa223?21?127 64位的浮点数(代入64)就变成这样： mantissa252?21?1023 规格化 大于0小于最大值 任何值

(1+mantissa2numMantissaBits)?2exponent+1?2numExponentBits?1. 32位的浮点数(代入32)就变成这样： (1+mantissa223)?2exponent?127 64位的浮点数(代入64)就变成这样： (1+mantissa252)?2exponent?1023 Infinity 最大值 (32位的浮点数就是：0xFF 64位的double就是： 0x7FF 0 Inf的结果来自一些不恰当的计算比如除0. NaN 最大值 (32位的浮点数就是：0xFF 64位的double就是： 0x7FF  大于0 NaN 意思是 "not a number" 出现这个结果通常是计算得到的结果不是真正的数字，比如对负数开方.

representation conversion valueMantissa valueExponent 32-bit 非规格化数 The floating point equation is value=mantissa223?21?127
Factor a 223 out of the exponent and cancel out the denominator
value=mantissa223?223?21?127?23
value=mantissa?21?127?23 mantissa ?149 32-bit 规格化数 The floating point equation is value=(1+mantissa223)?2exponent?127
Factor a 223 out of the exponent and cancel out the denominator
value=(1+mantissa223)?223?2exponent?127?23
value=(223+mantissa)?2exponent?127?23 223+mantissa exponent?150 64-bit 非规格化数 The floating point equation is value=mantissa252?21?1023
Factor a 252 out of the exponent and cancel out the denominator
value=mantissa252?252?21?1023?52
value=mantissa?21?1023?52 mantissa ?1074 64-bit 规格化数 The floating point equation is value=(1+mantissa252)?2exponent?1023
Factor a 252 out of the exponent and cancel out the denominator
value=(1+mantissa252)?252?2exponent?1023?52
value=(252+mantissa)?2exponent?1023?52 252+mantissa exponent?1075

## Big Integers大数

12.3就可以用除数123和被除数10来表示了，因为123=123/10.

64位的double类型时，我们要用到超过1000位的整数来处理一些值。

? 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 constintc_maxDigits = 256;   // input number to convert from: value = valueMantissa * 2^valueExponent doublevalue = 122.5; // in memory 0x405EA00000000000                       // (sign=0x0, exponent=0x405, mantissa=0xEA00000000000) uint64_t valueMantissa = 8620171161763840ull; int32_t valueExponent = -46;   // output decimal representation to convert to char  decimalDigits[c_maxDigits];// buffer to put the decimal representation in int   numDigits = 0;              // this will be set to the number of digits in the buffer int   firstDigitExponent = 0;     // this will be set to the the base 10 exponent of the first                                    //  digit in the buffer      // Compute the first digit‘s exponent in the equation digit*10^exponent // (e.g. 122.5 would compute a 2 because its first digit is in the hundreds place)  firstDigitExponent = (int)floor(log10(value) );   // We represent value as (valueNumerator / valueDenominator) tBigInt valueNumerator; tBigInt valueDenominator; if(exponent > 0) {     // value = (mantissa * 2^exponent) / (1)     valueNumerator = BigInt_ShiftLeft(valueMantissa, valueExponent);     valueDenominator = 1; } else {     // value = (mantissa) / (2^(-exponent))     valueNumerator = mantissa;     valueDenominator = BigInt_ShiftLeft(1, -valueExponent); }   // Scale the input value such that the first digit is in the ones place // (e.g. 122.5 would become 1.225). // value = value / pow(10, firstDigitExponent) if(digitExponent > 0) {     // The exponent is positive creating a division so we multiply up the denominator.     valueDenominator = BigInt_Multiply( valueDenominator, BigInt_Pow(10,digitExponent) ); } elseif(digitExponent < 0) {     // The exponent is negative creating a multiplication so we multiply up the numerator.     valueNumerator = BigInt_Multiply( valueNumerator, BigInt_Pow(10,digitExponent) ); }   // while there is a non-zero value to print and we have room in the buffer while(valueNumerator > 0.0 && numDigits < c_maxDigits) {     // Output the current digit.     doubledigit = BigInt_DivideRoundDown(valueNumerator,valueDenominator);     decimalDigits[numDigits] = ‘0‘+ (char)digit;// convert to an ASCII charater     ++numDigits;             // Compute the remainder by subtracting the current digit     // (e.g. 1.225 would becom 0.225)     valueNumerator = BigInt_Subtract( BigInt_Multiply(digit,valueDenominator) );               // Scale the next digit into the ones place.     // (e.g. 0.225 would becom 2.25)     valueNumerator = BigInt_Multiply(valueNumerator,10); }

## Logarithms对数

（不好意思，CSDN的的公式很难表示，我先翻译到这里。sorry)

firstDigitExponent=?log10value?

log10value=log2valuelog210

log210=log1010log102=1log102

log10value=log2value?log102

firstDigitExponent=?log2value?log102?

log2value=log2(valueMantissa)+valueExponent

Because our mantissa is represented as a binary integer, the base-2 logarithm is close to the index of its highest set bit. Explicitly, the highest bit index is:

highestMantissaBitIndex=?log2(valueMantissa)?

Given that valueExponent is an integer, we can define:

?log2value?=?log2(valueMantissa)?+valueExponent

?log2value?=highestMantissaBitIndex+valueExponent

Given that highestMantissaBitIndex is also an integer, we can write the following inequality.

log2(value)?1<highestMantissaBitIndex+valueExponentlog2(value)

Scaling by the positive value log102, we get:

log2(value)?log102?log102<log102?(highestMantissaBitIndex+valueExponent)log2(value)?log102

log10(value)?log102<log102?(highestMantissaBitIndex+valueExponent)log10(value)

Because we are using floating-point numbers, we need to account for drift when we multiply by log102. This could put us slightly above our upper bounds so instead we will subtract a small epsilon (e.g. 0.00001) such that:

log10(value)?log102?epsilon<log102?(highestMantissaBitIndex+valueExponent)?epsilon<log10(value)

Let‘s now consider the following estimate.

log10(value)log102?(highestMantissaBitIndex+valueExponent)?epsilon

The inequality implies that our estimate will either be within an error of log102+epsilon below the correct value. Evaluating log102 we get approximately 0.301. This works out great because our error is less than one (even with the epsilon added to it) and what we really want to compute is ?log10(value)?. Thus, taking the floor of our approximation will always evaluate to the correct result or one less than the correct result. This image should help visualize the different cases as the logarithm slides between the neighboring integers and how the error range will overlap the lower integer boundary when the logarithm is low enough.

We can detect and account for the incorrect estimate case after we divide our value by 10firstDigitExponent. If our resulting numerator is still more than 10 times the denominator, we estimated incorrectly and need to scale up the denominator by 10 and increment our approximated exponent.

When we get to writing the actual Dragon4 algorithm, we‘ll find that to simplify the number of operations, we actually compute 1+?log10v?. To do this we end up taking the ceiling of the approximation instead of the floor. This aligns with the method from the Burger and Dybvig paper. Also, just like in the paper, we are going to do a few less high-precision multiplies when we generate an inaccurate estimate. However, what the paper fails to take advantage of, is computing the estimation such that it tries to be inaccurate more often than accurate. This will let us take the faster code path whenever possible. Right now, our error plus the epsilon is somewhere between 0.301 and 0.302 below the correct value. We then take the ceiling of that approximation so any error below 1.0 produces results that are off by one or correct. As the error approaches 1.0, we will get more results that are off by one. By adding some additional error less than 1.0?log102?epsilon, we can shift our error bounds lower. Explicitly, I subtract 0.69 from the estimate which moves the error range from about (-0.301, -0.69] to about (-0.991, -0.69].

The final inequality is:

?log10(value)?<?log102?(highestMantissaBitIndex+valueExponent)?0.69??log10(value)?+1

## Minimal Output and Rounding

As it stands, our prototype algorithm will print a number until all precision is exhausted (i.e. when the remainder reaches zero). Because there are discrete deltas between each representable number, we don‘t need to print that many digits to uniquely identify a number. Thus, we are going to add support for printing the minimal number of digits needed to identify a number from its peers. For example, let‘s consider the following three adjacent 32-bit floating-point values.

binary exact decimal representation v?1 0x4123c28e 10.2349987030029296875 v0 0x4123c28f 10.23499965667724609375 v1 0x4123c290 10.2350006103515625

There is a margin of 0.00000095367431640625 between each of these numbers. For any value, all numbers within the exclusive range (value?margin2,value+margin2) can uniquely identify it. For the above value, v0, we can uniquely identify it with the shorter number v0s=10.235 because it is closer to v0 than it is to to either v?1 or v1.

Solving this problem is what lets the Dragon4 algorithm output "pretty" representations of floating point numbers. To do it, we can track the margin as a big integer with the same denominator as the value. Every time we advance to the next digit by multiplying the numerator by 10, we will also multiply the margin by 10. Once the next digit to print can be safely rounded up or down without leaving the bounds of the margin, we are far enough from the neighboring values to print a rounded digit and stop.

In this image, we can see how the algorithm decides to terminate on the shortest decimal representation (the chosen numbers don‘t correspond to any real world floating-point format and are just for illustration).

• A: We start off with an input value that has the exact representation 7.625. The margin between each value is 0.25 which means that the set of values within 0.125 units below or above the input value can be safely rounded to it. This plus or minus 0.125 range has been highlighted by the box around the end of the value bar.
• B: We divide our value by 10 to get the first printed digit, 7. We then set the value to the division‘s remainder: 0.625. The low end of the margin box is greater than zero and the high end is less than one. This means that if we were to stop printing now and round the current digit to 7 or 8, we would print a number closer to another value than it is to 7.625. Because neither the rounded-down or rounded-up digit is contained by our margin box, we continue.
• C: To get our next output digit in the ones place, we multiply our value by 10 giving us 6.25. We also multiply our margin by 10 giving us a distance of 2.5 between numbers. The safe rounding distance of plus or minus half the margin is 1.25.
• D: We divide our value by 10 to get the second printed digit, 6. We then set the value to the division‘s remainder: 0.25. The low end of the margin box is now less than zero and the high end is greater than one. This means we can safely stop printing and output the current digit as a 6 or 7 because they will both uniquely identify our number. By checking if the remainder of 0.25 is below or above 0.5 we can decide to round down or up (we‘ll discuss ties shortly). In this case, we round down and the shortest output for 7.625 is "7.6".

## Unequal Margins

Unfortunately, the IEEE floating-point format adds an extra wrinkle we need to deal with. As the exponent increases, the margin between adjacent numbers also increases. This is why floating-point values get less accurate as they get farther from zero. For normalized values, every time the mantissa rolls over to zero, the exponent increases by one. When this happens, the higher margin is twice the size of the lower margin and we need to account for this special case when testing the margin bounds in code.

In understanding this better, I find it helpful to make a small floating-point format where we can examine every value. The tables below show every positive floating-point value for a 6-bit format that contains 1 sign bit, 3 exponent bits and 2 mantissa bits.

### Positive Denormalized Numbers

value=mantissa22?21+1?23?1=mantissa4?21?3

Here we can see our format‘s minimal margin of 0.0625 between each number. No smaller delta can be represented.

binary exponent mantissa equation exact value minimal value 0 000 00 0 0 (0/4)*2^(-2) 0 0 0 000 01 0 1 (1/4)*2^(-2) 0.0625 0.06 0 000 10 0 2 (2/4)*2^(-2) 0.125 0.1 0 000 11 0 3 (3/4)*2^(-2) 0.1875 0.2

### Positive Normalized Numbers

value=(1+mantissa22)?2exponent+1?23?1=(1+mantissa4)?2exponent?3

Here we can see our format‘s lowest exponent has the minimal margin of 0.0625 between each number. Each time the exponent increases, so does our margin by a factor of two and for our largest values the margin is 2.0.

binary exponent mantissa equation exact value minimal value 0 001 00 1 0 (1 + 0/4)*2^(-2) 0.25 0.25 0 001 01 1 1 (1 + 1/4)*2^(-2) 0.3125 0.3 0 001 10 1 2 (1 + 2/4)*2^(-2) 0.375 0.4 0 001 11 1 3 (1 + 3/4)*2^(-2) 0.4375 0.44 0 010 00 2 0 (1 + 0/4)*2^(-1) 0.5 0.5 0 010 01 2 1 (1 + 1/4)*2^(-1) 0.625 0.6 0 010 10 2 2 (1 + 2/4)*2^(-1) 0.75 0.8 0 010 11 2 3 (1 + 3/4)*2^(-1) 0.875 0.9 0 011 00 3 0 (1 + 0/4)*2^(0) 1 1 0 011 01 3 1 (1 + 1/4)*2^(0) 1.25 1.2 0 011 10 3 2 (1 + 2/4)*2^(0) 1.5 1.5 0 011 11 3 3 (1 + 3/4)*2^(0) 1.75 1.8 0 100 00 4 0 (1 + 0/4)*2^(1) 2 2 0 100 01 4 1 (1 + 1/4)*2^(1) 2.5 2.5 0 100 10 4 2 (1 + 2/4)*2^(1) 3 3 0 100 11 4 3 (1 + 3/4)*2^(1) 3.5 3.5 0 101 00 5 0 (1 + 0/4)*2^(2) 4 4 0 101 01 5 1 (1 + 1/4)*2^(2) 5 5 0 101 10 5 2 (1 + 2/4)*2^(2) 6 6 0 101 11 5 3 (1 + 3/4)*2^(2) 7 7 0 110 00 6 0 (1 + 0/4)*2^(3) 8 8 0 110 01 6 1 (1 + 1/4)*2^(3) 10 10 0 110 10 6 2 (1 + 2/4)*2^(3) 12 11 0 110 11 6 3 (1 + 3/4)*2^(3) 14 12

### Positive Infinity

binary exponent mantissa value 0 111 00 7 0 Inf

### Positive NaNs

binary exponent mantissa value 0 111 01 7 1 NaN 0 111 10 7 2 NaN 0 111 11 7 3 NaN

As far as I can tell, the original Dragon4 paper by Steele and White actually computes the low margin wrong for the lowest normalized value because it divides it in half. It should actually remain the same because the transition from denormalized to normalized numbers does not change the exponent. This error is not present in the follow-up paper by Burger and Dybvig.

## Breaking Ties

When deciding how to round-off the final printed digit, it is possible for the exact value to lie directly between the lower digit and the higher digit. In this case, a rounding rule needs to be decided. You could round towards zero, away from zero, towards negative infinity, etc. In my implementation, I always round to the even digit because it is similar to the rule used by IEEE floating-point operations.

In the reverse algorithm of converting decimal strings into binary floating-point, you can hit a similar case where the string representation lies exactly between to representable numbers. Once again, a rounding rule must be chosen. If you know this rule in both algorithms, there are times when you can print out a unique representation of a number with one less digit. For example, consider how the digit printing algorithm tests if half the high-margin encompasses a digit to round-up to. I tested if the remaining value was greater than one. If we knew how the reverse algorithm would handle ties, we might be able to test greater than or equal to one. David Gay‘s dtoa implementation seems to do this; therefor it should always be run in conjunction with his strtod implementation.

Personally, I prefer printing extra digits until the rounded number is inclusively contained within the margin. This way the output is guaranteed to uniquely identify the number regardless of the tie breaking algorithm chosen by a strtod parser. Unfortunately, this is only guaranteed for a proper strtod implementation. The implementation that ships with Microsoft Visual C++ 2010 Express sometimes has an off-by-one error when rounding to the closest value. For example, the input "6.439804741657803e-031" should output the 64-bit double 0x39AA1F79C0000000, but instead outputs 0x39AA1F79BFFFFFFF. I can‘t say for certain where Microsoft‘s bug comes from, but I can say that the problem of reading a decimal string is more difficult than printing a decimal string so it‘s more likely to have bugs. Hopefully I can find the time to try my hands at it in the near future.

## Coding Dragon4

Next, we will program an optimized implementation of the Dragon4 algorithm. Click here to read part 2..