diff --git a/test/owens_t_T7.hpp b/test/owens_t_T7.hpp index 31a51abd2..fd688e639 100644 --- a/test/owens_t_T7.hpp +++ b/test/owens_t_T7.hpp @@ -95,22 +95,22 @@ namespace boost const RealType hs = h*h; const RealType as = a*a; - RealType u = one_div_two_pi() * a * exp(-half()*hs*(static_cast(1)+as)); + //RealType u = one_div_two_pi() * a * exp(-half()*hs*(static_cast(1)+as)); RealType v = c2[0]; - RealType val = u*v; + RealType val = v; RealType last_val = val+1; // last_val must not be the same as val - unsigned k = 0; + int k = 0; std::vector memory; - memory.push_back(u*v); + memory.push_back(v); - while((val != last_val) || (k(2*k+1)); @@ -133,12 +133,17 @@ namespace boost std::sort(ndx_4_sorted_data.begin(), ndx_4_sorted_data.end(), boost::bind(owens_t_sort_proxy::size_type, RealType>, _1, _2, &memory[0])); - val = static_cast(0); - for(unsigned i = 0; i != memory.size(); i++) + val = memory[ndx_4_sorted_data[0]]; + for(unsigned i = 1; i != memory.size(); i++) { val+=memory[ndx_4_sorted_data[i]]; } + // split the exponential to avoid values that go below the minimum floating pt value + val *= exp(-half()*hs*as); + val *= exp(-half()*hs); + val *= one_div_two_pi() * a; + return val; } // RealType compute_owens_t_T7(const RealType h, const RealType a) diff --git a/test/test_owens_t.cpp b/test/test_owens_t.cpp index 8709da55b..56137f482 100644 --- a/test/test_owens_t.cpp +++ b/test/test_owens_t.cpp @@ -159,8 +159,8 @@ void check_against_T7(RealType) using namespace std; // ADL of std names. // apply log scale because points near zero are more interesting - for(RealType a = static_cast(-10.0l); a < static_cast(3l); a+= static_cast(0.1l)) - for(RealType h = static_cast(-10.0l); h < static_cast(4l); h+= static_cast(0.1l)) + for(RealType a = static_cast(-10.0l); a < static_cast(3l); a+= static_cast(0.2l)) + for(RealType h = static_cast(-10.0l); h < static_cast(4l); h+= static_cast(0.2l)) { const RealType expa = exp(a); const RealType exph = exp(h);