Losing precision with dc

Question

I want to use dc to handle some base 16 numbers with hexadecimal points but I'm running into precision problems. For example, below I'm multiplying F423F.FD by 100, both hex. The expected answer is F423FFD, instead it is giving F423FFA.E1, close but not accurate enough even after rounding.

$ dc
16 d i o F423F.FD 100 * p
F423FFA.E1

I read that dc was an unlimited precision calculator, and this is not a large number by any means. Is there something I'm doing wrong?

Thank for your answers. Given the issues with dc, I bit the bullet and wrote my own parser for real numbers in other bases. If anyone is interested in the code, I can post it here.

Answer 1

Expressed as decimal (using dc to convert), this corresponds to 999999.98 (rounded down) × 256, i.e. 255999994.88, which is F423FFA.E1 in hexadecimal.

So the difference comes from dc’s rounding behaviour: instead of calculating 256 × (999999 + 253 ÷ 256), which would give 255999997, it rounds 253 ÷ 256 down and multiplies the result.

dc is an arbitrary precision calculator, which means it can calculate to any precision you want, but you have to tell it what that is. By default, its precision is 0, meaning division produces integer values only, and multiplication uses the number of digits in the input. To set the precision, use k (and bear in mind that precision is always expressed in decimal digits, regardless of the input or output radix):

10 k
16 d i o
F423FFD 100 / p
F423F.FD0000000
100 * p
F423FFD.000000000

(8 digits’ precision would be sufficient since that’s what you need to represent 1 ÷ 256 in decimal.)

Answer 2

Note that just printing the original number shows it is rounded:

$ dc <<<'16 d i o F423F.FD p'
F423F.FA

You can get around it by adding lots of trailing zeroes for more precision:

$ dc <<<'16 d i o F423F.FD000000 100 * p'
F423FFD.0000000

Answer 3

The Issue

The problem is the way in which dc (and bc) understand numeric constants.
For example, the value (in hex) 0.3 (divided by 1) gets transformed into a value close to 0.2

$ dc <<<"20k 16 d i o 0.3 1 / p"
.199999999999999999999999999

In fact, the plain constant 0.3 also gets changed:

$ dc <<<"20 k 16 d i o     0.3     p"
.1

It seems that it is in an odd way, but it is not (more later).
Adding more zeros makes the answer approach the correct value:

$ dc <<<"20 k 16 d i o     0.30     p"
.2E

$ dc <<<"20 k 16 d i o     0.300     p"
.2FD

$ dc <<<"20 k 16 d i o     0.3000     p"
.3000

The last value is exact and will stay exact no matter how may more zeros are added.

$ dc <<<"20 k 16 d i o     0.30000000     p"
.3000000

The problem is also present in bc:

$ bc <<< "scale=20; obase=16; ibase=16;    0.3 / 1"
.19999999999999999

$ bc <<< "scale=20; obase=16; ibase=16;    0.30 / 1"
.2E147AE147AE147AE

$ bc <<< "scale=20; obase=16; ibase=16;    0.300 / 1"
.2FDF3B645A1CAC083

$ bc <<< "scale=20; obase=16; ibase=16;    0.3000 / 1"
.30000000000000000

One digit per bit?

The very non intuitive fact for floating point numbers is that the number of digits required (after the dot) is equal to the number of binary bits (also after the dot). A binary number 0.101 is exactly equal to 0.625 in decimal. The binary number 0.0001110001 is (exactly) equal to 0.1103515625 (ten decimal digits)

$ bc <<<'scale=30;obase=10;ibase=2; 0.101/1; 0.0001110001/1'; echo ".1234567890"
.625000000000000000000000000000
.110351562500000000000000000000
.1234567890

Also, for a floating point number like 2^(-10), which in binary has only one (set) bit:

$ bc <<<"scale=20; a=2^-10; obase=2;a; obase=10; a"
.0000000001000000000000000000000000000000000000000000000000000000000
.00097656250000000000

Has the same number of binary digits .0000000001 (10) as decimal digits.0009765625 (10). That may not be the case in other bases but base 10 is the internal representation of numbers in both dc and bc and therefore is the only base that we really need to care about.

The math proof is at the end of this answer.

bc scale

The number of digits after the dot could be counted with the built-in function scale() form bc:

$ bc <<<'obase=16;ibase=16; a=0.FD; scale(a); a; a*100'
2
.FA
FA.E1

As shown, 2 digits is insufficient to represent the constant 0.FD.

And, also, just counting the number of characters used after the dot is a very incorrect way to report (and use) the scale of the number. The scale of a number (in any base) should calculate the number of bits needed.

Binary digits in a hex float.

As it is known, each hex digit use 4 bits. Therefore, each hex digit after the decimal dot require 4 binary digits, which due to the (odd?) fact above also require 4 decimal digits.

Therefore, a number like 0.FD will require 8 decimal digits to be represented correctly:

$ bc <<<'obase=10;ibase=16;a=0.FD000000; scale(a);a;a*100'
8
.98828125
253.00000000

Add zeros

The math is straightforward (for hex numbers):

• Count the number of hex digits (h) after the dot.
• Multiply h by 4.
• Add h×4 - h = h × (4-1) = h × 3 = 3×h zeros.

In shell code (for sh):

a=F423F.FD
h=${a##*.}
h=${#h}
a=$a$(printf '%0*d' $((3*h)) 0)
echo "$a"

echo "obase=16;ibase=16;$a*100" | bc

echo "20 k 16 d i o $a 100 * p" | dc

Which will print (correctly both in dc and bc):

$  sh ./script
F423F.FD000000
F423FFD.0000000
F423FFD.0000000

Internally, bc (or dc) could make the number of digits required match the number calculated above (3*h) to convert hex floats to the internal decimal representation. Or some other function for other bases (assuming the number of digits is finite in relation to base 10 (internal of bc and dc) in such other base). Like 2ⁱ (2,4,8,16,...) and 5,10.

posix

The posix specification states that (for bc, which dc is based on):

Internal computations shall be conducted as if in decimal, regardless of the input and output bases, to the specified number of decimal digits.

But "… the specified number of decimal digits." could be understood to be " … the needed number of decimal digits to represent the numeric constant" (as described above) without affecting the "decimal internal computations"

Because:

bc <<<'scale=50;obase=16;ibase=16; a=0.FD; a+1'
1.FA

bc is not really using 50 ("the specified number of decimal digits") as set above.

Only if divided it is converted (still incorrectly as it uses an scale of 2 to read the constant 0.FD before expanding it to 50 digits):

$ bc <<<'scale=50;obase=16;ibase=16; a=0.FD/1; a'
.FAE147AE147AE147AE147AE147AE147AE147AE147A

However, this is exact:

$ bc <<<'scale=50;obase=16;ibase=16; a=0.FD000000/1; a'
.FD0000000000000000000000000000000000000000

Again, reading numeric strings (constants) should use the correct number of bits.

---

Math proof

In two steps:

A binary fraction can be written as a/2ⁿ

A binary fraction is a finite sum of negative powers of two.

For example:

= 0.00110101101 = 
= 0. 0     0      1     1      0      1     0      1      1     0       1

= 0 + 0×2^-1 + 0×2^-2 + 1×2^-3 + 1×2^-4 + 0×2^-5 + 1×2^-6 + 0×2^-7 + 1×2^-8 + 1×2^-9 + 0×2^-10 + 1×2^-11

= 2^-3 + 2^-4 + 2^-6 + 2^-8 + 2^-9 + 2^-11 = (with zeros removed)

In a binary fraction of n bits, the last bit has a value of 2^-n, or 1/2ⁿ. In this example: 2^-11 or 1/2¹¹.

= 1/2³ + 1/2⁴ + 1/2⁶ + 1/2⁸ + 1/2⁹ + 1/2¹¹ = (with inverse)

In general, the denominator could become 2ⁿ with a positive numerator exponent of two. All terms can then be combined into a single value a/2ⁿ. For this example:

= 2⁸/2¹¹ + 2⁷/2¹¹ + 2⁵/2¹¹ + 2³/2¹¹ + 2²/2¹¹ + 1/2¹¹ = (expressed with 2¹¹)

= (2⁸ + 2⁷ + 2⁵ + 2³ + 2² + 1 ) / 2¹¹ = (extracting common factor)

= (256 + 128 + 32 + 8 + 4 + 1) / 2¹¹ = (converted to value)

= 429 / 2¹¹

Every Binary Fraction Can Be Expressed As b/10ⁿ

Multiply a/2ⁿ by 5ⁿ /5ⁿ, getting (a×5ⁿ)/(2ⁿ×5ⁿ) = (a×5ⁿ)/10ⁿ = b/10ⁿ, where b = a×5ⁿ. It has n digits.

For the example, we have:

(429·5¹¹)/10¹¹ = 20947265625 / 10¹¹ = 0.20947265625

It has been shown that every binary fraction is a decimal fraction with the same number of digits.