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