はじめに

『べき乗数を法としたモジュラ逆数はニュートン法で求められる』ということを知ったので、コードに落としてみました。

モジュラ逆数とは

モジュラ逆数についてよく知っている人はこの章を読み飛ばして結構です。
モジュラ逆数とは、魔法の数です…という説明では何も分からないと思います。とりあえず下記のコードを実行してみてください。

use std::num::Wrapping;

fn main() {
    let c = Wrapping::<u32>(0xaaaa_aaab);
    println!("c = 0x{:x}", c);
    for i in 0..20 {
        let i = Wrapping::<u32>(i);
        println!("c * {} = {}", i, c * i);
    }
}

補足ですが、Wrappingは、算術オーバーフローが発生したときのエラーを無視するための記述です。
結果も置いておきます。

c = 0xaaaaaaab
c * 0 = 0
c * 1 = 2863311531
c * 2 = 1431655766
c * 3 = 1
c * 4 = 2863311532
c * 5 = 1431655767
c * 6 = 2
c * 7 = 2863311533
c * 8 = 1431655768
c * 9 = 3
c * 10 = 2863311534
c * 11 = 1431655769
c * 12 = 4
c * 13 = 2863311535
c * 14 = 1431655770
c * 15 = 5
c * 16 = 2863311536
c * 17 = 1431655771
c * 18 = 6
c * 19 = 2863311537

乗算しか使っていないのに、3で割る処理が実現していることが分かりますね。
0xaaaaaaabが法$2^{32}$における3のモジュラ逆数だからです。
法nにおけるaのモジュラ逆数というのは、nで割った余りを考えたときに、aの逆数のような働きをする数を指しています。
現代的なコンピュータは、オーバーフローが発生したときに$2^{32},2^{64}$といった2のべき乗で割った余りを返せるため、モジュラ逆数をかけることで除算が実現できるのです。
こことかここで応用例を見ることができるでしょう。

ニュートン法

以下の漸化式は、実数$a$の逆数をニュートン法で求める式です。
$$x_{n+1} = x_n(2 - ax_n)$$
誤差を$x_n = d_n + a^{-1}$と置いてみるこんな感じです。

d_{n+1} + a^{-1} = (d_n + a^{-1})(1 - ad_n) \\
d_{n+1} = -a{d_n}^2

従って、実数の場合には初期条件の誤差が$|x_0 - a^{-1}| < a^{-1}$であれば収束するのですが、モジュラ逆数の場合は少々事情が違います。
$d_0\in 偶数$のとき同様に漸化式を適用すると、$d_1 \in 4の倍数$、$d_2 \in 2^{4}の倍数$、$d_3 \in 2^{8}の倍数$…となって一見すると誤差が増えていくように見えます。
ところが、法$2^{32}$の場合は$d_5$以降が突如$0$になって収束します。
また、$d_0\in 奇数$だったり、法が2のべき乗でなかったり、$a^{-1}$が存在しない場合などは収束しないこともあります。
ともあれ、法$2^{32}$の場合は以下のコードでモジュラ逆数が計算できます。

use std::num::Wrapping;

fn modular_inv(a: Wrapping<u32>) -> Option<Wrapping<u32>> {
    if a.0 % 2 == 0 {
        return None;
    }
    let mut x = a;
    for _ in 0..4 {
        x = x * (Wrapping(2) - a * x);
    }
    Some(x)
}

まとめ

法nが$2^{32}$のとき、ニュートン法でモジュラ逆数を求めるプログラムを示しました。
べき乗以外のときにも効率的なアルゴリズムはあるんでしょうか?

モジュラ逆数をニュートン法で求める

はじめに

モジュラ逆数とは

ニュートン法

まとめ

参考