|
These three fixed-point types, coded as C++ template classes, are intended as drop-in replacements for int or double.
August 01, 2005
URL:http://www.drdobbs.com/cpp/fixed-point-arithmetic-types-for-c/184401992
Computers cannot exactly represent the value of fractions such as 8/9, where the divisor is not a multiple of 2 [1]. In fact, the value of the integer expression 8/9 is zero, even though the mathematical quantity is much closer to 1. For computations involving fractions, most developers use floating-point arithmetic, which approximates the value of 8/9 as the decimal fraction 0.88889, a more sensible and accurate choice than the unintuitive result of the integer expression.
The floating-point numeric representation is complex internally. This makes floating-point math slower than integer math, even when floating-point instructions are built into the processor. Fortunately, there are other ways to represent fractional numbers.
A fixed-point number (for instance, the number "123.456") has a radix point with a fixed number of digits to the right, representing a fractional part. We say "radix point" instead of decimal point because there are fixed-point numeric representations for both decimal and binary numbers. In a decimal fixed-point number, the digits to the right of the radix point are decimal digits. In binary fixed-point, the fractional digits are binary bits.
Fixed-point numeric types are implemented using integer arithmetic operations that execute faster than corresponding floating-point operations. The range of a fixed-point type is limited by the range of the integer that represents it internally, whereas floating-point numbers can represent extremely large numbers [2]. For a large class of everyday problems that don't need the astronomical values representable by floating-point numbers, this trade-off is quite reasonable.
I coded C++ template classes implementing three fixed-point types. They are intended as drop-in replacements for int or double. Their precision and range can be tuned to suit different applications. On my Pentium 4, compiled by Visual C++ 7.1 with optimization turned on, the fixed-point types I present here are over twice as fast as equivalent scalar floating-point code. On embedded processors and integer digital signal processors (DSPs) with no floating-point hardware, I expect them to be many times faster than floating-point function libraries.
The fixed-point classes directly implement four arithmetic assignment operators (+=, -=, *=, /=), and two comparison operators (==, <). A trick from the Boost numerics library [3] uses a helper class to define five more arithmetical operators (+, -, *, /, unary -) and the remaining comparison operators (!=, >=, <=, >). The fixed-point classes derive from the helper template class, which takes the fixed-point class as a template argument. At first glance, this technique looks impossibly circular. It works because the compiler doesn't check the template functions until they are instantiated, and both classes have complete definitions at this point. Listing One illustrates how the helper template implements operator+() and operator>=() in terms of the derived class's operator+=() and operator<(), respectively.
// Listing1.h: helper base class
template <class D> class Helper
{
friend D operator+(D const& lhs, D const& rhs)
{ D tmp = lhs; tmp += rhs; return tmp; }
friend bool operator>=(D const& lhs, D const& rhs)
{ return !(lhs < rhs); }
};
class D : public Helper<D>
{
int rep_;
public:
D(int i=0) : rep_(i) {/*empty*/}
D& operator+=(D const& rhs)
{ rep_ += rhs.rep_; return *this; }
bool operator<(D const& rhs) const
{ return rep_ < rhs.rep_; }
};
Chip makers have spent much effort and silicon (perhaps 20 percent of the Pentium 4 die) to improve the performance of floating-point computations. Beating the Pentium's floating-point performance in my fixed-point types was thus a nontrivial challenge. I chose C++ template classes so the compiler would inline method calls and reduce template parameters to compile-time constants.
I benchmarked the fixed-point implementations against floating-point code, and checked the compiler-generated assembly-language output to see if unexpected artifacts were confounding my measurements.
Because different compilers implement different optimizations, users of these fixed-point classes may need to tweak them for best performance on other compilers or non-Intel processors. The Boost numerics library [3] contains an alternate implementation for compilers that don't perform the returned value optimization.
The rounding integer representation produces a more satisfying integer result by rounding the result to the nearest integer, rather than truncating it toward zero. Thus, the value of the expression 8/9 in rounding integer representation is 1, not 0. The speed penalty (on my Pentium 4 using VC7) is only 15 percent versus ordinary integer calculations involving division. Rounding integers have the same range as ordinary integers.
The rounding integer implementation considers an integer as a binary fixed-point number with no bits to the right of the radix point. Integer addition, subtraction, and multiplication produce exact results. They are implemented using the arithmetic operations for the underlying integral type.
Integer division truncates toward zero. To round the result of integer division, the code takes advantage of the C++ modulus (%) operator. The value of a%b is the integer remainder of the expression a/b. If the remainder is greater than or equal to half the divisor, the quotient is rounded up by adding 1. Rather than divide the divisor by 2 (reintroducing the undesirable truncation), I double the remainder (using multiplication, which has an exact result) and compare it to the divisor.
assert(n >= 0 && d > 0); q = n / d; r = n % d; if ((r<<1) >= d) q += 1;
The integer division instruction on most processors produces both a quotient and a remainder. The Visual C++ optimizer efficiently reuses the already-computed remainder instead of recomputing it. On more modest hardware, or less-optimizing compilers, the Standard C runtime function ldiv() can be used to produce a structure containing both quotient and remainder. ldiv() is not an intrinsic function in VC7, so it has a high performance cost.
The signs of the operands present a further challenge in rounding. The signs of the quotient and remainder depend on the signs of the divisor and dividend. Handling all four cases of operand signs generates some code bloat, but not much speed penalty. I have verified this code on Intel-architecture processors. Other processors implement the sign of the remainder differently, so users of this template should verify its correctness with their processor.
I investigated a branch-free rounding computation to improve performance, but there was no observable improvement on the PC. The branch-free code is compact, but was difficult to construct. Listing Two shows the division code with rounding.
// Listing2.h: scaled multiplication and rounded division
template <typename T> inline T scaled_multiplication(T a, T b, int s)
{
T p = rounded_division(a * b, (T)s);
return p;
}// end rounded_multiplication
template <typename T> inline T rounded_division(T n, T d)
{
T q = n / d;
T r = n % d;
# if defined BRANCH_FREE
q += sign(n^d) & neg(abs(d) - (abs(r)<<1) - 1);
# elif defined NO_ROUNDING
# else
if (n < 0)
if (d < 0)
{// num -, den -, quot +, rem -
if ((r << 1) < d)
++q;
}
else // d >= 0)
{// num -, den +, quot -, rem -
if ((-r << 1) > d)
--q;
}
else // n >= 0
if (d < 0)
{// num +, den -, quot -, rem +
if ((-r << 1) < d)
--q;
}
else // d >= 0
{// num +, den +, quot +, rem +
if ((r << 1) > d)
++q;
}
# endif
return q;
} // end rounded_division
In the decimal fixed-point template class, a template parameter N gives the number of decimal digits to the right of the decimal point. The decimal fixed-point type precisely represents all decimal fractions (all hundredths, thousandths, and so on, depending on the template parameter). It is an excellent choice for computations involving currency, because it exactly represents every penny. Whole integer values are represented internally by multiplying the integer by a scale factor that is the Nth power of 10; 102 or 100 for dollars-and-cents calculations. Thus, if N equals 2, the value 1.23 has an exact representation as 123.
Addition and subtraction are implemented as ordinary integer operations on the scaled values.
You multiply decimal fixed-point numbers using integer multiplication to obtain an intermediate product. The intermediate product has twice as many fractional digits as the operands, so if the operands must have the full 32 bits of significance, the intermediate result requires 64 bits. An implementation can use a long long variable for the intermediate product or accept overflow over a wider range of operands. The intermediate product is divided by the scale factor to produce the proper number of fractional digits in the final result. The scaled result is rounded up if the remainder is greater than or equal to half the scale factor.
For division, you multiply the dividend by the scale factor before dividing. This leaves the quotient with the proper number of decimal digits of precision. As with multiplication, the implementation must choose between using a long long intermediate dividend, or accepting overflow on a greater range of dividends. The quotient is rounded up if the remainder is greater than or equal to half the scale factor.
Listing Three shows multiplication and division code for decimal fixed-point numbers. I implemented scaling the intermediate product in terms of rounding division to avoid duplicating the tricky rounding code.
// Listing3.h: scaled multiply and divide
fixed_decimal& operator*=(fixed_decimal const& that)
{// multiplication with rounding
this->rep_ = scaled_multiplication(this->rep_, that.rep_, scale_);
return *this;
}
fixed_decimal& operator/=(fixed_decimal const& that)
{// division with rounding
this->rep_ = rounded_division(this->rep_ * scale_, that.rep_);
return *this;
}
The scale factor is the Nth power of 10. I obtain this number as a compile-time constant using a recursive template with a specialization to stop the recursion. This common template idiom is reproduced in Listing Four.
template <int N> struct exp10
{ enum { value = 10 * exp10<N-1>::value }; };
template <> struct exp10<0>
{ enum { value = 1 }; };
In the binary fixed-point template class, a template parameter N gives the number of binary bits to the right of the radix point. Whole integer values are represented internally by multiplying the integer by a scale factor that is the Nth power of 2; 27 or 128 gives about 2 decimal digits of precision. Thus, the value 1.23 is represented approximately as 1.23*128=157. Note that 157/128 is 1.2265625, which rounds to 1.23 but is not exactly 1.23. The multiplication involved in scaling integers may be efficiently implemented as a left shift, which is faster than integer multiplication on most processors.
Addition and subtraction are implemented as ordinary integer operations on the scaled values.
Multiplication and division work the same as for decimal fixed-point numbers, except the scale factor can, in principle, be more efficiently applied as a left shift. The presented implementation uses the same scaled multiplication function as for the decimal fixed-point type. This leaves it to the optimizer to discover the shortcut multiplication.
An int has 31 bits of magnitude in VC7, so a fixed_binary<7> has 31-7=24 bits of magnitude, providing an integer range of +/-16,777,215. This may not be enough to count Bill Gates's fortune, but it's plenty for many applications.
The fixed-point classes presented here provide useful basic functionality, but much more could be done. I have not yet investigated whether a template specialization for unsigned integral types would improve performance. There could be a function to convert among fixed-point types of differing precision. Multiplication (and division) could be distributed into two 32-bit multiplies (divides), to increase the range over which the product (quotient) can be computed without overflow and without 64-bit arithmetic. This would come at the cost of a second multiply (divide). Overflow could be handled either by throwing an exception, or more intriguingly, by saturation—the replacement of an overflowed result by the largest representable value. I imagine a traits class to parameterize all this functionality.
I have used fixed-point arithmetic to improve the precision of computations for calibration and interpolation in analog to digital conversion, for computing latency in networks, and for collecting performance statistics in servers. Fixed-point math also has well-known uses in graphics and digital-signal processing. All these applications have common properties—a need for speed, a desire for accuracy, limited range of the computed values, and computations involving multiplication and division. For these applications, there is no need to sacrifice speed for precision because fixed-point arithmetic gives both.
[1] Rational numbers such as 8/9 can be represented exactly as a separated numerator and denominator. Doing arithmetic using this representation requires repeated computation of the greatest common divisor, which is not an operation built into processors or optimized by compilers. I didn't code this implementation because I didn't expect a performance advantage over floating point.
[2] The values exactly representable by a fixed-point type are spaced evenly through its range, and include each integer value. This isn't true for floating-point types. Fully half the exactly represented floating-point values are in the interval (-1,1). Exactly represented floating-point values become increasingly sparse as their magnitude increases.
[3] Boost operators, by Dave Abrahams and Jeremy Siek. (http://www.boost.org/libs/ utility/operators.htm).
[4] http://www.guntheroth.com/ contains an updated copy of this paper and the fixed-point templates and support functions.
Kurt Guntheroth holds an MS in computer science from the University of Washington and has been developing commercial software for 20 years. He can be contacted at http:// www.guntheroth.com/.
Terms of Service | Privacy Statement | Copyright © 2024 UBM Tech, All rights reserved.