Square root by hand
Create a program that will calculate n digits of the square root of a number.
The program should continue forever (or until the number of digits is specified) calculating and outputting each decimal digit in succession. The program should be a "spigot algorithm" generating the digits of the number sequentially from left to right providing increasing precision as the algorithm proceeds.
C#
<lang csharp>using System; using static System.Math; using static System.Console; using BI = System.Numerics.BigInteger;
class Program {
static void Main(string[] args) {
BI i, j, k, d; i = 2; int n = -1; int n0 = -1;
j = (BI)Floor(Sqrt((double)i)); k = j; d = j;
DateTime st = DateTime.Now;
if (args.Length > 0) int.TryParse(args[0], out n);
if (n > 0) n0 = n; else n = 1;
do {
Write(d); i = (i - k * d) * 100; k = 20 * j;
for (d = 1; d <= 10; d++)
if ((k + d) * d > i) { d -= 1; break; }
j = j * 10 + d; k += d; if (n0 > 0) n--;
} while (n > 0);
if (n0 > 0) WriteLine("\nTime taken for {0} digits: {1}", n0, DateTime.Now - st); }
}</lang>
- Output:
14142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276415727350138462309122970249248360558507372126441214970999358314132226659275055927557999505011527820605714701095599716059702745345968620147285174186408891986095523292304843087143214508397626036279952514079896872533965463318088296406206152583523950547457502877599617298355752203375318570113543746034084988471603868999706990048150305440277903164542478230684929369186215805784631115966687130130156185689872372Time taken for 500 digits: 00:00:00.0092331
Go
<lang go>package main
import (
"fmt" "math/big"
)
var one = big.NewInt(1) var ten = big.NewInt(10) var twenty = big.NewInt(20) var hundred = big.NewInt(100)
func sqrt(n int64, limit int) {
i := big.NewInt(n)
j := new(big.Int).Sqrt(i)
k := new(big.Int).Set(j)
d := new(big.Int).Set(j)
t := new(big.Int)
digits := 0
for digits < limit {
fmt.Print(d)
t.Mul(k, d)
i.Sub(i, t)
i.Mul(i, hundred)
k.Mul(j, twenty)
d.Set(one)
for d.Cmp(ten) <= 0 {
t.Add(k, d)
t.Mul(t, d)
if t.Cmp(i) > 0 {
d.Sub(d, one)
break
}
d.Add(d, one)
}
j.Mul(j, ten)
j.Add(j, d)
k.Add(k, d)
digits = digits + 1
}
fmt.Println()
}
func main() {
sqrt(2, 480) // enough for demo purposes
}</lang>
- Output:
141421356237309504880168872420969807856967187537694807317667973799073247846210703885038753432764157273501384623091229702492483605585073721264412149709993583141322266592750559275579995050115278206057147010955997160597027453459686201472851741864088919860955232923048430871432145083976260362799525140798968725339654633180882964062061525835239505474575028775996172983557522033753185701135437460340849884716038689997069900481503054402779031645424782306849293691862158057846311159666871
Julia
Uses channels to iterate the spigot flow. <lang julia>function sqrt_spigot(number, limit=10000, bufsize=32)
spigot = Channel{Int8}(bufsize)
evaldigits(arr) = evalpoly(BigInt(10), arr)
""" Mark off pairs of digits, starting from the decimal point, working left. """
function markoff(n)
d = digits(n)
pairs, len = Vector{BigInt}[], length(d)
if isodd(len)
push!(pairs, [pop!(d)])
len -= 1
end
for i in len-1:-2:1
push!(pairs, [d[i], d[i+1]])
end
return pairs
end
""" look at first digit(s) and find largest i such that i^2 < that number """
function firststep(pairs)
i = evaldigits(popfirst!(pairs))
firstdigit = BigInt(findlast(x -> x * x <= i, 0:9) - 1)
put!(spigot, Int8(firstdigit))
return pairs, [firstdigit], i - firstdigit * firstdigit
end
"""
What is the largest number that we can put in the units and also multiply times
the divisor such that the result is still be less than or equal to what we have?
"""
function nextstep(pairs, founddigits, remain)
divisor = evaldigits(founddigits) * 2
remwithnext = remain * 100 + evaldigits(popfirst!(pairs))
d = BigInt(findlast(x -> x * (divisor * 10 + x) <= remwithnext, 0:9) - 1)
remain = remwithnext - (divisor * 10 + d) * d
pushfirst!(founddigits, d)
put!(spigot, Int8(d))
return pairs, founddigits, remain
end
""" start the process of adding digits to the channel """
function longhand_sqrt(n)
p = markoff(n)
pairs, founddigits, remain = firststep(p)
for _ in 1:limit
if remain == 0 # rational root, all zeros hereafter
put!(spigot, 0)
else
if isempty(pairs) # more zeros for part right of decimal point
push!(pairs, [0, 0], [0, 0], [0, 0], [0, 0])
end
pairs, founddigits, remain = nextstep(pairs, founddigits, remain)
end
end
end
@async(longhand_sqrt(number))
# return the channel from which to take! digits. return spigot
end
const spigot = sqrt_spigot(2)
for i in 1:500
print(take!(spigot)) i % 50 == 0 && println()
end
</lang>
- Output:
14142135623730950488016887242096980785696718753769 48073176679737990732478462107038850387534327641572 73501384623091229702492483605585073721264412149709 99358314132226659275055927557999505011527820605714 70109559971605970274534596862014728517418640889198 60955232923048430871432145083976260362799525140798 96872533965463318088296406206152583523950547457502 87759961729835575220337531857011354374603408498847 16038689997069900481503054402779031645424782306849 29369186215805784631115966687130130156185689872372
Phix
...whereas this is a spigot algorithm! <lang Phix>-- based on https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Decimal_(base_10) function bcd_sub(string a,b)
-- returns "a"-"b", coping with different lengths
-- (assumes a>=b, which it always will be here,
--- protected as it is by the bcd_le(b,a) call.)
integer c = 0, d = length(a)-length(b)
if d<0 then a = repeat('0',-d)&a
elsif d>0 then b = repeat('0', d)&b end if
for i=length(a) to 1 by -1 do
d = a[i]-b[i]-c
c = d<0
a[i] = d+c*10+'0'
end for
a = trim_head(a,"0") -- (note: "" equ "0")
return a
end function
function bcd_xp20x(string p, integer x)
-- returns x*(p*20+x)
integer c = 0, d, m = 1
p &= x+'0'
for i=length(p) to 1 by -1 do
d = (p[i]-'0')*m*x+c
p[i] = remainder(d,10)+'0'
c = floor(d/10)
m = 2
end for
if c then
p = (remainder(c,10)+'0')&p
c = floor(c/10)
if c then ?9/0 end if -- loop rqd?
end if
return p
end function
function bcd_le(string a,b)
-- returns a<=b numerically, taking care of different lengths
integer d = length(a)-length(b)
if d<0 then
a = repeat('0',-d)&a
elsif d>0 then
b = repeat('0',d)&b
end if
return a<=b
end function
function spigot_sqrt(string s, integer maxlen=50)
-- returns the square root of a positive string number to any precision
if find('-',s) or s="" then ?9/0 end if
integer dot = find('.',s)
if dot=0 then dot = length(s)+1 end if
if remainder(dot,2)=0 then s = "0"&s end if
dot += 1
string res = "", p = "", c = ""
integer i = 1
while true do -- (until (i>length && carry=0) or > maxlen)
if (i<=length(s) and s[i]='.')
or (i >length(s) and dot) then
res &= "."
dot = 0
i += 1
end if
c &= iff(i<=length(s)?s[i]:'0') &
iff(i<length(s)?s[i+1]:'0')
for x=9 to 0 by -1 do
string y = bcd_xp20x(p,x)
if bcd_le(y,c) then
c = bcd_sub(c,y)
res &= x+'0'
p &= x+'0'
exit
end if
if x=0 then ?9/0 end if -- (sanity check)
end for
i += 2
if (c="" and i>length(s)) or length(res)>maxlen then exit end if
end while
return res
end function ?spigot_sqrt("152.2756") ?spigot_sqrt("15241.383936") string r = spigot_sqrt("2",500) puts(1,join_by(r,1,100,""))</lang>
- Output:
(the final "2" was re-joined up by hand)
"12.34" "123.456" 1.41421356237309504880168872420969807856967187537694807317667973799073247846210703885038753432764157 2735013846230912297024924836055850737212644121497099935831413222665927505592755799950501152782060571 4701095599716059702745345968620147285174186408891986095523292304843087143214508397626036279952514079 8968725339654633180882964062061525835239505474575028775996172983557522033753185701135437460340849884 71603868999706990048150305440277903164542478230684929369186215805784631115966687130130156185689872372
REXX
This REXX version handles non-negative numbers less than unity, and may suppress superfluous trailing zeros.
It also handles the placing of a decimal point (if needed). <lang rexx>/*REXX program computes the square root by the old "by pen─n'─paper" (hand) method.*/ signal on halt /*handle the case of user interrupt. */ parse arg xx digs . /*obtain optional arguments from the CL*/ if xx== | xx=="," then xx= 2 /*Not specified? Then use the default.*/ if digs== | digs=="," then digs= 400 /* " " " " " " */ numeric digits digs + digs % 2 /*ensure enough decimal digits for calc*/ call sqrtHand xx, digs /*invoke the function for sqrt by hand.*/ halt: say /*pgm comes here for exact sqrt or HALT*/ exit 0 /*stick a fork in it, we're all done. */ /*──────────────────────────────────────────────────────────────────────────────────────*/ iSqrt: procedure; parse arg z; q= 1; r= 0; do while q<=z; q= q*4; end
do while q>1; q= q%4; _= z-r-q; r= r%2; if _>=0 then do; z= _; r= r+q; end; end
return r /*R is the integer square root of Z. */
/*──────────────────────────────────────────────────────────────────────────────────────*/ sqrtHand: parse arg x,places; parse value iSqrt(x) with j 1 k 1 ? /*j, k, ? ≡ iSqrt(x) */
if ?==0 then ?= /*handle the case of sqrt < 1. */
if j*j=x then do; say j; return; end /*have we found the exact sqrt?*/
L= length(?) /*L: used to place dec. point.*/
if L==0 then call charout , . /*handle dec. point for i < 1. */
do spit=1 for places
call charout , ?
if L>0 then do; call charout , .; L= 0; end /*process dec. point*/
if ?== then ?= 0 /*ensure the ? is a valid dig.*/
x= (x - k*?) * 100; ?= 1
k= j * 20
do while ?<=10
if (k + ?)*? > x then do; ?= ? - 1; leave; end
else ?= ? + 1
end /*while ?≤10*/
j= ? + j*10
k= ? + k
end /*spit*/
return</lang>
- output when using the default inputs:
1.41421356237309504880168872420969807856967187537694807317667973799073247846210703885038753432764157 2735013846230912297024924836055850737212644121497099935831413222665927505592755799950501152782060571 4701095599716059702745345968620147285174186408891986095523292304843087143214508397626036279952514079 8968725339654633180882964062061525835239505474575028775996172983557522033753185701135437460340849884 7160386899970699004815030544027790316454247823068492936918621580578463111596668713013015618568987237 2
- output when using the inputs of: .9 80
.9486832980505137995996680633298155601158665417975650480572514558377783315917714
- output when using the inputs of: 625
25
Visual Basic .NET
This is "spigot like", but not a true spigot, just an implementation of the "by hand" method of computing the square root, in this case, of two.<lang vbnet>Imports System.Math, System.Console, BI = System.Numerics.BigInteger
Module Module1
Sub Main(ByVal args As String())
Dim i, j, k, d As BI : i = 2
j = CType(Floor(Sqrt(CDbl(i))), BI) : k = j : d = j
Dim n As Integer = -1, n0 As Integer = -1,
st As DateTime = DateTime.Now
If args.Length > 0 Then Integer.TryParse(args(0), n)
If n > 0 Then n0 = n Else n = 1
Do
Write(d) : i = (i - k * d) * 100 : k = 20 * j
For d = 1 To 10
If (k + d) * d > i Then d -= 1 : Exit For
Next
j = j * 10 + d : k += d : If n0 > 0 Then n = n - 1
Loop While n > 0
If n0 > 0 Then WriteLine (VbLf & "Time taken for {0} digits: {1}", n0, DateTime.Now - st)
End Sub
End Module</lang>
- Output:
Execute without any command line parameters for it to run until it crashes (due to BigInteger variables eating up available memory). Output with command line parameter of 500:
14142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276415727350138462309122970249248360558507372126441214970999358314132226659275055927557999505011527820605714701095599716059702745345968620147285174186408891986095523292304843087143214508397626036279952514079896872533965463318088296406206152583523950547457502877599617298355752203375318570113543746034084988471603868999706990048150305440277903164542478230684929369186215805784631115966687130130156185689872372 Time taken for 500 digits: 00:00:00.0263710
Wren
<lang ecmascript>import "/big" for BigInt
var sqrt = Fn.new { |n, limit|
var i = BigInt.new(n)
var j = i.isqrt
var k = j
var d = j
var digits = 0
while (digits < limit) {
System.write(d)
i = (i - k*d) * 100
k = j * 20
d = BigInt.one
while (d <= 10) {
if ((k + d)*d > i) {
d = d.dec
break
}
d = d.inc
}
j = j*10 + d
k = k + d
digits = digits + 1
}
System.print()
}
sqrt.call(2, 480) // enough for demo purposes</lang>
- Output:
141421356237309504880168872420969807856967187537694807317667973799073247846210703885038753432764157273501384623091229702492483605585073721264412149709993583141322266592750559275579995050115278206057147010955997160597027453459686201472851741864088919860955232923048430871432145083976260362799525140798968725339654633180882964062061525835239505474575028775996172983557522033753185701135437460340849884716038689997069900481503054402779031645424782306849293691862158057846311159666871