#1
  1. No Profile Picture
    Contributing User
    Devshed Newbie (0 - 499 posts)

    Join Date
    Jul 2011
    Posts
    313
    Rep Power
    0

    Integer square root algorithm


    Hi,

    i'm searching for efficient algorithm which computes integer square root of 64 bit (say, signed) integer number (long long).

    That is: ISQ(n) = max k (integer) such that k*k <= n.

    My approach is stupid.
    I first compute k as integer part of sqrt(n) function (from <math.h>) and then try to increase or decrease k in a loop.

    Thanks in advance for any suggestions.
    Last edited by leszek31417; August 8th, 2011 at 04:15 AM.
  2. #2
  3. Contributed User
    Devshed Specialist (4000 - 4499 posts)

    Join Date
    Jun 2005
    Posts
    4,413
    Rep Power
    1871
    If you dance barefoot on the broken glass of undefined behaviour, you've got to expect the occasional cut.
    If at first you don't succeed, try writing your phone number on the exam paper
  4. #3
  5. No Profile Picture
    Contributing User
    Devshed Newbie (0 - 499 posts)

    Join Date
    Jul 2011
    Posts
    313
    Rep Power
    0
    Thank you.
    Especially the second citation is something that i should read.
  6. #4
  7. No Profile Picture
    Contributing User
    Devshed Newbie (0 - 499 posts)

    Join Date
    Feb 2009
    Posts
    191
    Rep Power
    50
    This little Python function is not optimal, but I think is nonetheless pretty fast:

    Code:
    def isqrt(n):
        # construct an initial guess value
        d = len(str(n))	# count of decimal digits
        x = '3'
        d = d / 2
        while d > 0:
            d = d - 1
            x = x + '1'
        x = int(x)
    
        # perform Newton-Raphson iteration
        prev = 0
        while x != prev:
            prev = x
            q = n / x
            x = (x + q) / 2
    
        return x
    Python doesn't have an integer square root, or at least didn't back when I wrote this function.

    The first block makes a crude approximation to the square root; the second repeats Newton-Raphson until it converges. For extra credit, prove that the while loop will always terminate (or, prove that it can fail to terminate)
  8. #5
  9. No Profile Picture
    Contributing User
    Devshed Newbie (0 - 499 posts)

    Join Date
    Jul 2011
    Posts
    313
    Rep Power
    0
    Thank you again !

    I'll play with this algorithm.

    How many steps the second loop (Newton-Raphson iteration) needs?
  10. #6
  11. No Profile Picture
    Contributing User
    Devshed Newbie (0 - 499 posts)

    Join Date
    Feb 2009
    Posts
    191
    Rep Power
    50
    How many steps the second loop (Newton-Raphson iteration) needs?
    Newton-Raphson for square root is a wonderful algorithm, ranking in my mind with Euclid's. It is both very simple, and extremely efficient. When my crude technique is used for the initial "guess" value, the maximum number of iterations required to reach the nearest integer square root of n is well approximated by log2(log2(n)). [log2 means "log base 2", or roughly the number of bits.]

    Because the program I posted uses one more iteration to detect convergence, a good estimate is 1 + log2(log2(n)). Note well that in some cases, the number of iterations will be much lower than this worst-case figure.

    Numbers up to 10,000 bits require no more than about 13 iterations; up to 100,000 bits, no more than about 16 iterations.

    As I hinted in the previous post, the Python script I gave doesn't always terminate: in integer math, the iteration can stably alternate between two values, in which case the test for same value as before will never evaluate true. Here's a better version:

    Code:
    def isqrt(n):
        # restrict to safe arguments
        if n < 1: return 0
    
        # construct an initial guess value
        d = len(str(n))	# count of decimal digits
        x = '1'
        d = d / 2
        while d > 0:
            d = d - 1
            x = x + '0'
        x = int(x)
    
        # perform Newton-Raphson iteration
        prev = 0
        prv2 = 0
        while x != prev and x != prv2:
            prv2 = prev
            prev = x
            q = n / x
            x = (x + q) / 2
    
        return x
    It can be improved by computing a more refined initial value for Newton-Raphson, which could reduce the worst case by probably one or two iterations.
  12. #7
  13. No Profile Picture
    Contributing User
    Devshed Newbie (0 - 499 posts)

    Join Date
    Jul 2011
    Posts
    313
    Rep Power
    0
    Thank you very much Mr. Mah$us !

    I played with this algorithm a little.

    Originally Posted by mah$us
    ... is well approximated by log2(log2(n)) ...
    Where i can find the proofs of the above?

    Here is my C version (without good initial guess).

    PHP Code:
    int isqrt int n ){ int p,r;
       
       
    n// initial guess

       
    do{
           
          
    r;
          
    = ( n/)/2// integer division!
          
       
    }while( );
       
       return 
    p;

  14. #8
  15. No Profile Picture
    Contributing User
    Devshed Newbie (0 - 499 posts)

    Join Date
    Feb 2009
    Posts
    191
    Rep Power
    50
    1. Though the original post specified 64 bits, and integers this size are easily handled in modern C implementations, computation related to cryptography sometimes requires taking square roots of much bigger numbers (more than 1000 bits).

    Anyone who wants to write a square-root program in C for large integers can use a large-integer library (there are several available). Python natively handles large integers, and so is very convenient for this kind of computation.

    2. leszek, the loop termination in the C program is valid, if and only if the successive approximations decrease monotonically. I suggest that you prove (or disprove) this monotonicity.

    3. The proof for the number of iterations is not difficult, and rather fun, so I'll leave it as an exercise. Hint: if y is the true square root, define the most recent approximation as y (1 + epsilon), where epsilon = 2^-n. Speaking a little loosely, we can say that this approximation is "good to n bits" (that is, the first n bits are correct). Using the iteration formula to compute the next approximation -- that approximation is good to how many bits? Once you know this, it is not difficult to compute the number of iterations required.
  16. #9
  17. No Profile Picture
    Contributing User
    Devshed Newbie (0 - 499 posts)

    Join Date
    Jul 2011
    Posts
    313
    Rep Power
    0

    second version


    PHP Code:
    // second version

    int isqrt int n ){ int p,r;

       
    1; while ( <= n/2*p;
       
       
    p// initial guess: isqrt(n) < r
       
       
    do{
      
          
    r;

          
    = ( n/)/2;
          
       }while( 
    r<);
       
       return 
    p;

    It works very fast! Now i believe in log2(log2(n)) convergence!

    I am now working on proof of correctness.
  18. #10
  19. No Profile Picture
    Contributing User
    Devshed Newbie (0 - 499 posts)

    Join Date
    Jul 2011
    Posts
    313
    Rep Power
    0
    Let k be the integer root of n (i.e., k is the number we are looking for).
    Let r' = ( r + n/r )/2 (integer divisions and r>0).

    Lemma. r'<r if and only if r>k.
    Proof. Easy.

    So the sequence of iterations is strictly decreasing, provided it stays ABOVE k.

    This proves the termination of algorithm: if we start with initial value > k,
    then there must be a step such that r' <= k < r .
    The next iteration value (say r'') is greater than or equal to r' (from the "only if" part),
    so the loop terminates.

    Claim. r' = k.
    Proof. A little bit more involved. As main tool use the well known fact
    that geometric mean <= arithmetic mean.

    Mr. Mah$us, thank you very much for bringing to my attention this beautiful algorithm !
    Last edited by leszek31417; August 17th, 2011 at 02:32 AM. Reason: english language ...

IMN logo majestic logo threadwatch logo seochat tools logo