secp256k1における平方剰余計算

secp256k1における平方剰余計算を考えてみます。

secp256k1の法素数
$q = {2}^{256} - {2}^{32} - {2}^{9} - {2}^{8} - {2}^{7} - {2}^{6} - {2}^{4} - 1$
はこの形をみてわかる通り、 $q \mod 4 = 3$ なので、
この記事mod p での平方剰余を計算する(p mod 4 = 3の場合) - Pebble Coding
で示したように、解は(存在する場合は)
$x = {a} ^ { \frac {q+1} {4} }$
と書けます。
コンピュータで計算する場合は、次のように計算します。
まず、pythonで値を見てみます。
q+1は4で割り切れるので、python3での整数の割り算には//が使用できます。

>>> q = 2 ** 256 - 2 ** 32 - 2 ** 9 - 2** 8 - 2 ** 7 - 2 ** 6 - 2 ** 4 -1
>>> print(q)
115792089237316195423570985008687907853269984665640564039457584007908834671663
>>> print((q+1)//4)
28948022309329048855892746252171976963317496166410141009864396001977208667916

この値の２進表現(バイナリ表現)を見てみます。

>>> print(format((q+1)//4, 'b'))
11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111110111111111111111111111100001100

とても綺麗な形をしており、1が連続する箇所が3箇所しかありません。
最初の1は223個、次の1は22個、最後の1は2個が連続で続いています。
やることはaをこの値分だけベキ乗することです。
問題を少し単純化してみてみましょう。

この２進数が仮に
10010 だったとしましょう。
計算すべきは ${a}^{10010}$ です。
${a}^{1001}$ を2乗すると ${a}^{10010}$ になります。
上位へビットシフトしたら2倍になるので、こうなります。
手順は、まず ${a}^{1}$ を計算し、3回ビットシフトつまり ${a}^{1}$ の2乗を3回行い、それにaを掛け、その値を1回ビットシフト、つまりその値を2乗すると ${a}^{10010}$ が得られます。
この考えを用いて今度は、 ${a}^{1が223個並んだもの}$ を計算したいのですが、効率よく計算するには次のようにします。

表現を簡単にするため、 ${a}^{2進表現}$ のことを以後、x111001...001のように表現します。
まず、x1を2乗するとx10となります、これにaを掛けるとx11となります。
これを2乗しx110となり、それにaを掛け、x111とします。
これの2乗を3回して、x111_000としますが、先ほど計算したx111と掛けx111_111となります。
これの2乗を3回して、x111_111_000となりますが、先ほど計算したx111を掛け、x111_111_111となります。
これの2乗を2回して、x11_111_111_100となりますが、先ほど計算したx11を掛け、x11_111_111_111となります。
これの2乗を11回して,...
これの2乗を22回して,...
これの2乗を44回して,...
これの2乗を88回して,...
...
こんな感じで、1が223個並んだ値が計算でき、このロジックを利用して、最終的に、
${a}^{11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111110111111111111111111111100001100}$ を計算します。

コンピュータで解があるかどうかはこうして計算したxを使って $x^{2} mod q$ を計算してaに等しければ解あり、等しくなければ解なしと判断します。

参考:

secp256k1/field_impl.h at master · bitcoin-core/secp256k1 · GitHub