题引
题解
题解显而易见就是对数列的前k项进行前缀乘法。需要注意到的是乘法是会爆longlong的。
这时你的直观做法就是使用快速幂的思想规避溢出情况。
形如这样的代码:
ll mult(ll a,ll b) {
a %= mod;b %= mod;
ll ans = 0;
while(b) {
if(b&1) ans = (ans + a) % mod;
a = (a + a) % mod;
b >>= 1;
}
return ans;
}
时间复杂度 : O(logn)
的确我们解决了乘法的溢出情况。但是为此我们多了一个log。(所以说他是慢速乘法)
不幸的是对于本题O(nlogn)并不能过掉此题。常数过大了。
万幸的是存在一种O(1)的快速乘法算法,能在保证乘法结果正确的同时,降低时间复杂度。
引自2009国家集训队论文:
骆可强:《论程序底层优化的一些方法与技巧》
ll fast_mult(ll a,ll b) {
return (a*b - (ll)((long double)a/mod*b)*mod+mod)%mod;
}
上述也称O(1)快速乘。
但是本题的出题人DLS并不希望O(1)快速乘通过此题。也就是说本题的正解并不是上诉的做法。以上O(1)快速乘法仅仅只能勉强卡过本题。
据出题人DLS一句话题解 正确的做法为 : Montgomery modular multiplication
网上也没几篇像样的Blog,Wiki读一晚上也没感悟出来。只能先留着吧。
Update
一篇关于黑科技的Blog Montgomery modular multiplication
附上Dls代码:
#include <bits/stdc++.h>
using namespace std;
#define rep(i,a,n) for (int i=a;i<n;i++)
#define per(i,a,n) for (int i=n-1;i>=a;i--)
#define pb push_back
#define mp make_pair
#define all(x) (x).begin(),(x).end()
#define fi first
#define se second
#define SZ(x) ((int)(x).size())
typedef vector<int> VI;
typedef long long ll;
typedef pair<int,int> PII;
const ll mod=1000000007;
ll powmod(ll a,ll b) {ll res=1;a%=mod; assert(b>=0); for(;b;b>>=1){if(b&1)res=res*a%mod;a=a*a%mod;}return res;}
ll gcd(ll a,ll b) { return b?gcd(b,a%b):a;}
// head
typedef unsigned long long u64;
typedef __int128_t i128;
typedef __uint128_t u128;
int _,k;
u64 A0,A1,M0,M1,C,M;
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; rep(i,0,5) 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;
u64 pmod(u64 a,u64 b,u64 p) {
u64 d=(u64)floor(a*(long double)b/p+0.5);
ll ret=a*b-d*p;
if (ret<0) ret+=p;
return ret;
}
void bruteforce() {
u64 ans=1;
for (int i=0;i<=k;i++) {
ans=pmod(ans,A0,M);
u64 A2=pmod(M0,A1,M)+pmod(M1,A0,M)+C;
while (A2>=M) A2-=M;
A0=A1; A1=A2;
}
printf("%llu\n",ans);
}
int main() {
for (scanf("%d",&_);_;_--) {
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);
for (int i=0;i<=k;i++) {
ans=ans*a0;
a2=m0*a1+m1*a0+c;
a0=a1; a1=a2;
}
printf("%llu\n",ans.get());
}
}