Skip to content

Commit

Permalink
implement limits of F8 and F9 for xi -> 0
Browse files Browse the repository at this point in the history
  • Loading branch information
Alexander Voigt authored and Alexander Voigt committed Feb 2, 2017
1 parent 60ea367 commit dce90c2
Show file tree
Hide file tree
Showing 2 changed files with 51 additions and 0 deletions.
27 changes: 27 additions & 0 deletions src/threshold_loop_functions.cpp
Expand Up @@ -193,6 +193,15 @@ static double F8_1_x2(double x1, double x2)
/pow6(-1. + sqr(x2));
}

/// F8(x1,x2) in the limit x1 -> 0
static double F8_0_x2(double x1, double x2)
{
const double lx22 = std::log(sqr(x2));

return -2. + (2.*sqr(x1)*lx22)/(-1. + sqr(x2)) +
(2.*sqr(x2)*lx22)/(-1. + sqr(x2));
}

// F8(x1,x2) in the limit x1 -> x2
static double F8_x1_x2(double x1, double x2)
{
Expand Down Expand Up @@ -245,6 +254,12 @@ double F8(double x1, double x2)
return F8_1_x2(x2, x1);
}

if (is_equal(x1, 0., 0.0001))
return F8_0_x2(x1, x2);

if (is_equal(x2, 0., 0.0001))
return F8_0_x2(x2, x1);

if (is_equal(x1, x2, 0.00001))
return F8_x1_x2(x1, x2);

Expand Down Expand Up @@ -303,6 +318,12 @@ static double F9_1_x2(double x1, double x2)
/pow6(-1. + sqr(x2));
}

/// F9(x1,x2) in the limit x1 -> 0
static double F9_0_x2(double x1, double x2)
{
return (2.*std::log(sqr(x2)))/(-1. + sqr(x2));
}

/// F9(x1,x2) in the limit x1 -> x2
static double F9_x1_x2(double x1, double x2)
{
Expand Down Expand Up @@ -345,6 +366,12 @@ double F9(double x1, double x2)
return F9_1_x2(x2, x1);
}

if (is_equal(x1, 0., 0.0001))
return F9_0_x2(x1, x2);

if (is_equal(x2, 0., 0.0001))
return F9_0_x2(x2, x1);

if (is_equal(x1, x2, 0.00001))
return F9_x1_x2(x1, x2);

Expand Down
24 changes: 24 additions & 0 deletions test/test_threshold_loop_functions.cpp
Expand Up @@ -428,6 +428,20 @@ BOOST_AUTO_TEST_CASE(test_F8)
x = 0.001; BOOST_CHECK_SMALL(F8(1,x) - F8_bare(1.00001,x), 1e-2);
x = 0.0001; BOOST_CHECK_SMALL(F8(1,x) - F8_bare(1.00001,x), 1e-2);

x = 0.1; BOOST_CHECK_CLOSE(F8(x,0), F8_bare(x,0.00001), 1e-2);
x = 0.02; BOOST_CHECK_CLOSE(F8(x,0), F8_bare(x,0.00001), 1e-2);
x = 0.01; BOOST_CHECK_CLOSE(F8(x,0), F8_bare(x,0.00001), 1e-2);
x = 0.001; BOOST_CHECK_CLOSE(F8(x,0), F8_bare(x,0.00001), 1e-2);
x = 0.0001; BOOST_CHECK_CLOSE(F8(x,0), F8_bare(x,0.00001), 1e-2);

x = 0.1; BOOST_CHECK_CLOSE(F8(0,x), F8_bare(0.00001,x), 1e-2);
x = 0.02; BOOST_CHECK_CLOSE(F8(0,x), F8_bare(0.00001,x), 1e-2);
x = 0.01; BOOST_CHECK_CLOSE(F8(0,x), F8_bare(0.00001,x), 1e-2);
x = 0.001; BOOST_CHECK_CLOSE(F8(0,x), F8_bare(0.00001,x), 1e-2);
x = 0.0001; BOOST_CHECK_CLOSE(F8(0,x), F8_bare(0.00001,x), 1e-2);

x = 0.0001; BOOST_CHECK_CLOSE(F8(0,0), F8_bare(x,x+0.00001), 1e-2);

BOOST_CHECK(!std::isnan(F8(0,0)));
BOOST_CHECK(!std::isnan(F8(0,1)));
BOOST_CHECK(!std::isnan(F8(1,0)));
Expand Down Expand Up @@ -477,6 +491,16 @@ BOOST_AUTO_TEST_CASE(test_F9)
x = 0.001; BOOST_CHECK_CLOSE(F9(1,x), F9_bare(1.00001,x), 1e-2);
x = 0.0001; BOOST_CHECK_CLOSE(F9(1,x), F9_bare(1.00001,x), 1e-2);

x = 0.1; BOOST_CHECK_CLOSE(F9(x,0), F9_bare(x,0.00001), 1e-2);
x = 0.02; BOOST_CHECK_CLOSE(F9(x,0), F9_bare(x,0.00001), 1e-2);
x = 0.01; BOOST_CHECK_CLOSE(F9(x,0), F9_bare(x,0.00001), 1e-2);
x = 0.001; BOOST_CHECK_CLOSE(F9(x,0), F9_bare(x,0.00001), 1e-2);

x = 0.1; BOOST_CHECK_CLOSE(F9(0,x), F9_bare(0.00001,x), 1e-2);
x = 0.02; BOOST_CHECK_CLOSE(F9(0,x), F9_bare(0.00001,x), 1e-2);
x = 0.01; BOOST_CHECK_CLOSE(F9(0,x), F9_bare(0.00001,x), 1e-2);
x = 0.001; BOOST_CHECK_CLOSE(F9(0,x), F9_bare(0.00001,x), 1e-2);

BOOST_CHECK(!std::isnan(F9(0,1)));
BOOST_CHECK(!std::isnan(F9(1,0)));
BOOST_CHECK(!std::isnan(F9(1,1)));
Expand Down

0 comments on commit dce90c2

Please sign in to comment.