[Quantlib-dev] gaussianorthogonalpolynomial.cpp : double to bool conversion and possible problem on x64 targets

[Kakhkhor A. <kab...@gm...>]

I compiled QuantLib with ICC for x64 target. Everything went well
except a glitch with Null class (I reported it in my previous e-mail)
and some issues in gaussianorthogonalpolynomial.cpp. Current
implementations of GaussJacobiPolynomial::alpha and
GaussJacobiPolynomial::beta have several issues.

1) On x64 targets Size is 64bit and conversion to double may generate
a warning. For some strange reason ICC refuses to compile this line
Real denom = (2.0*i+alpha_+beta_)*(2.0*i+alpha_+beta_+2);
ICC bails out with an internal error. I don't know why it happens (ICC
bug?), but casting i from Size into int fixes the problem.
2) Current implementation checks if num and denom variables are zero
and then tries to use l'Hospital rule. But denom==0 can be true only
if alpha_=0 and beta_=0 and/or i=0. In that case num=0 and l'Hospital
rule is never used. In anyway, direct conversion of floating point
numbers into bool is a fishy business at best. We should consider
using fuzzy comparison instead of direct comparison, something like
this:

//-----------------------------------------------------------------------------------
    Real GaussJacobiPolynomial::alpha(Size iIn) const {
        int i = static_cast<int>(iIn);
        Real num = beta_*beta_ - alpha_*alpha_;
        Real denom = (2.0*i+alpha_+beta_)*(2.0*i+alpha_+beta_+2);

        if (std::fabs(denom)<QL_EPSILON) {
            if (std::fabs(num)<QL_EPSILON) {
                QL_FAIL("can't compute a_k for jacobi integration\n");
            }
            else {
                // l'Hospital
                num  = 2*beta_;
                denom= 2*(2.0*i+alpha_+beta_+1);
                QL_ASSERT(std::fabs(num)>QL_EPSILON, "can't compute
a_k for jacobi integration\n");
            }
        }

        return num / denom;
    }

    Real GaussJacobiPolynomial::beta(Size iIn) const {
        int i = static_cast<int>(iIn);
        Real num = 4.0*i*(i+alpha_)*(i+beta_)*(i+alpha_+beta_);
        Real denom = (2.0*i+alpha_+beta_)*(2.0*i+alpha_+beta_)
                   * ((2.0*i+alpha_+beta_)*(2.0*i+alpha_+beta_)-1);

        if (std::fabs(denom)<QL_EPSILON) {
            if (std::fabs(num)<QL_EPSILON) {
                QL_FAIL("can't compute b_k for jacobi integration\n");
            } else {
                // l'Hospital
                num  = 4.0*i*(i+beta_)* (2.0*i+2*alpha_+beta_);
                denom= 2.0*(2.0*i+alpha_+beta_);
                denom*=denom-1;
                QL_ASSERT(std::fabs(num)>QL_EPSILON, "can't compute
b_k for jacobi integration\n");
            }
        }
        return num / denom;
    }
//-----------------------------------------------------------------------------------

QL_EPSILON can be replaced with something more suitable, but at least
it should be OK in most of the cases.


Regards,
Kakhkhor Abdijalilov.