はじめに(Introduction)
SHA-256の実装の中で、ハッシュ値の初期値や固定値は、素数の平方根、立方根だったのですが
通常の精度では、期待の値が取れなかったので、実装しようと考えました。
調べていくと、開平法、開立法といった仕組みで算出可能でした。
それらの原理と実装を行います。
開平法(extraction of square root)
$\sqrt{x}$を各桁の和で表します。$0\le a_i \le 9$
また、$m$の位より大きい位を整数にしたものを$p_m$とします。
\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$の場合は次のようになります。
\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$から以下の不等式が成り立ちます。
\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}$を求めます。
\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$を桁毎の和とした時に下のようになるとします。
\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桁整数部になります。
整数部だけを考えると以下のようになります。
\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*}
最上位から最下位へ向けて計算していきます。
\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$から始めます。
\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$となります。
\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$となります。
\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$とします。
\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$から以下の不等式が成り立ちます。
\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*}
開平法と同様に以下の手順で求めることができます。
\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)
\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)
\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桁の文字列を返します。
平方根
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;
}
立方根
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にあります。
結果は、ここで見る事ができます。