1

What I'm asking about is not a duplicate of this very popular question. For a randomly chosen input, some quick tests can be done and if they fail to say "not a square", some computation of the square root has to be done (I myself tried a solution, too).

When the numbers to be tested come from a simple sequence, the situation differs, as it's possible to use the previous (approximate) square root. For a trivial sequence it's trivial as well, e.g.,

long sqrt = 1;
for (long i=1; i<limit; ++i) {
   if (sqrt*sqrt == i) {
       handleSquare(i);
       ++sqrt;
   }
}

My question is what can be done for more complicated sequences like

x[i] = start + i*i;

or

x[i] = start - i*i*i;

I'm thinking about Newton's method, but I can't see how to make it fast (as division is a pretty expensive operation).

Community
  • 1
  • 1
maaartinus
  • 44,714
  • 32
  • 161
  • 320
  • I think that your "trivial" example is a little *too* trivial. Think of it from this angle: the square root of a number that is a perfect square is rational; that is, there is no decimal component to the number. – Makoto Apr 27 '15 at 02:43
  • @Makoto Right, it's too trivial as everybody would instead loop over `sqrt`. But forget it, the idea is that *keeping some information* between iterations helps. – maaartinus Apr 27 '15 at 02:52

1 Answers1

0

To what kind of sequence would you like to apply your algorithm ? Below is a solution that should work well when x[i] diverges but not too fast.

For instance if

x[i] = a*i^p + o(i^p)

and that i is large enough then you will have

x[i+1]-x[i] ~ p * a * i^{p-1}.

If y[i] denotes the largest integer such that

y[i]^2 <= x[i]

then you have

y[i] ~ sqrt(a) i^{p/2} 

and

y[i+1]-y[i] ~ 1/(2 y[i]) * (x[i+1]-x[i]) ~ p/2 * sqrt(a) i^{p/2-1}

So you can take this as a guess for y[i+1] and then update to the correct value, which should save you some iterations.

In general you can always use the formula

y[i+1]-y[i] ~ 1/(2 y[i]) * (x[i+1]-x[i])

as a guess but this will be useful only if x[i+1]-x[i] is small with respect to y[i]^2 -- ie, with respect to x[i]. It may also be worth improving a bit on the formula using the (exact) second order expansion of

y[i+1]^2 = y[i]^2 + 2y[i](y[i+1]-y[i]) + (y[i+1]-y[i])^2

in order to improve the guess of y[i+1].

Note that this won't work well if x[i] remains bounded when i is large or if x diverges exponentially fast.

vib
  • 2,254
  • 2
  • 18
  • 36