Quadrat- und Kubikwurzeln berechnen vor 30 Jahren und heute ;-)

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]