Montgomery 算法流程

下面默认对于模数\(m\)取模

Montgomery Reduction算法思想简介

在计算取模运算的过程中，将每一个元素\(T\)都乘上一个特定的值\(R(R>m,\gcd(R,m)=1)\)

用特殊的方法处理相乘时除掉一个\(R\)的过程，而避免取模运算

在是用的模数位常量时，编译器会自动加入类似的优化，因此实际上这个算法对于动态模数的情形更为适用

(你不一定写得过STL)

\[\ \]

编程上的应用简介

对于\(P\)为奇数的情况，取\(R=2^{32}\)，用自然溢出来加速运算

我们还需要令\(m' = -m^{-1} \mod R\)，有结论

对于某一个数\(T,0 \leq T < mR\)，若令\(U = Tm’ \mod R\)，则 \(\frac{T+Um}{R}\)为整数，且 \(\frac{T+Um}{R}=TR^{-1} \mod m\)

那么我们在计算\(\frac{T}{R}\)时，实际上只需要计算\(\frac{T+Um}{R}\)，可以预处理\(m'\)，溢出计算\(Tm'\)，位运算左移计算\(\frac{T+Um}{R}\)

实际使用时的实现，可以用一个类实现以下方法

在实现时需要尤其注意不要出现溢出

1.预处理\(m'\)

\((R-\lfloor \frac{R}{m}\rfloor )\cdot (R\mod m)\)

using u32=unsigned;
using i32=int;
using u64=unsigned long long;
using i64=long long;
// inv=m'
u32 m;
u32 getinv(){
	u32 inv=m;
	for(int i=0;i<4;++i) inv*=2-inv*m;
}

2.reduce方法

u32 reduce(u64 x) {
    u32 y = u32(x >> 32) - u32((u64(u32(x)*inv)*m) >> 32);
    // 先取u32(x)得到x mod R ，然后再转成u64进行乘法
    return i32(y) < 0 ? y + m : y;
}

3.普通数转Montgomery Reduction

我们要计算\(x\rightarrow xR=x\cdot 2^{32}\)，但是如果直接用取模就失去了意义。。。

方法是快速计算\(x\cdot R^2\)，然后reduce一次

u32 R2=-u64(m)%m;
u32 intToMont(i32 x){
    return reduce(u64(x)*R2);
}

\[\ \]

4.Montomery运算

u32 Add(u32 x,u32 y) {
    x+=y-m;
    return i32(x<0)?x+m:x;
}
u32 Dec(u32 x,u32 y){
    x-=y;
    return i32(x<0)?x+m:x;
}
u32 Mul(u32 x,u32 y){
    return reduce(u64(x)*y);
}

\[\ \]

5.Montomery Reduction转普通数

i32 get(u32 x){
    return reduce(x);
}

以上代码中的\(i32(x)<0\)均可替换为\(x>=P\)

封装之后，得到板子一号

using u32=uint32_t;
using i32=int32_t;
using u64=uint64_t;
using i64=int64_t;

static u32 m,inv,r2,P;
u32 getinv(){
    u32 inv=m;
    for(int i=0;i<4;++i) inv*=2-inv*m;
    return inv;
}
struct Mont{
private :
    u32 x;
public :
    static u32 reduce(u64 x){ 
        u32 y=(x+u64(u32(x)*inv)*m)>>32;
        return y>=m?y-m:y;
    }
    Mont(){ ; }
    Mont(i32 x):x(reduce(u64(x)*r2)) { }
    Mont& operator += (const Mont &rhs) { return x+=rhs.x-m,x>=m&&(x+=m),*this; }
    Mont& operator -= (const Mont &rhs) { return x-=rhs.x,x>=m&&(x+=m),*this; }
    Mont& operator *= (const Mont &rhs) { return x=reduce(u64(x)*rhs.x),*this; }
    friend Mont operator + (Mont x,const Mont &y) { return x+=y; }
    friend Mont operator - (Mont x,const Mont &y) { return x-=y; }
    friend Mont operator * (Mont x,const Mont &y) { return x*=y; }
    i32 get(){ return reduce(x); }
};
void Init(int m) { 
    ::m=m;
    inv=-getinv();
    r2=-u64(m)%m;
}