I was thinking of a soft-fp format on 8-bit MCUs with 40-bit mantissa \$m\$, radix 256, and a 7-bit exponent \$p\$, as in \$m \cdot 256^{p - 67} = m \cdot 2^{8 p - 536}\$ with a separate sign bit.
While it would mean the actual precision would vary in steps, cyclically, between 32 to 40 bits, it would eliminate bit shifts. It would seriously speed up additions and subtractions, but also multiplications. Calculating
c +=
a×
b, if we label the five bytes 0 (least significant) to 4 (most significant), then
c4:c3 += a4×b4 c3:c2 += a4×b3 + a3×b4 c2:c1 += a4×b2 + a3×b3 + a2×b4 c1:c0 += a4×b1 + a3×b2 + a2×b3 + a1×b4 c0:u4 += a4×b0 + a3×b1 + a2×b2 + a1×b3 + a0×b4 u4:u3 += a3×b0 + a2×b1 + a1×b2 + a0×b3 u3:u2 += a2×b0 + a1×b1 + a0×b2 u2:u1 += a1×b0 + a0×b1 u1:u0 += a0×b0While
c3 will be nonzero for nonzero products,
c4 may become zero (whenever
a4×b4 < 256). Then,
u4:u3 products may need to be calculated, to obtain the eight bits "shifted in" (while decrementing the exponent). Typically 15 multiplications and 30 byte-wise additions are needed (not counting carry rippling, which isn't that common).
The problem isn't
speed, but the amount of machine instructions this kind of code generates, in my opinion. Most 8-bitters don't have that much Flash. Even if you implement the above using loops, you end up generating many instructions to compute the addresses. One could write the multiply-add in assembly by hand, but even then it is quite a few instructions. It ends up
bulky, not slow.
On 32-bit MCUs with a fast 32×32-bit multiplication returning the high or low 32 bits of the product, software floating-point is usually not worth it, because you can do multi-word products like above. For example, if you have a fixed-point format with 32 fractional bits, 31 integral bits, and a sign bit, then
c =
a×
b is
c3 = HI(a1×b1) c2 = HI(a1×b0) + HI(a0×b1) + LO(a1×b1) c1 = HI(a0×b0) + LO(a0×b1) + LO(a1×b0) c0 = LO(a0×b0)of which only
c2:c1 forms the result
c. We only need the high bit of
c0 for half-ULP correct rounding, and that is the same as AND operation between the most significant bits of
a0 and
b0. (When that bit and the least significant bit of
c1 are both set, increment
c1, and carry over to
c2 and
c3 if overflow. Otherwise it is ignored.) If
c3 is nonzero, the product overflows.
Thus, for modular arithmetic, you only need 6 32-bit unsigned multiplications (three low, three high), some additions, and a couple of bit checks for half-ULP correct rounding. Addition and subtraction is just two 32-bit operations. Division can be implemented either bitwise, or via 32-bit/32-bit division-and-remainder.
Extending to 96 bits,
c5 = HI(a2×b2) c4 = HI(a2×b1) + HI(a1×b2) + LO(a2×b2) c3 = HI(a2×b0) + HI(a1×b1) + HI(a0×b2) + LO(a2×b1) + LO(a1×b2) c2 = HI(a1×b0) + HI(a0×b1) + LO(a2×b0) + LO(a1×b1) + LO(a0×b2) c1 = HI(a0×b0) + LO(a0×b1) + LO(a1×b0) c0 = LO(a0×b0)For 32 fractional bits, you use
c1,
c2, and
c3. Again, only the most significant bit of
c0 matters, and if
c4 or
c5 are nonzero, the calculation overflows. For modular arithmetic, you only need 13 multiplications.
For 64 fractional bits, you use
c2,
c3, and
c4. For half-ULP correct rounding, you need to calculate the most significant bit of
c1, but it again only matters when the least significant bit of
c2 is also set.
For branchless ripple carry, after adding to e.g.
c2, you need to add zero with carry to
c3 (and if using 64 fractional bits, to
c4), so there are quite a few additions as well, but they tend to be very fast.
13 multiplications can sound like a lot, but 96 bits is a crapton of precision, easily beating even IEEE 754 double-precision floating point, unless you for some strange reason need the range. I rarely do, even in atomic simulations.
From above, we can easily extend into
arbitrary precision numbers, simply by making the number of integral and fractional words –
"limbs" – variable, in both directions.
Arbitrary precision numbers are not floating-point, because the decimal point is always where the two branches (integral and fractional) meet.
Combining both, you can represent any real number as a continuous sequence of (say 32-bit unsigned) words, i.e.
$$v = \left(-1\right)^S \, \sum_{k=0}^{N-1} m_k \, 2^{32 ( k + B) }$$
where \$S\$ is the sign bit, \$B\$ is the position of the least significant bit (in units of 32 bits), \$N\$ is the number of words in the value, and \$m_0\$ through \$m_{N-1}\$ contain the bit representation of the value.
Arbitrary precision fixed- and floating-point types are generally called
multiprecision numbers.
These are scary powerful: no limits on what you can do, really. For example, you can easily use these – say, as provided by the
GNU MPFR (multiple-precision floating-point library) or
GNU MPC (multi-precision library), or
any of the many alternatives to calculate exactly-correct IEEE Binary32 (float), Binary64 (double), Binary128(quad), or whatever precision you need. In fact, most Linux/BSD/Unix machines have
bc installed, which is a command-line arbitrary-precision calculator. For example, if you want the first 101 decimal digits of pi, you use the arctan(1) function:
echo 'scale=100; 4*a(1)' | bc -l(Yes,
bc shortens arctan() to
a(), sin to
s() cosine to
c(), natural logarithm to
l(), exponential to
e(), when the math library is enabled via the
-l option.)
Do let me know if you see any bugs in above. I've been buggy lately!
Home
Help
Search
Login
Register
