課題: EDSACのプログラム技法

• EDSAC - Wikipedia

  • パラメトロン計算機: EDSACのプログラム技法

  • From EDSAC to Raspberry Pi: 75 years of “computers that work” | University of Cambridge

• Methods of computing square roots - Wikipedia#A_two-variable_iterative_method

  • High order algorithms to find square roots and inverses

  • Pythagoras’ Constant : \(\sqrt{2}\)

  • (PDF) A Note on an Iterative Method for Root Extraction

\(f(x) = \frac{1}{x^2} - A\) とする. ニュートン法により

(1)\[\begin{eqnarray} x_{n+1} &=& x_{n} - \frac{{x_{n}}^{-2}-A}{-2{x_{n}}^{-3}} \\ &=& x_{n} + \frac{1}{2} x_{n} \left( 1 - A {x_{n}}^2 \right) \\ &=& x_{n} - \frac{1}{2} x_{n} c_{n} \\ &=& x_{n} \left( 1-\frac{c_n}{2} \right) \end{eqnarray}\]

ここで \(c_{n} = A{x_{n}}^2-1\) と置き換えた.

この式を \({x_{n}}^2 = \frac{1}{A} (1+c_{n})\) と変形し,上式の両辺を自乗して代入すると

\[(1+c_{n+1}) = (1+c_{n}) \left( 1-\frac{c_{n}}{2}\right)^{2}\]

(2)\[\begin{eqnarray} c_{n+1} &=& (1+c_{n})(1-c_{n}+\frac{{c_{n}}^2}{4}) - 1 \\ &=& \frac{1}{4}{c_{n}}^2 \left( c_{n} - 3 \right) \end{eqnarray}\]

この二重漸化式は, 除算を行う必要がないため (~~\(1/2\) や \(1/4\)は有限小数の乗算~~ビットシフト演算), 高精度で高速な演算を期待できる.

この式は, 次のように一般化できる:

(3)\[\begin{eqnarray} x_{n+1} &=& x_{n} \left( 1 - \frac{c_{n}}{n}\right) \\ 1 + c_{n+1} &=& (1 + c_{n})\left( 1 - \frac{c_{n}}{n} \right)^{n} \end{eqnarray}\]

\(x_{0} = a\), \(c_{0} = b-1\) かつ \(c_{n} \rightarrow 0\) のとき

\[x_{n} \rightarrow \frac{a}{b^{\frac{1}{n}}}\]

x, c = 2, 2-1

while True:
    x = x - x * c / 2
    c = c**2 * (c - 3) / 4
    print(x, c)
    if abs(c) < 0.00000001:
        break

1.0 -0.5
1.25 -0.21875
1.38671875 -0.03850555419921875
1.413416936993599 -0.001126281109816385
1.4142128893918142 -9.517390282169533e-07
1.4142135623726146 -6.793555988965238e-13

%%ai vioras
次のプログラムについて

```python
x, c = 2, 2-1

while True:
    x = x - x * c / 2
    c = c**2 * (c - 3) / 4
    print(x, c)
    if abs(c) < 0.00000001:
        break
```

このプログラムは、xとcという二つの変数に対して特定の計算を繰り返し行うものです。ここでの計算は以下の通りです:

1. xの新しい値は、現在のxからxとcの積の半分を引いたものになります。

2. cの新しい値は、現在のcの2乗と(c - 3)の積を、4で割ったものになります。

3. その後、新しいxとcを出力します。

この計算を続けるかどうかは、cの絶対値が0.00000001より小さくなった場合に停止します。

初期値はx = 2およびc = 2 - 1 = 1となっています。

ただし、特定の数学的な背景や目的があるかどうかはこのコードだけからは明確には判断できません。たとえばこれは特定の数列を生成するためのものであったり、特定の制約下での最適化問題を解くためのものである可能性もあります。それは具体的な応用コンテキストや、このコードが書かれた背景によります。

%%ai vioras
`c`は0に収束し、`a`は$2$の平方根を近似します。

ありがとうございます、その説明が非常に有用です。このコードは、実際には平方根の近似値を計算するためのもので、特に$2$の平方根を計算しています。

まず初期値としてx=2とc=1を設定します。そしてループ内では、以下の操作を繰り返します:

1. xからx * c / 2を引くことで、xの新しい値を計算します。

2. cの新しい値は、cの2乗と(c - 3)の積を4で割ったものに更新します。

そして、cの絶対値が十分に小さくなる（具体的には 0.00000001より小さくなる）までこれを繰り返します。

このプロセスの結果、cは0に収束し、xは$2$の平方根に収束します。これは、xの更新式が$2$の平方根を逐次近似する方法となっているからです。ただし、この近似方法自体がどのように導出されたか、またはどのような理論的背景に基づいているかについては、このコードだけからは判断できません。

%%ai vioras
このアルゴリズムは、EDSACのために開発されたニュートン法に基づく二変数の漸化式です。

なるほど、それは興味深い情報です。このアルゴリズムは、EDSAC (Electronic Delay Storage Automatic Calculator) のために開発されたニュートン法に基づいています。ニュートン法は、方程式の根を求めるための反復的なアルゴリズムで、初期近似値から始めてその値を逐次改良していきます。

この具体的なコードは、二変数の漸化式を用いて、2の平方根の近似値を求めています。xとcの更新ルールは、ニュートン法の更新ステップを表現しています。

ニュートン法の基本的なアイデアは、関数 $f(x)$ の近似解を見つけるために、初期点 $x_0$ から始めて、タンジェントラインを使って次の近似値 $x_{n+1}$ を求めるというものです。具体的には、次の式で更新を行います：

$$ x_{n+1} = x_n - \frac{f(x_n)}{f’(x_n)} $$

ここで、$f’(x)$ は $f(x)$ の導関数です。このプロセスを反復すると、$x_n$ は $f(x) = 0$ の解に収束します。

このコードでは、更新ステップを x = x - x * c / 2 と c = c**2 * (c - 3) / 4 で表現しており、これらはニュートン法に基づく近似ステップを表しています。そして、c が十分に小さくなるまでこの更新を繰り返すことで、x は2の平方根に収束します。

%%ai vioras
$f(x) = \frac\\{1\\}\\{x^2\\} - A$です。

理解しました。この場合、関数 $f(x)$ は次のようになります：

$$ f(x) = \frac{1}{x^2} - A $$

この関数の根を求めることは、$ \frac{1}{x^2} = A $ つまり $x = \sqrt{1 / A}$ を求めることに相当します。この例では $A=2$ なので、求める数値は $2$ の平方根となります。

ニュートン法を用いると、次のような更新式が得られます。まず、$f(x)$ の導関数 $f’(x)$ を求めます：

$$ f’(x) = -2 \cdot \frac{1}{x^3} $$

次に、ニュートン法の更新式を適用します：

$$ x_{n+1} = x_n - \frac{f(x_n)}{f’(x_n)} $$

これを $f(x)$ および $f’(x)$ に代入すると、次のような更新式が得られます：

$$ x_{n+1} = x_n - \frac{\frac{1}{x_n^2} - A}{-2 \cdot \frac{1}{x_n^3}} $$

これを整理すると、次のような形になります：

$$ x_{n+1} = x_n \left(1 + \frac{A x_n^2 - 1}{2} \right) $$

これがx = x - x * c / 2と一致するためには、$c$は次のように定義されている必要があります：

$$ c = A x_n^2 - 1 $$

そして、c の更新ルール c = c**2 * (c - 3) / 4 を満たすためには、$c$ が次のような漸化式で更新されている必要があります：

$$ c_{n+1} = (A x_{n+1}^2 - 1) = c_n^2 \left(c_n - 3 \right) / 4 $$

この関係が成り立つかどうかは、その適用される具体的な数学的背景によります。

%%ai vioras
$c_0$によって収束は異なります。

はい、その通りです。ニュートン法のような反復的なアルゴリズムでは、初期値の選択が結果の収束に大きな影響を与えます。

具体的には、初期値$c_0$が異なると、その後の値$c_n$の進行と収束が異なる可能性があります。ニュートン法は局所的な方法であり、初期値が解の近くになければ適切に収束しない可能性があります。また、初期値が特定の”不運な”値であると、アルゴリズムは収束しないか、あるいは非常に遅く収束する可能性があります。

そのため、このアルゴリズムを使用する場合、適切な初期値を選択することが重要です。通常、問題の具体的な知識を用いて初期値を選択します。このコードでは、初期値として$c_0 = 2 - 1 = 1$が選択されています。

%%ai vioras
EDSACが平方根の求解に何故このアルゴリズムを必要としたか考察してください。

EDSAC (Electronic Delay Storage Automatic Calculator) は、1940年代に開発された最初期の一般目的のデジタルコンピュータの一つです。その当時、コンピュータの計算能力は現代と比べて非常に限られており、複雑な計算を効率的に行うためのアルゴリズムが必要とされていました。

平方根の計算は、多くの科学的、工学的な問題で必要とされる基本的な演算であり、その効率的な計算方法は重要です。しかし、当時のコンピュータでは乗除算の計算コストが高かったため、それを避ける工夫が必要でした。

ニュートン法はこの問題に対する一つの解答で、反復的に解を改良していく方法で、初期値から始めてその値を逐次改良していきます。その更新ステップは比較的単純で、乗除算を必要としません。

このコードのアルゴリズムは、2の平方根を求めるニュートン法の一種で、xとcの二つの変数を用いて反復計算を行っています。xの更新はx = x - x * c / 2という単純な形をしており、これは乗算と減算のみで計算できます。また、cの更新もc = c**2 * (c - 3) / 4という形で、乗算と減算のみで計算できます。

したがって、このアルゴリズムは、当時のコンピュータの限られた計算能力を考慮に入れた上で、平方根の計算を効率的に行うための工夫された方法と言えるでしょう。

%%ai vioras
除算に関する考察が足りません。

すみません、その点は見落としていました。確かに除算は現代のコンピュータにおいても、加算や乗算に比べてコストの高い演算です。1940年代の初期のコンピュータ、特にEDSACでは、除算は非常に高コストで、時間と資源を大量に消費する可能性がありました。

ニュートン法は、一般的には除算を含む更新式 $x_{n+1} = x_n - \frac{f(x_n)}{f’(x_n)}$ に基づいていますが、このアルゴリズムは除算を避けるように工夫されています。

x の更新 x = x - x * c / 2 では、除算は明示的には行われておらず、代わりに乗算が用いられています。c の更新 c = c**2 * (c - 3) / 4 では、定数 4 での除算が行われていますが、これは固定の除数であるため、シフト演算による最適化等が可能で、高コストな除算を避けることができます。

したがって、このアルゴリズムは除算のコストを抑えるように設計されていたと言えるでしょう。これは、初期のコンピュータの限られた計算能力を考慮した、効率的なアルゴリズム設計の一例です。

%%ai vioras
EDSACは、$2$と$4$の除算をビットシフト演算で行っていました。

その通りです。ビットシフト演算は、2のn乗での乗算や除算を行う際に、非常に効率的な手法です。

一般的に、2進数表現においては、左にビットを1つシフトすることは2倍（乗算）を、右にビットを1つシフトすることは半分（除算）を意味します。例えば、2進数の 1001 (10進数で9) を1ビット右にシフトすると 100 (10進数で4) になり、これは9を2で除した値と等しくなります。

EDSACはこの特性を利用し、2や4での除算をビットシフト演算で行っていたとのことです。したがって、このコードにおける c = c**2 * (c - 3) / 4 の操作は、高コストな除算を行っているように見えますが、実際には効率的なビットシフト演算を行っていたと言えます。このような工夫により、計算資源が限られていた当時のコンピュータでも効率的に平方根を計算することが可能になっていたというわけです。

Note

以下、補足です

位取り記法の基数による乗算と除算

%%ai vioras
位取り記法の基数による乗算と除算が簡単な理由

位取り記法の基数による乗算と除算が簡単な理由は、以下のような要点に集約されます：

1. 基数の特性: 基数の性質により、乗算と除算は互いに逆の操作となります。たとえば、10進法では、乗算は10の累乗によるスケーリングと解釈でき、除算はその逆の操作です。同じように、2進法では、左にビットをシフトすることは乗算に、右にシフトすることは除算に相当します。

2. 位取り記法: 位取り記法は、数字の位置によってその値が変わるという特性があります。これにより、乗算や除算を行う際に、各桁ごとに独立して計算を行うことが可能になります。これは計算を簡単にし、特に大きな数に対する計算を可能にします。

3. 再利用性: 一度学習した乗算や除算の手順は、任意の基数に対して再利用することができます。たとえば、10進法で学んだ乗算の手順は、2進法や16進法でもそのまま利用することが可能です。

4. 計算の分割: 乗算や除算は、複数の小さな計算に分割できます。このため、全体の計算を一度に行う代わりに、各部分を個別に計算し、最終的に結果を組み合わせることが可能になります。これにより、計算の複雑さを大幅に減らすことができます。

これらの理由から、位取り記法の基数による乗算と除算は、比較的簡単な手順で行うことができます。

%%ai vioras
４番目の計算の分割について具体例を示してください

もちろんです。4番目の「計算の分割」について、10進法を例に具体的な説明をします。

10進法における乗算は、各桁ごとに行われ、それぞれの結果が加算されることで最終的な結果が得られます。例えば、123 × 45を計算する場合、以下のように分割して計算を行います。

   123 
×  45 
------
  615 (これは123 × 5の結果)
+4920 (これは123 × 4の結果、一つ左にシフトされている)
------
 5535 (これが最終的な結果)

上記の計算では、123という数を1の位の3、10の位の2、100の位の1という3つの数に分割し、それぞれ45と乗算しています。それぞれの乗算の結果を得た後、それらの結果を合計することで最終的な答えを得ています。このように、位取り記法を用いることで大きな数の乗算を、より小さな数の乗算に分割し、計算をシンプルにすることができます。

同様に、除算もある数を基数の累乗で表現し、それぞれの部分で除算を行い、最後に結果を組み合わせることで計算を行うことができます。これにより、大きな数の除算を、より小さな数の除算に分割し、計算を簡単にすることができます。