Re: gaussianorthogonalpolynomial.cpp : double to bool conversion and possible problem on x64 targets

Posted by Kim Kuen Tang on
URL: http://quantlib.414.s1.nabble.com/gaussianorthogonalpolynomial-cpp-double-to-bool-conversion-and-possible-problem-on-x64-targets-tp13296p13297.html


Hi Kakhkhor,

Kakhkhor Abdijalilov schrieb:
> 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.
many polynomials functions are already implemented in boost math
toolkit. Can you try using these functions?

Kim

> 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
>
>  


------------------------------------------------------------------------------
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