# Sqrt Challenge

April 28 2012 at 4:07 PM

Write a program that accurately calculates the square root of any positive integer less than 1 million to 100 significant digits.

 Respond to this message

Moderator

# This seems to work.

May 5 2012, 4:06 PM

It's roughly based on:

http://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Decimal_.28base_10.29

I do it in a decimal base so that I don't have to convert a huge binary number to decimal afterwards. I suspect that this could be optimized quite a bit more in assembly, perhaps using BCD functionality.

I tried to optimize the way y is calculated. In doing so, I've confused myself a bit, so that I have a bit of a struggle following it.

P.S.

It doesn't round, so I have a target of 101 digits so that you know which way to round the 100th digit.

--------------------------------------------------------

For clarity, here is a simplified version that works for c <= 99, but overflows quickly.

--------------------------------------------------------

'public domain
DEFLNG A-Z
DIM y(0 TO 1)
c = 2 ' must be <= 99
'won't show decimal point
'will overflow for non-termating results.

CLS
DO
x = 1
w = 1
o = 0
y(o) = 0
DO
z = z + 1
y(w) = y(o) + z
IF y(w) > c THEN
z = z - 1
x = x - 1
EXIT DO
END IF
z = z + 1
IF x = 9 THEN EXIT DO
x = x + 1
o = w
w = w XOR 1
LOOP
w = x AND 1
PRINT HEX\$(x);
c = c - y(w)
c = c * 100
IF c = 0 THEN EXIT DO
z = z * 10
LOOP
PRINT
END

--------------------------------------------------------

And hopefully this will somewhat clarify what I'm doing with z:

--------------------------------------------------------

DEFLNG A-Z

p = 3 'arbitrary test number

CLS
z = 20 * p
FOR x = 1 TO 9
z = z + 1
y = y + z
PRINT "x:"; x, "y:"; y, "x*(20*p+x):"; x * (20 * p + x)
z = z + 1
NEXT
END

--------------------------------------------------------

Here is the program:

--------------------------------------------------------

'public domain, 2012 may, michael calkins
' http://www.network54.com/Forum/202193/message/1335654428/Sqrt+Challenge
' basically using:
' http://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Decimal_.28base_10.29
' but with attempts at optimization.

CONST cubound = &H1FF
CONST zubound = &H7F
CONST yubound = &H7F

CONST pnt = cubound - 5
'pnt should be even. subtract enough from cubound to make room for however
'many digits init will have left of the decimal place.

CONST target = 101
'if you inrease the target too much, you might need to increase the array
'ubounds.

CONST yl = 0
'these will be arrays of decimal digits stored in little endian
'make these next 3 arrays _BYTEs if using QB64.
DIM c(0 TO cubound) AS INTEGER
DIM z(0 TO zubound) AS INTEGER
DIM y(0 TO yubound, 0 TO 1) AS INTEGER

'the l and u values are the indecies of the low and high digits
DIM cl AS INTEGER
DIM cu AS INTEGER
DIM zl AS INTEGER
DIM zu AS INTEGER
'yl is always 0
DIM yu(0 TO 1) AS INTEGER
DIM init AS LONG
DIM x AS INTEGER
DIM w AS INTEGER        'selects which y()
DIM o AS INTEGER        'w xor 1 (where needed)
DIM i AS INTEGER        'temporary variables
DIM i1 AS INTEGER
DIM ac AS INTEGER       'temp variable (accumulator)
DIM t AS INTEGER        'return value for pseudo-function

CLS
DO
PRINT
DO
INPUT "Enter an integer in the range of 1 to 999999"; init
LOOP WHILE init < 0 OR init > 999999
IF 0 = init THEN EXIT DO
PRINT LTRIM\$(STR\$(init ^ .5#))

'initialize arrays and variables:

FOR i = 0 TO pnt - 1
c(i) = 0
NEXT
FOR i = 0 TO zubound
z(i) = 0
NEXT

'this could be modified to support non-integer inits.

i = pnt - 1
DO
i = i + 1
c(i) = init MOD 10
init = init \ 10
LOOP WHILE init
cu = i
cl = (cu OR 1) - 1
zl = UBOUND(z) - 1
zu = UBOUND(z) - 1

'correlation between variables and the wikipedia article:

' x is x

' c is c

' y(w) is y at the point where it is compared with c
' there are two y() arrays so that y(o) acts as a cache for the previous
' value of y(w).

' z is an optimization related to 20 * p.
' at the start of each outer loop, z is 20 * p.
' inside the inner loop, it is manipulated to assist with calculating y(w).
' after the end of the inner loop, z is 20 * p + 2 * x.

DO
x = 1
w = 1
o = 0
yu(0) = yl: y(yl, 0) = 0       ' y(o) = 0
DO
GOSUB incz            ' z = z + 1
GOSUB ywaddyoz        ' y(w) = y(o) + z
GOSUB gtywc
IF t THEN             ' if y(w) > c then
GOSUB decz           ' z = z - 1
x = x - 1
EXIT DO
END IF
GOSUB incz            ' z = z + 1
IF x = 9 THEN EXIT DO
x = x + 1
o = w
w = w XOR 1
LOOP
w = x AND 1            ' o is not and does not need to be updated
PRINT CHR\$(x + &H30);
IF cl = pnt THEN PRINT ".";
GOSUB subcyw           ' c = c - y(w)
IF c(cu) = 0 THEN cu = cu - 1          'drop a leading zero (not necessary)
cl = cl - 2            ' c = c * 100 + next two digits
FOR i = cl TO cu
IF c(i) THEN EXIT FOR
NEXT
' if c = 0 or found target digits after decimal point then exit do
IF i = cu + 1 OR cl = pnt - 2 * target THEN EXIT DO
zl = zl - 1            ' z = z * 10
LOOP

LOOP
SYSTEM

incz:                   ' z = z + 1
IF z(zl) < 9 THEN
z(zl) = z(zl) + 1
ELSE
z(zl) = 0
IF zl = zu THEN
zu = zu + 1
z(zu) = 1
ELSE
z(zl + 1) = z(zl + 1) + 1
' the carry will not propagate more than 1 place.
' this is because the right 2 digits of the initial number will not be
' above 80, and we will not be incrementing more than 18 times.
END IF
END IF
RETURN

decz:                   ' z = z - 1
IF z(zl) THEN
z(zl) = z(zl) - 1
ELSE
z(zl) = 0
IF zu > zl + 1 OR z(zl + 1) > 1 THEN
z(zl + 1) = z(zl + 1) - 1
' the borrow will not propagate more than 1 place.
ELSE
zu = zl
END IF
END IF
RETURN

subzx:                  ' z = z - x
'this sub is never called. I thought it was necessary, but it doesn't seem
'to be.
IF z(zl) >= x THEN
z(zl) = z(zl) - x
ELSE
z(zl) = z(zl) + 10 - x
IF zu > zl + 1 OR z(zu) > 1 THEN
z(zu) = z(zu) - 1
' the borrow will not propagate more than 1 place.
ELSE
zu = zl
END IF
END IF
RETURN

gtywc:                  ' t = y(w) > c
IF yu(w) > cu - cl THEN
t = -1: RETURN
ELSE
i = cu - cl
i1 = cu
DO
IF i <= yu(w) THEN
IF y(i, w) > c(i1) THEN
t = -1: RETURN
ELSEIF y(i, w) < c(i1) THEN
t = 0: RETURN
END IF
ELSE
IF c(i1) THEN t = 0: RETURN
END IF
IF i = yl THEN t = 0: RETURN
i = i - 1
i1 = i1 - 1
LOOP
END IF
RETURN

ywaddyoz:               ' y(w) = y(o) + z
i = yl
i1 = zl
ac = 0
DO
IF i <= yu(o) THEN ac = ac + y(i, o)
IF i1 <= zu THEN ac = ac + z(i1)
IF ac < 10 THEN
y(i, w) = ac
ac = 0
ELSE
y(i, w) = ac - 10
ac = 1
END IF
i = i + 1
i1 = i1 + 1
LOOP WHILE i <= yu(o) OR i1 <= zu OR ac
yu(w) = i - 1
RETURN

subcyw:                 ' c = c - y(w)
i = cl
i1 = yl
ac = 0
DO
IF i <= cu THEN ac = ac + c(i)
IF i1 <= yu(w) THEN ac = ac - y(i1, w)
IF ac >= 0 THEN
c(i) = ac
ac = 0
ELSE
c(i) = ac + 10
ac = -1
END IF
i = i + 1
i1 = i1 + 1
LOOP WHILE i <= cu OR i1 <= yu(w) OR ac
cu = i - 1
RETURN

printz:
PRINT
PRINT "z: ";
FOR i = zu TO zl STEP -1
PRINT HEX\$(z(i));
NEXT
PRINT
RETURN

printc:
PRINT
PRINT "c: ";
FOR i = cu TO cl STEP -1
PRINT HEX\$(c(i));
NEXT
PRINT
RETURN

printyw:
PRINT
PRINT "y(w): ";
FOR i = yu(w) TO yl STEP -1
PRINT HEX\$(y(i, w));
NEXT
PRINT
RETURN

--------------------------------------------------------

 This message has been edited by MCalkins on May 5, 2012 4:18 PM

 Respond to this message

Moderator

# a bug fix.

May 6 2012, 11:47 AM

----- change:

'this could be modified to support non-integer inits.

----- to:

'this could be modified to support non-integer inits. anything related to the
'init variable would need to be modified, as well as the part that checks if c
'is zero.

----- and change:

FOR i = cl TO cu
IF c(i) THEN EXIT FOR
NEXT
' if c = 0 or found target digits after decimal point then exit do
IF i = cu + 1 OR cl = pnt - 2 * target THEN EXIT DO

----- to:

IF cl <= pnt - 2 THEN
'if we have reached the decimal point.
FOR i = cl TO cu
IF c(i) THEN EXIT FOR
NEXT
' if c = 0 or found target digits after decimal point then exit do
' this assumes an integer init. to support non-integer inits, this would
' need to be modified.
IF i = cu + 1 OR cl = pnt - 2 * target THEN EXIT DO
END IF

----------

That should fix missing zeros for inits like 10000 or 14400, and incorrect results for inits like 10001.

In the case of 10001, for example, the old code would have taken the number two digits at a time: "01", "00", "01". The "01" would have generated no remainder, and when the code would have looked at "00", it would have thought that it was done.

----------

To support non-integer "init"s, you might allow the number to be input as a string, then copy it's digits into the "c" array. "pnt" should probably become a variable instead of a constant. Also, the "FOR i = cl TO cu" part would have to be changed to search the entire actual number, not just the currently visible part of the number. Some new variable might be used in place of "cl" on that line. You could just change "cl" to "0" in the "FOR" line, but using a new variable would seem better. As it is, the code is okay, because the "IF cl <= pnt - 2 THEN" line means that we will have reached the decimal point, and therefore, we will be dealing with the entire integer "init".

Regards,
Michael

P.S. Actually, as the wikipedia article suggests, the decimal point has no bearing on the actual digits. So, knowing the the position of the decimal point in the above program has two roles currently: to know when to print the "." to the screen, and to know when you have reached at least to the end of the initial number. Those are actually two distinct roles. For example, if the code that puts the "init" value into the "c" array knows where the last nonzero pair of digits is, that is what is important. So for example, in the case of an init like 10000 or 14400, there is no reason the algorithm should need to go all the way to the decimal point, except for the purpose of displaying the proper trailing zeros. But, in any case, the program must not stop before it has reached the last non-zero pair of digits from the "init", regardless of where the decimal point is.

 This message has been edited by MCalkins on May 6, 2012 12:01 PMThis message has been edited by MCalkins on May 6, 2012 12:00 PMThis message has been edited by MCalkins on May 6, 2012 11:52 AM

 Respond to this message
lawgin