Document number: P0037R5
Date: 2018-05-07
Reply-to: John McFarlane, fixed-point@john.mcfarlane.name
Audience: SG6, SG14, LEWG
This proposal introduces a system for performing fixed-point arithmetic using integral types.
Floating-point types are an exceedingly versatile and widely supported method of expressing real numbers on modern architectures.
However, there are certain situations where fixed-point arithmetic is preferable:
Integer types provide the basis for an efficient representation of binary fixed-point real numbers. However, laborious, error-prone steps are required to normalize the results of certain operations and to convert to and from fixed-point types.
A set of tools for defining and manipulating fixed-point types is proposed. These tools are designed to make work easier for those who traditionally use integers to perform low-level, high-performance fixed-point computation. They are composable such that a wide range of trade-offs between speed, accuracy and safety are supported.
This proposal is a pure library extension. It does not require changes to any standard classes or functions. It adds several new class and function templates to new header file, <fixed_point>. Some optional deduction guides, member functions and operator overloads rely on types proposed in [P0828] and [P1050].
The design is driven by the following aims in roughly descending order:
More generally, the aim of this proposal is to contain within a single API all the tools necessary to perform fixed-point arithmetic. The design facilitates a wide range of competing compile-time strategies for avoiding overflow and precision loss, but implements only the simplest by default. Similarly, orthogonal concerns such as run-time overflow detection and rounding modes are deferred to the underlying integer types used as storage.
Fixed-point numbers are specializations of
template <class Rep, int Exponent, int Radix>
class fixed_point;
where the template parameters are described as follows.
Rep Type Template ParameterThis parameter indicates the integral type used as storage. Fundamental integral types other than bool are ideal choices but any suitably integer-like type can be used.
Other than scale, the characteristics of fixed_point<Rep> are the characteristics of Rep including:
Exponent and Radix Non-Type Template ParametersThe radix (or base) of a fixed-point type is typically two to denote scaling by powers of two. In financial applications, accurate representation of decimal fractions requires a radix of ten. Thus while Radix can be any number greater than one, 2 is the default.
The exponent of a fixed-point type is the equivalent of the exponent field in a floating-point type and shifts the stored value by the requisite number of digits necessary to produce the desired range. The default value of Exponent is zero, giving fixed_point<T> the same range as T. By far the most common use of fixed-point is to store values with fractional digits. Thus, the exponent is typically a negative value.
The resolution of an instantiation of fixed_point is
pow(Radix, Exponent)
and the minimum and maximum values are
std::numeric_limits<Rep>::min() * pow(Radix, Exponent)
and
std::numeric_limits<Rep>::max() * pow(Radix, Exponent)
respectively.
Any usage that results in values of Exponent which lie outside the range, (INT_MIN / 2, INT_MAX / 2), may result in undefined behavior and/or overflow or underflow. This range of exponent values is far in excess of the largest built-in floating-point type and should be adequate for all intents and purposes.
While effort is made to ensure that significant digits are not lost during conversion, no effort is made to avoid rounding errors. Whatever would happen when converting to and from Rep largely applies to fixed_point objects also. For example:
fixed_point<int, -1>{.499}==0
...equates to true and is considered an acceptable rounding error.
Rep ValueIt is sometimes necessary to read from and write to the Rep value contained in a fixed_point<Rep> object. This is supported through numeric traits, to_rep and from_rep respectively. These traits are described in paper, [P0675].
constexpr auto a = from_rep<fixed_point<int, -8>>()(320);
static_assert(a == 1.25);
constexpr auto b = to_rep(a);
static_assert(b == 320); // 1.25*(1<<8)
The type of a fixed_point object can be deduced by an integer initializer:
auto a = fixed_point(0ul);
static_assert(is_same_v<decltype(a), fixed_point<unsigned long, 0>>);
It can also be deduced with an integral constant of type constant (described in [P0827]):
constexpr auto b = fixed_point(constant<0xFF00000000L>{});
static_assert(is_same_v<decltype(b), const fixed_point<int, 32>>);
static_assert(to_rep(b) == 0xFF);
For Exponent, the highest value which does not incur data loss is used. This minimizes the required range of the underlying integer value which reduces the likelihood of out-of-range errors during arithmetic operations. For Rep, a fundamental integer type of int width is preferred unless a wider type is required.
Any arithmetic, comparison, logic and bitwise operators that might be applied to integer types can also be applied to fixed-point types. A guiding principle of operator overloads is that they perform as little run-time computation as is practically possible.
With the exception of shift and comparison operators, binary operators can take any combination of:
numeric_limits is specialized.Assuming a binary operation, @, in the form
auto R = S @ T;
where S is of type fixed_point<RepS, ExponentS, 2> and T is a numeric type — possibly another fixed_point instantiation — then result, R, of the operation is determined as follows.
If T is a floating-point type, Float, then S is cast to Float and a floating-point operation takes place, e.g.:
auto a = fixed_point<long long>(3) + 4.f;
static_assert(is_same_v<decltype(a), decltype(3.f + 4.f)>);
If T is a constant of integer type, Integer, then:
a. If the operator is bitwise left shift (<<), then the result is fixed_point<RepS, ExponentS+T::value> with the shift operator applied.
b. If the operator is bitwise right shift (>>), then the result is fixed_point<RepS, ExponentS-T::value> with the shift operator applied.
c. Otherwise, T is converted to fixed_point(T{}), e.g.
auto b = fixed_point(200U) - constant<100L>{};
static_assert(is_same_v<decltype(b), decltype(fixed_point<unsigned>(200) - fixed_point<int>(100))>);
and proceeding rule #4 subsequently applies.
If T is an integer type, Integer, then:
a. If the operator is bitwise shift (<< or >>), then the result is type S with the shift operator applied.
b. Otherwise, T is cast to fixed_point<Integer, 0>, e.g.
auto c = fixed_point<>(5) * 6ul;
static_assert(is_same_v<decltype(c), decltype(fixed_point<>(5) * fixed_point<unsigned long>(6))>);
and proceeding rule #4 subsequently applies.
If T is type, fixed_point<RepT, ExponentT, 2>, then:
a. If the operator is multiplication (*), then the result is fixed_point<decltype(RepS*RepT), ExponentS+ExponentT>, e.g.:
constexpr auto d = fixed_point<uint8_t, -7>{1.25} * fixed_point<uint8_t, -3>{8};
static_assert(is_same_v<decltype(d), const fixed_point<int, -10>>);
static_assert(d == 10);
b. If the operator is division (/), then the result is fixed_point<decltype(RepS/RepT), ExponentS-ExponentT>, e.g.:
constexpr auto e = fixed_point<short, -5>{1.5} / fixed_point<short, -3>{2.5};
static_assert(is_same_v<decltype(e), const fixed_point<int, -2>>);
static_assert(e == .5);
c. If the operator is modulo (%), then the result is fixed_point<decltype(RepS%RepT), ExponentS>, e.g.:
constexpr auto f = fixed_point<short, -5>{1.5} % fixed_point<short, -3>{2.5};
static_assert(is_same_v<decltype(f), const fixed_point<int, -5>>);
static_assert(f == .25);
d. If the operator is addition (+) or subtraction (-), then the operand with the greater exponent is converted such that its exponent matches the other operands' exponent. Then the result is fixed_point<decltype(RepS@RepT), min(ExponentS,ExponentT)>, eg.:
constexpr auto g = fixed_point<int8_t, -2>{12.5} - fixed_point<short, 0>{8};
static_assert(is_same_v<decltype(g), const fixed_point<int, -2>>);
static_assert(g == 4.5);
e. If the operator is comparison (==, !=, <, >, <= or >=), then the operand with the greater exponent is converted such that its exponent matches the other operands' exponent. Then the result is decltype(RepS@RepT), eg.:
constexpr auto h = fixed_point<int8_t, -2>{12.5} <= fixed_point<short, 0>{8};
static_assert(is_same_v<decltype(h), const bool>);
static_assert(h == false);
(See section, Extended Comparison Range, for additional details.)
f. If the operator is bitwise or (|) or xor (^) then the same rules as addition (+) are applied.
g. If the operator is bitwise and (&) then the same rules as bitwise or (|) are applied except that the greater — not less — exponent is preferred.
Some details have been left out for brevity. Unary operators are supported. Some minor variations occur when S is not fixed_point and T is fixed_point. Rules for bit-shifting values where Radix!=2 do not necessarily involve a different result type. Other binary operations involving different radixes produce a return type which is optimized to contain the precise result with the minimum value stored in Rep and the minimum viable value of Radix.
The complete set of rules may appear to be large and complex. However, this mostly reflects the existing complexity in the behavior of arithmetic types. Relatively few design principles govern these rules:
fixed_point<T, 0> should follow the same behavior as T to the greatest extend practical, reflecting the facts that: a) all integers are fixed-point — rather than floating-point types and b) integer arithmetic generally provides the best efficiency and performance characteristics.safe_integer and elastic_integer types discussed in [P0554].The fixed_point division operator, /, performs the least work possible and, combined with the modulo operator, %, produces lossless results. However, it behaves very differently from floating-point division and is likely to be a source of surprises for some users.
In particular, the choice of quotient type can have a dramatic effect on precision. If, for example, the dividend and divisor have the same Exponent and Radix, then the quotient's Exponent will be zero and all fractional digits will be dropped. In contrast to floating-point division, the choice of Exponent cannot be tailored to the result.
For this reason, the fractional type [P1050] is provided in order to facilitate two important use cases.
Firstly a 'sane default' result type can be calculated automatically in line with the formula detailed in P0106 [P0106]. Here, a deduction guide does the work of determining that a division involving a dividend with 31 integer digits should result in a quotient with 31 fractional digits:
constexpr auto i = fixed_point{fractional{1, 3}};
static_assert(i == 0.333333333022892475128173828125L);
static_assert(is_same_v<decltype(i), const fixed_point<int64_t, -31>>);
Alternatively, the user can forgo CTAD and choose the template parameters explicitly:
constexpr auto j = fixed_point<int, -16>{fractional{1, 3}};
static_assert(j == 0.3333282470703125);
static_assert(is_same_v<decltype(j), const fixed_point<int, -16>>);
RepUsing built-in integral types as the default underlying representation minimizes certain costs:
However, this choice also brings with it many of the deficiencies of built-in types. For example:
long long int.The effort involved in addressing these deficiencies is non-trivial and on-going (for example [P0105]). As solutions are made available, it should become easier to define custom integral types which address concerns surrounding robustness and correctness. How to combine such numeric types is the topic of [P0554], Composition of Arithmetic Types.
Of particular note is the elastic_integer type detailed in [P0828]. When used in combination with fixed_point, the resultant composite type is able to avoid a large proportion of the out-of-range errors associated with fixed-point arithmetic while avoiding expensive run-time overflow checks.
For a type to be suitable as parameter, Rep, of fixed_point, it must meet the following requirements:
numeric_limitsmake_signed and make_unsignednum_digits and set_num_digits,to_rep, from_rep and from_value.Note that make_signed and make_unsigned cannot be specialized for custom types. Unless this rule can be relaxed, some equivalent mechanism must be introduced in order for custom types to be used with fixed_point<>. One possibility is the addition of numeric_limits<>::signed and numeric_limits<>::unsigned type aliases.
The following function, magnitude, calculates the magnitude of a 3-dimensional vector.
template<class Fp>
constexpr auto magnitude(Fp x, Fp y, Fp z)
{
return sqrt(x*x+y*y+z*z);
}
And here is a call to magnitude.
auto m = magnitude(
fixed_point<uint16_t, -12>(1),
fixed_point<uint16_t, -12>(4),
fixed_point<uint16_t, -12>(9));
// m === fixed_point<uint32_t, -24>{9.8994948863983154}
namespace std {
template <class Rep, int Exponent, int Radix> class fixed_point;
// for each unary arithmetic, comparison, logic and bitwise operator, @
template <class RhsRep, int RhsExponent, int RhsRadix>
constexpr auto operator@(
const fixed_point<RhsRep, RhsExponent, int RhsRadix> & rhs);
// for each binary arithmetic, comparison, logic and bitwise operator, @
template <class LhsRep, int LhsExponent, int LhsRadix, class RhsRep, int RhsExponent, int RhsRadix>
constexpr auto operator@(
const fixed_point<LhsRep, LhsExponent, int LhsRadix> & lhs,
const fixed_point<RhsRep, RhsExponent, int RhsRadix> & rhs);
template <class LhsRep, int LhsExponent, int LhsRadix, class RhsFloat, typename = _impl::enable_if_t<numeric_limits<RhsFloat>::is_iec559>>
constexpr auto operator@(
const fixed_point<LhsRep, LhsExponent, int LhsRadix> & lhs,
const RhsFloat & rhs);
template <class LhsFloat, class RhsRep, int RhsExponent, int RhsRadix, typename = _impl::enable_if_t<numeric_limits<LhsFloat>::is_iec559>>
constexpr auto operator@(
const LhsFloat & lhs,
const fixed_point<RhsRep, RhsExponent, int RhsRadix> & rhs);
template <class LhsRep, int LhsExponent, int LhsRadix, class RhsInteger, typename = _impl::enable_if_t<numeric_limits<RhsInteger>::is_integer>>
constexpr auto operator@(
const fixed_point<LhsRep, LhsExponent, int LhsRadix> & lhs,
const RhsInteger & rhs);
template <class LhsInteger, class RhsRep, int RhsExponent, int RhsRadix, typename = _impl::enable_if_t<numeric_limits<LhsInteger>::is_integer>>
constexpr auto operator@(
const LhsInteger & lhs,
const fixed_point<RhsRep, RhsExponent, int RhsRadix> & rhs);
template <class LhsRep, int LhsExponent, int LhsRadix, auto RhsValue>
constexpr auto operator@(
const fixed_point<LhsRep, LhsExponent, int LhsRadix> & lhs,
constant<RhsValue>);
template <auto LhsValue, class RhsRep, int RhsExponent, int RhsRadix>
constexpr auto operator@(
constant<LhsValue>,
const fixed_point<RhsRep, RhsExponent, int RhsRadix> & rhs);
// for each arithmetic, comparison, logic and bitwise compound assignment operator, @=
template <class LhsRep, int LhsExponent, int LhsRadix, class RhsRep, int RhsExponent, int RhsRadix>
constexpr auto operator@=(
fixed_point<LhsRep, LhsExponent, int LhsRadix> & lhs,
const fixed_point<RhsRep, RhsExponent, int RhsRadix> & rhs);
template <class LhsRep, int LhsExponent, int LhsRadix, class RhsFloat, typename = _impl::enable_if_t<numeric_limits<RhsFloat>::is_iec559>>
constexpr auto operator@=(
fixed_point<LhsRep, LhsExponent, int LhsRadix> & lhs,
const RhsFloat & rhs);
template <class LhsFloat, class RhsRep, int RhsExponent, int RhsRadix, typename = _impl::enable_if_t<numeric_limits<LhsFloat>::is_iec559>>
constexpr auto operator@=(
LhsFloat & lhs,
const fixed_point<RhsRep, RhsExponent, int RhsRadix> & rhs);
template <class LhsRep, int LhsExponent, int LhsRadix, class RhsInteger, typename = _impl::enable_if_t<numeric_limits<RhsInteger>::is_integer>>
constexpr auto operator@=(
fixed_point<LhsRep, LhsExponent, int LhsRadix> & lhs,
const RhsInteger & rhs);
template <class LhsInteger, class RhsRep, int RhsExponent, int RhsRadix, typename = _impl::enable_if_t<numeric_limits<LhsInteger>::is_integer>>
constexpr auto operator@=(
LhsInteger & lhs,
const fixed_point<RhsRep, RhsExponent, int RhsRadix> & rhs);
template <class LhsRep, int LhsExponent, int LhsRadix, auto RhsValue>
constexpr auto operator@=(
fixed_point<LhsRep, LhsExponent, int LhsRadix> & lhs,
constant<RhsValue>);
template <auto Value>
fixed_point(::cnl::constant<Value>)
-> /* ... */;
template <class Integer>
fixed_point(Integer)
-> fixed_point<Integer, 0>;
}
fixed_point<> Class Templatetemplate <class Rep = int, int Exponent = 0, int Radix = 0>
class fixed_point
{
public:
using rep = Rep;
using radix = Radix;
constexpr static int exponent;
constexpr fixed_point();
template<class FromRep, int FromExponent, int FromRadix>
constexpr fixed_point(fixed_point<FromRep, FromExponent> const&);
template<CNL_IMPL_CONSTANT_VALUE_TYPE Value>
constexpr fixed_point(constant<Value>);
template<class S, _impl::enable_if_t<numeric_limits<S>::is_integer, int> Dummy = 0>
constexpr fixed_point(S const&);
template<class S, _impl::enable_if_t<numeric_limits<S>::is_iec559, int> Dummy = 0>
constexpr fixed_point(S);
template<class Numerator, class Denominator>
constexpr fixed_point(const fractional<Numerator, Denominator>&);
template<class S, _impl::enable_if_t<numeric_limits<S>::is_integer, int> Dummy = 0>
constexpr fixed_point& operator=(S);
template<class S, _impl::enable_if_t<numeric_limits<S>::is_iec559, int> Dummy = 0>
constexpr fixed_point& operator=(S);
template<class FromRep, int FromExponent, int FromRadix>
constexpr fixed_point& operator=(fixed_point<FromRep, FromExponent> const&);
template<class S, _impl::enable_if_t<numeric_limits<S>::is_integer, int> Dummy = 0>
explicit constexpr operator S() const;
template<class S, _impl::enable_if_t<numeric_limits<S>::is_iec559, int> Dummy = 0>
explicit constexpr operator S() const;
};Because the aim is to provide an alternative to existing arithmetic types which are supported by the standard library, it is conceivable that a future proposal might specialize existing class templates and overload existing functions.
Possible candidates for overloading include the functions defined in <cmath> and a templated specialization of numeric_limits. A new type trait, is_fixed_point, would also be useful.
While fixed_point is intended to provide drop-in replacements to existing built-ins, it may be preferable to deviate slightly from the behavior of certain standard functions. For example, overloads of functions from <cmath> will be considerably less concise, efficient and versatile if they obey rules surrounding error cases. In particular, the guarantee of setting errno in the case of an error prevents a function from being defined as pure. This highlights a wider issue surrounding the adoption of the functional approach and compile-time computation that is beyond the scope of this document.
One suggested addition is a specialization of std::complex. This would take the form:
template<class Rep, int Exponent, int Radix>
class complex<fixed_point<Rep, Exponent, Radix>>;
This type's arithmetic operators would differ from existing specializations because fixed_point<> operators often return results of a different type to their operands. Hence signatures such as
template<class T>
complex<T> operator*( const complex<T>& lhs, const complex<T>& rhs);
would need to be replaced with:
template<class T>
auto operator*( const complex<T>& lhs, const complex<T>& rhs);
It is quite possible that the order of the three template parameters of fixed_point are the wrong way around and should be
template<int Exponent=0, Rep=int, int Radix=2>
class fixed_point;
or:
template<int Exponent=0, int Radix=2, Rep=int>
class fixed_point;
While being able to express fixed_point<Integer> reinforces the fact that integers are merely fixed-point types with radix point fixed at zero, the second parameter, Exponent is perhaps more important to specify.
Comparison operations between two fixed_point operands require that they both have the same exponent. When they do not, conversion takes place to ensure they do. Unfortunately, if the difference in exponents is too great, the conversion may cause an out-of-bounds condition.
However, where two operands have bits whose values are in ranges that do not overlap, it may not be necessary to perform a conversion which results in out-of-range results: a result that ensures they continue to not overlap may be sufficient. For example,
static_assert(fixed_point<uint8_t, 0>{0} < fixed_point<uint8_t, 128>{4.e38});
requires that the right-hand operand be converted to fixed_point<uint8_t, 0>. This will result in the underlying integer being scaled up by 1000 bits, resulting in undefined behavior and/or a flushed value. But in this case, it only needed to be scaled by 8 bits in order for none of its bit values to overlap with those of the left-hand operand.
The possibility of scaling by other than an exponent has been discussed previously [future]. It is likely that fixed_point is just one example of a more general numeric type, scaled:
// a numeric type
template<class Rep, class Scale>
class scaled;
// a member of concept, "scale"
template<int Exponent, int Radix>
class power;
// fixed_point as a numeric type
template<class Rep, int Exponent, int Radix>
using fixed_point = scaled<Rep, power<Exponent, Radix>>;
This approach separates two concerns currently combined in fixed_point:
In this dichotomy, scaled holds onto a member variable of type, Rep, while power is stateless and provides an API which converts to and from the Rep type through function parameters. power also provides arithmetic support.
Alternatives to power would include:
ratio (with additional arithmetic support)This design is a matter for future experimentation but it has immediate consequences for the API of fixed_point.
Does it make sense to allow binary operations which take, say, a base-2 and a base-10 number? The answer is relatively straight-forward when one considers that all base-2 numbers can be expressed using base-10 numbers. What about a base-2 and a base-3 number? At this point, we may need to convert them to a base-6 number to proceed.
Next, how do the exponents interact in situations when the radixes are different? For example, when adding fixed_point<int, -2, 2> and fixed_point<int, -1, 4>, is the result the former or the latter? They are computationally equivalent because they both represent units of 0.25.
The likely solution is to choose a result type with the minimum radix and then the minimum exponent necessary in order to be able to represent all possible values. However, this may result in a set of operators which are surprising to the user. Thus is it tempting to simply forbid inter-radix operations. (Note: a similar problem is faced by common_type(chrono::duration).)
Many examples of fixed-point support in C and C++ exist. While almost all of them aim for low run-time cost and expressive alternatives to raw integer manipulation, they vary greatly in detail and in terms of their interface.
One especially interesting dichotomy is between solutions which offer a discrete selection of fixed-point types and libraries which contain a continuous range of exponents through type parameterization.
One example of the former is found in proposal N1169 [N1169], the intent of which is to expose features found in certain embedded hardware. It introduces a succinct set of language-level fixed-point types and impose constraints on the number of integer or fractional digits each can possess.
As with all examples of discrete-type fixed-point support, the limited choice of exponents is a considerable restriction on the versatility and expressiveness of the API.
Nevertheless, it may be possible to harness performance gains provided by N1169 fixed-point types through explicit template specialization. This is likely to be a valuable proposition to potential users of the library who find themselves targeting platforms which support fixed-point arithmetic at the hardware level.
There are many other C++ libraries available which fall into the latter category of continuous-range fixed-point arithmetic [mizvekov] [schregle] [viboes]. In particular, an existing library proposal by Lawrence Crowl, P0106 [P0106] (formerly N3352), aims to achieve very similar goals through similar means and warrants closer comparison than N1169.
P0106 introduces four class templates covering the quadrant of signed versus unsigned and fractional versus integer numeric types. It is intended to replace built-in types in a wide variety of situations and accordingly, is highly compile-time configurable in terms of how rounding and overflow are handled. Parameters to these four class templates include the range in bits and - for fractional types - the resolution.
The fixed_point class template could probably - with a few caveats - be generated using the two fractional types, nonnegative and negatable, replacing the Rep parameter with the integer bit count of Rep, specifying fastest for the rounding mode and specifying undefined as the overflow mode.
However, fixed_point more closely and concisely caters to the needs of users who already use integer types and simply desire a more concise, less error-prone form. It more closely follows the four design aims of the library and - it can be argued - more closely follows the spirit of the standard in its pursuit of zero-cost abstraction.
Some aspects of the design of the P0106 API which back up these conclusion are that:
A more detailed comparison of the approaches taken in this paper and P0106 can be found in [P0554].
SG6: Lawrence Crowl
SG14: Guy Davidson, Michael Wong
Contributors: Ed Ainsley, Billy Baker, Lance Dyson, Marco Foco, Mathias Gaunard, Clément Grégoire, Nicolas Guillemot, Kurt Guntheroth, Matt Kinzelman, Joël Lamotte, Sean Middleditch, Paul Robinson, Patrice Roy, Peter Schregle, Ryhor Spivak
This paper revises P0037R4:
multiply and divide functionsRadix non-type template parameter to fixed_point
fixed_point to put Exponent firstfractional and elastic_integerfractional c'tor and assignment operatorP0037R4 revises P0037R3:
width and set_widthRep template parameterRep Value and Class Template Deductionmake_fixed and make_ufixed function templatesconstant operatorsadd and subtract function templatesoperator/ and divideAn in-development implementation of the fixed_point class template and its essential supporting functions and types is available [github].
Despite a focus on usable interface and direct translation from integer-based fixed-point operations, there is an overwhelming expectation that the source code result in minimal instructions and clock cycles. A few preliminary numbers are presented to give a very early idea of how the API might perform.
Some notes:
A few test functions were run, ranging from single arithmetic operations to basic geometric functions, performed against integer, floating-point and fixed-point types for comparison.
Figures were taken from a single CPU, OS and compiler, namely:
Fixed inputs were provided to each function, meaning that branch prediction rarely fails. Results may also not represent the full range of inputs.
Details of the test harness used can be found in the source project mentioned in Appendix 1;
Times are in nanoseconds;
Code has not yet been optimized for performance.
Where applicable various combinations of integer, floating-point and fixed-point types were tested with the following identifiers:
uint8_t, int8_t, uint16_t, int16_t, uint32_t, int32_t, uint64_t and int64_t built-in integer types;float, double and long double built-in floating-point types;Plus, minus, multiplication and division were tested in isolation using a number of different numeric types with the following results:
name cpu_time
add(float) 1.78011
add(double) 1.73966
add(long double) 3.46011
add(u4_4) 1.87726
add(s3_4) 1.85051
add(u8_8) 1.85417
add(s7_8) 1.82057
add(u16_16) 1.94194
add(s15_16) 1.93463
add(u32_32) 1.94674
add(s31_32) 1.94446
add(int8_t) 2.14857
add(uint8_t) 2.12571
add(int16_t) 1.9936
add(uint16_t) 1.88229
add(int32_t) 1.82126
add(uint32_t) 1.76
add(int64_t) 1.76
add(uint64_t) 1.83223
sub(float) 1.96617
sub(double) 1.98491
sub(long double) 3.55474
sub(u4_4) 1.77006
sub(s3_4) 1.72983
sub(u8_8) 1.72983
sub(s7_8) 1.72983
sub(u16_16) 1.73966
sub(s15_16) 1.85051
sub(u32_32) 1.88229
sub(s31_32) 1.87063
sub(int8_t) 1.76
sub(uint8_t) 1.74994
sub(int16_t) 1.82126
sub(uint16_t) 1.83794
sub(int32_t) 1.89074
sub(uint32_t) 1.85417
sub(int64_t) 1.83703
sub(uint64_t) 2.04914
mul(float) 1.9376
mul(double) 1.93097
mul(long double) 102.446
mul(u4_4) 2.46583
mul(s3_4) 2.09189
mul(u8_8) 2.08
mul(s7_8) 2.18697
mul(u16_16) 2.12571
mul(s15_16) 2.10789
mul(u32_32) 2.10789
mul(s31_32) 2.10789
mul(int8_t) 1.76
mul(uint8_t) 1.78011
mul(int16_t) 1.8432
mul(uint16_t) 1.76914
mul(int32_t) 1.78011
mul(uint32_t) 2.19086
mul(int64_t) 1.7696
mul(uint64_t) 1.79017
div(float) 5.12
div(double) 7.64343
div(long double) 8.304
div(u4_4) 3.82171
div(s3_4) 3.82171
div(u8_8) 3.84
div(s7_8) 3.8
div(u16_16) 9.152
div(s15_16) 11.232
div(u32_32) 30.8434
div(s31_32) 34
div(int8_t) 3.82171
div(uint8_t) 3.82171
div(int16_t) 3.8
div(uint16_t) 3.82171
div(int32_t) 3.82171
div(uint32_t) 3.81806
div(int64_t) 10.2286
div(uint64_t) 8.304
Among the slowest types are long double. It is likely that they are emulated in software. The next slowest operations are fixed-point multiply and divide operations - especially with 64-bit types. This is because values need to be promoted temporarily to double-width types. This is a known fixed-point technique which inevitably experiences slowdown where a 128-bit type is required on a 64-bit system.
Here is a section of the disassembly of the s15:16 multiply call:
30: mov %r14,%rax
mov %r15,%rax
movslq -0x28(%rbp),%rax
movslq -0x30(%rbp),%rcx
imul %rax,%rcx
shr $0x10,%rcx
mov %ecx,-0x38(%rbp)
mov %r12,%rax
4c: movzbl (%rbx),%eax
cmp $0x1,%eax
↓ jne 68
54: mov 0x8(%rbx),%rax
lea 0x1(%rax),%rcx
mov %rcx,0x8(%rbx)
cmp 0x38(%rbx),%rax
↑ jb 30
The two 32-bit numbers are multiplied together and the result shifted down - much as it would if raw int values were used. The efficiency of this operation varies with the exponent. An exponent of zero should mean no shift at all.
A fast sqrt implementation has not yet been tested with fixed_point. (The naive implementation takes over 300ns.) For this reason, a magnitude-squared function is measured, combining multiply and add operations:
template <class FP>
constexpr FP magnitude_squared(const FP & x, const FP & y, const FP & z)
{
return x * x + y * y + z * z;
}
Only real number formats are tested:
float 2.42606
double 2.08
long double 4.5056
s3_4 2.768
s7_8 2.77577
s15_16 2.752
s31_32 4.10331
Again, the size of the type seems to have the largest impact.
A similar operation includes a comparison and branch:
template <class Real>
bool circle_intersect_generic(Real x1, Real y1, Real r1, Real x2, Real y2, Real r2)
{
auto x_diff = x2 - x1;
auto y_diff = y2 - y1;
auto distance_squared = x_diff * x_diff + y_diff * y_diff;
auto touch_distance = r1 + r2;
auto touch_distance_squared = touch_distance * touch_distance;
return distance_squared <= touch_distance_squared;
}
float 3.46011
double 3.48
long double 6.4
s3_4 3.88
s7_8 4.5312
s15_16 3.82171
s31_32 5.92
Again, fixed-point and native performance are comparable.