(Please let me know what you think — take the poll at the end of this article.)

Why I Don’t Like ‘Binary Fraction’

Of all the existing terms, binary fraction is probably the most commonly used. I don’t like it because its analog, decimal fraction, is not clearly defined. I want to avoid a term that inherits this problem.

Decimal fraction is commonly defined as any number with an explicit or implicit power of ten denominator, either entirely fractional or not. For example, 254/1000, 0.254, 15/10, and 1.5 are decimal fractions. But what about 2/5? It can be written as 4/10 or 0.4, although as written its denominator is not a power of ten. And what about 1/3? Its equivalent form as a decimal is 0.3 — a repeating decimal. It can never be written with a power of ten denominator. To make things more confusing, 2/5 and 1/3 are fractions — fractions written in decimal numerals.

You could argue that decimal fraction includes 2/5 but excludes 1/3 and still have a reasonable definition. However, I’m looking for the equivalent of decimal, a term which includes 0.4 and 0.3, but excludes 2/5, 4/10, and 1/3.

‘Bicimal’

I discovered the term bicimal on the Web and in Google Books, but I don’t know its origin. I pronounce it “bye’ suh mull”, or as Merriam-Webster might express it, \ˈbī-sə-məl\. A bicimal is built with negative powers of two, whereas a decimal is built with negative powers of ten.

Like the term decimal, bicimal usually means a pure fractional value, like 0.11001. However, in some contexts, it could mean numbers with a whole and fractional part, like 101.11. In this case, nonnegative powers of two come into play — for the whole part.

Why I Like ‘Bicimal’

Ideally, there would be a base-independent term for the fractional part of a number. I invented the term fractional for this purpose. I’ve called a fractional decimal number a decimal fractional, and a fractional binary number a binary fractional. The purpose of this new term was to separate the form of a number — a number with a “point” in it — from its base. If 0.427 is a decimal, does that make 0.11001 a binary decimal? You can see why we need a better term.

One problem with my terminology is that I’ve created two terms (decimal fractional and binary fractional) when I really only needed to create one. Why not stick with decimal and invent a new term just for binary? Decimal is easier to say than decimal fractional, and everyone knows what it means. So what’s a good replacement for binary fractional? I’ve come to like the term bicimal.

At first glance, there’s not much to like about bicimal. It’s base-dependent, and it is a portmanteau for binary decimal. It is a poorly formed portmanteau at that. While the prefix ‘bi-’ is perfectly acceptable, its pairing with the suffix ‘-cimal’ seems ill-formed. Binimal seems linguistically the better choice; it swaps out the prefix ‘dec-’ for the prefix ‘bin-’ and retains the suffix ‘-imal’.

On the other hand, say bicimal and binimal outloud, over and over; I think you’ll find that bicimal sounds better, as do its associated terms: bicimal point, bicimal places, bicimal part, terminating bicimal, repeating bicimal, infinite bicimal, etc. And bicimal produces natural sounding phrases like “multiply bicimals”, “convert a decimal to a bicimal”, “convert a bicimal to a decimal”, “convert a bicimal to a fraction”, “convert a fraction to a bicimal”, etc.

I think bicimals will be immediately understood by newcomers. It evokes all the feelings and terminology and operations of decimals (for better or worse ). I don’t think binimals — or any of the alternative terms — has this property. So all things considered, I like bicimal the best.

What Do You Think?

Please leave a comment or take the poll (Polldaddy poll is available only if Javascript is enabled). If you like binary fraction, please tell me why.

By Rick Regan (Copyright © 2008-2012 Exploring Binary)

Bicimals

1.8254370818746402660437411213933955878019332885742187

As Geza noted, the problem is in the bigcomp() function, an optimization that kicks in for long decimal inputs. I traced his example through bigcomp() — I’ll show you what’s going on.

(Geza has reported the bug, and David Gay has already fixed it.)

strtod() Returns a Value That Is One ULP Too High

In binary, 1.8254370818746402660437411213933955878019332885742187 is a number with a non-terminating binary fraction portion; here are its first 175 significant bits:

1.11010011010011111101100000110111100011101010100000110111111111111
1111111111111111111111111111111111111111111111111111111111111111111
111111111111111111111111111111111111111110...

This is extremely close to — but less than — halfway between two double-precision numbers: it has a 0 at significant bit 54 (highlighted), followed by 120 1s. (A halfway case has a 1 at bit 54 followed by all 0s.)

Correctly rounded to the nearest double-precision value this is

1.1101001101001111110110000011011110001110101010000011,

or 0x1.d34fd8378ea83p+0 in hexadecimal “%a” notation.

In decimal, this equals

1.8254370818746401550214386588777415454387664794921875.

strtod() returns 0x1.d34fd8378ea84p+0, which is 1 binary ULP too high.

What strtod() Does

strtod() starts by computing its floating-point approximation as 0x1.d34fd8378ea83p+0 (which happens to be the correct result). Next, it uses big integer arithmetic to compare the approximation to the decimal input. With the bigcomp optimization enabled, which it is by default, only the first 18 digits of the decimal input are used for the comparison; this truncation, in a close case like this, means strtod() must call bigcomp() to determine if the approximation is correct.

bigcomp()

bigcomp() decides which of two adjacent floating-point numbers is closest to the decimal input string. It generates the decimal digits of the number halfway between these two values, comparing them digit-by-digit to the input string. In this case, the halfway value is

1.82543708187464026604374112139339558780193328857421875,

which is the same as the input but with a ‘5’ appended as the 54th digit. bigcomp() quits after examining the 53 digits of its input — as it should — but it fails to account for the implied 0 at digit 54. Mistakenly, it deems the input greater than the halfway value, when it should see that the implied 0 makes the input less than the halfway value.

The Fix

Here is the description of the fix, as logged by David Gay:

20120417
 dtoa.c:  augment a test in bigcomp() to correctly compare a long
input string of digits with a computed string after all of the
input string has otherwise been processed.  This matters only
rarely; a string for which it matters (reported by Herman Geza [sic])
is 1.8254370818746402660437411213933955878019332885742187 .

The fix changes one line of code; this line

if (b->x[0] || b->wds > 1)

is changed to this:

if (b->x[0] || b->wds > 1 || dig > 0)

In this example, dig = 5, so the ‘if’ branch is taken, resulting in the smaller of the two neighboring values to be chosen.

Alternate Fix

You can avoid the bug by disabling bigcomp, which you can do by “#defining” the variable “NO_STRTOD_BIGCOMP”.

By Rick Regan (Copyright © 2008-2012 Exploring Binary)

A Bug in the Bigcomp Function of David Gay’s strtod()

Example of Binary Division

The pencil-and-paper method of binary division is the same as the pencil-and-paper method of decimal division, except that binary numerals are manipulated instead. As it turns out though, binary division is simpler. There is no need to guess and then check intermediate quotients; they are either 0 are 1, and are easy to determine by sight.

Decimal Division

Pencil-and-paper division, also known as long division, is the hardest of the four arithmetic algorithms. Like the other algorithms, it requires you to solve smaller subproblems of the same type. But unlike the other algorithms, there is no limited set of “facts” that solve all possible subproblems. Solving these division subproblems requires estimation, guessing, and checking. In addition to these division subproblems, multiplication and subtraction are required as well.

Let’s review how decimal division is done, so that we can set the stage for how division is done in binary. Here is an example:

Example of Decimal Division

The algorithm is a series of steps, each step having these four substeps:

1. Divide: Divide the working portion of the dividend by the divisor. (The dividend is the number under the line.)
2. Multiply: Multiply the quotient (a single digit) by the divisor.
3. Subtract: Subtract the product from the working portion of the dividend. (I don’t write down minus signs — they’re implied.)
4. Bring down: Copy down the next digit of the dividend to form the new working portion.

Here’s the example again, step-by-step:

Steps of Decimal Division

• Step 0

  Does 88 go into 8? No, because it’s greater than 8. Does 88 go into 83? No, because it’s greater than 83. Does 88 go into 831? Yes, because it’s less than or equal to 831.

  (The first step of long division, as commonly practiced, combines several steps and their substeps into one. That’s why I call this step 0. Technically, 88 goes into 8 zero times, so we should write down a 0, multiply 88 by 0, subtract 0 from 8, and then bring down the 3. Next, we should write down a 0 because 88 goes into 83 zero times, multiply 88 by 0, subtract 0 from 83, and bring down the 1. We’re just eliminating a bunch of stuff that produces superfluous leading zeros.)

• Step 1
  1. Divide: Does 88 go into 831? (Yes, we already know that from step 0.) How many times does it go in? 9 times. (Normally you have to guess the answer and do the multiplication and subtraction to verify you guessed correctly; I’ll just show the correct answer, to keep things from getting too messy.)
  2. Multiply: 9 x 88 = 792. (If 9 were a guess, your first check is to see if 792 is less than 831. It is, so so far, so good.)
  3. Subtract: 831 – 792 = 39. (If 9 were a guess, your second check is to see if 39 is less than 88. It is, so your guess was correct. Actually, in this case, this has to check out, since we can’t go higher than 9.)
  4. Bring down: Bring down the 2 to make 392.
• Step 2
  1. Divide: Does 88 go into 392? Yes, 4 times.
  2. Multiply: 4 x 88 = 352.
  3. Subtract: 392 – 352 = 40.
  4. Bring down: Bring down the implied trailing 0 to make 400.
• Step 3
  1. Divide: Does 88 go into 400? Yes, 4 times.
  2. Multiply: 4 x 88 = 352.
  3. Subtract: 400 – 352 = 48.
  4. Bring down: Bring down the implied trailing 0 to make 480.
• Step 4
  1. Divide: Does 88 go into 480? Yes, 5 times.
  2. Multiply: 5 x 88 = 440.
  3. Subtract: 480 – 440 = 40.
  4. Bring down: Bring down the implied trailing 0 to make 400.
• Step 5

  Stop the presses! We tried to divide 400 by 88 before — two steps ago. That means we have a two-digit cycle (45) from here on out. The answer is 9.445.

The red digits are the carries that occur during the multiplication substeps (the multiplication is done as if the divisor — the bigger number — is on top, by convention). Each red digit is crossed out before the next multiplication. To avoid clutter, I have chosen not to mark the borrows that occur during subtraction.

I Picked The Hardest Type of Example

My example has a multi-digit divisor, and has an answer with a remainder that I wrote as a repeating decimal. I wanted one example that showed long division to its fullest. I could have picked a problem with a single-digit divisor (which would require no guessing, assuming you know the multiplication facts), or one that produced an integer quotient, or one that produced a quotient with a fractional part that terminated. I could have expressed the fractional part as an integer remainder, or in fraction form.

Other Cases

• If the divisor or dividend is negative, you can remove the signs and apply the appropriate sign to the answer at the end.
• If the divisor has a decimal point, shift the decimal point right until the divisor is an integer, and shift the dividend by the same number of places.
• If the divisor is greater than the dividend, just proceed with the algorithm as is. Trailing zeros will be brought down to form the appropriate subproblems.

Binary Division

Let’s return to the example of the introduction, 1011.11/11. Here it is broken down into steps, following the same algorithm I used for decimal numbers:

Steps of Binary Division

• Step 0

  Does 11 go into 1? No, because it’s greater than 1. Does 11 go into 10? No, because it’s greater than 10. Does 11 go into 101? Yes, because it’s less than or equal to 101. (Remember, these are binary numerals; pronounce them “one-one”, “one-zero”, “one-zero-one”, etc.)

• Step 1
  1. Divide: Does 11 go into 101? (Yes, we already know that from step 0.) How many times does it go in? One time. There is no guessing. It’s easy to see 11 is less than 101, so we know it goes in. And if it goes in, it goes in only once.
  2. Multiply: 1 x 11 = 11. (Remember how simple it is to “multiply” a binary number by a single digit — just copy the number down if that single digit is 1, or write down 0 if that single digit is 0.)
  3. Subtract: 101 – 11 = 10.
  4. Bring down: Bring down the 1 to make 101.
• Step 2
  1. Divide: Does 11 go into 101? Yes, 1 time.
  2. Multiply: 1 x 11 = 11.
  3. Subtract: 101 – 11 = 10.
  4. Bring down: Bring down the 1 to make 101.
• Step 3
  1. Divide: Does 11 go into 101? Yes, 1 time.
  2. Multiply: 1 x 11 = 11.
  3. Subtract: 101 – 11 = 10.
  4. Bring down: Bring down the 1 to make 101.
• Step 4
  1. Divide: Does 11 go into 101? Yes, 1 time.
  2. Multiply: 1 x 11 = 11.
  3. Subtract: 101 – 11 = 10.
  4. Bring down: Bring down the 0 to make 100.
• Step 5
  1. Divide: Does 11 go into 100? Yes, 1 time.
  2. Multiply: 1 x 11 = 11.
  3. Subtract: 100 – 11 = 1.
  4. Bring down: Bring down the 0 to make 10.
• Step 6
  1. Divide: Does 11 go into 10? No (write down a 0).
  2. Multiply: (We don’t need to record this step; we’re just going to get 0.)
  3. Subtract: (We don’t need to record this step; we’re just going to get 10.)
  4. Bring down: Bring down the 0 to make 100.
• Step 7

  We stop here, recognizing that we divided 100 by 11 two steps ago. This means we have a two-digit cycle (10) from here on out. The quotient is 11.1110.

Checking the Answer

When the answer has a repeating fractional part, checking it is not as straightforward as it is for the other arithmetic operations. What we can do is approximate the quotient to a finite number of places and then check that it comes close to the expected answer.

You can check the answer in a few ways. One way is by doing binary multiplication by hand: you verify that the approximated quotient (11.11101011, for example) multiplied by the divisor (11) equals the dividend (1011.11). (I’ll leave that as an exercise, but the answer is 1011.11000001, which is very close to 1011.11).

Another way to check is to convert the operands to decimal, do decimal division, and then convert the approximate decimal answer to binary. 1011.11 = 11.75, and 11 = 3. 11.75/3 = 3.916. Estimating that as 3.91666666666666667, for example, my binary converter says it equals 11.111010101010101010101010101010101010 when truncated to 36 places. That looks like it wants to be 11.1110, the answer we got using binary division.

You can also check the answer using my binary calculator. It says 1011.11/11 is 11.111010101010 (to 12 places, for example). Again, that looks like 11.1110.

If you want to verify the repeating part directly, you can use this conversion tool; here’s what to enter:

• Initial Base: 2
• Integral part: 11
• Non-Repeating Fractional Part: 11
• Repeating Fractional Part: 10
• New Base: 10

It gives the decimal answer we expect: 3.916.

You can also use this tool to convert in the opposite direction, verifying that 3.916 converts to 11.1110.

Discussion

Like the other arithmetic algorithms, I described the division algorithm in a base-independent way. I wanted to stress the mechanical procedure, not why it works (in either decimal or binary).

When you do binary long division, you might find yourself doing some of the substeps in your head in decimal (e.g., 101 – 11 is 5 – 3 = 2, which is 10 in binary).

Although binary division is easier than decimal division (because there’s no guessing and effectively no multiplication), you will find that always having the same number (the divisor) as the subtrahend will produce a pattern that will start mesmerizing you; it’s easy to get lost in that sea of 1s and 0s. (Be thankful my example only had a two-digit repeating cycle!)

If you play around with binary division you’ll see that it produces more repeating fractional numbers than decimal division does. For example, 2/5 = 0.4, but 10/101 = 0.0110.

Further Reading

There are many explanations of binary division on the Web; one that I like in particular, and that comes closest to what I’ve explained, is Dr. Math’s “Long Division in Binary.”

By Rick Regan (Copyright © 2008-2012 Exploring Binary)

Binary Division

Example of Binary Multiplication

The pencil-and-paper method of binary multiplication is just like the pencil-and-paper method of decimal multiplication; the same algorithm applies, except binary numerals are manipulated instead. The way it works out though, binary multiplication is much simpler. The multiplier contains only 0s and 1s, so each multiplication step produces either zeros or a copy of the multiplicand. So binary multiplication is not multiplication at all — it’s just repeated binary addition!

Decimal Multiplication

To multiply two multiple-digit decimal numbers, you first need to know how to multiply two single-digit decimal numbers. This requires the memorization of 100 facts, or 55 facts if you exclude the commutative or “turnaround” facts. These facts are usually represented in a “multiplication table,” also known as a “times table.” Example facts are 2 x 9 = 18, 9 x 7 = 63, and 1 x 6 = 6.

A multiplication problem is written with one number on top, called the multiplicand, and one number on the bottom, called the multiplier. The algorithm has two phases: the multiplication phase, where you produce what are called partial products, and the addition phase, where you add the partial products to get the result.

In the multiplication phase, the digits of the multiplier are stepped through one at a time, from right to left. Each digit of the multiplicand is then multiplied, in turn, by the current multiplier digit; taken together, these single-digit multiplications form a partial product. The answer to each single-digit multiplication comes from the multiplication table. Some of these answers are double-digit numbers, in which case the least significant digit is recorded and the most significant digit is carried over to be added to the result of the next single-digit multiplication.

For example, let’s multiply 3.87 and 5.3:

Example of Decimal Multiplication

There are two digits in the multiplier, so there are two partial products: 1161 and 19350. Each partial product has its own set of carries, which are crossed out before computation of the next partial product. Here is the multiplication phase, broken down into steps:

Steps of Decimal Multiplication (Multiplication Phase Only)

When the multiplication phase is done, the partial products are added, and the decimal point is placed appropriately. (If there were any minus signs, they would be taken into account at this point as well.) This gives the answer 20.511.

Binary Multiplication

Binary multiplication uses the same algorithm, but uses just three order-independent facts: 0 x 0 = 0, 1 x 0 = 0, and 1 x 1 = 1 (these work the same as in decimal). If you perform the multiplication phase with these facts, you’ll notice two things: there are never any carries, and the partial products will either be zeros or a shifted copy of the multiplicand.

Observing this, you’ll realize there’s no need for digit-by-digit multiplication, which means there’s no need to consult a times table — which means there’s no multiplication, period! Instead, you just write down 0 when the current digit of the multiplier is 0, and you write down the multiplicand when the current digit of the multiplier is 1.

In the introduction, I showed this example: 1011.01 x 110.1. I wrote it as if you followed the decimal algorithm to the letter. Here’s how it looks if you follow the simpler “write zero or multiplicand” algorithm (it’s the same result, but with blanks representing 0s; this matches better conceptually with what we are now doing):

Example of Binary Multiplication (Simplified View)

Here’s what the “multiplication” phase looks like, step-by-step:

Steps of Binary Multiplication (Multiplication Phase Only)

Each step is the placement of an entire partial product, unlike in decimal, where each step is a single-digit multiplication (and possible addition of a carry).

In the addition phase, the partial products are added using binary addition, and then the radix point is placed appropriately. This gives the answer 1001001.001.

Checking the Answer

You can check the answer by converting the operands to decimal, doing decimal multiplication, and then converting the decimal answer to binary. 1011.01 = 11.25, and 110.1 = 6.5. Their product is 73.125, which is 1001001.001, the answer we got using binary multiplication.

You can also check the answer using my binary calculator.

Discussion

Computers don’t multiply in exactly this way, but they do exploit the simplified view of binary multiplication that I’ve described.

By Rick Regan (Copyright © 2008-2012 Exploring Binary)

Binary Multiplication

A Binary Multiplication Problem Expressed With One Hundred Acorn Tops.

(Notice it’s a little asymmetric. The zero in the upper right has only 15 acorn tops, not 16 like the rest; I stole one for the times sign.)

There are two ways to solve this problem. The first way is to convert the operands to decimal and then multiply using decimal arithmetic:

1. Binary 1010 = 1 x 2³ + 0 x 2² + 1 x 2¹ + 0 x 2⁰ = 2³ + 2¹ = 8 + 2 = 10 decimal.
2. 10 x 10 = 100.

The second way is to multiply using binary arithmetic and then convert the answer to decimal:

1.       1010
       x 1010
       ------
         0000
        10100
       000000
      1010000
      -------
      1100100
   

2. Binary 1100100 = 1 x 2⁶ + 1 x 2⁵ + 0 x 2⁴ + 0 x 2³ + 1 x 2² + 0 x 2¹ + 0 x 2⁰ = 2⁶ + 2⁵ + 2² = 64 + 32 + 4 = 100 decimal.

The first way is easier, but it’s not as much fun .

By Rick Regan (Copyright © 2008-2012 Exploring Binary)

One Hundred Acorns in Binary

Example of Binary Subtraction

The pencil-and-paper method of binary subtraction is just like the pencil-and-paper method of decimal subtraction you learned in elementary school. Instead of manipulating decimal numerals, however, you manipulate binary numerals, according to a basic set of rules or “facts.”

Decimal Subtraction

For decimal subtraction, the basic facts are things like 5 – 1 = 4, 9 – 8 = 1, and 18 – 9 = 9. In each case, the answer is a single-digit, nonnegative integer. Most of the facts are “single-digit minus single-digit” problems, but some are “double-digit minus single-digit” problems (the double digits are the numbers 10 through 18). The latter represent cases of borrowing, which is the process by which negative answers are prevented.

Here’s an example of decimal subtraction:

Example of Decimal Subtraction

After the points are aligned, subtraction proceeds from right to left. Red marks indicate borrowing. If a non-zero digit is borrowed from, it is crossed out, one is subtracted from it, and the decremented digit is written above it; a 1 is then placed next to the digit in the borrowing position, making it a two-digit number. If a zero digit is borrowed from, the borrow “cascades” until a non-zero digit is found.

Here’s the example again, step-by-step:

Steps of Decimal Subtraction

• Step 1: 5 – 0 = 5.
• Step 2: Borrow to make 12 – 5 = 7.
• Step 3: Borrow to make 15 – 7 = 8.
• Step 4: Borrow to make 10 – 1 = 9.

Some people refer to this as the “American method” (although this is just one variation of it — see Salman Khan’s video, for example). Whatever your method is though, you can apply it to binary numbers.

Binary Subtraction

For binary subtraction, there are four facts instead of one hundred:

• 0 – 0 = 0
• 1 – 0 = 1
• 1 – 1 = 0
• 10 – 1 = 1

The first three are the same as in decimal. The fourth fact is the only new one; it is the borrow case. It applies when the “top” digit in a column is 0 and the “bottom” digit is 1. (Remember: in binary, 10 is pronounced “one-zero” or “two.”)

Now let’s subtract 1011.11 from 10101.101, following the same algorithm I used for decimal numbers:

Steps of Binary Subtraction

• Step 1: 1 – 0 = 1.
• Step 2: Borrow to make 10 – 1 = 1.
• Step 3: Borrow to make 10 – 1 = 1.
• Step 4: Cascaded borrow to make 10 – 1 = 1.
• Step 5: 1 – 1 = 0.
• Step 6: 0 – 0 = 0.
• Step 7: Borrow to make 10 – 1 = 1.

Since there are lots of 0s in binary numbers, there can be lots of borrows — and lots of messy looking cross-outs.

Checking the Answer

You can check the answer in a few ways. One way is to add the result (1001.111) to the subtrahend (1011.11), and check that that answer matches the minuend (10101.101):

Check Binary Subtraction Using Binary Addition

Another way is to convert the operands to decimal, do decimal subtraction, and then convert the decimal answer to binary. 10101.101 = 21.625 and 1011.11 = 11.75, and 21.625 – 11.75 = 9.875. 9.875 = 1001.111, the answer we got using binary subtraction.

You can also check the answer using my binary calculator.

Subtracting a Bigger Number From a Smaller Number

To subtract a bigger number from a smaller number, just swap the numbers, do the subtraction, and negate the result.

Discussion

Notice that I didn’t discuss the number base when describing the algorithm; it is base-independent. Nonetheless, I could have talked about powers of ten and powers of two, and how the process can be visualized by regrouping. My goal was to explain just the mechanized algorithm (presumably you do decimal subtraction mechanically, no longer thinking about why it works).

Subtracting Using Complements

Computers don’t subtract this way; they subtract by adding complements. It’s more efficient. You can do subtraction by complements with pencil-and-paper, but you won’t find it more efficient. (In decimal, you would use nine’s complement or ten’s complement; in binary, you would use ones’ complement or two’s complement.)

By Rick Regan (Copyright © 2008-2012 Exploring Binary)

Binary Subtraction

“In every speech I give, I talk to people who are still running Apple IIs, and they say those machines are still running after this many years.”

So I got it out of the attic and powered it up. The dozen or so dead keys notwithstanding, it still works — after 30 years!

Some BASIC Commands I Tried On My Keyboard-Challenged But Otherwise Still Working Apple II

Getting The Monitor To Work

I was expecting not to have the cable or connector I needed to hook the computer up to a TV (as I recall, there was a splitter that attached from the video cable to the TV’s VHF inputs). Also, I was wondering if I had a TV with the right inputs. It turns out I had what I needed.

The video cable was already attached to the computer, inside the case to an RF modulator. Hooking that to the TV’s video input did not work. After looking online for documentation (I no longer have the manuals), I found a copy of the “Apple II Basic Programming Manual”; on page 5 I saw this:

“If you have a color (or black and white) monitor, just connect the appropriate cable from the jack marked “VIDEO OUT” (on the rear of the Apple II) to the input of the monitor.

If you have an ordinary TV, you will have to install an RF modulator.”

I never used a monitor back in the day — just a TV. So all I did was move the cable to the rear input — and I got video!

Getting Into BASIC

All that appeared on the screen was the “Apple ][” logo. My disk drive was attached, and it was spinning (there was no floppy disk in it — I don’t have those either). It occurred to me to remove the floppy drive controller card so that the computer would boot into ROM BASIC. I powered up again and, bingo!

A Partially Working Keyboard

Unfortunately, many keys did not work (among others: a, b, i, o, y, 4, -, =, and ,). I opened the bottom of the case to see if there was anything obviously wrong, but I didn’t see anything (surprisingly though, there was no dust inside the case).

An Example BASIC Program

I wanted to make some binary related program using the available keys. Luckily, I had double quotes and the question mark, so I could print stuff. I also had exponentiation (‘^’), addition, and division, which were enough to make the small program pictured above (the picture shows our old TV, before I realized that with correct attachment, it would work on our newer TV):

• Line 10: Prints “Rick Regan exploringbinary.com 2/5/12”, or the closest thing to it with missing keys and all caps.
• Line 20: Prints 2³², which displays in scientific notation, correctly rounded to eight significant digits, as 4.2949673E+09.
• Line 30: Prints 2^-32, which displays in scientific notation, correctly rounded to nine significant digits, as 2.32830644E-10.
• Line 50 (remember my ‘4’ key doesn’t work): Prints 2⁰ + 2² + 2³ + 2⁵ + 2⁷ + 2¹² + 2¹⁵, which corresponds to the binary representation of 37037.
• Line 60: Prints 2.2250738585072011/10³⁰⁸ (which is how I had to enter 2.2250738585072011e-308 without a minus sign). It overflows, but at least it doesn’t hang!

Is the Apple II Really a Computer?

My kids watched me as I did this. “Where’s the mouse?” “How do you click on things?” “Does it play Angry Birds?” The problem was I kept calling it a computer, but it wasn’t like any computer they’d ever seen.

By Rick Regan (Copyright © 2008-2012 Exploring Binary)

Exploring Binary On My Apple II

Example of Binary Addition

Binary Arithmetic

Binary arithmetic is of interest because that’s how computers do math. When implemented in computers, many things must be taken into account: format (fixed-point, floating-point, etc.), word size (8-bit, 16-bit, 32-bit, etc.), sign representation (sign-magnitude, ones’ complement, two’s complement, etc.), overflow (when numbers are too big), and underflow (when numbers are too small). Furthermore, there is the design of the circuits to perform the algorithms (full adders, half adders, ripple carry adders, carry-lookahead adders, restoring dividers, non-restoring dividers, etc.). These are not inherent properties of binary arithmetic — they’re just implementation issues.

The truth is, in its purest form, binary arithmetic is very simple. The same pencil-and-paper algorithms you learned in grade school work for binary numbers. All you do differently is apply the “facts” for binary numerals, the (small) set of rules for manipulating 0s and 1s. And binary numbers on paper are written as you’d expect: without leading zeros, and with a minus sign (‘-’) if negative.

In my series of articles, I will explain the pure forms of the four basic operations of binary arithmetic: addition, subtraction, multiplication, and division. In this first article, I will discuss binary addition.

Decimal Addition

To add two multiple-digit decimal numbers, you first need to know how to add two single-digit decimal numbers. This requires the memorization of 100 facts, or 55 facts if you exclude the commutative or “turnaround” facts. Also, because of carries, you need to know ten additional facts: 10 + 0 = 10, 10 + 1 = 11, … , 10 + 9 = 19. The latter apply when there’s a carry (always 1) and the “top” digit is 9.

For example, let’s add 19.7 and 12.8:

Example of Decimal Addition

Let’s think about what a carry does. It turns our neat “add two single-digit numbers per column” problem into “add three single-digit numbers per column.” But addition is a binary (in a different sense) operation, meaning it operates on two numbers at a time. So what we do is add the carry to the “top” digit, and then add that result to the “bottom” digit. The first addition may result in a double digit answer (10), which means we’d have to use one of the ten extra facts to do the second addition.

In our example, the carry is added to 9, giving 10, and then 10 is added to 2, giving 12. The second addition used the extra “10 + 2 = 12″ fact.

Decimal Addition Showing Intermediate Additions

Binary Addition

Binary addition works the same way as decimal addition, except it uses a different — and much smaller — set of facts. There are only four single-digit facts, or three if you exclude the commutative fact. To handle carries, you need to know two additional facts, 10 + 0 = 10 and 10 + 1 = 11. (Resist the urge to pronounce 10 as ten and 11 as eleven; either call them one-zero and one-one, or two and three.) The latter apply when there’s a carry (always 1) and the “top” digit is 1.

Of the five order-independent facts, two are the same as in decimal (0 + 0 = 0 and 1 + 0 = 1), and one is trivial since it just adds zero (10 + 0 = 10). This leaves only two facts to memorize: 1 + 1 = 10 and 10 + 1 = 11!

To add two binary numbers, proceed as in decimal addition:

• If one or both numbers has a fractional part, line up the radix points.
• Proceeding from right to left, add the digits in each “column,” according to the facts table.
• If the result has two-digits, write down the least significant digit; carry the most significant digit to the next column.

Let’s return to the example in the introduction, 1011.01 + 11.011, this time showing it step-by-step:

Steps of Binary Addition

Eventually you won’t need to write down the intermediate additions, and you might even start doing ‘1 + 1 + 1’ as one step.

Checking the Answer

You can verify the answer by converting the operands to decimal, doing decimal addition, and then converting the decimal answer to binary (of course you can do it that way period, avoiding binary arithmetic entirely!). 1011.01 = 11.25, and 11.011 = 3.375. They sum to 14.625, which is 1110.101, the answer we got using binary addition.

You can also check the answer using my binary calculator.

Discussion

We didn’t need to know the place value of columns, or that a carry represents a power of the base; the algorithm is base-independent. We did however need a base-dependent set of facts, which allowed us to manipulate binary numerals according to the basic rules of binary addition.

By Rick Regan (Copyright © 2008-2012 Exploring Binary)

Binary Addition

In this article, I’ll explain how bigcomp works, and when it applies. Also, I’ll talk briefly about its performance; my informal testing shows that, under the default setting, bigcomp actually worsens performance for some inputs.

When Bigcomp Applies

Bigcomp is enabled by default in strtod() (it can be disabled by “#defining” the variable “NO_STRTOD_BIGCOMP”). It applies to decimal strings longer than STRTOD_DIGLIM significant digits, which by default is “#defined” to 40. For decimal inputs dStr of 40 significant digits or less, the normal approximation check is done, using all significant digits of dStr for its big integer representation d; for inputs longer than 40 digits, only the first 18 digits are used for d. This truncation of d leads to smaller big integers, but it can lead to cases where the approximation check cannot decide between adjacent floating-point values. When that happens, the bigcomp procedure kicks in, through a call to the bigcomp() function.

Specifically, bigcomp() is called when the approximation check cannot decide between the approximation b and its successor, b + u. (As with all my articles in this series, I will only be discussing how strtod() works for round-to-nearest/round-half-to-even rounding.) This happens when the truncated value of d falls between b and b + h; that is, between b and one-half ULP above b.

(In this article, I distinguish between the input dStr and its big integer representation d, which may or may not be a truncation of dStr.)

Overview of the Algorithm

bigcomp() compares the significant digits of dStr to a string representation of the significant digits of b + h. It goes left-to-right, digit by significant digit, stopping when the two strings differ. (In the bulk of cases, the difference occurs between the 16th and 18th digits, although in the worst case the difference can occur at the 54th digit; it depends on how close the decimal number is to b + h.) If dStr is less than b + h, bigcomp() returns b; if dStr is greater than b + h, bigcomp() returns b + u; if dStr equals b + h, either b or b + u is chosen, according to round-half-to-even rounding.

bigcomp() generates the significant digits of b + h using essentially the same algorithm as the dtoa() function, the floating-point to decimal conversion routine co-located in David Gay’s dtoa.c. It uses big integer arithmetic, including big integer division.

bigcomp() starts by creating a fraction that represents the significant digits of b + h (not b + h itself). In this form, it generates a digit by multiplying the numerator of the fraction by ten, and then using integer division to form a quotient and remainder. The quotient represents the next digit; the remainder represents the remaining digits, which are accessed one at a time by repeating this process. Digits are generated until there is a mismatch between dStr and b + h.

The Algorithm, By Example

I’ll break the algorithm into five steps, and illustrate it with two examples. I set STRTOD_DIGLIM to 19 so I could generate shorter examples (both are 20 digits).

Example 1: Check 1.3694713649464322631e-11

In this example, we want to check if dStr = 1.3694713649464322631e-11 is closer to b = 0x1.e1d703bb5749cp-37 or b + u = 0x1.e1d703bb5749dp-37. We do this by comparing dStr to b + h = 0x1.e1d703bb5749c8p-37. (I describe b + h using a hexadecimal floating-point constant, even though it won’t be represented as an IEEE floating-point number. This number has 54 bits, with bit 54 being a 1. It is b with a hexadecimal ‘8’ tacked on at the end.)

1. Represent b + h as an integer times a power of two

We start by representing b as a 53-bit integer times a power of two, as we do in the basic approximation check:

b = bInt x 2bExp
  = 8476617176609948 x 2-89

Since bInt is 53-bit integer, u = 2^bExp, and h is 2^bExp-1. Therefore,

b + h = bInt x 2bExp + 2bExp-1
      = (bInt x 2 + 1) x 2bExp-1
      = (8476617176609948 x 2 + 1) x 2-90
      = 16953234353219897 x 2-90

In code: We create bInt and bExp from b’s underlying floating-point representation: we treat its significand as a 53-bit integer, and subtract 52 from its power of two exponent. We create a big integer from bInt, and then we make it the integer portion of b + h by shifting it left and ‘ORing’ a 1 into its least significant bit. We maintain the exponent of the power of two as a native integer, setting it to bExp-1.

2. Scale b + h so its first significant digit is just to the left of the decimal point

If you were to compute b + h from the expression above, you’d see its value is

0.0000000000136947136494643226044416541235590733258456475063269408565
24743139743804931640625

(We’re not going to do that; I just wanted to show you where we’re headed. By the way, I computed it using PARI/GP.)

What we want to do is scale b + h so that it looks like

1.3694713649464322604441654123559073325845647506326940856524743139743
804931640625

In other words, we want to multiply it by a power of ten so that, when viewed as a decimal number, its significant digits look normalized; 10¹¹ does the trick, effectively shifting b + h left by 11 decimal places:

Scaled b + h = 16953234353219897 x 2-90 x 1011

This “normalization” sets up the digit generation algorithm.

In code: A full-blown floating-point to decimal conversion would need to use logarithms to compute the power of ten; but in this case, we can obtain the power of ten exponent directly from dStr.

3. Combine powers of two

If we decompose the power of ten into a power of two and a power of five, we get

Scaled b + h = 16953234353219897 x 2-90 x 211 x 511

We can now combine powers of two to reduce the size of the resulting big integers.

Scaled b + h = 16953234353219897 x 2-79 x 511

In code: The code is simple: just add the exponents of the power of two from b + h and the power of two from the power of ten scaling factor.

4. Express scaled b + h as a fraction

The negative power of two, 2^-79, makes this a fraction:

Scaled b + h = (16953234353219897 x 511)/279
             = 827794646153315283203125/604462909807314587353088

In code: Compute the powers as big integers, using the absolute values of the exponents; then, multiply the terms out, as appropriate.

5. Generate and check digits

Scaled b + h has a numerator num and denominator den:

num = 827794646153315283203125
den = 604462909807314587353088

The significant digits are generated from these integers using the “repeated multiplication by ten” algorithm; this is what’s done at each iteration (dig stands for ‘digit’, and rem stands for ‘remainder’):

1. dig = num/den (just the integer part of the quotient)
2. rem = num – dig x den
3. num = rem x 10

In other words, at each iteration, a digit is peeled off and the rest of the number is shifted left one decimal place. This step is repeated until a corresponding digit from dStr and b + h don’t match; here are the first three iterations:

1. 827794646153315283203125/den  = 1 rem 223331736346000695850037
2. 2233317363460006958500370/den = 3 rem 419928634038063196441106
3. 4199286340380631964411060/den = 6 rem 572508881536744440292532
4. ...

This stops after step 19; that is, dStr and b + h differ at the 19th digit. Here are the digits of dStr and b + h aligned, to show where they differ:

Digits of dStr  = 13694713649464322631
Digits of b + h = 1369471364946432260444165412355907332584564750632
6940856524743139743804931640625

Since 3 is greater than 0, dStr is closer to b + u than it is to b, so b + u is the answer.

Example 2: Check 9.3170532238714134438e+16

In this example, we want to check if dStr = 9.3170532238714134438e+16 is closer to b = 0x1.4b021afd9f651p+56 or b + u = 0x1.4b021afd9f652p+56. We do this by comparing dStr to b + h = 0x1.4b021afd9f6518p+56.

1. Represent b + h as an integer times a power of two

b = bInt x 2bExp
  = 5823158264919633 x 24

b + h = (bInt x 2 + 1) x 2bExp-1
      = (5823158264919633 x 2 + 1) x 23
      = 11646316529839267 x 23

2. Scale b + h so its first significant digit is just to the left of the decimal point

If you were to compute b + h from the expression above, you’d see its value is

93170532238714136

We want to shift it right 16 decimal places so it looks like

9.3170532238714136

Scaled b + h = 11646316529839267 x 23 x 10-16

3. Combine powers of two

If we decompose the power of ten into a power of two and a power of five, we get

Scaled b + h = 11646316529839267 x 23 x 2-16 x 5-16

Combining powers of two we get

Scaled b + h = 11646316529839267 x 2-13 x 5-16

4. Express scaled b + h as a fraction

The negative powers of two and five make this a fraction:

Scaled b + h = 11646316529839267/(213 x 516)
             = 11646316529839267/1250000000000000

5. Generate and check digits

Scaled b + h has a numerator num and denominator den:

num = 11646316529839267
den = 1250000000000000

Here are the first three iterations of the repeated multiplication by ten algorithm:

1. 11646316529839267/den = 9 rem 396316529839267
2. 3963165298392670/den  = 3 rem 213165298392670
3. 2131652983926700/den  = 1 rem 881652983926700
4. ...

This stops after step 17; that is, dStr and b + h differ at the 17th digit. Here are the digits of dStr and b + h aligned, to show where they differ:

Digits of dStr  = 93170532238714134438
Digits of b + h = 93170532238714136

Since 4 is less than 6, dStr is closer to b than it is to b + u, so b is the answer.

Discussion

strtod() has an extra step, which I’ve not described. It scales the numerator and denominator by a common power of two, which is necessary for its ‘quorem()’ function to compute the quotient and remainder properly. The power of two exponent is determined by a function called ‘dshift().’

If dStr represents a fractional value, step 2 (scale by power of ten) and step 3 (combine powers of two) are technically unnecessary. Without them, step 5 (generate digits) would generate leading zeros, which could be ignored. (Of course, the pre-scaling is more elegant, as it reduces everything to one form.)

The 40-Digit Cutoff Is Too Low

bigcomp() does a lot of work with big integers, which made me wonder why it performed better than the basic approximation check. I ran some tests, and found that it doesn’t perform better — not at the default STRTOD_DIGLIM value of 40! My tests indicate that bigcomp() only performs better starting at somewhere between 80 and 83 digits. (I don’t know why 40 was chosen as the default.) bigcomp() outperforms the basic check for really long inputs; for example, the difference in performance is great at 500 digits, and substantial at 1000 digits.

Let’s look at a 1000-digit input to see the difference in the big integers involved:

dStr = 5.39393324717760434968487301635560570642059608762708229108908
24537754965408622951477281124864324524316631524646435056101989215441
13243584680626358335990970896315873327586087827333915843584231417344
18113410301502801068079096892041478634628161830670245841932386407800
12745803879948980583567722824049594060350336856210415679895612943123
50098271611832693583464221850190462604938332826437767915224673375308
74570398083131101089815004834148500912320646940469604931721018872877
13566715387976682236017863055187831755129573243931855002598816011501
73130862094578091926124182802567489337171671008075773833668718543916
59331687129078431660248055444327878892028873933579030052796118406453
59419797286387142281725295340184928613290361074947899276429780250293
85401861377815243939607988608126354001619571808095966555370152710852
57007714188421395421027623995592231545518458559072989336658831420578
08271261804364513313549623924640140003927894134813865873451090057841
21193123838848522389829029095139051820951352412966707477e+16

Without bigcomp, using the basic approximation check, these are the variables involved:

dS = 539393324717760434968487301635560570642059608762708229108908245
37754965408622951477281124864324524316631524646435056101989215441132
43584680626358335990970896315873327586087827333915843584231417344181
13410301502801068079096892041478634628161830670245841932386407800127
45803879948980583567722824049594060350336856210415679895612943123500
98271611832693583464221850190462604938332826437767915224673375308745
70398083131101089815004834148500912320646940469604931721018872877135
66715387976682236017863055187831755129573243931855002598816011501731
30862094578091926124182802567489337171671008075773833668718543916593
31687129078431660248055444327878892028873933579030052796118406453594
19797286387142281725295340184928613290361074947899276429780250293854
01861377815243939607988608126354001619571808095966555370152710852570
07714188421395421027623995592231545518458559072989336658831420578082
71261804364513313549623924640140003927894134813865873451090057841211
93123838848522389829029095139051820951352412966707477

bS = 539393324717760400000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000

hS = 400000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000

|dS-bS| = 3496848730163556057064205960876270822910890824537754965408
62295147728112486432452431663152464643505610198921544113243584680626
35833599097089631587332758608782733391584358423141734418113410301502
80106807909689204147863462816183067024584193238640780012745803879948
98058356772282404959406035033685621041567989561294312350098271611832
69358346422185019046260493833282643776791522467337530874570398083131
10108981500483414850091232064694046960493172101887287713566715387976
68223601786305518783175512957324393185500259881601150173130862094578
09192612418280256748933717167100807577383366871854391659331687129078
43166024805544432787889202887393357903005279611840645359419797286387
14228172529534018492861329036107494789927642978025029385401861377815
24393960798860812635400161957180809596655537015271085257007714188421
39542102762399559223154551845855907298933665883142057808271261804364
51331354962392464014000392789413481386587345109005784121193123838848
522389829029095139051820951352412966707477

|dS-bS| is less than hS, so b is chosen.

With bigcomp, there are two sets of big integers involved: one set for the basic approximation check (on the truncated value), and one set for bigcomp. Here are the big integers used for the basic approximation check:

dS = 539393324717760434
bS = 539393324717760400
hS = 40
|dS-bS| = 34

Here are the big integers used for bigcomp (not shown is the extra factor of 2⁸ in the numerator and denominator due to dshift()):

num = 13484833117944011
den = 2500000000000000

bigcomp() finds that digit 17 of dStr (3) is less than digit 17 of b + h (4), so it chooses b.

Acknowledgement

Thanks to Mark Dickinson for his commentary on bigcomp, both in Python’s version of David Gay’s dtoa.c, and in private communication.

By Rick Regan (Copyright © 2008-2012 Exploring Binary)

Bigcomp: Deciding Truncated, Near Halfway Conversions

WordPress.com Stats Reports One Million Views

This site went live on November 14, 2008; that works out to an average of over 920 views per day. My decimal/binary converter (an HTML form with an accompanying article) accounts for nearly half of those views (445,700+). Not counting my converter, these are my top ten articles:

1. Java Hangs When Converting 2.2250738585072012e-308 (134,500+ views)
2. PHP Hangs On Numeric Value 2.2250738585072011e-308 (93,600+ views)
3. Binary Dates in 2010 and 2011 (45,600+ views)
4. Ten Ways to Check if an Integer Is a Power Of Two in C (37,400+ views)
5. Converting Floating-Point Numbers to Binary Strings in C (15,200+ views)
6. Why “Volatile” Fixes the 2.2250738585072011e-308 Bug (12,200+ views)
7. FPUpdater Fixes the Java 2.2250738585072012e-308 Bug (10,600+ views)
8. How I Taught Third Graders Binary Numbers (9,800+ views)
9. The Answer is One (Unless You Use Floating-Point) (9,700+ views)
10. A Closer Look at the Java 2.2250738585072012e-308 Bug (9,200+ views)

Five of my top ten articles are related to the PHP and Java floating-point conversion bugs; together they total 260,200+ views.

By Rick Regan (Copyright © 2008-2012 Exploring Binary)

Exploring Binary: One Million Views

Overview

The goal of David Gay’s strtod() is to perform efficient, correctly rounded decimal to floating point conversion. To that end, it minimizes the overhead due to arbitrary-precision arithmetic, which is unavoidable for some conversions. In particular, strtod() avoids division of big integers, as would be required by a full-blown big integer algorithm. Instead of dividing, it multiplies, adds, subtracts, and compares.

The key to the algorithm is the floating-point approximation. It can be compared to the original decimal input to see if it’s too small, correct, or too big (it will be within 10 ULPs or so of correct). The approximation can be adjusted accordingly.

(My discussion will deal only with normal double-precision numbers; strtod() handles subnormal numbers as well.)

Correctly Rounded Conversions

For a floating-point number b to be considered the correctly rounded conversion of a decimal number d, it must be the closest floating-point number to d. That definition sounds simple enough, but there’s a problem: the meaning of “closest” can depend on the value of b. In any case, the tests for closeness are all based on a measure called a ULP — a binary (for our purposes) unit in the last place. A ULP is the gap between adjacent floating-point numbers.

I’ll present the tests as three separate cases, although in code they must be considered together.

Case 1: d is less than one-half ULP away from b

If d is less than one-half of a ULP from b, then b is the correctly rounded conversion of d (this test is ambiguous when b is a power of two — I’ll explain below). Mathematically, this is stated as

b – h < d < b + h ,

where h = 1/2·u, and where u equals one ULP.

Here’s what the relationship looks like:

Correct rounding when d is less than one-half ULP away from b

Equivalently, this relationship can be written as

-h < d – b < h ,

or

|d – b| < h .

Case 2: d is one-half ULP away from b

If d is one-half of a ULP away from b, then b is the correctly rounded conversion of d only if b’s least significant bit is 0 (this test is ambiguous when b is a power of two — I’ll explain below). This is called round-half-to-even rounding.

The relationship is written as

|d – b| = h

Correct rounding when d is one-half ULP away from b

Case 3: b is a power of two

If b is a power of two, two different ULP values come into play, making our half ULP test ambiguous. One ULP below b is half as much as one ULP above b, since a ULP doubles at each increasing power of two boundary. This means that if b is a power of two, it is the correctly rounded conversion of d only if d is no more than one-half ULP above b or no more than one-quarter of that same (upper) ULP below b; here is the inequality:

-h/2 ≤ d – b ≤ h .

This is what it looks like:

Correct rounding when b is a power of two

We use the “≤” signs in the inequality because b is the correctly rounded result even if d = b – h/2 or d = b + h. This is because the significands of (normal) powers of two are zero.

Comparing Using Big Integers

To explain how the algorithm works, I will use this example: d = 1.7864e-45 and b = 0x1.465a72e467d88p-149 (the latter is expressed as a hexadecimal floating-point constant). Here’s how d and b are related:

1.7864 x 10^-45 In Relation to Its Binary Floating-Point Approximation, 0x1.465a72e467d88p-149

The diagram shows that b is the closest floating-point number to d — but how do we check that programmatically, using binary arithmetic? We need to apply the above test: |d – b| < h. (I know a priori this is not a special case — I wanted to keep my example simple.) But there’s a problem: how do we represent d? Answer: as a big integer. If we also represent b and h as big integers, such that they are proportional to d, we can apply the test given by the inequality.

I’ve broken the process into four steps:

Step 1: Represent the significant digits of d and b as integers

The first step is to represent d as an integer times a power of ten, and b as an integer times a power of two:

d = dInt x 10^dExp
b = bInt x 2^bExp

d was already put in this form when the floating point approximation b was computed (but for the purposes of the comparison, all digits of d are used, not just the first 16). b can be put in this form simply, by accessing its underlying floating-point representation: we “unnormalize” its significand into an integer by decreasing its power of two exponent. We represent the significant digits as a 53-bit integer, which means we subtract 52 from its power of two exponent.

For our example, d = 1.7864e-45 is represented as 17864 x 10^-49, and b = 0x1.465a72e467d88p-149, which in binary scientific notation is

1.0100011001011010011100101110010001100111110110001 x 2^-149 ,

is represented (in decimal numerals) as 5741268244528520 x 2^-201.

Step 2: Represent h as a power of two

Since bInt is 53-bit integer, u = 2^bExp, and h is half that: 2^bExp-1. For consistency in my discussion, I’ll substitute the variable hExp for bExp-1, giving

h = 2^hExp .

Step 3: Scale the inequality

To complete the conversion to integers, we need to get rid of any negative powers; we do this by scaling both sides of the inequality. Let’s look at the inequality again:

|d – b| < h .

Here are the three values as we expressed them above:

d = dInt x 10^dExp
b = bInt x 2^bExp
h = 2^hExp

For efficiency, the scaling process will factor out common powers of two from the three powers; to facilitate this, we rewrite 10^dExp as 2^dExp x 5^dExp; now we have

d = dInt x 2^dExp x 5^dExp
b = bInt x 2^bExp
h = 2^hExp

If any of the exponents dExp, bExp, or hExp are negative, we make them nonnegative by scaling. To reflect any scaling, let’s rename d, b, and h to dS, bS, and hS, respectively. The scaled inequality is now written

|dS – bS| < hS ,

and the initial values of dS, bS, and hS are

dS = dInt x 2^dExp x 5^dExp
bS = bInt x 2^bExp
hS = 2^hExp

To show how scaling works, let’s continue with our example above; here are the initial values of the variables:

dS = 17864 x 2^-49 x 5^-49
bS = 5741268244528520 x 2^-201
hS = 2^-202

In this example, all three exponents — dExp, bExp, and hExp — are negative; we scale the values with a common factor S that makes them all nonnegative. The least common factor that does this is S = 2²⁰² x 5⁴⁹; our variables then become

dS = 17864 x 2¹⁵³
bS = 5741268244528520 x 2¹ x 5⁴⁹
hS = 5⁴⁹

Multiplied out, these are

dS = 203970822259994138521801764465966248930731085529088
bS = 203970822259994122305215569213032722473144531250000
hS = 17763568394002504646778106689453125

Step 4: Apply the test

Now all we have to do is check the inequality, |dS – bS| < hS:

|dS – bS| = 16216586195252933526457586554279088, which is less than hS = 17763568394002504646778106689453125. This proves that b is the correct conversion of d.

(In the diagram for this example above, I’ve drawn the relationship of d and b proportionally: (dS – bS)/hS is approximately 0.91, which means that d is about 0.91h above b.)

Computing the Scaling Exponents In Code

In code, the scaling is done by operating on native integer variables representing the exponents of the powers of two and five in each of the three scaled values. After the exponents are created, common factors of two are removed to reduce the size of the resulting big integers. Here is how I coded it in C (this is essentially how it’s done in David Gay’s strtod() and Java’s FloatingDecimal class; I’ve renamed the variables and structured it a bit differently to make it clearer):

//Initialize scaling exponents
dS_Exp2 = 0;
dS_Exp5 = 0;
bS_Exp2 = 0;
bS_Exp5 = 0;
hS_Exp2 = 0;
hS_Exp5 = 0;

//Adjust for decimal exponent
if (dExp >= 0) {
   dS_Exp2 += dExp;
   dS_Exp5 += dExp;
   }
else {
   bS_Exp2 -= dExp;
   bS_Exp5 -= dExp;
   hS_Exp2 -= dExp;
   hS_Exp5 -= dExp;
   }

//Adjust for binary exponent
if (bExp >= 0) {
    bS_Exp2 += bExp;
   }
else {
   dS_Exp2 -= bExp;
   hS_Exp2 -= bExp;
   }

//Adjust for half ulp exponent
if (hExp >= 0) {
   hS_Exp2 += hExp;
   }
else {
   dS_Exp2 -= hExp;
   bS_Exp2 -= hExp;
   }

//Remove common power of two factor from all three scaled values
common_Exp2 = min(dS_Exp2, min(bS_Exp2, hS_Exp2));
dS_Exp2 -= common_Exp2;
bS_Exp2 -= common_Exp2;
hS_Exp2 -= common_Exp2;

Here is how the code executes for our example (I’ll show the state of the exponent variables after adjusting for each exponent):

//After initializing scaling exponents
dS_Exp2 = 0;
dS_Exp5 = 0;
bS_Exp2 = 0;
bS_Exp5 = 0;
hS_Exp2 = 0;
hS_Exp5 = 0;

//After adjusting for decimal exponent (dExp = -49)
dS_Exp2 = 0;
dS_Exp5 = 0;
bS_Exp2 = 49;
bS_Exp5 = 49;
hS_Exp2 = 49;
hS_Exp5 = 49;

//After adjusting for binary exponent (bExp = -201)
dS_Exp2 = 201;
dS_Exp5 = 0;
bS_Exp2 = 49;
bS_Exp5 = 49;
hS_Exp2 = 250;
hS_Exp5 = 49;

//After adjusting for half ulp exponent (hExp = -202)
dS_Exp2 = 403;
dS_Exp5 = 0;
bS_Exp2 = 251;
bS_Exp5 = 49;
hS_Exp2 = 250;
hS_Exp5 = 49;

//After removing common power of two factor of 2^250
dS_Exp2 = 153;
dS_Exp5 = 0;
bS_Exp2 = 1;
bS_Exp5 = 49;
hS_Exp2 = 0;
hS_Exp5 = 49;

After this code executes, the powers are multiplied out as big integers and used to create big integers dS, bS, and hS.

When d, b, and h are already integers

When all three exponents are nonnegative, d, b, and h start out as integer values. In this case, they are scaled down, by their common power of two factor.

Testing the Inequality In Code

In strtod(), the checking centers around a comparison of |d – b| (actually, |b – d|, which is equivalent) to h, which results in three outcomes: |d – b| < h, |d – b| = h, or |d – b| > h. The special cases are worked into the first two outcomes. (Java’s FloatingDecimal class does similar checks, but it is structured differently.)

Regarding the required check at a power of two boundary (comparing d – b to -h/2), strtod() scales d – b by two — to ensure that all values remain integers. This gives the equivalent comparison of 2(d – b) to -h. (In many cases, h is divisible by two; strtod() could just as easily divide h by two instead of multiplying d – b by two. Java’s FloatingDecimal class does the former.)

Adjusting b

strtod() makes theses checks in a loop, and adjusts b as necessary (that is beyond the scope of this article).

strtod()’s Optimization Called “Bigcomp”

strtod() has an optimization called “bigcomp”, which limits the big integer value of d to 40 decimal digits. For decimal inputs of 40 significant digits or less, all 40 are used for d; for decimal inputs greater than 40 digits, only the first 18 are used. Compensating for this truncation of d leads to a more complicated test procedure. (This optimization is on by default; you have to #define variable “NO_STRTOD_BIGCOMP” to turn it off.)

References

• David Gay’s dtoa.c (see the strtod() function).
• David Gay’s Technical Report “Correctly Rounded Binary-Decimal and Decimal-Binary Conversions”.
• William D. Clinger’s paper “How to Read Floating Point Numbers Accurately” (see ‘algorithmR’).
• FloatingDecimal.java (see the doubleValue() method).
• Python’s version of David Gay’s dtoa.c.

By Rick Regan (Copyright © 2008-2012 Exploring Binary)

Using Integers to Check a Floating-Point Approximation

Overview of the Approximation Calculation

Conceptually, the approximation calculation involves two floating-point values: the significant digits of the decimal number, expressed as an integer, and a corresponding power of ten. For example, for the decimal string 163.118762e+109, the calculation is 163118762 x 10¹⁰³; for the decimal string 8.453127e-67, the calculation is 8453127 x 10^-73. If there are more than 16 significant digits, only the first 16 are used. The decimal string 6.2187331579177550499956283e+100, for example, would be treated as 6.218733157917755e+100, making the calculation 6218733157917755 x 10⁸⁵.

The powers of ten can be very small or very large. For double-precision floating-point (I’m considering normal numbers only), they range from 10^-323 (e.g., for DBL_MIN = 2.2250738585072014e-308) to 10³⁰⁸ (e.g., for 1e+308). The goal is to compute them efficiently.

Right-to-Left Binary Exponentiation

An efficient way to compute powers of ten is using right-to-left binary exponentiation. Conceptually, the process is this: given a power bⁿ to compute, write n as a sum of powers of two; then, using the laws of exponents, rewrite bⁿ as a product of b raised to each power of two. For example:

10¹⁰³ = 10^{1 + 2 + 4 + 32 + 64} = 10¹ * 10² * 10⁴ * 10³² * 10⁶⁴ .

Evaluating this expression requires only ten multiplications — six to compute the powers of ten 10² through 10⁶⁴ (by successive squaring), and four to multiply the five required powers together. If 10¹⁰³ were computed naively, it would require 102 multiplications.

As an example for computing negative exponents, consider

10^-73 = 10^{-(1 + 8 + 64)} = 10^-1 * 10^-8 * 10^-64 .

This requires eight multiplications — six to compute the powers of ten 10^-2 through 10^-64, and two to multiply the three required powers together. The naive calculation would do 72 multiplications.

Required Set of Powers

For the range of exponents of double-precision values, the required set of “binary powers” (powers of ten raised to powers of two) for right-to-left binary exponentiation are

• {10^-1, 10^-2, 10^-4, 10^-8, 10^-16, 10^-32, 10^-64, 10^-128, 10^-256} for negative exponents, and
• {10¹, 10², 10⁴, 10⁸, 10¹⁶, 10³², 10⁶⁴, 10¹²⁸, 10²⁵⁶} for positive exponents.

Instead of computing these powers for each conversion, they can be computed once and stored in a table. This reduces the number of multiplications further. In the case of our examples, 10¹⁰³ would now require only four multiplications, and 10^-73 would require only two.

The Approximation Calculation in strtod()

The approximation calculation in strtod() is based on the calculations I just described, but there are a few differences. The significant digits and power of ten are not built in two variables; one variable is initialized with the significant digits, and then the power of ten is factored in incrementally. This allows for a two-phased calculation: the first phase is a fast-path like calculation that handles any component power of ten between 10^-15 and 10¹⁵; the second phase uses binary exponentiation to finish the computation of the power of ten, for component binary powers 10^-16 and smaller and 10¹⁶ and greater.

In the first phase of the calculation, only positive powers of ten, 10¹ to 10¹⁵, are used. When negative powers of ten 10^-1 to 10^-15 are called for, the calculation is changed to the equivalent calculation that divides the significant digits by the corresponding positive power of ten. A separate table contains each of the powers of ten from 10¹ to 10¹⁵, which are looked up directly. The one division or multiplication replaces the up to three multiplications that might be needed for binary exponentiation (10¹ * 10² * 10⁴ * 10⁸). This not only is potentially more efficient, but reduces the number of possibly inexact operations (operations in which the result must be rounded).

Here are two examples: 163.118762e+109, which is parsed as 163118762 x 10¹⁰³, is calculated as 163118762 * 10⁷ * 10³² * 10⁶⁴; 8.453127e-67, which is parsed as 8453127 x 10^-73, is calculated as 8453127/10⁹ * 10^-64.

The Calculation In Code

I’ve taken strtod()‘s approximation code and distilled it into my own C program. The main difference between my code and strtod() is that strtod() handles overflow (for values near DBL_MAX) and underflow (for values that are near DBL_MIN or are subnormal). I’ve also renamed a few variables and made a few other small changes to make the code more readable.

double strtod_approx(char* decimal)
{
 int exponent, sign, i;
 long long sigDigits16;

 static const double fastPosTens[] = {
   1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
   1e10, 1e11, 1e12, 1e13, 1e14, 1e15};
 static const double binaryPosTens[] = {
   1e16, 1e32, 1e64, 1e128, 1e256};
 static const double binaryNegTens[] = {
   1e-16, 1e-32, 1e-64, 1e-128, 1e-256};

 double d;

 parseDecimalApprox(decimal,&sign,&sigDigits16,&exponent);

 d = sigDigits16;
 if (sign)
   d = -d;

 if (exponent > 0)
   {
    if (i = exponent & 0xF)
      d *= fastPosTens[i];
    if (exponent >>= 4)
      for(i = 0; exponent > 0; i++, exponent >>= 1)
        if (exponent & 1)
          d *= binaryPosTens[i];
   }
 else
   if (exponent < 0)
     {
      exponent = -exponent;
      if (i = exponent & 0xF)
        d /= fastPosTens[i];
      if (exponent >>= 4)
        for(i = 0; exponent > 0; i++, exponent >>= 1)
          if (exponent & 1)
            d *= binaryNegTens[i];
     }

 return d;
}

Notes

• This function calls an auxiliary function I wrote called parseDecimalApprox(), which I’ve not shown. It returns the sign of the number, the value of its first 16 significant digits, and its power of ten exponent.
• There are three tables of powers of ten, all precomputed at compile time:
  • fastPosTens: This holds the 15 exactly representable powers of ten needed for the first phase of the calculation (1e0 is not used — it’s just a placeholder).
  • binaryPosTens: This holds the five positive powers of ten used for binary exponentiation.
  • binaryNegTens: This holds the five negative powers of ten used for binary exponentiation.
• If the exponent is zero, no power of ten is factored in (it would be 10⁰ = 1).

Realized in code, you can see why the algorithm to compute the power of ten is called right-to-left binary exponentiation. Converting the exponent to a sum of powers of two is the same as converting it to binary; in code, the exponent is a binary integer. The bits of the exponent are accessed from right to left — least significant to most significant — using the right shift operator; a 1-bit means the corresponding power of two is part of the exponent. The exponent will be a maximum of 9 bits long.

Analysis

There are at most five floating-point multiplications and divisions in this calculation: three multiplications plus one multiplication or division for the worst case powers of ten, and one multiplication to combine the significant digits and the power of ten. If the significant digits start out inexact — either because they are truncated to 16 digits, or they are 16 digits (truncated or otherwise) but represent a value greater than 2⁵³ – 1 — that’s a maximum of six inexact operations.

Using my code above, I measured approximations as far as 10 ULPs off; 1.00431469722921494e-140 is one example. (This is consistent with my analysis of David Gay’s code, where the maximum difference I found was 11 ULPs. Why the difference? It could be due to the code that handles underflow and overflow, or it could be because I didn’t run the same exact tests in both cases.)

Contrast this with my quick and dirty decimal to floating-point conversion program, which does a multiplication or division for every digit. Using that program, I found an example that was 14 ULPs off.

Discussion

A Program That Converts Decimal Values Uses Decimal Values?

The code converts decimal numbers to floating-point, and yet it itself needs decimal numbers converted to floating-point — isn’t that circular? Not really. The compiler’s implementation would have to be different, or at least modified so that the decimal constants are specified in some other way — say as hexadecimal floating-point constants.

This raises another issue: what if the compiler converts these decimal literals incorrectly? After all, this is part of a larger conversion routine whose goal is correct conversion. My guess is that it does not matter. The worst that can happen is that the approximation will be a little less accurate; it will ultimately be corrected — by the arbitrary-precision based reconciliation process in strtod(). (For what it’s worth, I checked — Visual C++ got all of the conversions right.)

Java’s FloatingDecimal Class

The same approximation calculation is done in the doubleValue() method of Java’s FloatingDecimal Class, which is modeled on David Gay’s strtod().

Using Division For the Negative Exponent Case

For the negative exponent case of the binary exponentiation calculation, I tried using division (by positive powers of ten) instead of multiplication (by negative powers of ten). The division-based calculation ran a little slower, which was expected. But interestingly, the approximations were a little less accurate, on average. I don’t know why.

By Rick Regan (Copyright © 2008-2012 Exploring Binary)

strtod()’s Initial Decimal to Floating-Point Approximation