34
$\begingroup$

Is there an efficient method to determine if a very large number (300 - 600 digits) is a perfect square without finding its square root. I tried finding the square root but I realized that even for perfect squares, it wasn't efficient (I used newton's approximation method and it was coded in JAVA). I read somewhere about using modular arithmetic but I'm not quite sure how to do that. Can anyone give me a hint on how to proceed?

Thanks in advance.

Improve this question
$\endgroup$
4
  • 2
    $\begingroup$ You can efficiently compute the Jacobi symbol of $n$ relative to $q$ for any $q$ that is relatively prime to $n$; if the result is $-1$, then you know that $n$ is not a square, but if the result is $1$ then the test is inconclusive. Is this the kind of thing you might have read? $\endgroup$ – Arturo Magidin May 25 '11 at 19:12
  • 3
    $\begingroup$ You need to define "efficient" algorithm. Are you seeking to minimize the maximum time to determine if it is a square, or the expected time, or something else? There are some quick observations that could tell you it is not a square, e.g. if the last digit is 2,3,7 or 8 then it is not a square. Using such methods you can reduce the expected time but probably not the maximum time. $\endgroup$ – Philip Gibbs May 25 '11 at 21:30
  • $\begingroup$ I know this is an old question, but I seem to have found a near log-time method, linked to in my answer below! This is much better than Newton. $\endgroup$ – user253804 Jul 13 '15 at 23:43
  • $\begingroup$ does this not discuss the criteria. I'm confused $\endgroup$ – ldgorman Apr 20 '17 at 14:06

10 Answers 10

Active Oldest Votes
22
$\begingroup$

Faster than binary search is to use an integer version of Newton's method: for $\sqrt{A}$ and an initial guess $x_0$ of the right order of magnitude, let $x_{n+1} = \left \lfloor \frac{x_n^2 + A}{2 x_n} \right \rfloor$. I tried a 1200-digit number for $A$, with $x_0$ a 600-digit number, and $x_{10}$ was already correct. In Maple 15 on a not-very-fast computer, the CPU time was given as 0.

Improve this answer
$\endgroup$
12
  • 3
    $\begingroup$ Are you sure it's faster? At each step you need to perform one division. $\endgroup$ – Yuval Filmus May 25 '11 at 21:23
  • 1
    $\begingroup$ Once you get two consecutive equal values, you must repeat from there. The recursion only looks back one step! $\endgroup$ – David E Speyer May 25 '11 at 22:00
  • 3
    $\begingroup$ @Yuval Filmus: Yes, but that's pretty fast for numbers of this size. Maple uses the GMP library. An integer division of a 1200-digit number by a 600-digit number might take about .0001 second. @obinna: I actually tested for 2 consecutive equal values. Note that $x_n = x_{n-1}$ iff $x_n^2 \le A < x_n^2 + 2x_n$. You actually run into an infinite loop with $x_n$ alternating between $\lfloor \sqrt{A} \rfloor$ and $\lceil \sqrt{A} \rceil$ if $A+1$ happens to be a square, so I should really have checked for that, but I didn't bother. $\endgroup$ – Robert Israel May 25 '11 at 22:05
  • 3
    $\begingroup$ When I wrote this in GMP (a C/C++ library Gnu Multiple Precision) I found it safer to stop when, say, $$ | x_{n} - x_{n-1} | \leq 5 $$ and then adjust what you have by increments of $\pm 1.$. The problem is that is is quite possible to get a period 2 cycle repeating forever. $\endgroup$ – Will Jagy May 26 '11 at 2:00
  • 5
    $\begingroup$ @YuvalFilmus With the right algorithms, this approach is asymptotically faster; binary search for an $n$-digit square root takes $\Theta(n^2)$ bit operations, but it can be proven that this Newton method can be made to converge in $\Theta(\log n)$ operations, and division is known to be asymptotically as fast as multiplication, so this method is $\Theta(M(n)\log n)$ where $M(n)$ is the time to multiply two $n$-bit numbers; it's known that fast algorithms can get $M(n)$ down to $O(n\log n\log\log n)$. $\endgroup$ – Steven Stadnicki Nov 12 '13 at 9:11
11
$\begingroup$

I agree with Yuval Filmus that you should use binary search, but in case you're still curious about the approach using modular arithmetic: a number $a$ is a square $\bmod p$ for $p$ a prime not dividing $a$ if and only if the Jacobi symbol $\left( \frac{a}{p} \right)$ is equal to $1$. You can efficiently compute the Jacobi symbol using quadratic reciprocity, and if you get an answer of $1$ for $n$ primes, then $a$ is square with probability about $1 - \frac{1}{2^n}$.

(Alternately, $a$ is a square $\bmod p$ if and only if $a^{ \frac{p-1}{2} } \equiv 1 \bmod p$. I don't know how the efficiency of computing this compares to the efficiency of computing the Jacobi symbol.)

Improve this answer
$\endgroup$
1
  • 6
    $\begingroup$ +1. Also, searching for small prime factors is worthwhile -- a number divisible by 2 but not 4, 3 but not 9, etc., can't be a square. Further, some residue classes like 2 mod 3 are impossible. It's often worthwhile to check mod a precomputed number and use a lookup table to see if the number is a possible square mod that number. Continuing along this path leads to pseudosquares and the research of Sorenson and such. $\endgroup$ – Charles May 25 '11 at 19:28
8
$\begingroup$

The square root can be found using binary search. The resulting algorithm should run very fast on a 600 digit number, my guess is under a second.

You can implement the binary search without repeated squaring - each step you're only shifting and adding. That's why it's so fast. But even if you were squaring at each step, it would still be very quick and certainly feasible.

Any reasonable bignum package will contain a function computing the square root, so you don't even need to code the binary search yourself.

Improve this answer
$\endgroup$
8
  • $\begingroup$ Try mpz_sqrt() if you use gmp. $\endgroup$ – Yuval Filmus May 25 '11 at 19:14
  • $\begingroup$ JAVA's Big Integer package has no square root function. What do you mean by the phrase "shifting and adding". I mean how can I determine whether the test case is greater than or less than the number without squaring it? $\endgroup$ – Obinna Okechukwu May 25 '11 at 19:23
  • 1
    $\begingroup$ @obinna: until you get within a factor of $2$ or so you can compare your guess with the actual number just by comparing the floors of their base-$2$ logarithms (particularly easy to do on a computer). "Shifting and adding" refers to the process for generating new guesses: see www2.lv.psu.edu/ojj/courses/cmpsc-201/numerical/bisection.html . $\endgroup$ – Qiaochu Yuan May 25 '11 at 19:32
  • $\begingroup$ @Qiaochu Yuan: If you have convenient access to the floor of the base-2 logarithm you can get within a factor of 2 directly with just shifts: 1 << (L >> 1). $\endgroup$ – Charles May 25 '11 at 19:45
  • 3
    $\begingroup$ @obinna: I'm afraid your comment reveals that you haven't really thought about the problem. If your number $n$ is even, just count the number $b$ of trailing binary zeroes. If $b$ is odd, then $n$ is certainly not a square. If $b$ is even, then $n$ is a square if and only if the (odd!) number $n/(2^b)$ is. So you gain nothing from the knowledge that $n$ is odd, because the even case can trivially be reduced to the odd case. $\endgroup$ – TonyK May 25 '11 at 21:06
5
$\begingroup$

There is direct analogue of algorithm mentioned here http://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Digit-by-digit_calculation for arbitrary big integer. For binary representation it is only needed subtractions, comparisons and shifts. So, despite it finding integer part of square root precisely, I doubt it is possible to find something faster - the complexity level of the algorithm is comparable with integer division of two numbers. I suspect that Yuval Filmus also mentioned something similar in note about “shifting and adding” - here instead of that is shifting and subtracting.

Improve this answer
$\endgroup$
1
  • $\begingroup$ "For binary representation" only mean internal representation of big integer on computer. Example of the algorithm in wikipedia need only small modification to work with arbitrary size integers. $\endgroup$ – Alex 'qubeat' May 26 '11 at 9:31
5
$\begingroup$

I'm late to the party. But, 2-adic Newton is pretty fast. Here's a 64-bit version:

/* If x is square, return +sqrt(x).
 * Otherwise, return -1.
 */
long long isqrt(unsigned long long x) {
    /* Precondition: make x odd if possible. */
    int sh = __builtin_ctzll(x);
    x >>= (sh&~1);

    /* Early escape. */
    if (x&6) return -1;

    /* 2-adic Newton */
    int i;
    const int ITERATIONS = 5; // log log x - 1
    unsigned long long z = (3-x)>>1, y=x*z;
    for (i=0; i<ITERATIONS; i++) {
        unsigned long long w = (3 - z*y) >> 1;
        y *= w;
        z *= w;
    }
    assert(x==0 || (y*z == 1 && x*z == y));

    /* Get the positive square root.  Also the top bit
     * might be wrong. */
    if (y & (1ull<<62)) y = -y;
    y &= ~(1ull<<63);

    /* Is it the square root in Z? */
    if (y >= 1ull<<32) return -1;

    /* Yup. */
    return y<<(sh/2);
}

The function listed above takes 31 Haswell cycles when inlined. So it's about half as fast as the double-precision sqrtsd() call. It computes the square root of its argument, which is not what OP asked for. But it takes advantage of the fact that you only care about integer square roots. As a result, it will be much faster on extended-precision integers than an approximate square root algorithm.

This works for longer integer types, with a few changes.

  • It's a fixed-precision algorithm. If you use mpn_t for this, it will allocate new limbs instead of wrapping which is not what you want. You will need to explicitly clamp to the lowest n limbs.
  • You need log log x - 1 ITERATIONS. So if x is 2048 bits, you need 10 instead of 5. That is, it's Newton's method so each iteration doubles the precision.
  • You can compute the first 4 iterations with longs if you want.
  • You need to replace the 32 with half the number of bits of your integer type, likewise 62 and 63 with 2 less and 1 less.
$\endgroup$
6
  • 1
    $\begingroup$ This is brilliant! Thank you so much for posting this. On my system, this test is actually 5 times faster than using 80-bit hardware square root, e.g., uint64_t y = (uint64_t)sqrtl((long double)x); bool is_square = (y * y == x); across the whole range from 0 all the way up to 2^64-1. I'm amazed!!! $\endgroup$ – Todd Lehman Sep 9 '14 at 7:57
  • $\begingroup$ Question: From what is this derived? You say it's "2-adic Newton," which I'm not familiar with. I'll go read about it, but in particular I'm wondering why strange things like 3-x, which underflows an unsigned value, and the x&6 != 0 test. It's like magic and I'd like to understand what's really going on under the hood. $\endgroup$ – Todd Lehman Sep 9 '14 at 8:00
  • $\begingroup$ By the way, what happens when x==0 on input? My compiler doesn't have __builtin_ctzll(), so instead I did int sh = 0; while (x % 4 == 0) { x /= 4; sh++; } and for the return part did return y << sh. The while loop gets stuck forever if x == 0, so I check that specially before doing anything else. Just curious what happens in your original version using the trailing zero count. Does it right-shift x by 64 and then fall through with x being 0? $\endgroup$ – Todd Lehman Sep 9 '14 at 8:05
  • $\begingroup$ Also curious: Your comment with the if (x&6) return -1; statement says "Early escape" — implying an optimization shortcut. But I noticed that if I comment it out, it actually returns incorrect result for numbers like 80, 128, 192, 320, 512, 768, 1280, 2048, etc. So it's not really so much an early-out test, but rather is fundamental and crucial to the algorithm, correct? $\endgroup$ – Todd Lehman Sep 9 '14 at 8:24
  • 1
    $\begingroup$ @ToddLehman: It's Newton's algorithm, computing in the $p$-adic field $\mathbb{Q}_2$. Since you're working with $2$-adic integers, working with 64-bits of precision is identical to doing arithmetic modulo $2^{64}$. I.e. addition, subtraction, and multiplication of uint64_t does exactly the right thing. $\endgroup$ – user14972 Aug 26 '18 at 10:02
5
$\begingroup$

Since there seems to be some confusion about this, here is the almost trivial pseudo-code to test for integer square roots in a language that has bigints, and integer division without remainder that I'll write div (as / might be confused with exact rational division):

boolean is_square(bigint n)
{ bigint a=n // or any value whose square is known to be at least n
; while a*a>n do a=div(a+div(n,a),2) od
; return a*a==n
}

No cycles are possible, no fuss (but negative numbers as input may cause it to loop forever, so that should be avoided). The outermost division is just a bit-shift, so it should be rather inexpensive. The clue to the correctness of this simple procedure it that, when the initial guess $a$ has $a^2>n$, then the exact value $b=\frac12(a+\frac na)\in\Bbb Q$ can be shown to also have $b^2>n$; if $\sqrt n$ is integer, then rounding down $b$ to an integer either gives exactly $\sqrt n$, or another too high estimate. If rounding down gives a too small value, then one can be sure that $\sqrt n$ is not integer. One also clearly has $b<a$, so the sequence (with rounding down at each step) is strictly decreasing and termination is assured.

For very large numbers the given code wastes a lot of time to get a ballpark-figure for the square root. This could be avoided by using a better initial guess for $a$, using for instance a logarithm, or basically anything that tells you the (approximate) number of bits in $n$ (there might well be a built-in function for that). But I think any pure integer arithmetic method to do this will suffer from the same problem.

As a numerical example, for a $999$-digit perfect square as input it took $1667$ iterations to find its square root. But just repeated halving would take about $\log_2(10^{499})\approx 1658$ iterations to get down to the proper order of magnitude, so that is what the algorithm is essentially doing most of the time. If instead of $n\approx10^{998}$ one takes an initial guess of $10^{500}$ (still about $10$ times larger than the square root) it takes only $13$ iterations to find the square root. The last iteration step hits the square root on the head from a distance of more than $10^{142}$ (that is, the error of the approximation goes from that gigantic number directly to $0$, which is a witness of how much Newton's method improves the approximation here). By contrast, a binary search procedure even with a fairly good initial guess like $10^{500}$ would need as number of iterations at least $1658$, namely the number binary digits of the square root (since it adds one binary digit precision each iteration). For the record, the input used in the experiment was the square of $10^{499}+345674632452435$ (you can guess the type of "random generation" procedure I used for it).

$\endgroup$
2
  • $\begingroup$ Apparently what I did was make an initial guess of $1,$ and not test $(n-a^2)$ during update loop. With your pseudocode letters, the new $a$ is $$ \left\lfloor \frac{a + \left\lfloor \frac{n}{a} \right\rfloor}{2} \right\rfloor. $$ In the very special case $n = b^2 + 2 b = (b+1)^2 -1,$ when $a=b$ the update is $a=b+1,$ hence cycle, as the next such update gives $a=b$ again. I certainly made no special effort to give a seed above the correct root. So, your way is better, but that's not what I actually did. $\endgroup$ – Will Jagy Nov 12 '13 at 20:24
  • 1
    $\begingroup$ In Pari/GP one can get the initial guess for $a$ by a=1 << (#binary(n) \2) (giving $a=2^{\lfloor \log_2(n)\rfloor /2} $ )which is near the squareroot and does need nearly nothing. For n with 4000-digits I got roughly 13 iterations (with cycle detection) and on average 40 msec, for 2000 digits on average 10 msec with 11 to 12 iterations. (Standard notebook of 2012) $\endgroup$ – Gottfried Helms Jan 19 '18 at 3:38
3
$\begingroup$

In my answer here I've described a probabilistic test for squaredness that can achieve an arbitrarily small error with time complexity very close to $\mathcal O$(log $n$) for an integer $n$. This is much better than Newton et. al.

It's the same style of probabilistic test used to determine primality for huge numbers, except it's based on properties of quadratic residues. It's sort of a generalization of the approach used in the GMP library.

$\endgroup$
2
$\begingroup$

One of the ways to quickly eliminate several possibilities upfront is to use the test whether the sum of the digits of the given number adds up to either 1,4,7,9. The number not satisfying this can't be a perfect square. If it does, then you can carry on with one of the efficient methods discussed in the responses to check whether the number is perfect square.

$\endgroup$
5
  • $\begingroup$ Just to clarify, if $n^2 \equiv p (\mathrm{mod}\ q)$, $p$ is called a quadratic residue modulo $q$. What Atharva is referring to is the fact that the only residues modulo $9$ are $0, 1, 4, 7,$ and $9$. The reason this is useful is that for a number $k$, $k \equiv m (\textrm{mod}\ 9)$ where $m$ is the sum of the digits of $k$. So if $k$ is a perfect square, then $m$, repeatedly reduced modulo $9$, will be one of the five quadratic residues. If you do this process and it is not one of these five numbers, $k$ is not a perfect square. $\endgroup$ – Michael Albanese Sep 22 '12 at 4:05
  • $\begingroup$ The reason why you need to include $0$ as one of the possible residues is that your potential perfect square might be divisible by $9$. For example, $81$ is a perfect square and the sum of its digits modulo $9$ is 0. $\endgroup$ – Michael Albanese Sep 22 '12 at 4:07
  • $\begingroup$ Thanks Michael for suggesting the improvement to my answer in your comments. $\endgroup$ – Atharva Patel Sep 24 '12 at 10:00
  • 1
    $\begingroup$ This is a long shot after such a time , but how does this work for a number such as 45. As the digit sum is 9 so either 9 mod 9 = 0 ; or 45-9 = 36 and this is a multiple of 9 so the residual is also 0. So by the rules 45 is a perfect square. Clearly I have really misunderstood somewhere. Can you explain please? $\endgroup$ – porphyrin Jan 27 '17 at 15:43
  • $\begingroup$ @porphyrin this is for eliminating numbers which cannot be perfect squares, says nothing about whether they are or not $\endgroup$ – Sridhar Thiagarajan Feb 17 '19 at 7:24
2
$\begingroup$

An algorithm based on subtractions only.

Let $n$ be the number, $a=5n$ and $b=5$. While $a\geq b$, replace $a$ with $a-b$ and add $10$ to $b$. When you get $a<b$, then $n$ is a square if and only if $a=0$.

$\endgroup$
1
  • 1
    $\begingroup$ This answers lacks explanation, says nothing about time complexity of the algorithm, and introduces an useless factor 5. The reason why it works is that the loop preserves the following invariant: $a_i = 5(n - i^2) \wedge b_i = 5(1+2i)$ where $i$ is the number of previous iterations. Summing the sequence of odd integers gives the sequence of square numbers, this algorithm is simply a disguised way of doing that until $n$ is reached. So you would compute for $\mathcal{O}(\sqrt{n})$ steps. This is exponential time, terribly inefficient! The OP was concerned with 600-digit numbers. $\endgroup$ – Maëlan Sep 28 '19 at 12:58
-1
$\begingroup$

This is my way of finding square roots larger than $100$. Take $1089$ for example. You first find the closest number (without going over) that has a square root divisible by $10$. That would be $900$; $30 \times 30$. The tens digit is a $3$. Next subtract that number leaving me with $189$. Now take the square root of $900$. Take the $10$s digit, multiply from $20$, and add $1$. That would be $61$. Subtract that number from $189; 128$. Now subtract $63, 65, 67$, etc until you reach $0$. Including $61$, count the amount of numbers that you subtracted that is the units digit. For this problem it would be $189-61-63-65=0$, so the units is $3$. Your final answer is $\sqrt{1089}=33$.

$\endgroup$
1
  • 2
    $\begingroup$ As noted by Ross Millikan in a comment to this duplicate answer, this is just the standard digit-by-digit square root algorithm. The OP wanted a way without computing the square root. $\endgroup$ – robjohn♦ Apr 23 '13 at 23:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.