Hundredsoft

前回、モンゴメリ乗算 サンプルコードを書きましたが、よく使われる冪剰余についてもC++サンプルを載せておきます。

// モンゴメリ冪剰余演算 (a^b) mod nを求める
// gcd(n,r) == 1 となること。(nは奇数であれば良い)
//
#ifndef MONTGOMERY_INCLUDED
#define MONTGOMERY_INCLUDED

#include "Bigint.h"

class montgomery{
public:
    montgomery(const Bigint& n) : m_n(), m_r2(), m_nb(0), m_n2(){
        m_n = n;
        // n を法とする
        // Bit長を調べる(nb)
        m_nb = n.length();

        // r: nより大きな２のべき
        Bigint r(1);
        r <<= m_nb;

        // r2 = r^2 mod n
        m_r2 = r * r;
        m_r2 %= n;

        // n*n2 ≡ -1 mod r
        m_n2 = 0; /* 求めるN' */
        Bigint t = 0;
        Bigint vi = 1;
        int ni = m_nb;
        while (ni--){ /* Rのトップビットを除いたビット数分繰り返す */
            if ((t & 1) == 0){
                /* ゼロになっているビットがあったら、N'のその部分を1にする（NはRと互いに素なので必ず奇数）*/
                t += n; /* 掛け算だが、二進数一桁の掛け算なので実質は足し算 */
                m_n2 += vi; /* N'のその部分を1にする */
            }
            t >>= 1; /* 必ず端数が出るが切り捨てる */
            vi <<= 1; /* Rは2の冪なので、絶対端数は出ない */
        }
    }

    ~montgomery(){}

    Bigint reduction(const Bigint& t){
        // tのモンゴメリリダクション
        Bigint c = t * m_n2;
        c.hibitmask(m_nb); // mod Rの意味,(m_nb)bit以上を0クリア
        c *= m_n;
        c += t;
        c >>= m_nb; // 1/r の意味
        if (c >= m_n) c -= m_n;
        return c;
    }

    // 冪剰余 a^b mod n, バイナリ法の下位桁から計算
    Bigint exp(const Bigint& a, const Bigint& b){
        Bigint p = reduction(a * m_r2);
        Bigint x = reduction(m_r2);
        Bigint y(b);
        while (y > 0){
            if (y & 1){
                x = reduction(x * p);
            }
            p = reduction(p * p);
            y >>= 1;
        }
        return reduction(x);
    }

protected:
    montgomery() : m_n(), m_r2(), m_nb(0), m_n2(){}

private:
    Bigint m_n;   // 法(m_n)
    Bigint m_r2;  // rを(m_n)より大きな２のべきとして、(m_r2)=(r)^2 mod n
    Bigint m_n2;  // (m_n)*(m_n2) ≡ -1 mod r を満たす(m_n2)
    int    m_nb;  // (m_n)の有効Bit数, 例えば(m_n)=5 のとき(m_nb)=3
};

#endif

wikipediaの説明にあるようにバイナリ法を用いています。
初めに法(n)とするコンストラクタを呼び出します。
冪剰余は、exp()メソッドで行います。


このアルゴリズムを利用した16進・２進・10進電卓です。
整数だけですがおよそ100万桁まで計算できます。


(注)Bigintクラスは、多倍長整数を格納する仮想のクラスであり、実在するものではありません。

---

Tags: プログラムメモ