愚直な平方根の計算

ある講演会資料作成のため、コンピューター技術の歴史を振り返るときにチューリングマシンを復習する機会がありました。
一桁づつの計算でも手順を工夫することで複雑な計算が可能ということが示され、その後のコンピューター開発に影響を与えました。

素晴らしい発想力ですね。その想像力のほんのカケラだけでもあやかりたいものです。

さて私ができることといえば・・一桁づつということにヒントを得て、BCD配列で平方根を求めるプログラムを書きました。
配列同士で足し算・掛け算を行い、平方根を求めます。

C:\temp>python hsqrt-1.1.py 2 50
------------------------------------------------------------
10桁までの計算結果
1.414213562
------------------------------------------------------------
20桁までの計算結果
1.4142135623730950488
------------------------------------------------------------
20桁までの計算結果
1.4142135623730950488
------------------------------------------------------------
40桁までの計算結果
1.414213562373095048801688724209698078569
------------------------------------------------------------
50桁までの計算結果
1.4142135623730950488016887242096980785696718753769
------------------------------
---------------結果---------------
1.4142135623730950488016887242096980785696718753769
----------------------------------

計算開始時刻:2018-06-13 17:43:11.532379
計算終了時刻:2018-06-13 17:43:15.250831
計算時間: 0:00:03.718452

(Intel Core-i7920 @2.67GHz)

20桁目の表示が2回あるのは表示判定部分の高速化(手抜き)のためです・・


RaspberryPi(BCM2835)だと、同じ50桁で

計算開始時刻:2018-06-13 08:45:38.218279
計算終了時刻:2018-06-13 08:46:27.735490
計算時間: 0:00:49.517211

う~ん、愚直すぎて遅い・・・
しかし、このようにして無理数が計算されていくのは、なんとも不思議な感じがします。


<2018.06.16追記>
インデントが失われておりましたので再掲載します。失礼いたしました。


# -*- coding: utf-8 -*-

# hsqrt-1.1.py
#
# BCDで平方根を求めるプログラム
# Writtend by T'z Factory Co.,Ltd
# Ver 1.0  2018/05/25
# Ver 1.1  2018/06/04  Python2系に対応するため、copy属性をlist関数に書き換え。表示を一部sys.stdout関数に書き換え

#==========================================
# インポートモジュール
#==========================================
try:
  import datetime, logging, sys, time

except ImportError:
  print('必要なモジュールをインポートできませんでした。環境を確認してください。')
  sys.exit()


#==========================================
#変数定義
#==========================================

SEP =  '---------------'

#==========================================
# 上位桁から有効数字までの'0'を削除する
#==========================================
def zerosupress(a):
  i = len(a)
  j = 0
  while (i > 0 and j < 1):
    if a[i-1] == '0':
      a.pop()
    else:
      j += 1
    i -= 1
  i = len(a)
  if a[i-1] == '.':
    a.append('0')
  return(a)

#==========================================
# 配列同士の足し算
#==========================================
def plus(a1, a2):
  a = list(a1)
  b = list(a2)

# 小数点以下の桁数を合わせる
  p = a.index('.') - b.index('.')
  if p > 0:
    for i in range(p):
      b.insert(0, '0')
  if p < 0:
    for i in range(-p):
      a.insert(0, '0')

# 整数部分の桁数を合わせる
  p = (len(a) - a.index('.')) - (len(b) - b.index('.'))
  if p > 0:
    for i in range(p):
      b.append('0')
  if p < 0:
    for i in range(-p):
      a.append('0')
  l_a = len(a)

# 和格納用配列初期化
  c = []
  for i in range(l_a + 1):
    c.append('0')

# キャリー初期化
  cy = 0

# 計算
  for i in range(l_a):
    if a[i] == '.':
      c[i] = '.'
    else:
      t = int(a[i]) + int(b[i]) + cy
      if t > 9:
        cy = 1
        t = int(str(t)[1])
      else:
        cy = 0
      c[i] = str(t)

# 最後のキャリーを足し込む
  if cy == 1:
    c[i + 1] = str(cy)
  return(c)

#==========================================
# 配列同士の掛け算
#==========================================
def multi(a1, a2):
  a = list(a1)
  b = list(a2)

# ピリオドを一旦最下位に移動(整数化)
  p_a = a.index('.')
  p_b = b.index('.')
  a.pop(p_a)
  b.pop(p_b)
  a.insert(0, '.')
  b.insert(0, '.')
  l_a = len(a)
  l_b = len(b)
  l_d = l_a - l_b
  l_p = l_a + l_b

# 桁数を合わせる
  if l_d > 0:
    for i in range(l_d):
      b.append('0')
  if l_d < 0:
    for i in range(-l_d):
      a.append('0')

# 積格納用配列初期化
  d = []
  for i in range(l_p):
    d.append('0')
  d.insert(0, '.')

# 計算
  for i in range(1, l_a):

# 和格納用配列初期化
    c = []
    for j in range(1, l_a):
      c.append('0')
    c.insert(0, '.')
    for k in range(int(b[i])):
      c = plus(c, a)
    c.pop(0)
    for k in range(1, i):
      c.insert(0,'0')
    c.insert(0, '.')
    d = plus(d, c)
  d.pop(0)
  d.insert(p_a + p_b, '.')
  return(zerosupress(d))

#==========================================
# 配列同士の比較
#==========================================
# 戻り値
# a1 > a2 -- 2
# a1 < a2 -- 1
# a1 = a2 -- 0
# エラー --- -1
#==========================================
def compare(a1, a2):
  a = list(a1)
  b = list(a2)

# 小数点以下の桁数を合わせる
  p = a.index('.') - b.index('.')
  if p > 0:
    for i in range(p):
      b.insert(0, '0')
  if p < 0:
    for i in range(-p):
      a.insert(0, '0')

# 整数部分の桁数を合わせる
  p = (len(a) - a.index('.')) - (len(b) - b.index('.'))
  if p > 0:
    for i in range(p):
      b.append('0')
  if p < 0:
    for i in range(-p):
      a.append('0')
  l_a = len(a)
  r = -1

# 比較
  for i in range(l_a):
    j = l_a - i -1
    if a[j] > b[j]:
      r = 2
      break
    if a[j] < b[j]:
      r = 1
      break
    if a[j] == b[j]:
      r = 0
  return(r)

#==========================================
# 計算結果を表示
#==========================================
def display(a):
  if len(a) == 0:
    return(False)
  l_a = len(a)
  for i in range(l_a):
    j = l_a - i - 1
    if j > 0:
      b = a[j]
    if j == 0 and a[j] == '.':
      b = ' '
    else:
      b = a[j]
    sys.stdout.write(b)
  print('\r')

#==========================================
# 平方根の計算
#==========================================
# x -- 平方根を求める値
# c -- 求める平方根の桁数
#==========================================
def sqrt(x, c):

# 平方根の初期値と保存用
  z = ['.', '0']
  z1 = []

# 加算値
  p = ['.', '1']
  i = j = 0
  r = 9
  while i < c:
    z2 = multi(z, z)
    r = compare(x, z2)

# 平方根が得られた
    if r == 0:
      return(z)

# 二乗値が目標値を超えた
    if r == 1:
      z = z1
      if p[len(p)-1] == '1':
        p[len(p)-1] = '0'
        p.insert(0, '1')
      else:
        p[0] = '0'
        p.insert(0, '1')
      i = len(z) - 1
      if i % 10 == 0:
        print(SEP + SEP)
        print(str(i) +'桁までの計算結果')
        display(z)
        print(SEP + SEP)

# 目標値より少ない
    if r == 2:
      z1 = list(z)
      z = zerosupress(plus(z, p))
  return(z)

#==========================================
# メイン
#==========================================
if __name__ == '__main__':
  try:
    xi = list(str(sys.argv[1]))
    c = int(sys.argv[2])
    x_len = len(xi)
    if c == 0:
      raise exception
  except:
    print('計算回数を入力してください')
    print('Usage: sqrt.py xxx count')
    sys.exit()
    
  x = []
  y = []
  p = 0
  for i in range(x_len):
    if xi[i] == '.':
      p += 1
    x.insert(0, xi[i])

# ピリオドが2個以上ならエラー
  if p > 1:
    print('入力した数字が誤っています')
    sys.exit()

# ピリオドが含まれていない場合は最終桁にピリオドを追加する
  if p == 0:
    x.insert(0, '.')

# 計算開始時刻を記録
  print(SEP + SEP)
  st = datetime.datetime.now()

# 平方根を計算
  y = sqrt(x, c)

# 計算終了時刻を記録
  et = datetime.datetime.now()

# 表示
  print('\r')
  print(SEP + '結果' + SEP)
  display(y)
  print(SEP + '----' + SEP)

# 終了処理

  print('\r')
  print('\r')
  sys.stdout.write('計算開始時刻:')
  print(str(st))
  sys.stdout.write('計算終了時刻:')
  print(str(et))
  sys.stdout.write('計算時間:     ')
  print(str(et-st))

  sys.exit()


2018年06月13日