Search code examples
c++mathsqrt

How to prevent overflow in sqrt?


I have the code

unsigned long long f(unsigned long long x, unsigned long long y){
   return sqrt( ( (5*x*x + 1) * (5*y*y +1) - 1 ) / 4 );         
}

but if x or y is too big, the thing overflows even though the output is supposed to be small. is there a way around this?


Solution

  • Split the argument of the square root up algebraically, maybe?:

    return sqrt((3*x*x+1)/2) * sqrt((6*y*y-5)/2);
    

    Or split it up further based on your needs.

    If x is big enough you can ignore the +1 and make the first term:

    sqrt((3*x*x)/2) = fabs(x) * sqrt(3.0/2.0)
    

    and similarly with y for the second term, making it

    sqrt((6*y*y)/2) = fabs(y) * sqrt(3.0);
    

    EDIT: After OP edited his question to be:

    return sqrt(((3*x*x+1)*(6*y*y-5)-1)/4);  
    

    You can in fact split things up. You just have to be a little bit more careful. Bottom line is that if x is really big, then the +1 can be ignored. If y is really big then the -5 can be ignored. If both (3*x*x+1) and (6*y*y-5) are positive and either is really big, then the -1 can be ignored. You can use these tips and some additional surrounding logic to break this down a bit further. Like this:

     if(fabs(x) > BIGNUMBER && fabs(y) > BIGNUMBER)
     {
         return fabs(x) * fabs(y) * sqrt(18.0/4.0);
     }
     if(fabs(x) > BIGNUMBER && fabs(y) > 1.0)  // x big and y term positive
     {
         return fabs(x) * sqrt(6*y*y-5) * sqrt(3.0/2.0);
     }
     if(fabs(y) > BIGNUMBER)  // x term positive and y big
     {
         return sqrt(3*x*x+1) * fabs(y) * sqrt(6.0/2.0);
     }
     return sqrt(((3*x*x+1)*(6*y*y-5)-1)/4); 
    

    You can optimize this, but this is just meant to show the point.