Comparing sqrt(n) with the rational p/q

Question

You are given an integer n and a rational p/q (p and q are integers).

How do you compare sqrt(n) and p/q?

Solution 1: sqrt(n) <= (double) p / q
Should work, but calls sqrt which is slower than just using multiplication/division.

Solution 2: (double) n * q * q <= p * p
Better, but I can't help thinking that because we are using floats, we might get an incorrect answer if p/q is very close to sqrt(n). Moreover, it requires converting integers to floats, which is (marginally) slower than just working with integers.

Solution 3: n*q*q <= p*p
Even better, but one runs into trouble if p and q get big because of overflow (typically, if p or q >= 2^32 when working with 64 bits integers).

Solution 4: Use solution 3 with a bignum library / in a programming language that has unbound integers.

Solution 5: (q / p) * n <= p / q
Successfully avoids any overflow problems, but I am not sure that this is correct in all cases, because of integer division...

So... I would happily go with solution 2 or 4, but I was wondering if anyone has clever tricks to solve this problem or maybe a proof (or counter example) that solution 5 works (or not).

Answer 1

As I commented, a simple and elegant solution is to use bignum, especially if builtin, or easily available in the chosen language. It will work without restriction on n,p,q.

I will develop here an alternate solution based on IEEE floating point when:

• n,p,q are all representable exactly with the given floating point precision (e.g. are within 24 or 53 bits for single or double IEEE 754)
• a fused multiply add is available.

I will note f the float type, and f(x) the conversion of value x to f, presumably rounded to nearest floating point, tie to even.

fsqrt(x) will denote the floating point approximation of exact square root.

let f x = fsqrt(f(n)) and f y = f(p) / f(q).

By IEEE 754 property, both x and y are nearest floating point to exact result, and n=f(n), p=f(p), q=f(q) from our preliminary conditions.

Thus if x < y then problem is solved sqrt(n) < p/q.

And if x > y then problem is solved too sqrt(n) > p/q.

Else if x == y we can't tell immediately...

Let's note the residues f r = fma(x,x,-n) and f s = fma(y,q,-p).

We have r = x*x - n and s = y*q - p exactly. Thus s/q = y - p/q (the exact operations, not the floating point ones).

Now we can compare the residual errors. (p/q)^2 = y^2-2*y*s/q+ (s/q)^2. How does it compare to n = x^2 - r?

n-(p/q)^2 = 2*y*s/q - r - (s/q)^2.

We thus have an approximation of the difference d, at 1st order: f d = 2*y*s/f(q) - r. So here is a C like prototype:

int sqrt_compare(i n,i p,i q)
/* answer -1 if sqrt(n)<p/q, 0 if sqrt(n)==p/q, +1 if sqrt(n)>p/q */
/* n,p,q are presumed representable in f exactly */
{
  f x=sqrt((f) n);
  f y=(f) p / (f) q;
  if(x<y) return -1;
  if(x>y) return +1;
  f r=fma(x,x,-(f) n);
  f s=fma(y,(f) q,-(f) p);
  f d=y*s/(f) q - r;
  if(d<0) return -1;
  if(d>0) return +1;
  if(r==0 && s==0) return 0; /* both exact and equal */
  return -1; /* due to 2nd order */
}

As you can see, it's relatively short, should be efficient, but is hard to decipher, so at least from this POV, I would not qualify this solution as better than trivial bignum.

Answer 2

You might consider solution 3 with integers 2x the size,

n * uint2n_t{q} * q <= uint2n_t{p} * p

This overflows if n * q * q overflows, but in that case you return false anyway.

uint2n_t nqq;
bool overflow = __builtin_mul_overflow(uint2n_t{n} * q, q, &nqq);
(!overflow) && (uint2n_t{n} * q * q <= uint2n_t{p} * p);