I have the following code (I'm sorry for the lengthiness):
double primeValue( const func1D &func,
const double lowerBound, const double upperBound,
const double pole )
{
// check bounds
if( lowerBound >= upperBound )
throw runtime_error( "lowerBound must be smaller than upperBound!" );
// C++0x way of writing: fullFunc(x) = func(x)/(x-a)
func1D fullFunc =
bind( divides<double>(), // division of
bind(func, _1), // f(x), with _1 = x
bind(minus<double>(), _1, pole) ); // by x-a, with _1 = x
// pole not in domain
if( pole<lowerBound || pole>upperBound)
{
cout << "Case 1" << endl;
return integrateSimpson( fullFunc, 1000, lowerBound, upperBound );
}
// pole closer to upper bound
else if( upperBound-pole < pole-lowerBound )
{
cout << "Case 2" << endl;
// C++0x way of writing g(x) := [f(x)-f(2a-x)]/(x-a)
func1D specialFirstFunc =
bind( std::divides<double>(), // division of
bind(minus<double>(), // numerator:
bind(func, _1), // f(x) minus
bind(func, bind(minus<double>(), 2.*pole, _1))), //f(2a-x)
bind(minus<double>(), _1, pole) ); // denominator: x-a
const double trickyPart = integrateSimpson( specialFirstFunc, 1000, pole+.000001, upperBound );
const double normalPart = integrateSimpson( fullFunc, 1000, lowerBound, 2.*pole-upperBound );
cout << "Tricky part: " << trickyPart << endl;
cout << "Normal part: " << normalPart << endl;
return trickyPart + normalPart;
}
else // pole closer to lower bound
{
cout << "Case 3" << endl;
// C++0x way of writing g(x) := [f(x)-f(2a-x)]/(x-a)
func1D specialFirstFunc =
bind( std::divides<double>(), // division of
bind(minus<double>(), // numerator:
bind(func, _1), // f(x) minus
bind(func, bind(minus<double>(), 2.*pole, _1))), //f(2a-x)
bind(minus<double>(), _1, pole) ); // denominator: x-a
const double trickyPart = integrateSimpson( specialFirstFunc, 1000, lowerBound, pole-.00001 );
const double normalPart = integrateSimpson( fullFunc, 1000, 2.*pole-lowerBound, upperBound );
cout << "Tricky part: " << trickyPart << endl;
cout << "Normal part: " << normalPart << endl;
return trickyPart + normalPart;
}
}
It integrates functions over the real axis that contain a singularity (pole) using the principal values concept from the math domain of Complex Analysis. The bind
and function
parts modify the original function f(x) to
(f(x)-f(2*pole-x))/(x-a)
It even gives he correct result for my simple test case function. Additional details I can provide if requested:
typedef std::function<double (double)> func1D;
double integrateSimpson( func1D, const size_t nSteps, const double lowerBound, const double upperBound);
The latter integrates using the simple Simpson integration rule. Code can be provided, but isn't very relevant to the problem at hand.
This compiles fine with GCC 4.4+ (tested with 4.4.5 and 4.5.2 prerelease, CFLAGS="-O2 -std=c++0x -pedantic -Wall -Wextra"), but produces internal header errors (C2664) on MSVC 2010. (I can provide error output if needed, there are no references at all to my code (!)).
Have I found a bug in MSVC?
Thanks!