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