Unsere Rechner können sehr viele Rechenoperationen schon auf der CPU erledigen, wenn es darum geht, mit den Standardtypen (z.b. double, dem gängigen 8-Byte-Fließkommaformat) zu rechnen. Das war gegen Anfang der 80er Jahre noch nicht so, da konnten typische CPUs nur gerade mit 8-Bit-Ganzzahlen addieren, subtrahieren und ein paar Bit-Operationen ausführen und doch ließ sich daraus alles aufbauen. Langzahladdition und -subtraktion waren trivial, Multiplikation und Division etwas schwieriger, aber noch gut machbar. Aber für die Quadratwurzeln musste man schon etwas genauer überlegen. Das mitgelieferte Basic konnte das für 5-Byte-Fließkommazahlen, aber das Verfahren war einfach nur dumm und langsam, den es wurde der Logarithmus genommen, halbiert und dann die Exponentialfuntion. Mathematisch völlig richtig, aber eben langsam und mit unnötigen Rundungsfehlern behaftet. Wer sich ein bisschen mit dem Thema auskennt, landet recht schnell bei dem Newton-Verfahren. Man schätzt die Quadratwurzel, z.B. von 2, dividiert 2 durch den Schätzwert und nimmt dann den Mittelwert von Schätzung und Quotient für die nächste Iteration. Zu offensichtlich, um es nicht so zu machen und in der Praxis auch ganz brauchbar. Es gibt aber ein Verfahren, wie man früher ähnlich wie die Division auch die Quadratwurzeln „handschriftlich“ berechnen konnte und dieses Verfahren lässt sich sehr gut verwenden, um eine Quadratwurzelberechnung zu schreiben. Hier ist eine Ruby-Implementierung, die von Ganzzahlen ausgeht:
Zu einer ganzen Zahl 
werden ganze Zahlen 
und 
gesucht so dass 
und 
gilt.
def check_is_int(x, name="x")
raise TypeError, "#{name}=#{x.inspect} must be Integer" unless x.kind_of? Integer
end
#
# helper method for internal use: checks if word_len is a reasonable
# size for splitting a number into parts
#
def check_word_len(word_len, name="word_len")
raise TypeError, "#{name} must be a positive number <= 1024" unless (word_len.kind_of? Fixnum) && word_len > 0 && word_len <= 1024
word_len
end
#
# split number (Integer) x into parts of word_len bits each such
# that the concatenation of these parts as bit patterns is x
# (the opposite of merge_from_words)
#
def split_to_words(x, word_len = 32)
check_word_len(word_len)
check_is_int(x, "x")
m = x.abs
s = (x <=> 0)
bit_pattern = (1 << word_len) - 1
words = []
while (m != 0 || words.length == 0) do
w = m & bit_pattern
m = m >> word_len
words.unshift(w)
end
if (s < 0) then
words[0] = -words[0]
end
words
end
#
# calculate the an integer s >= 0 and a remainder r >= 0 such that
# x = s**2 + r and s**2 <= x < (s+1)**2
# the wordwise algorithm is used, which works well for relatively
# large values of x. n defines the word size to be used for the
# algorithm. It is good to use half of the machine word, but the
# algorithm would also work for other values.
#
def sqrtw_with_remainder(x, n = 16)
check_is_int(x, "x")
check_is_int(n, "n")
n2 = n<<1
n1 = n+1
check_word_len(n2, "2*n")
s = (x <=> 0)
if (s == 0) then
return [0, 0]
elsif (s < 0)
a = sqrtw_with_remainder(-x)
return [ Complex(0, a[0]), a[1]]
end
xwords = split_to_words(x, n2)
if (xwords.length == 1) then
return sqrtb_with_remainder(xwords[0])
end
xi = (xwords[0] << n2) + xwords[1]
a = sqrtb_with_remainder(xi)
yi = a[0]
if (xwords.length <= 2) then
return a
end
xi -= yi*yi
2.upto(xwords.length-1) do |i|
xi = (xi << n2) + xwords[i]
d0 = (yi << n1)
q = (xi / d0).to_i
q0 = q
j = 0
was_negative = false
while (true) do
d = d0 + q
r = xi - (q * d)
break if (0 <= r && (r < d || was_negative))
if (r < 0) then
was_negative = true
q = q-1
else
q = q+1
end
j += 1
if (j > 10) then
break
end
end
xi = r
yi = (yi << n) + q
end
return [ yi, xi ]
end
#
# calculate an integer s >= 0 and a remainder r >= 0 such that
# x = s**2 + r and s**2 <= x < (s+1)**2
# the bitwise algorithm is used, which works well for relatively
# small values of x.
#
def sqrtb_with_remainder(x)
check_is_int(x, "x")
s = (x <=> 0)
if (s == 0) then
return [0, 0]
elsif (s < 0)
a = sqrtb_with_remainder(-x)
return [ Complex(0, a[0]), a[1]]
end
xwords = split_to_words(x, 2)
xi = xwords[0] - 1
yi = 1
1.upto(xwords.length-1) do |i|
xi = (xi << 2) + xwords[i]
d0 = (yi << 2) + 1
r = xi - d0
b = 0
if (r >= 0) then
b = 1
xi = r
end
yi = (yi << 1) + b
end
return [yi, xi]
end
Daraus lässt sich dann relativ einfach eine Quadratwurzelfunktion mit entsprechender Nachkommstellenberechnung bilden.
Witzerweise lässt sich das Verfahren auch für Kubikwurzeln sinngemäß anwenden.
Hier ein kleiner Test:
irb(main):138:0> sqrtw_with_remainder(2) => [1, 1] irb(main):139:0> sqrtw_with_remainder(200) => [14, 4] irb(main):140:0> sqrtw_with_remainder(20000) => [141, 119] irb(main):141:0> sqrtw_with_remainder(2000000) => [1414, 604] irb(main):142:0> sqrtw_with_remainder(200000000000000000000000000000000000000000000000000000000000) => [447213595499957939281834733746, 228299936041363866321288807484] irb(main):143:0> sqrtw_with_remainder(2000000000000000000000000000000000000000000000000000000000000) => [1414213562373095048801688724209, 1974464361663955412145937324319] irb(main):144:0> sqrtw_with_remainder(3000000000000000000000000000000000000000000000000000000000000) => [1732050807568877293527446341505, 3021967735564464902990914334975] irb(main):145:0> sqrtw_with_remainder(3) => [1, 2] irb(main):146:0> sqrtw_with_remainder(300) => [17, 11] irb(main):147:0> sqrtw_with_remainder(30000) => [173, 71] irb(main):148:0> sqrtw_with_remainder(3000000) => [1732, 176]
