開平法、開立法

# はじめに(Introduction)

[SHA-256の実装]()の中で、ハッシュ値の初期値や固定値は、素数の平方根、立方根だったのですが
通常の精度では、期待の値が取れなかったので、実装しようと考えました。
調べていくと、開平法、開立法といった仕組みで算出可能でした。
それらの原理と実装を行います。

# 開平法(extraction of square root)

$\sqrt{x}$を各桁の和で表します。$0\le a_i \le 9$
また、$m$の位より大きい位を整数にしたものを$p_m$とします。


```math
\begin{align*}
\sqrt{x} &= \sum_{k\le n} 10^ka_k \\
&= 10^na_n+\cdots+10a_1+a_0+10^{-1}a_{-1}+\cdots \\
\\
p_m &=\sum_{k=m+1}^{n}10^{k-m-1}a_k \\
&= 10^{n-m-1}a_n+\cdots +10a_{m+2}+a_{m+1} &\\
\end{align*}
```

例えば、$x=2,m=-2$の場合は次のようになります。

```math
\begin{align*}
\sqrt{x} &= 1+10^{-1}\cdot 4+10^{-2}\cdot 1+10^{-3}\cdot 4+\cdots \\
p_{-2} &= 10\cdot 1+4 \\
\end{align*}
```

$\sqrt{x},p_m$から以下の不等式が成り立ちます。

```math
\begin{align*}
10^{m+1}p_m + 10^ma_m &\le \sqrt{x} \\
(10^{m+1}p_m + 10^ma_m)^2 &\le (\sqrt{x})^2 \\
10^{2m+2}\cdot p_m^2 + 2\cdot 10^{2m+1}\cdot p_m\cdot a_m + 10^{2m}\cdot a_m^2 &\le x \\
(2\cdot 10\cdot p_m + a_m)a_m &\le 10^{-2m}\cdot x - 10^{2}\cdot p_m^2 \\
\end{align*}
```

この不等式を満たす最大の$a_m$を見つけて行きます。
ここで、左辺$(2\cdot 10\cdot p_m + a_m)a_m$は必ず整数となるので、
右辺$10^{-2m}\cdot x - 10^{2}\cdot p_m^2$の整数部分のみを考えれば良いわけです。

左辺を$r_m=(2\cdot 10\cdot p_m + a_m)a_m$
右辺を$y_m=10^{-2m}\cdot x - 10^{2}\cdot p_m^2\ $と置きます。
右辺$y_m$の整数部分を$z_m$とした時、$z_{m+1}$を求めます。

```math
\begin{align*}
y_m &= 10^{-2m}\cdot x - 10^{2}\cdot p_m^2 \\
&= 10^2\{10^{-2(m+1)}\cdot x - (10\cdot p_{m+1}+a_{m+1})^2\} \\
&= 10^2\{10^{-2(m+1)}\cdot x - 10^2\cdot p_{m+1}^2\}-10^2(2\cdot 10\cdot p_{m+1}+a_{m+1}) a_{m+1} \\
&= 10^2\cdot y_{m+1} -10^2\cdot r_{m+1} \\
\end{align*}
```

$x$を桁毎の和とした時に下のようになるとします。

```math
\begin{align*}
&x = 10^nx_n+\cdots +10\cdot x_1+x_0+10^{-1}x_{-1}+\cdots
\end{align*}
```

$z_m$は$y_m$の整数部なので、$10^2\cdot y_{m+1}$の整数部は、$10^2\cdot z_{m+1}+10\cdot x_{2m+1} + x_{2m} $となります。
小数部分が$10^2$により2桁整数部になります。
整数部だけを考えると以下のようになります。

```math
\begin{align*}
y_m &= 10^2\cdot y_{m+1} -10^2\cdot r_{m+1} \\
z_m &= 10^2\cdot z_{m+1}+10\cdot x_{2m+1} + x_{2m} -10^2\cdot r_{m+1} \\
 &= 10^2(z_{m+1}-r_{m+1})+10\cdot x_{2m+1} + x_{2m}  \\
\end{align*}
```

最上位から最下位へ向けて計算していきます。

```math
\begin{align*}
z_m &= 10^2(z_{m+1}-r_{m+1})+10\cdot x_{2m+1} + x_{2m} \\
r_m &= (2\cdot 10\cdot p_m+a_m)a_m \\
& r_m \le z_m \text{ を満たす最大の }a_m\text{ を求めます。}\\
p_{m-1} &= 10\cdot p_m + a_m \\
\end{align*}
```

### 例

$\sqrt{2}$を考えます。
$z_1=0,r_1=0,p_0=0$から始めます。


```math
\begin{align*}
z_0&=10^2(z_1-r_1)+10\cdot 0 + 2 \\
r_0&=(2\cdot 10\cdot p_0+a_0)a_0 = a_0\cdot a_0 \\
&r_0\le z_0\text{ を満たす最大の }a_0\text{ を求めます。} &\\
\\
&z_0=2,a_0=1,r_0=1 \\
&p_{-1}=10\cdot p_0 + a_0 = 1 \\
\end{align*}
```

$10^0$の位は$a_0=1$となります。

```math
\begin{align*}
z_{-1}&=10^2(z_0-r_0)+10\cdot 0 + 0 = 100 \\
r_{-1}&=(2\cdot 10\cdot p_{-1} + a_{-1})a_{-1} =(20+ a_{-1}) a_{-1} \\
&r_1\le z_1\text{ を満たす最大の}a_{-1}\text{ を求めます。} \\
\\
&z_{-1}=100,a_{-1}=4,r_{-1}=96 \\
&p_{-2}=10\cdot p_{-1} + a_{-1} = 14 \\
\end{align*}
```

$10^{-1}$の位は$a_{-1}=4$となります。

```math
\begin{align*}
z_{-2}&=10^2(z_{-1}-r_{-1})+10\cdot 0 + 0 = 400 \\
r_{-2}&=(2\cdot 10\cdot p_{-2} + a_{-2})a_{-2} =(280+ a_{-2}) a_{-2} \\
&r_2\le z_2\text{ を満たす最大の}a_{-2}\text{ を求めます。} \\
\\
&z_{-2}=400,a_{-2}=1,r_{-2}=281 \\
&p_{-3}=10\cdot p_{-2} + a_{-2} = 141 \\
\end{align*}
```

$10^{-2}$の位は$a_{-2}=1$となります。

このように再帰的に求める事ができます。

# 開立法(extraction of cubic root)

開平法と同様に考えます。
$\sqrt[3]{x}$を各桁の和で表します。$0\le a_i \le 9$
また、$m$の位より大きい位を整数にしたものを$p_m$とします。


```math
\begin{align*}
\sqrt[3]{x} &= \sum_{k\le n} 10^ka_k \\
&= 10^na_n+\cdots+10a_1+a_0+10^{-1}a_{-1}+\cdots \\
\\
p_m &=\sum_{k=m+1}^{n}10^{k-m-1}a_k \\
&= 10^{n-m-1}a_n+\cdots +10a_{m+2}+a_{m+1} &\\
\end{align*}
```

$\sqrt[3]{x},p_m$から以下の不等式が成り立ちます。

```math
\begin{align*}
10^{m+1}p_m + 10^ma_m &\le \sqrt[3]{x} \\
(10^{m+1}p_m + 10^ma_m)^3 &\le (\sqrt{x})^3 \\
10^{3m+3}\cdot p_m^3+3\cdot 10^{3m+2}\cdot p_m^2\cdot a_m +3\cdot 10^{3m+1}\cdot p_m\cdot a_m^2 + 10^{3m}\cdot a_m^3 &\le x \\
(3\cdot 10^2\cdot p_m^2+3\cdot 10\cdot a_m+a_m^2)a_m &\le 10^{-3m}\cdot x - 10^{3}\cdot p_m^3 \\
\end{align*}
```

開平法と同様に以下の手順で求めることができます。

```math
\begin{align*}
z_m &= 10^3(z_{m+1}-r_{m+1})+10^2\cdot x_{3m+2} + 10\cdot x_{3m+1} + x_{3m} \\
r_m &= (3\cdot 10^2\cdot p_m^2+3\cdot 10\cdot a_m+a_m^2)a_m \\
& r_m \le z_m \text{ を満たす最大の }a_m\text{ を求めます。}\\
p_{m-1} &= 10\cdot p_m + a_m \\
\end{align*}
```

# 16進数(hex)

10進数で考えていた為、10を基数として考えていました。
16進数の場合、16を基数とすれば良いだけです。
基数を$b$とした時、開平法と開立法は以下のとおりになります。

### 開平法(extraction of square root)

```math
\begin{align*}
z_m &= b^2(z_{m+1}-r_{m+1})+b\cdot x_{2m+1} + x_{2m} \\
r_m &= (2\cdot b\cdot p_m+a_m)a_m \\
& r_m \le z_m \text{ を満たす最大の }a_m\text{ を求めます。}\\
p_{m-1} &= b\cdot p_m + a_m \\
\end{align*}
```

### 開立法(extraction of cubic root)

```math
\begin{align*}
z_m &= b^3(z_{m+1}-r_{m+1})+b^2\cdot x_{3m+2}+b\cdot x_{3m+1}+x_{3m} \\
r_m &= (3\cdot b^2\cdot p_m^2+3\cdot b\cdot a_m+a_m^2)a_m \\
& r_m \le z_m \text{ を満たす最大の }a_m\text{ を求めます。}\\
p_{m-1} &= b\cdot p_m + a_m \\
\end{align*}
```

# 実装(Programing)

javascriptで実装してみます。
$x$は整数、小数点以下16桁の文字列を返します。

#### 平方根

```js:sqrt
function sqrt(v) {
    let n = BigInt(v);
    let x = [];
    let b = 16n;
    while (n != 0n) {
        let r = n % b;
        n = (n - r) / b;
        x.push(r);
    }
    if (x.length % 2 == 1) {
        x.push(0n);
    }
    x.reverse();
    let z = 0n;
    let r = 0n;
    let p = 0n;
    let s = "";
    let d = 16;
    for (let i = 0; i < x.length + d * 2; i += 2) {
        let xx = 0n;
        if (i < x.length) {
            xx = b * x[i] + x[i + 1];
        }
        if (i == x.length) {
            s += ".";
        }
        let a = b - 1n;
        // z = b^2(z-r)+xx
        z = b * b * (z - r) + xx;
        while (true) {
            // r = (2*b*p+a)a
            r = (2n * b * p + a) * a;
            if (r <= z) {
                p = b * p + a;
                s += a.toString(Number(b));
                break;
            }
            a--;
        }
    }
    return s;
}
```

#### 立方根

```js:cbrt
function cbrt(v) {
    let n = BigInt(v);
    let x = [];
    let b = 16n;
    while (n != 0n) {
        let r = n % b;
        n = (n - r) / b;
        x.push(r);
    }
    while (x.length % 3 != 0) {
        x.push(0n);
    }
    x.reverse();
    let z = 0n;
    let r = 0n;
    let p = 0n;
    let s = "";
    let d = 16;
    for (let i = 0; i < x.length + d * 3; i += 3) {
        let xxx = 0n;
        if (i < x.length) {
            xxx = b * b * x[i] + b * x[i + 1] + x[i + 2];
        }
        if (i == x.length) {
            s += ".";
        }
        let a = b - 1n;
        z = b * b * b * (z - r) + xxx;
        while (true) {
            r = (3n * b * b * p * p + 3n * b * p * a + a * a) * a;
            if (r <= z) {
                p = b * p + a;
                s += a.toString(Number(b));
                break;
            }
            a--;
        }
    }
    return s;
}
```

# さいごに(Conclusions)

ソースコードは、[Github](https://github.com/tnakagawa/education/tree/master/root)にあります。
結果は、[ここ](https://tnakagawa.github.io/education/root/root.html)で見る事ができます。

# 参考文献(Reference)

[開平法 - Wikipedia](https://ja.wikipedia.org/wiki/%E9%96%8B%E5%B9%B3%E6%B3%95)
[開立法 - Wikipedia](https://ja.wikipedia.org/wiki/%E9%96%8B%E7%AB%8B%E6%B3%95)