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.
'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
'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 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 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 PM This message has been edited by MCalkins on May 6, 2012 12:00 PM This message has been edited by MCalkins on May 6, 2012 11:52 AM
I recall the method you used as being taught in school as a manual way of finding square roots. This was before the era of hand-held calculators. I doubt that it's still being taught in school today. Anyway, you did a remarkable job of converting the method to computer code. I compared your results to those of mine and found no discrepancies other than the exact number of digits displayed, but that's of no consequence.
I used the continued fraction method which was also covered in the Wikipedia article you cited and also in this one:
http://mathworld.wolfram.com/PellEquation.html
The math was a little over my head, but I was able to adapt the general description to code. The independent functions, add2$ and subtract2$ are not part of the solution. They are used to add and subtract the big integers generated by the algorithm.
DIM p(1000) AS SINGLE
DIM q(1000) AS SINGLE
DIM a(1000) AS SINGLE
DIM h(1000) AS STRING
DIM k(1000) AS STRING
DO
n = n + 1
p(n) = a(n - 1) * q(n - 1) - p(n - 1)
q(n) = (tnum - p(n) ^ 2) / q(n - 1)
a(n) = INT((a(0) + p(n)) / q(n))
FOR i = 1 TO a(n)
h(n) = add2(h(n), h(n - 1))
k(n) = add2(k(n), k(n - 1))
NEXT i
h(n) = add2(h(n), h(n - 2))
k(n) = add2(k(n), k(n - 2))
LOOP UNTIL LEN(k(n)) > 51
temp$ = k(n)
DO
m = 0
DO UNTIL (LEN(h(n)) < LEN(temp$)) OR (LEN(h(n)) = LEN(temp$) AND h(n) < temp$)
h(n) = subtract2(h(n), k(n))
m = m + 1
LOOP
ans$ = ans$ + LTRIM$(STR$(m))
h(n) = h(n) + "0"
LOOP UNTIL LEN(ans$) = 100
dp = LEN(STR$(a(0))) - 1
PRINT LEFT$(ans$, dp) + "." + MID$(ans$, dp + 1)
END
FUNCTION add2$ (a1$, a2$)
alen = LEN(a1$) - LEN(a2$)
IF alen > 0 THEN
a2$ = STRING$(alen, "0") + a2$
ELSE
a1$ = STRING$(ABS(alen), "0") + a1$
END IF
FOR b = LEN(a1$) TO 1 STEP -1
j = VAL(MID$(a1$, b, 1)) + VAL(MID$(a2$, b, 1)) + cry
cry = j \ 10
j$ = LTRIM$(STR$(j MOD 10))
sum$ = j$ + sum$
NEXT b
IF cry THEN sum$ = LTRIM$(STR$(cry)) + sum$
add2$ = sum$
END FUNCTION
FUNCTION subtract2$ (a1$, a2$)
alen = LEN(a1$) - LEN(a2$)
IF alen >= 0 THEN
a2$ = STRING$(alen, "0") + a2$
ELSE
a1$ = STRING$(ABS(alen), "0") + a1$
END IF
FOR b = LEN(a1$) TO 1 STEP -1
d1 = VAL(MID$(a1$, b, 1))
d2 = VAL(MID$(a2$, b, 1))
d = d1 - d2
IF d > 0 THEN j = d1 - d2 - bor: bor = 0
IF d < 0 THEN j = d1 + 10 - d2 - bor: bor = 1
IF d = 0 AND bor = 1 THEN j = d1 + 10 - d2 - bor: bor = 1
IF d = 0 AND bor = 0 THEN j = d1 - d2
j$ = LTRIM$(STR$(j))
dif$ = j$ + dif$
NEXT b
DO
j$ = LEFT$(dif$, 1)
IF j$ = "0" AND LEN(dif$) > 1 THEN dif$ = MID$(dif$, 2)
LOOP UNTIL j$ <> "0" OR LEN(dif$) = 1
subtract2$ = dif$
END FUNCTION