别问我数学问题, 数学原理我一点不懂, 我只是照着网上Montgomery模乘的论文写下来的.
======================

#include <iostream>
using namespace std;

unsigned long ModInv(unsigned long A, unsigned long P)
{
    unsigned long u = A, v = P;
    unsigned long x1 = 1, x2 = 0, flag=0;

    unsigned long q, r, x;
    while (u != 1)
    {
        q = v/u; r = v%u; x = q*x1 + x2;

        v = u; u = r; x2 = x1; x1 = x;

        flag ^= 1;
    }

    if (flag == 1)
    {
        P -= x1;
        return P;
    }

    return x1%P;
};

//    偏移数.
unsigned deg_num(unsigned long A)
{
    unsigned i = 0;
    while (A != 1)
    {
        A >>= 1;
        i += 1;
    }
    return i;
};
unsigned long deg_bit(unsigned long A)
{
    unsigned long i = 1;
    while (A != 1)
    {
        A >>= 1;
        i <<= 1;
    }
    return i;
};

// return R^(-1) * R * R mod P
unsigned long GetendR(unsigned long R, unsigned long P)
{
    unsigned long long t;

    // t = R^(-1)
    t = ModInv(R, P);

    t *= R;
    t %= P;
    t *= R;
    t %= P;
    return (unsigned long)t;
};

// Montgomery模乘: 对R的"模"和"除"都可以通过"移位"和"比特与"来进行, 这里只是演示过程.
// 而这恰恰是提高效率的部分, 省去模和除.
unsigned long Mont(unsigned long A, unsigned long B, unsigned long P1, unsigned long P, unsigned long R)
{
    unsigned long long t1, t2;

    // t1 = A * B;
    t1 = A;
    t1 *= B;

    // t2 = (t1 * P1) % R;
    t2 = (t1 * P1) % R;

    // P1 是 P mod R的负逆, 所以减.
    t2 = R - t2;

    // t1 = (t1 + t2 * P) / R;
    t1 += (t2 * P);
    t1 /= R;

    if (t1 >= P) t1 -= P;
    return (unsigned long)t1;
};

int main()
{
    cout << "Hello, world" << endl;

////////////////////////////////////////

    unsigned long A, B, P = 65521;
    unsigned long n, R, P1; // P1 = -P^(-1) mod R

    for (unsigned i=0; i<40; i++)
    {
        A = rand()*rand();
        B = rand()*rand();
        A %= P;
        B %= P;

        // 先求出正确值.
        {
            unsigned long long t = A;
            t *= B;
            t %= P;
            cout << "\n(A*B)  mod P : " << t << endl;
        }

        ////////////////

        // P 与 R必须互素, 鉴于P是素数, 就没有必要计算 GCD(P, R) == 1;
        // 2^(n-1) <= P <= 2^n;  R = 2^n;

        // P = 65521 ==  1111111111110001 (素数)
        // R = 65536 == 10000000000000000
        n = deg_num(P) + 1;
        R = 1;
        R <<= n;

        // P1 = "-" P^(-1) mod R; P1 是 P mod R的逆, 而且是负数.
        P1 = ModInv(P, R);

        {
            // t = (A * B * R^(-1)) mod P;
            unsigned long long t = Mont(A, B, P1, P, R);

            // 由于 t = (A * B * R^(-1)) mod P;
            // 所以 t 转换为 (A * B) mod P; (脱去蒙哥马利表示)
            // A * B mod P = t * R * R * R^(-1) mod P;
//            t *= GetendR(R, P);

            // ??? 乘以一个R不可以吗? 为什么要乘以R * R * R^(-1)
            // A * B mod P = t * R mod P;
            t *= R;

            t %= P;
            cout << "Mont() mod P : " << t << endl;
        }
    }

    return 0;
};

[课程]Android-CTF解题方法汇总！