libs/multiprecision/example/constexpr_float_arithmetic_examples.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 <iostream>
   7 #include <boost/math/constants/constants.hpp>
   8 #ifdef BOOST_HAS_FLOAT128
   9 #include <boost/multiprecision/float128.hpp>
  10 #endif
  11
  12 //[constexpr_circle
  13
  14 template <class T>
  15 inline constexpr T circumference(T radius)
  16 {
  17    return 2 * boost::math::constants::pi<T>() * radius;
  18 }
  19
  20 template <class T>
  21 inline constexpr T area(T radius)
  22 {
  23    return boost::math::constants::pi<T>() * radius * radius;
  24 }
  25 //]
  26
  27 template <class T, unsigned Order>
  28 struct const_polynomial
  29 {
  30  public:
  31    T data[Order + 1];
  32
  33  public:
  34    constexpr const_polynomial(T val = 0) : data{val} {}
  35    constexpr const_polynomial(const std::initializer_list<T>& init) : data{}
  36    {
  37       if (init.size() > Order + 1)
  38          throw std::range_error("Too many initializers in list");
  39       for (unsigned i = 0; i < init.size(); ++i)
  40          data[i] = init.begin()[i];
  41    }
  42    constexpr T& operator[](std::size_t N)
  43    {
  44       return data[N];
  45    }
  46    constexpr const T& operator[](std::size_t N) const
  47    {
  48       return data[N];
  49    }
  50    template <class U>
  51    constexpr T operator()(U val)const
  52    {
  53       T result = data[Order];
  54       for (unsigned i = Order; i > 0; --i)
  55       {
  56          result *= val;
  57          result += data[i - 1];
  58       }
  59       return result;
  60    }
  61    constexpr const_polynomial<T, Order - 1> derivative() const
  62    {
  63       const_polynomial<T, Order - 1> result;
  64       for (unsigned i = 1; i <= Order; ++i)
  65       {
  66          result[i - 1] = (*this)[i] * i;
  67       }
  68       return result;
  69    }
  70    constexpr const_polynomial operator-()
  71    {
  72       const_polynomial t(*this);
  73       for (unsigned i = 0; i <= Order; ++i)
  74          t[i] = -t[i];
  75       return t;
  76    }
  77    template <class U>
  78    constexpr const_polynomial& operator*=(U val)
  79    {
  80       for (unsigned i = 0; i <= Order; ++i)
  81          data[i] = data[i] * val;
  82       return *this;
  83    }
  84    template <class U>
  85    constexpr const_polynomial& operator/=(U val)
  86    {
  87       for (unsigned i = 0; i <= Order; ++i)
  88          data[i] = data[i] / val;
  89       return *this;
  90    }
  91    template <class U>
  92    constexpr const_polynomial& operator+=(U val)
  93    {
  94       data[0] += val;
  95       return *this;
  96    }
  97    template <class U>
  98    constexpr const_polynomial& operator-=(U val)
  99    {
 100       data[0] -= val;
 101       return *this;
 102    }
 103 };
 104
 105 template <class T, unsigned Order1, unsigned Order2>
 106 inline constexpr const_polynomial<T, (Order1 > Order2 ? Order1 : Order2)> operator+(const const_polynomial<T, Order1>& a, const const_polynomial<T, Order2>& b)
 107 {
 108    if
 109       constexpr(Order1 > Order2)
 110       {
 111          const_polynomial<T, Order1> result(a);
 112          for (unsigned i = 0; i <= Order2; ++i)
 113             result[i] += b[i];
 114          return result;
 115       }
 116    else
 117    {
 118       const_polynomial<T, Order2> result(b);
 119       for (unsigned i = 0; i <= Order1; ++i)
 120          result[i] += a[i];
 121       return result;
 122    }
 123 }
 124 template <class T, unsigned Order1, unsigned Order2>
 125 inline constexpr const_polynomial<T, (Order1 > Order2 ? Order1 : Order2)> operator-(const const_polynomial<T, Order1>& a, const const_polynomial<T, Order2>& b)
 126 {
 127    if
 128       constexpr(Order1 > Order2)
 129       {
 130          const_polynomial<T, Order1> result(a);
 131          for (unsigned i = 0; i <= Order2; ++i)
 132             result[i] -= b[i];
 133          return result;
 134       }
 135    else
 136    {
 137       const_polynomial<T, Order2> result(b);
 138       for (unsigned i = 0; i <= Order1; ++i)
 139          result[i] = a[i] - b[i];
 140       return result;
 141    }
 142 }
 143 template <class T, unsigned Order1, unsigned Order2>
 144 inline constexpr const_polynomial<T, Order1 + Order2> operator*(const const_polynomial<T, Order1>& a, const const_polynomial<T, Order2>& b)
 145 {
 146    const_polynomial<T, Order1 + Order2> result;
 147    for (unsigned i = 0; i <= Order1; ++i)
 148    {
 149       for (unsigned j = 0; j <= Order2; ++j)
 150       {
 151          result[i + j] += a[i] * b[j];
 152       }
 153    }
 154    return result;
 155 }
 156 template <class T, unsigned Order, class U>
 157 inline constexpr const_polynomial<T, Order> operator*(const const_polynomial<T, Order>& a, const U& b)
 158 {
 159    const_polynomial<T, Order> result(a);
 160    for (unsigned i = 0; i <= Order; ++i)
 161    {
 162       result[i] *= b;
 163    }
 164    return result;
 165 }
 166 template <class U, class T, unsigned Order>
 167 inline constexpr const_polynomial<T, Order> operator*(const U& b, const const_polynomial<T, Order>& a)
 168 {
 169    const_polynomial<T, Order> result(a);
 170    for (unsigned i = 0; i <= Order; ++i)
 171    {
 172       result[i] *= b;
 173    }
 174    return result;
 175 }
 176 template <class T, unsigned Order, class U>
 177 inline constexpr const_polynomial<T, Order> operator/(const const_polynomial<T, Order>& a, const U& b)
 178 {
 179    const_polynomial<T, Order> result;
 180    for (unsigned i = 0; i <= Order; ++i)
 181    {
 182       result[i] /= b;
 183    }
 184    return result;
 185 }
 186
 187 //[hermite_example
 188 template <class T, unsigned Order>
 189 class hermite_polynomial
 190 {
 191    const_polynomial<T, Order> m_data;
 192
 193  public:
 194    constexpr hermite_polynomial() : m_data(hermite_polynomial<T, Order - 1>().data() * const_polynomial<T, 1>{0, 2} - hermite_polynomial<T, Order - 1>().data().derivative())
 195    {
 196    }
 197    constexpr const const_polynomial<T, Order>& data() const
 198    {
 199       return m_data;
 200    }
 201    constexpr const T& operator[](std::size_t N)const
 202    {
 203       return m_data[N];
 204    }
 205    template <class U>
 206    constexpr T operator()(U val)const
 207    {
 208       return m_data(val);
 209    }
 210 };
 211 //]
 212 //[hermite_example2
 213 template <class T>
 214 class hermite_polynomial<T, 0>
 215 {
 216    const_polynomial<T, 0> m_data;
 217
 218  public:
 219    constexpr hermite_polynomial() : m_data{1} {}
 220    constexpr const const_polynomial<T, 0>& data() const
 221    {
 222       return m_data;
 223    }
 224    constexpr const T& operator[](std::size_t N) const
 225    {
 226       return m_data[N];
 227    }
 228    template <class U>
 229    constexpr T operator()(U val)
 230    {
 231       return m_data(val);
 232    }
 233 };
 234
 235 template <class T>
 236 class hermite_polynomial<T, 1>
 237 {
 238    const_polynomial<T, 1> m_data;
 239
 240  public:
 241    constexpr hermite_polynomial() : m_data{0, 2} {}
 242    constexpr const const_polynomial<T, 1>& data() const
 243    {
 244       return m_data;
 245    }
 246    constexpr const T& operator[](std::size_t N) const
 247    {
 248       return m_data[N];
 249    }
 250    template <class U>
 251    constexpr T operator()(U val)
 252    {
 253       return m_data(val);
 254    }
 255 };
 256 //]
 257
 258 void test_double()
 259 {
 260    constexpr double radius = 2.25;
 261    constexpr double c      = circumference(radius);
 262    constexpr double a      = area(radius);
 263
 264    std::cout << "Circumference = " << c << std::endl;
 265    std::cout << "Area = " << a << std::endl;
 266
 267    constexpr const_polynomial<double, 2> pa = {3, 4};
 268    constexpr const_polynomial<double, 2> pb = {5, 6};
 269    static_assert(pa[0] == 3);
 270    static_assert(pa[1] == 4);
 271    constexpr auto pc = pa * 2;
 272    static_assert(pc[0] == 6);
 273    static_assert(pc[1] == 8);
 274    constexpr auto pd = 3 * pa;
 275    static_assert(pd[0] == 3 * 3);
 276    static_assert(pd[1] == 4 * 3);
 277    constexpr auto pe = pa + pb;
 278    static_assert(pe[0] == 3 + 5);
 279    static_assert(pe[1] == 4 + 6);
 280    constexpr auto pf = pa - pb;
 281    static_assert(pf[0] == 3 - 5);
 282    static_assert(pf[1] == 4 - 6);
 283    constexpr auto pg = pa * pb;
 284    static_assert(pg[0] == 15);
 285    static_assert(pg[1] == 38);
 286    static_assert(pg[2] == 24);
 287
 288    constexpr hermite_polynomial<double, 2> h1;
 289    static_assert(h1[0] == -2);
 290    static_assert(h1[1] == 0);
 291    static_assert(h1[2] == 4);
 292
 293    constexpr hermite_polynomial<double, 3> h3;
 294    static_assert(h3[0] == 0);
 295    static_assert(h3[1] == -12);
 296    static_assert(h3[2] == 0);
 297    static_assert(h3[3] == 8);
 298
 299    constexpr hermite_polynomial<double, 9> h9;
 300    static_assert(h9[0] == 0);
 301    static_assert(h9[1] == 30240);
 302    static_assert(h9[2] == 0);
 303    static_assert(h9[3] == -80640);
 304    static_assert(h9[4] == 0);
 305    static_assert(h9[5] == 48384);
 306    static_assert(h9[6] == 0);
 307    static_assert(h9[7] == -9216);
 308    static_assert(h9[8] == 0);
 309    static_assert(h9[9] == 512);
 310
 311    static_assert(h9(0.5) == 6481);
 312
 313 }
 314
 315 void test_float128()
 316 {
 317 #ifdef BOOST_HAS_FLOAT128
 318    //[constexpr_circle_usage
 319
 320    using boost::multiprecision::float128;
 321
 322    constexpr float128 radius = 2.25;
 323    constexpr float128 c      = circumference(radius);
 324    constexpr float128 a      = area(radius);
 325
 326    std::cout << "Circumference = " << c << std::endl;
 327    std::cout << "Area = " << a << std::endl;
 328
 329    //]
 330
 331    constexpr hermite_polynomial<float128, 2> h1;
 332    static_assert(h1[0] == -2);
 333    static_assert(h1[1] == 0);
 334    static_assert(h1[2] == 4);
 335
 336    constexpr hermite_polynomial<float128, 3> h3;
 337    static_assert(h3[0] == 0);
 338    static_assert(h3[1] == -12);
 339    static_assert(h3[2] == 0);
 340    static_assert(h3[3] == 8);
 341
 342    //[hermite_example3
 343    constexpr hermite_polynomial<float128, 9> h9;
 344    //
 345    // Verify that the polynomial's coefficients match the known values:
 346    //
 347    static_assert(h9[0] == 0);
 348    static_assert(h9[1] == 30240);
 349    static_assert(h9[2] == 0);
 350    static_assert(h9[3] == -80640);
 351    static_assert(h9[4] == 0);
 352    static_assert(h9[5] == 48384);
 353    static_assert(h9[6] == 0);
 354    static_assert(h9[7] == -9216);
 355    static_assert(h9[8] == 0);
 356    static_assert(h9[9] == 512);
 357    //
 358    // Define an abscissa value to evaluate at:
 359    //
 360    constexpr float128 abscissa(0.5);
 361    //
 362    // Evaluate H_9(0.5) using all constexpr arithmetic:
 363    //
 364    static_assert(h9(abscissa) == 6481);
 365    //]
 366 #endif
 367 }
 368
 369 int main()
 370 {
 371    test_double();
 372    test_float128();
 373    std::cout << "Done!" << std::endl;
 374 }