CNL  2.0.2 (development)
Compositional Numeric Library
math.h
Go to the documentation of this file.
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 #if !defined(CNL_IMPL_SCALED_INTEGER_MATH_H)
11 #define CNL_IMPL_SCALED_INTEGER_MATH_H
12 
13 #include "definition.h"
14 #include "rep_of.h"
15 #include "tag_of.h"
16 
18 namespace cnl {
19 
21  // implementation-specific definitions
22 
23  namespace _impl {
24  namespace fp {
25 
26  template<class ScaledInteger>
27  [[nodiscard]] constexpr auto rounding_conversion(double d)
28  {
29  using one_longer = scaled_integer<
30  set_digits_t<rep_of_t<ScaledInteger>, digits_v<ScaledInteger> + 1>,
31  power<tag_of_t<ScaledInteger>::exponent - 1>>;
32  return from_rep<ScaledInteger>(
33  static_cast<rep_of_t<ScaledInteger>>((_impl::to_rep(one_longer{d}) + 1) >> 1));
34  }
35 
36  template<class ScaledInteger>
37  using unsigned_rep = typename std::make_unsigned<rep_of_t<ScaledInteger>>::type;
38 
39  template<class Input>
40  using make_largest_ufraction =
41  scaled_integer<unsigned_rep<Input>, power<-digits_v<unsigned_rep<Input>>>>;
42 
43  static_assert(
44  std::is_same<
45  make_largest_ufraction<scaled_integer<int32_t, power<-15>>>,
46  scaled_integer<uint32_t, power<-32>>>::value);
47 
48  // TODO: template magic to get the coefficients automatically
49  // from the number of bits of precision
50  // Define the coefficients as [[nodiscard]] constexpr,
51  // to make sure they're converted to fp
52  // at compile time
53  template<class CoeffType>
54  struct poly_coeffs {
55  static constexpr CoeffType a1{rounding_conversion<CoeffType>(0.6931471860838825)};
56  static constexpr CoeffType a2{rounding_conversion<CoeffType>(0.2402263846181129)};
57  static constexpr CoeffType a3{rounding_conversion<CoeffType>(0.055505126858894846)};
58  static constexpr CoeffType a4{rounding_conversion<CoeffType>(0.009614017013719252)};
59  static constexpr CoeffType a5{
60  rounding_conversion<CoeffType>(0.0013422634797558564)};
61  static constexpr CoeffType a6{
62  rounding_conversion<CoeffType>(0.00014352314226313836)};
63  static constexpr CoeffType a7{
64  rounding_conversion<CoeffType>(0.000021498763160402416)};
65  };
66 
67  template<typename A, typename B>
68  [[nodiscard]] constexpr auto safe_multiply(A const& a, B const& b) requires(digits_v<decltype(a * b)> <= digits_v<A> + digits_v<B>)
69  {
70  return set_digits_t<A, digits_v<A> + digits_v<B>>{a}
71  * set_digits_t<B, digits_v<A> + digits_v<B>>{b};
72  }
73 
74  template<typename A, typename B>
75  [[nodiscard]] constexpr auto safe_multiply(A const& a, B const& b) requires(digits_v<A> + digits_v<B> <= digits_v<decltype(a * b)>)
76  {
77  return a * b;
78  }
79 
80  template<class Rep, int Exponent>
81  [[nodiscard]] inline constexpr auto evaluate_polynomial(
82  scaled_integer<Rep, power<Exponent>> xf)
83  {
84  using fp = scaled_integer<Rep, power<Exponent>>;
85 
86  // Use a polynomial min-max approximation to generate the exponential of
87  // the fraction part. Note that the constant 1 of the polynomial is added later,
88  // this gives us one more bit of precision here for free
89  using coeffs = poly_coeffs<fp>;
90  return fp{safe_multiply(
91  xf,
92  (coeffs::a1
93  + fp{safe_multiply(
94  xf,
95  (coeffs::a2
96  + fp{safe_multiply(
97  xf,
98  (coeffs::a3
99  + fp{safe_multiply(
100  xf,
101  (coeffs::a4
102  + fp{safe_multiply(
103  xf,
104  (coeffs::a5
105  + fp{safe_multiply(
106  xf,
107  (coeffs::a6
108  + fp{safe_multiply(
109  coeffs::a7,
110  xf)}))}))}))}))}))}))};
111  }
112 
113  // Computes 2^x - 1 for a number x between 0 and 1, strictly less than 1
114  // If the exponent is not negative, there is no fraction part,
115  // so this is always zero
116  template<class Rep, int Exponent>
117  [[nodiscard]] inline constexpr auto exp2m1_0to1(scaled_integer<Rep, power<Exponent>>) requires(Exponent >= 0)
118  {
119  // Cannot construct from 0, since that would be a shift by more than width of type!
120  return from_rep<make_largest_ufraction<scaled_integer<Rep, power<Exponent>>>>(0);
121  }
122 
123  // for a positive exponent, some work needs to be done
124  template<class Rep, int Exponent>
125  [[nodiscard]] inline constexpr auto exp2m1_0to1(scaled_integer<Rep, power<Exponent>> x) requires(Exponent < 0)
126  {
127  // Build the type with the same number of bits, all fraction,
128  // and unsigned. That should be enough to exactly hold enough bits
129  // to guarantee bit-accurate results
130  using im = make_largest_ufraction<scaled_integer<Rep, power<Exponent>>>;
131  // The intermediate value type
132 
133  return evaluate_polynomial(im{scaled_integer<rep_of_t<im>, power<Exponent>>{x}});
134  }
135 
136  template<class Rep, int Exponent, int Radix>
137  requires(-digits_v<Rep> < Exponent)
138  [[nodiscard]] constexpr auto fractional(
139  scaled_integer<Rep, power<Exponent, Radix>> const& x,
140  Rep const& floored)
141  {
142  return x - floored;
143  }
144 
145  template<class Rep, int Exponent, int Radix>
146  requires(-digits_v<Rep> >= Exponent)
147  [[nodiscard]] constexpr auto fractional(scaled_integer<Rep, power<Exponent, Radix>> const& x, Rep const&)
148  {
149  return x;
150  }
151 
152  template<class Intermediate, typename Rep, int Exponent>
153  [[nodiscard]] constexpr auto exp2(
154  scaled_integer<Rep, power<Exponent>> const& x, Rep const& floored)
155  {
156  return floored <= Exponent
157  ? rep_of_t<Intermediate>{1} // return immediately if the shift would
158  // result in all bits being shifted out
159  // Do the shifts manually. Once the branch with shift operators is
160  // merged, could use those
161  : (_impl::to_rep(exp2m1_0to1<Rep, Exponent>(fractional(
162  x,
163  floored))) // Calculate the exponent of the fraction part
164  >> (-tag_of_t<Intermediate>::exponent + Exponent
165  - floored)) // shift it to the right place
166  + (Rep{1}
167  << (floored
168  - Exponent)); // The constant term must be one, to
169  // make integer powers correct
170  }
171  }
172  }
173 
182  template<class Rep, int Exponent>
183  [[nodiscard]] constexpr auto
184  exp2(scaled_integer<Rep, power<Exponent>> x) -> scaled_integer<Rep, power<Exponent>>
185  {
186  using out_type = scaled_integer<Rep, power<Exponent>>;
187  // The input type
188  using im = _impl::fp::make_largest_ufraction<out_type>;
189 
190  // Calculate the final result by shifting the fraction part around.
191  // Remember to add the 1 which is left out to get 1 bit more resolution
192  return _impl::from_rep<out_type>(_impl::fp::exp2<im>(x, static_cast<Rep>(floor(x))));
193  }
194 
195 }
196 
197 #endif /* CNL_IMPL_SCALED_INTEGER_MATH_H */
std::is_same
cnl::digits_v
constexpr int digits_v
provide number of binary digits of the given type
Definition: digits.h:31
cnl::exp2
constexpr auto exp2(scaled_integer< Rep, power< Exponent >> x) -> scaled_integer< Rep, power< Exponent >>
Calculates exp2(x), i.e. 2^xAccurate to 1LSB for up to 32 bit underlying representation.
Definition: math.h:184
cnl::power::exponent
constexpr static int exponent
value of template parameter, Exponent
Definition: definition.h:32
cnl
compositional numeric library
Definition: abort.h:15
cnl::scaled_integer
_impl::wrapper< Rep, Scale > scaled_integer
literal real number approximation that uses fixed-point arithmetic
Definition: definition.h:52
cnl::power
tag representing the scaling of an integer by a fixed factor
Definition: declaration.h:13
std::make_unsigned