6

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

R2B2
  • 1,413
  • 1
  • 11
  • 17
  • 1
    Floating-point is not necessarily slower than integer. See https://stackoverflow.com/questions/2550281/floating-point-vs-integer-calculations-on-modern-hardware – llllllllll Feb 03 '18 at 11:41
  • 2
    Counter-example to 5: `n = 4; p = 2; q = 1;` Expected outcome: `sqrt(n) == p/q;` Actual outcome: `(1 / 2) * 4` with integer division is `4`, while `2/1` is `2`, yielding `4 <= 2` which is false. – Niet the Dark Absol Feb 03 '18 at 11:42
  • What is the criteria by which you reckon one solution is *better* than another ? Or what are the criteria ? And are you, as your question hints, concerned only with integers representable in 32 or 64 bits ? – High Performance Mark Feb 03 '18 at 12:09
  • @HighPerformanceMark Criteria would be efficiency and elegance. Concerning the integers: the problem only really arises with 64 bits in my settings. If I have integers of 32 bits or less, I can use 64 bits to do my computations and I do not have overflow problems. If I have integers of more than 64 bits, again overflows would not be a problem (since I would be using a bignum library). But the problem could still be translated to any size: when your inputs have _n_ bits, but you cannot use more than _2*n-1_ bits for your computations. – R2B2 Feb 03 '18 at 13:05
  • Without prior knowledge of a bound on n,p,q, the answer is 4, use bignum – aka.nice Feb 03 '18 at 16:38

2 Answers2

3

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.

aka.nice
  • 8,465
  • 1
  • 25
  • 39
1

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);
Veedrac
  • 50,755
  • 12
  • 100
  • 159