Finding square roots using long division.

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