Finding square roots using long division.

[jahuja73]

Below is code for long division based finding square root.
But, it fails for decimals, as well as non-square integers.
I mean it gives the same result for 1225, and 1235, i.e. 35.

INFINITY_ = 9999999
 
def sqrtByLongDivision(n,l): 
    i = 0
    udigit, j = 0, 0 # Loop counters 
    cur_divisor = 0
    quotient_units_digit = 0
    cur_quotient = 0
    cur_dividend = 0
    cur_remainder = INFINITY_
    #print 'length of input:', len(n)
    a = [0]*l
 
    # Dividing the number into segments 
    while (n > 0): 
        a[i] = int( n % 100)
        n = int(n // 100)
        i += 1
 
    # Last index of the array of segments 
    i -= 1

    for j in range(i, -1, -1): 
        cur_remainder = INFINITY_
        cur_dividend = cur_dividend * 100 + a[j] 
 
        for udigit in range(10): 
            if (cur_dividend  - (cur_divisor * 10 + udigit) * udigit >= 0):
            #if ( (cur_remainder >(cur_dividend - ((cur_divisor * 10 + udigit) * udigit))) and (cur_dividend - ((cur_divisor * 10 + udigit) * udigit) >= 0)):
                # Calculating the remainder 
                cur_remainder = cur_dividend - ((cur_divisor * 10 + udigit)* udigit) 
                # Updating the units digit of the quotient 
                quotient_units_digit = udigit 
                if cur_remainder == 0:
                    print 'Exiting: cur_remainder ==0'
                    break
                print ('New rem. :', cur_remainder)
 
            else:
                print(' else break')
                break
 
        # Adding units digit to the quotient 
        cur_quotient = cur_quotient * 10 + quotient_units_digit 
 
        # New divisor is two times quotient 
        cur_divisor = cur_quotient * 2
 
        # Including the remainder in new dividend 
        cur_dividend = cur_remainder 
 
    return cur_quotient 
 
# Driver code 
x = 1235.789
# Find length of integer
y=x
length =0
while (y > 0): 
    y //=100
    length+= 1
print 'length :', length
 
print(sqrtByLongDivision(x, length))

Also, hope once can find the square root for non-square integers, can also find for decimals.

[Serafim]

The long division method, such as you have programmed it gives you the integer part of the square root, so it is correct that you get 35 both for 1225 and for 1235 as you stop counting when you have used the full length of the array ("a" in your code).

If you want to continue beyond this you must add new groups of zeroes to "a", one group for the number of decimals you want

Eg. the square root of 3:
With your method you get 1.
Next step, decide the number of decimals, say 3.
Extend "a" to [03, 00, 00, 00]
and count with the same method

        1.732
 ---------------
        3.000000
   1    1
         200
  27     189
          1100
 343      1029
            7900
3462        6924
             176

You can experiment with your own code and use 1225000000 => 35000, while 1235000000 => 35142. It is just a matter of where to place the decimal point...

[jahuja73]

If you want to continue beyond this you must add new groups of zeroes to "a", one group for the number of decimals you want

Eg. the square root of 3:
With your method you get 1.
Next step, decide the number of decimals, say 3.
Extend "a" to [03, 00, 00, 00]
and count with the same method

Thanks, seems easy for integers, but still clueless for floats as 1235.6789.
The coding is too difficult as have jobs:
1. In the driver code, need know length of float.
Need break in units of two based on that..

If make modifications in the Driver code as follows:

#Driver code 
x = 1235.789
# Find length of integer
y=x
length =0
"""
while (y > 0): 
    y //=100
    length+= 1
"""
while (y > 0): 
    if (y>=1):
        print 'integer part'
        y = y //100
        print 'y :', y
        length += 1
    elif (y <1) :
        print 'decimal part'
        y = y*100 
        y /= 100
        print 'y :', y
        length += 1
 

Then, get:

Output:

integer part
y : 12.0
integer part
y : 0.0
length : 2

Unable to find a way out.

[Serafim]

It is actually pretty easy. No need to break it up into an integer and a decimal part.

Say you want 10 decimals for the square root of 1235.789 (taken from your program).

Then multiply the value with 1020 = 123578900000000000000000

Running your function on it gives 351537906917. Divide by 1010 = 35.1537906917
Correct answer!

If you want many decimals, then don't calculate the answer as an integer, just keep the figures calculated in a list [3, 5, 1, 5, 3, 7, 9, 0, 6, 9, 1, 7]. The length of the list minus the number of decimals tells you where to place the decimal point.

[Serafim]

I have made an implementation that needs some explanation.

1. I convert the input to a string as I want to simply move the decimal point out of the way. Stepwise floating point multiplication gives raise to round-off errors that severely affects the result.
   
2. Two methods are not shown as they are tedious in their details but simple to implement.
   • "prepare" takes two parameters, "n" is the number to be "square rooted" an "d" is the number of desired decimals. The return value is a list with double figure integers (shortened if there is an intial zero) and the calculated number of decimals. E.g. prepare(123.456, 10) returns a=[1,23,45,60,0,0,0,0,0,0,0,0], d=10, while prepare(1.234, 0) returns a=[1,23,40] and d=2
     
   • "normalize" takes a string representation of a result and removes initial and trailing zeroes as well as the decimal point if there are no decimals
     
3. The implementation has a flaw due to the nature of the long division method that I didn't bother to correct. There might be a round-off error in the last digit of a result with decimals so all such results ought to be rounded with omission of the last digit. The long division method is devised for integer arithmetics and is accurate when dealing with integers but as the method collects one digit at the time there is no telling if the first omitted digit in a potentially infinite number of decimals would be less than 5 or not.
   

The implementation of the main method (note that I collect the figures in a string):

def sqrtLongDiv(n, d=0):
    a, d = prepare(n, d) # See the comment in the explanation
    b = ''               # Here we collect the result
    drange = range(10)[::-1]
    divisor = 0
    dividend = 0
    remainder = 0
    quotient = 0
    unit = 0
    dec_point = len(a)-d # the position of the decimal point
    for i in range(len(a)):
        if i == dec_point: # Insert the decimal point here
            b += '.'
        dividend = dividend * 100 + a[i] # Calculate the new dividend
        for dig in drange:  # Testing with 9,8,..,0 gives simple stop condition
            if (divisor * 10 + dig) * dig <= dividend:
                # collect the resulting digit
                b += str(dig)

                # calculate the reminder
                remainder = dividend - ((divisor * 10 + dig) * dig)

                # Remember the newly calculated digit
                unit = dig
                # and break out as we have found a correct digit
                break

        # Calculate the new intermediary quotient
        quotient = quotient * 10 + unit
        # as well as the new intermediary divisor
        divisor = quotient * 2
        # and let the calculated remainder become the base for the next dividend
        dividend = remainder
    return normalize(b)

Below there is a test run (I slightly changed the return value for the sake of the show).

Output:
Just integers and no decimals, sqrtLongDiv(n):
0      => 0
1      => 1
2      => 1
3      => 1
4      => 2
5      => 2
6      => 2
1225   => 35
1226   => 35
1227   => 35
1228   => 35
1234   => 35
1235   => 35

With up to 9 decimals, sqrtLongDiv(n, 9):
12.25  => 3.5
12.255 => 3.500714212
12.26  => 3.50142828
12.265 => 3.502142201
12.27  => 3.502855977
12.275 => 3.503569608
12.28  => 3.504283093
12.285 => 3.504996433

With up to 50 decimals, sqrtLongDiv(n, 50):
0      => 0
1      => 1
2      => 1.41421356237309504880168872420969807856967187537694
3      => 1.73205080756887729352744634150587236694280525381038
4      => 2
5      => 2.23606797749978969640917366873127623544061835961152
6      => 2.44948974278317809819728407470589139196594748065667
7      => 2.64575131106459059050161575363926042571025918308245
1225   => 35
1226   => 35.0142828000231928211878077093306242214076547909407
1227   => 35.02855977627398828194295855627241569536941311316346
1228   => 35.04283093587046371749192546958435590451070490205804
0.5    => 0.70710678118654752440084436210484903928483593768847
1.5    => 1.22474487139158904909864203735294569598297374032833
2.5    => 1.5811388300841896659994467722163592668597775696626

With 250 decimals, sqrtLongDiv(n, 250):
1235 => 35.1425667816111576639565564906336937149320439161531722366224396398535743816887293934179865891453178536207436977322785579237218729545867849156703998844975163669968675086712825842393426162278364026764080153029875549199784450044423207555532767716521039431

All results are verified by WolframAlpha (some are as predicted
not correctly rounded-off in the last decimal)

[jahuja73]

I have made an implementation that needs some explanation.

1. I convert the input to a string as I want to simply move the decimal point out of the way. Stepwise floating point multiplication gives raise to round-off errors that severely affects the result.
   
2. Two methods are not shown as they are tedious in their details but simple to implement.
   • "prepare" takes two parameters, "n" is the number to be "square rooted" an "d" is the number of desired decimals. The return value is a list with double figure integers (shortened if there is an intial zero) and the calculated number of decimals. E.g. prepare(123.456, 10) returns a=[1,23,45,60,0,0,0,0,0,0,0,0], d=10, while prepare(1.234, 0) returns a=[1,23,40] and d=2
     
   • "normalize" takes a string representation of a result and removes initial and trailing zeroes as well as the decimal point if there are no decimals
     
3. The implementation has a flaw due to the nature of the long division method that I didn't bother to correct. There might be a round-off error in the last digit of a result with decimals so all such results ought to be rounded with omission of the last digit. The long division method is devised for integer arithmetics and is accurate when dealing with integers but as the method collects one digit at the time there is no telling if the first omitted digit in a potentially infinite number of decimals would be less than 5 or not.
   

1. Though simple, request the two modules: 'prepare', 'normalize' for the sake of consistency.

2. It seems the parameter 'd' to prepare is user defined. The only way to custom (input based) definition seems to convert input to string, and check for decimal. Say, '.' in string (if correct). Then, based on number of digits after decimal, can decide value of 'd'. Cannot though think of a criteria.

3. The possible flaw referred by you should be due to limited number of decimal digits? If yes, then the accuracy is affected. If so, then rounding is to be done by a seperate module?

[Serafim]

Here goes:

def prepare(n, d):
    a = []
    m = str(n)
    int_part = ''
    dec_part = ''
    if '.' in m:
        int_part, dec_part = m.split('.')
    else:
        int_part = m
    if len(int_part) % 2 != 0:
        int_part = '0' + int_part
    for i in range(0, len(int_part), 2):
        a.append(int(int_part[i:i+2]))
    if len(dec_part) % 2 != 0:
        dec_part += '0'
    for i in range(0, len(dec_part), 2):
        a.append(int(dec_part[i:i+2]))
    min_dec = len(dec_part) // 2
    d = max (d, min_dec)
    for i in range(min_dec, d):
        a.append(0)
    return (a, d)

def normalize(s):
    if '.' in s:
        while s[-1] == '0':
            s = s[:-1]
        if s[-1] == '.':
            s = s[:-1]
    while len(s) > 1 and s[0] == '0' and s[1] != '.':
        s = s[1:]
    return s

[Serafim]

I didn't answer your last question (3):
I think normalize is the place. Add one digit to the decimals and then shorten the result. It's a messy business as the result is actually a string. Maybe consider the two last digits as a decimal number and round it of as an integer. Then replacing the two last digits in the result string by that integer.
E.g. the last two digits = '37' => 3.7 => rounded to 4 and replace '37' by '4'.

Maybe there are other ways.

[jahuja73]

I didn't answer your last question (3):
I think normalize is the place. Add one digit to the decimals and then shorten the result. It's a messy business as the result is actually a string. Maybe consider the two last digits as a decimal number and round it of as an integer. Then replacing the two last digits in the result string by that integer.
E.g. the last two digits = '37' => 3.7 => rounded to 4 and replace '37' by '4'.

Maybe there are other ways.

Very bold answer as never got this detail printed anywhere on SO. Seems they want reader to figure it out. But came only near to this idea.
But, still (after your answer) feel to comprehend fully need code or at least elaborate algorithm.
Thanks in anticipation.

[Serafim]

Sorry, I am off to my next project. Try implementing it yourself. It is not that big an effort. Use my hint in previous message and the round() function. That would do the trick.

By the way, I found the code that you have copied "your" implementation from. Seems like you have taken someone elses code and want us to explain why it doesn't work the way you want it to work. I think I explained why. My implementation is based on the fact that squareroot(100*x) = 10 * squareroot(x).