1d8407755Smaya/**************************************************************************
2d8407755Smaya *
3d8407755Smaya * (C) Copyright VMware, Inc 2010.
4d8407755Smaya * (C) Copyright John Maddock 2006.
5d8407755Smaya * Use, modification and distribution are subject to the
6d8407755Smaya * Boost Software License, Version 1.0. (See accompanying file
7d8407755Smaya * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
8d8407755Smaya *
9d8407755Smaya **************************************************************************/
10d8407755Smaya
11d8407755Smaya
12d8407755Smaya/*
13d8407755Smaya * This file allows to compute the minimax polynomial coefficients we use
14d8407755Smaya * for fast exp2/log2.
15d8407755Smaya *
16d8407755Smaya * How to use this source:
17d8407755Smaya *
18d8407755Smaya * - Download and build the NTL library from
19d8407755Smaya *   http://shoup.net/ntl/download.html , or install libntl-dev package if on
20d8407755Smaya *   Debian.
21d8407755Smaya *
22d8407755Smaya * - Download boost source code matching to your distro.
23d8407755Smaya *
24d8407755Smaya * - Goto libs/math/minimax and replace f.cpp with this file.
25d8407755Smaya *
26d8407755Smaya * - Build as
27d8407755Smaya *
28d8407755Smaya *   g++ -o minimax -I /path/to/ntl/include main.cpp f.cpp /path/to/ntl/src/ntl.a
29d8407755Smaya *
30d8407755Smaya * - Run as
31d8407755Smaya *
32d8407755Smaya *    ./minimax
33d8407755Smaya *
34d8407755Smaya * - For example, to compute exp2 5th order polynomial between [0, 1] do:
35d8407755Smaya *
36d8407755Smaya *    variant 0
37d8407755Smaya *    range 0 1
38d8407755Smaya *    order 5 0
39d8407755Smaya *    step 200
40d8407755Smaya *    info
41d8407755Smaya *
42d8407755Smaya *  and take the coefficients from the P = { ... } array.
43d8407755Smaya *
44d8407755Smaya * - To compute log2 4th order polynomial between [0, 1/9] do:
45d8407755Smaya *
46d8407755Smaya *    variant 1
47d8407755Smaya *    range 0 0.111111112
48d8407755Smaya *    order 4 0
49d8407755Smaya *    step 200
50d8407755Smaya *    info
51d8407755Smaya *
52d8407755Smaya * - For more info see
53d8407755Smaya * http://www.boost.org/doc/libs/1_47_0/libs/math/doc/sf_and_dist/html/math_toolkit/toolkit/internals2/minimax.html
54d8407755Smaya */
55d8407755Smaya
56d8407755Smaya#define L22
57d8407755Smaya#include <boost/math/bindings/rr.hpp>
58d8407755Smaya#include <boost/math/tools/polynomial.hpp>
59d8407755Smaya
60d8407755Smaya#include <cmath>
61d8407755Smaya
62d8407755Smayaboost::math::ntl::RR exp2(const boost::math::ntl::RR& x)
63d8407755Smaya{
64d8407755Smaya      return exp(x*log(2.0));
65d8407755Smaya}
66d8407755Smaya
67d8407755Smayaboost::math::ntl::RR log2(const boost::math::ntl::RR& x)
68d8407755Smaya{
69d8407755Smaya      return log(x)/log(2.0);
70d8407755Smaya}
71d8407755Smaya
72d8407755Smayaboost::math::ntl::RR f(const boost::math::ntl::RR& x, int variant)
73d8407755Smaya{
74d8407755Smaya   switch(variant)
75d8407755Smaya   {
76d8407755Smaya   case 0:
77d8407755Smaya      return exp2(x);
78d8407755Smaya
79d8407755Smaya   case 1:
80d8407755Smaya      return log2((1.0 + sqrt(x))/(1.0 - sqrt(x)))/sqrt(x);
81d8407755Smaya   }
82d8407755Smaya
83d8407755Smaya   return 0;
84d8407755Smaya}
85d8407755Smaya
86d8407755Smaya
87d8407755Smayavoid show_extra(
88d8407755Smaya   const boost::math::tools::polynomial<boost::math::ntl::RR>& n,
89d8407755Smaya   const boost::math::tools::polynomial<boost::math::ntl::RR>& d,
90d8407755Smaya   const boost::math::ntl::RR& x_offset,
91d8407755Smaya   const boost::math::ntl::RR& y_offset,
92d8407755Smaya   int variant)
93d8407755Smaya{
94d8407755Smaya   switch(variant)
95d8407755Smaya   {
96d8407755Smaya   default:
97d8407755Smaya      // do nothing here...
98d8407755Smaya      ;
99d8407755Smaya   }
100d8407755Smaya}
101d8407755Smaya
102