How to prevent overflow in sqrt?

Question

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?

Congdon

How exact of an answer do you need? – Aesthesia

Rename your function...that is terribly unobvious what it is supposed to do. And what templatetypedef said. – Capper

@Brendan: Has it occurred to you that this might be part of some propriety code, and that using the real name for f might violate the OP's NDA? For example purposes, f is fine. – Mercaptopurine

@Brendan: ...Or that it simply serves as an example and, in this case, the name of the function is completely irrelevant? – Monomolecular

sqrt won't overflow. Surely you must mean the computations leading up to it. – Tundra

@NathanErnst yes. I am just wondering if there is a programmatical way to deal with this or if I have to mess with the algebra – Congdon

Have you considered posting this on the math stackexchange as well? – Orthodox

Are there any constraints on the relative sizes of x and y? Or is it possible that one is very large and the other very small (in absolute value)? – Mckinney

@DanielFischer X tends to be small, Y tends to be very large. I've actually changed the numbers in my example a bit but the general principle applies. – Congdon

In that case, just pull y/2 out of the square root. Mathematically, sqrt(e*y*y) == abs(y)*sqrt(e), so sqrt(((5*x*x+1)*(5*y*y+1)-1)/4) becomes fabs(y)*0.5*sqrt((5*x*x+1)*(5 + 1/(y*y)) - 1/(y*y)). For most uses that's precise enough and avoids overflow. If that's not precise enough, you'll probably need a different datatype anyway. – Mckinney

@DanielFischer Doesn't 1/(ull*ull) auto to 0 anyway? – Congdon

Ah, missed that, thought they were doubles. Cast to double for the division then, 1/((double)y*y). – Mckinney

@DanielFischer Wow, I think that worked... is there a way to check if the result is a whole number? Or can I only do that by returning double? Will numbers < 10^10 be able to be stored in double? – Congdon

double can store integers < 2^53, roughly 10^15 exactly. But if you want to know whether the expression under the square root is a perfect square, i.e. the square root is an integer, you have to be aware of floating point errors. The computed result may be an integer even though the exact result isn't, and if intermediate results are too large to be represented exactly, it might turn out the other way. What exactly are you trying to do? There may be better methods. – Mckinney

@DanielFischer My function returns the area of a particular form of triangle, but I need to make sure it is integral – Congdon

let us continue this discussion in chat – Mckinney

@AgainstASicilian What's up? Back to the chatroom? – Mckinney

@DanielFischer If you have the time/willingness, of course – Congdon

Answer 1

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.

Boat

Note: Due to integer division, sqrt(3/2) == sqrt(1) == 1. – Nuli

@Robᵩ Good point, I switched to a more symbolic notation, but I'll make this C/C++ correct. – Boat

@AgainstASicilian same concept applies. – Boat

I can't split it into two products because I have a -1 – Congdon

@AgainstASicilian you can probably still do some stuff to break this up algebraically. Just might take more work to do it. – Boat

Unfortunately, I've tried with Wolfram. If I had been able to do it I wouldn't have posted the question! :) – Congdon

@AgainstASicilian See my added note it provides from additional ideas. – Boat

Typical programmer imposing solution on the hapless product owner :-) – Tragicomedy

Answer 2

I believe that at large x or y, the -1 can be neglected w.r.t (x*x)*(y*y) ... and since your function returns long long, the floating point precision won't matter.

You can check whether x or y are big or not, and accordingly, you can choose to neglect the -1 and do as Chris or Rob say.

Rosen