views:

56

answers:

1

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!

A: 

Why not just use a lambda? All of the binding stuff has been deprecated for this kind of purpose.

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)
    auto fullFunc = [=](double d) {
        return func(d) / (d - pole);
    };

    // 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)
       auto specialFirstFunc = [=](double x) -> double {
           double numerator = func(x) - func(2*pole - x);
           return numerator / (x - pole);
       };
        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)
       auto specialFirstFunc = [=](double x) -> double {
           double numerator = func(x) - func(2*pole - x);
           return numerator / (x - pole);
       };
         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;
    }
}

As for whether or not you actually found a bug in MSVC, I don't know, but your solution is definitely way unnecessary - this code is far cleaner and easier to maintain.

DeadMG
Well, as most linux distro's and myself still use GCC 4.4, lambda's aren't an option. I don't think `std::bind` has been deprecated already... But the lambda's are way shorter than demonstrated to me previously: http://stackoverflow.com/questions/4008369/combining-stdfunction-objects
rubenvb
@rubenvb: Lambdas were specifically introduced into the language because bind sucked, among other reasons. Just slap an #ifdef on it. Most likely, MSVC's libraries are somewhat out of date with what C++0x now has as the Standard - for example, _1 is defined as std::placeholders::_1 instead of just std::_1. I checked that link - he explicitly set reference or value for every variable he captured, whereas I set mine to capture all by value.
DeadMG
@DeadMG: Well, the lambda workaround sounded nice, but I think the internal implementation uses `std::bind` anyways (it gives the same error :s). I checked the `#ifdef` branches used, and it's ccompiling the lambda's using `std::bind` internally...
rubenvb
@ruben: I tested on MSVC10 and can assure you that the code posted above compiles in MSVC10.
DeadMG
@DeadMG: OK, I'll test it again as soon as I have time.
rubenvb
I had another nested bind function code bit, which I have replaced with its lambda equivalent and now the code compiles fine under MSVC 10. Thanks again!
rubenvb