Posted by
Kakhkhor Abdijalilov on
Aug 10, 2010; 2:49pm
URL: http://quantlib.414.s1.nabble.com/gaussianorthogonalpolynomial-cpp-double-to-bool-conversion-and-possible-problem-on-x64-targets-tp13296.html
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.
------------------------------------------------------------------------------
This SF.net email is sponsored by
Make an app they can't live without
Enter the BlackBerry Developer Challenge
http://p.sf.net/sfu/RIM-dev2dev
_______________________________________________
QuantLib-dev mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/quantlib-dev