libs/multiprecision/test/constexpr_test_cpp_int_6.cpp

   1 //  (C) Copyright John Maddock 2019.
   2 //  Use, modification and distribution are subject to the
   3 //  Boost Software License, Version 1.0. (See accompanying file
   4 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   5
   6 #include "boost/multiprecision/cpp_int.hpp"
   7 #include "test.hpp"
   8
   9 template <class T, unsigned Order>
  10 struct const_polynomial
  11 {
  12  public:
  13    T data[Order + 1];
  14
  15  public:
  16    constexpr const_polynomial(T val = 0) : data{val} {}
  17    constexpr const_polynomial(const const_polynomial&) = default;
  18    constexpr const_polynomial(const std::initializer_list<T>& init) : data{}
  19    {
  20       if (init.size() > Order + 1)
  21          throw std::range_error("Too many initializers in list");
  22       for (unsigned i = 0; i < init.size(); ++i)
  23          data[i] = init.begin()[i];
  24    }
  25    constexpr T& operator[](std::size_t N)
  26    {
  27       return data[N];
  28    }
  29    constexpr const T& operator[](std::size_t N) const
  30    {
  31       return data[N];
  32    }
  33    template <class U>
  34    constexpr T operator()(U val) const
  35    {
  36       T result = data[Order];
  37       for (unsigned i = Order; i > 0; --i)
  38       {
  39          result *= val;
  40          result += data[i - 1];
  41       }
  42       return result;
  43    }
  44    constexpr const_polynomial<T, Order - 1> derivative() const
  45    {
  46       const_polynomial<T, Order - 1> result;
  47       for (unsigned i = 1; i <= Order; ++i)
  48       {
  49          result[i - 1] = (*this)[i] * i;
  50       }
  51       return result;
  52    }
  53    constexpr const_polynomial operator-()
  54    {
  55       const_polynomial t(*this);
  56       for (unsigned i = 0; i <= Order; ++i)
  57          t[i] = -t[i];
  58       return t;
  59    }
  60    template <class U>
  61    constexpr const_polynomial& operator*=(U val)
  62    {
  63       for (unsigned i = 0; i <= Order; ++i)
  64          data[i] = data[i] * val;
  65       return *this;
  66    }
  67    template <class U>
  68    constexpr const_polynomial& operator/=(U val)
  69    {
  70       for (unsigned i = 0; i <= Order; ++i)
  71          data[i] = data[i] / val;
  72       return *this;
  73    }
  74    template <class U>
  75    constexpr const_polynomial& operator+=(U val)
  76    {
  77       data[0] += val;
  78       return *this;
  79    }
  80    template <class U>
  81    constexpr const_polynomial& operator-=(U val)
  82    {
  83       data[0] -= val;
  84       return *this;
  85    }
  86 };
  87
  88 template <class T, unsigned Order1, unsigned Order2>
  89 inline constexpr const_polynomial<T, (Order1 > Order2 ? Order1 : Order2)> operator+(const const_polynomial<T, Order1>& a, const const_polynomial<T, Order2>& b)
  90 {
  91    if
  92       constexpr(Order1 > Order2)
  93       {
  94          const_polynomial<T, Order1> result(a);
  95          for (unsigned i = 0; i <= Order2; ++i)
  96             result[i] += b[i];
  97          return result;
  98       }
  99    else
 100    {
 101       const_polynomial<T, Order2> result(b);
 102       for (unsigned i = 0; i <= Order1; ++i)
 103          result[i] += a[i];
 104       return result;
 105    }
 106 }
 107 template <class T, unsigned Order1, unsigned Order2>
 108 inline constexpr const_polynomial<T, (Order1 > Order2 ? Order1 : Order2)> operator-(const const_polynomial<T, Order1>& a, const const_polynomial<T, Order2>& b)
 109 {
 110    if
 111       constexpr(Order1 > Order2)
 112       {
 113          const_polynomial<T, Order1> result(a);
 114          for (unsigned i = 0; i <= Order2; ++i)
 115             result[i] -= b[i];
 116          return result;
 117       }
 118    else
 119    {
 120       const_polynomial<T, Order2> result(b);
 121       for (unsigned i = 0; i <= Order1; ++i)
 122          result[i] = a[i] - b[i];
 123       return result;
 124    }
 125 }
 126 template <class T, unsigned Order1, unsigned Order2>
 127 inline constexpr const_polynomial<T, Order1 + Order2> operator*(const const_polynomial<T, Order1>& a, const const_polynomial<T, Order2>& b)
 128 {
 129    const_polynomial<T, Order1 + Order2> result;
 130    for (unsigned i = 0; i <= Order1; ++i)
 131    {
 132       for (unsigned j = 0; j <= Order2; ++j)
 133       {
 134          result[i + j] += a[i] * b[j];
 135       }
 136    }
 137    return result;
 138 }
 139 template <class T, unsigned Order, class U>
 140 inline constexpr const_polynomial<T, Order> operator*(const const_polynomial<T, Order>& a, const U& b)
 141 {
 142    const_polynomial<T, Order> result(a);
 143    for (unsigned i = 0; i <= Order; ++i)
 144    {
 145       result[i] *= b;
 146    }
 147    return result;
 148 }
 149 template <class U, class T, unsigned Order>
 150 inline constexpr const_polynomial<T, Order> operator*(const U& b, const const_polynomial<T, Order>& a)
 151 {
 152    const_polynomial<T, Order> result(a);
 153    for (unsigned i = 0; i <= Order; ++i)
 154    {
 155       result[i] *= b;
 156    }
 157    return result;
 158 }
 159 template <class T, unsigned Order, class U>
 160 inline constexpr const_polynomial<T, Order> operator/(const const_polynomial<T, Order>& a, const U& b)
 161 {
 162    const_polynomial<T, Order> result;
 163    for (unsigned i = 0; i <= Order; ++i)
 164    {
 165       result[i] /= b;
 166    }
 167    return result;
 168 }
 169
 170 template <class T, unsigned Order>
 171 class hermite_polynomial
 172 {
 173    const_polynomial<T, Order> m_data;
 174
 175  public:
 176    constexpr hermite_polynomial() : m_data(hermite_polynomial<T, Order - 1>().data() * const_polynomial<T, 1>{0, 2} - hermite_polynomial<T, Order - 1>().data().derivative())
 177    {
 178    }
 179    constexpr const const_polynomial<T, Order>& data() const
 180    {
 181       return m_data;
 182    }
 183    constexpr const T& operator[](std::size_t N) const
 184    {
 185       return m_data[N];
 186    }
 187    template <class U>
 188    constexpr T operator()(U val) const
 189    {
 190       return m_data(val);
 191    }
 192 };
 193
 194 template <class T>
 195 class hermite_polynomial<T, 0>
 196 {
 197    const_polynomial<T, 0> m_data;
 198
 199  public:
 200    constexpr       hermite_polynomial() : m_data{1} {}
 201    constexpr const const_polynomial<T, 0>& data() const
 202    {
 203       return m_data;
 204    }
 205    constexpr const T& operator[](std::size_t N) const
 206    {
 207       return m_data[N];
 208    }
 209    template <class U>
 210    constexpr T operator()(U val)
 211    {
 212       return m_data(val);
 213    }
 214 };
 215
 216 template <class T>
 217 class hermite_polynomial<T, 1>
 218 {
 219    const_polynomial<T, 1> m_data;
 220
 221  public:
 222    constexpr       hermite_polynomial() : m_data{0, 2} {}
 223    constexpr const const_polynomial<T, 1>& data() const
 224    {
 225       return m_data;
 226    }
 227    constexpr const T& operator[](std::size_t N) const
 228    {
 229       return m_data[N];
 230    }
 231    template <class U>
 232    constexpr T operator()(U val)
 233    {
 234       return m_data(val);
 235    }
 236 };
 237
 238
 239
 240 int main()
 241 {
 242    using namespace boost::multiprecision::literals;
 243
 244    typedef boost::multiprecision::checked_int1024_t  int_backend;
 245
 246    // 8192 x^13 - 319488 x^11 + 4392960 x^9 - 26357760 x^7 + 69189120 x^5 - 69189120 x^3 + 17297280 x
 247    constexpr hermite_polynomial<int_backend, 13> h;
 248
 249    static_assert(h[0] == 0);
 250    static_assert(h[1] == 17297280);
 251    static_assert(h[2] == 0);
 252    static_assert(h[3] == -69189120);
 253    static_assert(h[4] == 0);
 254    static_assert(h[5] == 69189120);
 255    static_assert(h[6] == 0);
 256    static_assert(h[7] == -26357760);
 257    static_assert(h[8] == 0);
 258    static_assert(h[9] == 4392960);
 259    static_assert(h[10] == 0);
 260    static_assert(h[11] == -319488);
 261    static_assert(h[12] == 0);
 262    static_assert(h[13] == 8192);
 263
 264    return boost::report_errors();
 265 }
 266