1. 程式人生 > >Modular arithmetic and Montgomery form 實現快速模乘

Modular arithmetic and Montgomery form 實現快速模乘

ini a* targe -a init cat 數列 htm can

題目

電音之王


題解

求數列前n項相乘並取模


思路

①、這題的乘法是爆long long的,可以通過快速冪的思想去解決(按數位對其中的一個數進行剖分)。當然你的乘法會多出一個log的復雜度...

②、O(1)快速乘:一種O(1)復雜度求解整數相乘取模的思路(它對於64位的整型也是適用的):

  來自2009年國家集訓隊論文:駱可強:《論程序底層優化的一些方法與技巧》 (參考中附原文鏈接)

技術分享圖片
typedef long long ll;
#define MOL 123456789012345LL

inline ll mul_mod_ll(ll a,ll b)
{
    ll d 
= (ll)floor(a * (double)b / MOL + 0.5); ll ret = a * b - d * MOL; if(ret < 0) ret += MOL; return ret; }
View Code

③、正解dls一句話題解(當然是看不懂了...)

  參考中附一篇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()); } }
View Code

參考

論程序底層優化的一些方法與技巧

除算?剰余算の高速化

Modular arithmetic and Montgomery form 實現快速模乘