@@ -1,7 +1,7 @@ | |||||
// | // | ||||
// LolRemez - Remez algorithm implementation | // LolRemez - Remez algorithm implementation | ||||
// | // | ||||
// Copyright © 2005-2015 Sam Hocevar <sam@hocevar.net> | |||||
// Copyright © 2005—2015 Sam Hocevar <sam@hocevar.net> | |||||
// | // | ||||
// This program is free software. It comes without any warranty, to | // This program is free software. It comes without any warranty, to | ||||
// the extent permitted by applicable law. You can redistribute it | // the extent permitted by applicable law. You can redistribute it | ||||
@@ -317,21 +317,21 @@ private: | |||||
// r_mul <- "*" r_factor | // r_mul <- "*" r_factor | ||||
// r_div <- "/" r_factor | // r_div <- "/" r_factor | ||||
// term <- r_factor ( r_mul / r_div ) * | |||||
// r_term <- r_factor ( r_mul / r_div ) * | |||||
struct r_mul : ifapply<seq<_, one<'*'>, _, r_factor>, do_op<op::mul>> {}; | struct r_mul : ifapply<seq<_, one<'*'>, _, r_factor>, do_op<op::mul>> {}; | ||||
struct r_div : ifapply<seq<_, one<'/'>, _, r_factor>, do_op<op::div>> {}; | struct r_div : ifapply<seq<_, one<'/'>, _, r_factor>, do_op<op::div>> {}; | ||||
struct r_term : seq<r_factor, | struct r_term : seq<r_factor, | ||||
star<sor<r_mul, r_div>>> {}; | star<sor<r_mul, r_div>>> {}; | ||||
// add <- "+" r_term | |||||
// sub <- "-" r_term | |||||
// r_expr <- r_term ( add / sub ) * | |||||
// r_add <- "+" r_term | |||||
// r_sub <- "-" r_term | |||||
// r_expr <- r_term ( r_add / r_sub ) * | |||||
struct r_add : ifapply<seq<_, one<'+'>, _, r_term>, do_op<op::add>> {}; | struct r_add : ifapply<seq<_, one<'+'>, _, r_term>, do_op<op::add>> {}; | ||||
struct r_sub : ifapply<seq<_, one<'-'>, _, r_term>, do_op<op::sub>> {}; | struct r_sub : ifapply<seq<_, one<'-'>, _, r_term>, do_op<op::sub>> {}; | ||||
struct r_expr : seq<r_term, | struct r_expr : seq<r_term, | ||||
star<sor<r_add, r_sub>>> {}; | star<sor<r_add, r_sub>>> {}; | ||||
// stmt <- r_expr <end> | |||||
// r_stmt <- r_expr <end> | |||||
struct r_stmt : seq<_, r_expr, _, eof> {}; | struct r_stmt : seq<_, r_expr, _, eof> {}; | ||||
}; | }; | ||||
@@ -1,7 +1,7 @@ | |||||
// | // | ||||
// LolRemez - Remez algorithm implementation | // LolRemez - Remez algorithm implementation | ||||
// | // | ||||
// Copyright © 2005-2015 Sam Hocevar <sam@hocevar.net> | |||||
// Copyright © 2005—2015 Sam Hocevar <sam@hocevar.net> | |||||
// | // | ||||
// This program is free software. It comes without any warranty, to | // This program is free software. It comes without any warranty, to | ||||
// the extent permitted by applicable law. You can redistribute it | // the extent permitted by applicable law. You can redistribute it | ||||
@@ -1,7 +1,7 @@ | |||||
// | // | ||||
// LolRemez - Remez algorithm implementation | // LolRemez - Remez algorithm implementation | ||||
// | // | ||||
// Copyright © 2005-2015 Sam Hocevar <sam@hocevar.net> | |||||
// Copyright © 2005—2015 Sam Hocevar <sam@hocevar.net> | |||||
// | // | ||||
// This program is free software. It comes without any warranty, to | // This program is free software. It comes without any warranty, to | ||||
// the extent permitted by applicable law. You can redistribute it | // the extent permitted by applicable law. You can redistribute it | ||||
@@ -20,16 +20,16 @@ using namespace lol; | |||||
*/ | */ | ||||
template<typename T> | template<typename T> | ||||
struct LinearSystem : public array2d<T> | |||||
struct linear_system : public array2d<T> | |||||
{ | { | ||||
inline LinearSystem<T>(int cols) | |||||
inline linear_system<T>(int cols) | |||||
{ | { | ||||
ASSERT(cols > 0); | ASSERT(cols > 0); | ||||
this->SetSize(ivec2(cols)); | this->SetSize(ivec2(cols)); | ||||
} | } | ||||
void Init(T const &x) | |||||
void init(T const &x) | |||||
{ | { | ||||
int const n = this->GetSize().x; | int const n = this->GetSize().x; | ||||
@@ -39,12 +39,12 @@ struct LinearSystem : public array2d<T> | |||||
} | } | ||||
/* Naive matrix inversion */ | /* Naive matrix inversion */ | ||||
LinearSystem<T> inverse() const | |||||
linear_system<T> inverse() const | |||||
{ | { | ||||
int const n = this->GetSize().x; | int const n = this->GetSize().x; | ||||
LinearSystem a(*this), b(n); | |||||
linear_system a(*this), b(n); | |||||
b.Init((T)1); | |||||
b.init((T)1); | |||||
/* Inversion method: iterate through all columns and make sure | /* Inversion method: iterate through all columns and make sure | ||||
* all the terms are 1 on the diagonal and 0 everywhere else */ | * all the terms are 1 on the diagonal and 0 everywhere else */ | ||||
@@ -1,7 +1,7 @@ | |||||
// | // | ||||
// LolRemez - Remez algorithm implementation | // LolRemez - Remez algorithm implementation | ||||
// | // | ||||
// Copyright © 2005-2015 Sam Hocevar <sam@hocevar.net> | |||||
// Copyright © 2005—2015 Sam Hocevar <sam@hocevar.net> | |||||
// | // | ||||
// This program is free software. It comes without any warranty, to | // This program is free software. It comes without any warranty, to | ||||
// the extent permitted by applicable law. You can redistribute it | // the extent permitted by applicable law. You can redistribute it | ||||
@@ -91,9 +91,9 @@ void remez_solver::remez_init() | |||||
fxn.Push(eval_func(m_zeroes[i])); | fxn.Push(eval_func(m_zeroes[i])); | ||||
} | } | ||||
/* We build a matrix of Chebishev evaluations: row i contains the | |||||
/* We build a matrix of Chebyshev evaluations: row i contains the | |||||
* evaluations of x_i for polynomial order n = 0, 1, ... */ | * evaluations of x_i for polynomial order n = 0, 1, ... */ | ||||
LinearSystem<real> system(m_order + 1); | |||||
linear_system<real> system(m_order + 1); | |||||
for (int n = 0; n < m_order + 1; n++) | for (int n = 0; n < m_order + 1; n++) | ||||
{ | { | ||||
auto p = polynomial<real>::chebyshev(n); | auto p = polynomial<real>::chebyshev(n); | ||||
@@ -127,9 +127,9 @@ void remez_solver::remez_step() | |||||
for (int i = 0; i < m_order + 2; i++) | for (int i = 0; i < m_order + 2; i++) | ||||
fxn.Push(eval_func(m_control[i])); | fxn.Push(eval_func(m_control[i])); | ||||
/* We build a matrix of Chebishev evaluations: row i contains the | |||||
/* We build a matrix of Chebyshev evaluations: row i contains the | |||||
* evaluations of x_i for polynomial order n = 0, 1, ... */ | * evaluations of x_i for polynomial order n = 0, 1, ... */ | ||||
LinearSystem<real> system(m_order + 2); | |||||
linear_system<real> system(m_order + 2); | |||||
for (int n = 0; n < m_order + 1; n++) | for (int n = 0; n < m_order + 1; n++) | ||||
{ | { | ||||
auto p = polynomial<real>::chebyshev(n); | auto p = polynomial<real>::chebyshev(n); | ||||
@@ -173,7 +173,7 @@ void remez_solver::remez_step() | |||||
*/ | */ | ||||
void remez_solver::find_zeroes() | void remez_solver::find_zeroes() | ||||
{ | { | ||||
m_stats_cheby = m_stats_func = m_stats_weight = 0.f; | |||||
Timer t; | |||||
for (int i = 0; i < m_order + 1; i++) | for (int i = 0; i < m_order + 1; i++) | ||||
{ | { | ||||
@@ -210,8 +210,7 @@ void remez_solver::find_zeroes() | |||||
m_zeroes[i] = c.x; | m_zeroes[i] = c.x; | ||||
} | } | ||||
printf(" -:- timings for zeroes: estimate %f func %f weight %f\n", | |||||
m_stats_cheby, m_stats_func, m_stats_weight); | |||||
printf(" -:- timing for zeroes: %f ms\n", t.Get() * 1000.f); | |||||
} | } | ||||
/* | /* | ||||
@@ -225,14 +224,14 @@ void remez_solver::find_zeroes() | |||||
*/ | */ | ||||
void remez_solver::find_extrema() | void remez_solver::find_extrema() | ||||
{ | { | ||||
Timer t; | |||||
using std::printf; | using std::printf; | ||||
m_control[0] = -1; | m_control[0] = -1; | ||||
m_control[m_order + 1] = 1; | m_control[m_order + 1] = 1; | ||||
m_error = 0; | m_error = 0; | ||||
m_stats_cheby = m_stats_func = m_stats_weight = 0.f; | |||||
for (int i = 1; i < m_order + 1; i++) | for (int i = 1; i < m_order + 1; i++) | ||||
{ | { | ||||
struct { real x, err; } a, b, c, d; | struct { real x, err; } a, b, c, d; | ||||
@@ -277,8 +276,7 @@ void remez_solver::find_extrema() | |||||
m_control[i] = c.x; | m_control[i] = c.x; | ||||
} | } | ||||
printf(" -:- timings for extrema: estimate %f func %f weight %f\n", | |||||
m_stats_cheby, m_stats_func, m_stats_weight); | |||||
printf(" -:- timing for extrema: %f ms\n", t.Get() * 1000.f); | |||||
printf(" -:- error: "); | printf(" -:- error: "); | ||||
m_error.print(m_decimals); | m_error.print(m_decimals); | ||||
printf("\n"); | printf("\n"); | ||||
@@ -306,26 +304,17 @@ void remez_solver::print_poly() | |||||
real remez_solver::eval_estimate(real const &x) | real remez_solver::eval_estimate(real const &x) | ||||
{ | { | ||||
Timer t; | |||||
real ret = m_estimate.eval(x); | |||||
m_stats_cheby += t.Get(); | |||||
return ret; | |||||
return m_estimate.eval(x); | |||||
} | } | ||||
real remez_solver::eval_func(real const &x) | real remez_solver::eval_func(real const &x) | ||||
{ | { | ||||
Timer t; | |||||
real ret = m_func.eval(x * m_k2 + m_k1); | |||||
m_stats_func += t.Get(); | |||||
return ret; | |||||
return m_func.eval(x * m_k2 + m_k1); | |||||
} | } | ||||
real remez_solver::eval_weight(real const &x) | real remez_solver::eval_weight(real const &x) | ||||
{ | { | ||||
Timer t; | |||||
real ret = m_has_weight ? m_weight.eval(x * m_k2 + m_k1) : real(1); | |||||
m_stats_weight += t.Get(); | |||||
return ret; | |||||
return m_has_weight ? m_weight.eval(x * m_k2 + m_k1) : real(1); | |||||
} | } | ||||
real remez_solver::eval_error(real const &x) | real remez_solver::eval_error(real const &x) | ||||
@@ -1,7 +1,7 @@ | |||||
// | // | ||||
// LolRemez - Remez algorithm implementation | // LolRemez - Remez algorithm implementation | ||||
// | // | ||||
// Copyright © 2005-2015 Sam Hocevar <sam@hocevar.net> | |||||
// Copyright © 2005—2015 Sam Hocevar <sam@hocevar.net> | |||||
// | // | ||||
// This program is free software. It comes without any warranty, to | // This program is free software. It comes without any warranty, to | ||||
// the extent permitted by applicable law. You can redistribute it | // the extent permitted by applicable law. You can redistribute it | ||||
@@ -56,10 +56,5 @@ private: | |||||
lol::array<lol::real> m_control; | lol::array<lol::real> m_control; | ||||
lol::real m_k1, m_k2, m_epsilon, m_error; | lol::real m_k1, m_k2, m_epsilon, m_error; | ||||
/* Statistics */ | |||||
float m_stats_cheby; | |||||
float m_stats_func; | |||||
float m_stats_weight; | |||||
}; | }; | ||||