牛客国庆Day5B(蒙哥马利算法)

题目链接：https://ac.nowcoder.com/acm/contest/205/B?&headNav=www

原来窝被这算法虐过么= =！

这个算法其实是用来解决RSA中计算公钥的问题的，是个底层常数优化算法。。

主要解决两个问题：

• 蒙哥马利约减，即$t{p}^{-1}(modm)$

• 蒙哥马利乘模，即 $xy(modm)$

这个算法的目的就是避免取模和除法来降低常数
算法过程，直接看杜教板子比看网上的博客来得简单。。

首先预处理几个参数：

b为数制，这里所有的数都以b进制的蒙哥马利表示法表示
$p=b^k$且$p\ge m$
$inv=m^{-1}modp$
$r2=p^{2}modm$

取b进制是为了方便进行位运算，所以b一般取2，然后p取$2^{64}$的话取模的时候自然溢出就行了，比较方便
然后关于蒙哥马利表示法，就是x的蒙哥马利表示为$xp(modm)$
所有数的运算都在此表示法下进行

蒙哥马利约减：

若$t\le m^2$，可以用$(t-(t\ast invmodp)\ast m)/p$ 代替 $t{p}^{-1}modm$
证明：
因为$t-t\ast inv\ast m\equiv t+t(-m^{-1}\ast m)\equiv 0(modp)$
所以$p|(t+t\ast inv\ast m)$，即$p|(t+(t\ast invmodp)\ast m)$
故有$(t+(t\ast invmodp)\ast m)/p\equiv t{p}^{-1}+(t\ast invmodp)p^{-1}m\equiv t{p}^{-1}(modm)$
在数值上，由于
$(t-(t\ast invmodp\ast m)/p)\le m^2/p\le m$
$(t-(t\ast invmodp\ast m)/p)\ge -m$
所以只需判断正负即可

然后就可以用蒙哥马利约减来求出一个数的蒙哥马利表示了，直接对$x\ast r2$进行约减即可
从蒙哥马利表示法中还原出原数也只需做一个约减就可以了

蒙哥马利乘模：

令
$\hat{x}=xp(modm)$
$\hat{y}=yp(modm)$
那么在蒙哥马利表示法中，$\widehat{xy}=xyp$
即$\widehat{xy}\equiv \hat{x}\hat{y}/p(modm)$
直接对$x\ast y$进行一次约减就可以了

基本就这么多。。说那么多其实只要有板子就可以。。
这里就贴杜教的板子，窝的板子实在是太丑了。。。

#include <bits/stdc++.h>
using namespace std;

typedef long long ll;
typedef unsigned long long u64;
typedef __int128_t i128;
typedef __uint128_t u128;

struct Mod64 {
    Mod64() :n_(0) {}
    Mod64(u64 n) :n_(init(n)) {}
    static u64 init(u64 w) { return reduce(u128(w) * r2); }
    static void set_mod(u64 m) {
        mod = m; assert(mod & 1);
        inv = m; for (int i = 0; i < 5; ++i) inv *= 2 - inv * m;
        r2 = -u128(m) % m;
    }
    static u64 reduce(u128 x) {
        u64 y = u64(x >> 64) - u64((u128(u64(x)*inv)*mod) >> 64);
        return ll(y)<0 ? y + mod : y;
    }
    Mod64& operator += (Mod64 rhs) { n_ += rhs.n_ - mod; if (ll(n_)<0) n_ += mod; return *this; }
    Mod64 operator + (Mod64 rhs) const { return Mod64(*this) += rhs; }
    Mod64& operator -= (Mod64 rhs) { n_ -= rhs.n_; if (ll(n_)<0) n_ += mod; return *this; }
    Mod64 operator - (Mod64 rhs) const { return Mod64(*this) -= rhs; }
    Mod64& operator *= (Mod64 rhs) { n_ = reduce(u128(n_)*rhs.n_); return *this; }
    Mod64 operator * (Mod64 rhs) const { return Mod64(*this) *= rhs; }
    u64 get() const { return reduce(n_); }
    static u64 mod, inv, r2;
    u64 n_;
};

u64 Mod64::mod, Mod64::inv, Mod64::r2;


int t, k;
u64 A0, A1, M0, M1, C, M;

void Run()
{
    scanf("%d", &t);
    while (t--)
    {
        scanf("%llu%llu%llu%llu%llu%llu%d", &A0, &A1, &M0, &M1, &C, &M, &k);
        Mod64::set_mod(M);
        Mod64 a0(A0), a1(A1), m0(M0), m1(M1), c(C), ans(1), a2(0);
        ans *= a0; ans *= a1;
        for (int i = 2; i <= k; ++i)
        {
            a2 = a1;
            a1 = m0 * a1 + m1 * a0 + c;
            a0 = a2;
            ans *= a1;
        }
        printf("%llu\n", ans.get());
    }
}

int main()
{
    Run();
    return 0;
}