fixed_point  rev.2
Binary Fixed-Point Arithmetic Library in C++
fixed_point_math.h
1 
2 // Copyright Timo Alho 2016.
3 // Distributed under the Boost Software License, Version 1.0.
4 // (See accompanying file ../LICENSE_1_0.txt or copy at
5 // http://www.boost.org/LICENSE_1_0.txt)
6 
9 
10 #ifndef FIXED_POINT_MATH_H_
11 #define FIXED_POINT_MATH_H_
12 
13 #include <sg14/fixed_point>
14 
16 namespace sg14 {
17 
19  // implementation-specific definitions
20 
21  namespace _impl {
22  namespace fp {
23 
24  template<class FixedPoint>
25  constexpr FixedPoint rounding_conversion(double d) {
26  using one_longer = make_fixed<FixedPoint::integer_digits, FixedPoint::fractional_digits + 1>;
27  return FixedPoint::from_data(static_cast<typename FixedPoint::rep>((one_longer{ d }.data() + 1) >> 1));
28  }
29 
30  template<class FixedPoint>
31  using unsigned_rep = typename std::make_unsigned<typename FixedPoint::rep>::type;
32 
33  template<class Input>
34  using make_largest_ufraction = fixed_point<unsigned_rep<Input>, -std::numeric_limits<unsigned_rep<Input>>::digits>;
35 
36  static_assert(std::is_same<make_largest_ufraction<fixed_point<int32_t, -15>>, fixed_point<uint32_t, -32>>::value, "");
37 
38  //TODO: template magic to get the coefficients automatically
39  //from the number of bits of precision
40  //Define the coefficients as constexpr,
41  //to make sure they're converted to fp
42  //at compile time
43  template<class CoeffType>
44  struct poly_coeffs {
45  static constexpr CoeffType a1 { rounding_conversion<CoeffType>(0.6931471860838825) };
46  static constexpr CoeffType a2 { rounding_conversion<CoeffType>(0.2402263846181129) };
47  static constexpr CoeffType a3 { rounding_conversion<CoeffType>(
48  0.055505126858894846) };
49  static constexpr CoeffType a4 { rounding_conversion<CoeffType>(
50  0.009614017013719252) };
51  static constexpr CoeffType a5 { rounding_conversion<CoeffType>(
52  0.0013422634797558564) };
53  static constexpr CoeffType a6 { rounding_conversion<CoeffType>(
54  0.00014352314226313836) };
55  static constexpr CoeffType a7 { rounding_conversion<CoeffType>(
56  0.000021498763160402416) };
57  };
58 
59  template<class CoeffType>
60  constexpr CoeffType poly_coeffs<CoeffType>::a1;
61  template<class CoeffType>
62  constexpr CoeffType poly_coeffs<CoeffType>::a2;
63  template<class CoeffType>
64  constexpr CoeffType poly_coeffs<CoeffType>::a3;
65  template<class CoeffType>
66  constexpr CoeffType poly_coeffs<CoeffType>::a4;
67  template<class CoeffType>
68  constexpr CoeffType poly_coeffs<CoeffType>::a5;
69  template<class CoeffType>
70  constexpr CoeffType poly_coeffs<CoeffType>::a6;
71  template<class CoeffType>
72  constexpr CoeffType poly_coeffs<CoeffType>::a7;
73 
74  template<class Rep, int Exponent>
75  constexpr inline fixed_point<Rep, Exponent> evaluate_polynomial(
76  fixed_point<Rep, Exponent> xf) {
77  using fp = fixed_point<Rep, Exponent>;
78 
79  //Use a polynomial min-max approximation to generate the exponential of
80  //the fractional part. Note that the constant 1 of the polynomial is added later,
81  //this gives us one more bit of precision here for free
82  using coeffs = poly_coeffs<fp>;
83  return fp{multiply(xf, (coeffs::a1+fp{multiply(xf, (coeffs::a2+fp{multiply(xf, (coeffs::a3+fp{multiply(xf, (coeffs::a4
84  +fp{multiply(xf, (coeffs::a5+fp{multiply(xf, (coeffs::a6+fp{multiply(fp{coeffs::a7}, fp{xf})}))}))}))}))}))}))};
85  }
86 
87  //Computes 2^x - 1 for a number x between 0 and 1, strictly less than 1
88  //If the exponent is not negative, there is no fractional part,
89  //so this is always zero
90  template<class Rep, int Exponent, _impl::enable_if_t<(Exponent>=0), int> dummy = 0>
91  inline constexpr make_largest_ufraction<fixed_point<Rep, Exponent>> exp2m1_0to1(
92  fixed_point<Rep, Exponent>) {
93  return make_largest_ufraction<fixed_point<Rep, Exponent>>::from_data(
94  0); //Cannot construct from 0, since that would be a shift by more than width of type!
95  }
96  //for a positive exponent, some work needs to be done
97  template<class Rep, int Exponent, _impl::enable_if_t<(Exponent<0), int> dummy = 0>
98  constexpr inline make_largest_ufraction<fixed_point<Rep, Exponent>> exp2m1_0to1(
99  fixed_point<Rep, Exponent> x) {
100 
101  //Build the type with the same number of bits, all fractional,
102  //and unsigned. That should be enough to exactly hold enough bits
103  //to guarantee bit-accurate results
104  using im = make_largest_ufraction<fixed_point<Rep, Exponent>>;
105  //The intermediate value type
106 
107  return evaluate_polynomial(im{x}); //Important: convert the type once, to keep every multiply from costing a cast
108  }
109 
110  template<class Rep, int Exponent>
111  constexpr inline Rep floor(fixed_point<Rep, Exponent> x) {
112  return Rep { (x.data()) >> -Exponent };
113  }
114 
115  }
116  }
117 
126  template<class Rep, int Exponent>
128  using namespace _impl::fp;
129 
130  using out_type = fixed_point<Rep, Exponent>;
131  // The input type
132  using im = make_largest_ufraction<out_type>;
133 
134  //Calculate the final result by shifting the fractional part around.
135  //Remember to add the 1 which is left out to get 1 bit more resolution
136  return out_type::from_data(
137  floor(x) <= Exponent ?
138  typename im::rep{1}//return immediately if the shift would result in all bits being shifted out
139  :
140  //Do the shifts manually. Once the branch with shift operators is merged, could use those
141  (exp2m1_0to1<Rep, Exponent>(static_cast<out_type>(x - floor(x))).data()//Calculate the exponent of the fractional part
142  >> (-im::exponent + Exponent - floor(x)))//shift it to the right place
143  + (Rep { 1 } << (floor(x) - Exponent))); //The constant term must be one, to make integer powers correct
144  }
145 
146 }
147 
148 #endif /* FIXED_POINT_MATH_H_ */
literal real number approximation that uses fixed-point arithmetic
Definition: fixed_point_type.h:20
constexpr fixed_point< Rep, Exponent > exp2(fixed_point< Rep, Exponent > x)
Calculates exp2(x), i.e. 2^xAccurate to 1LSB for up to 32 bit underlying representation.
Definition: fixed_point_math.h:127
study group 14 of the C++ working group
Definition: const_integer.h:22