拉格朗日插值:
/// 注意mod，使用前须调用一次 polysum::init(int M); namespace polysum { #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--) typedef long long ll; const ll mod=1e9+7; /// 取模值 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;} const int D=101000; /// 最高次限制 ll a[D],f[D],g[D],p[D],p1[D],p2[D],b[D],h[D][2],C[D]; ll calcn(int d,ll *a,ll n) { if (n<=d) return a[n]; p1[0]=p2[0]=1; rep(i,0,d+1) { ll t=(n-i+mod)%mod; p1[i+1]=p1[i]*t%mod; } rep(i,0,d+1) { ll t=(n-d+i+mod)%mod; p2[i+1]=p2[i]*t%mod; } ll ans=0; rep(i,0,d+1) { ll t=g[i]*g[d-i]%mod*p1[i]%mod*p2[d-i]%mod*a[i]%mod; if ((d-i)&1) ans=(ans-t+mod)%mod; else ans=(ans+t)%mod; } return ans; } void init(int M) { /// M：最高次 f[0]=f[1]=g[0]=g[1]=1; rep(i,2,M+5) f[i]=f[i-1]*i%mod; g[M+4]=powmod(f[M+4],mod-2); per(i,1,M+4) g[i]=g[i+1]*(i+1)%mod; } ll polysum(ll n,ll *arr,ll m) { // a[0].. a[m] \sum_{i=0}^{n} a[i] for(int i = 0; i <= m; i++) a[i] = arr[i]; a[m+1]=calcn(m,a,m+1); rep(i,1,m+2) a[i]=(a[i-1]+a[i])%mod; return calcn(m+1,a,n); } ll qpolysum(ll R,ll n,ll *a,ll m) { // a[0].. a[m] \sum_{i=0}^{n-1} a[i]*R^i if (R==1) return polysum(n,a,m); a[m+1]=calcn(m,a,m+1); ll r=powmod(R,mod-2),p3=0,p4=0,c,ans; h[0][0]=0;h[0][1]=1; rep(i,1,m+2) { h[i][0]=(h[i-1][0]+a[i-1])*r%mod; h[i][1]=h[i-1][1]*r%mod; } rep(i,0,m+2) { ll t=g[i]*g[m+1-i]%mod; if (i&1) p3=((p3-h[i][0]*t)%mod+mod)%mod,p4=((p4-h[i][1]*t)%mod+mod)%mod; else p3=(p3+h[i][0]*t)%mod,p4=(p4+h[i][1]*t)%mod; } c=powmod(p4,mod-2)*(mod-p3)%mod; rep(i,0,m+2) h[i][0]=(h[i][0]+h[i][1]*c)%mod; rep(i,0,m+2) C[i]=h[i][0]; ans=(calcn(m,C,n)*powmod(R,n)-c)%mod; if (ans<0) ans+=mod; return ans; } }
例题:https://ac.nowcoder.com/acm/contest/139/F?&headNav=www
代码:
#pragma GCC optimize(2) #pragma GCC optimize(3) #pragma GCC optimize(4) #include<bits/stdc++.h> using namespace std; #define y1 y11 #define fi first #define se second #define pi acos(-1.0) #define LL long long #define ls rt<<1, l, m #define rs rt<<1|1, m+1, r //#define mp make_pair #define pb push_back #define ULL unsigned LL #define pll pair<LL, LL> #define pli pair<LL, int> #define pii pair<int, int> #define piii pair<pii, int> #define pdi pair<double, int> #define pdd pair<double, double> #define mem(a, b) memset(a, b, sizeof(a)) #define debug(x) cerr << #x << " = " << x << "\n"; #define fio ios::sync_with_stdio(false);cin.tie(0);cout.tie(0); //head const int N = 1e3 + 5; const int MOD = 1e9 + 7; int a[N], n; LL b[N]; /// 注意mod，使用前须调用一次 polysum::init(int M); namespace polysum { #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--) typedef long long ll; const ll mod=1e9+7; /// 取模值 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;} const int D=101000; /// 最高次限制 ll a[D],f[D],g[D],p[D],p1[D],p2[D],b[D],h[D][2],C[D]; ll calcn(int d,ll *a,ll n) { if (n<=d) return a[n]; p1[0]=p2[0]=1; rep(i,0,d+1) { ll t=(n-i+mod)%mod; p1[i+1]=p1[i]*t%mod; } rep(i,0,d+1) { ll t=(n-d+i+mod)%mod; p2[i+1]=p2[i]*t%mod; } ll ans=0; rep(i,0,d+1) { ll t=g[i]*g[d-i]%mod*p1[i]%mod*p2[d-i]%mod*a[i]%mod; if ((d-i)&1) ans=(ans-t+mod)%mod; else ans=(ans+t)%mod; } return ans; } void init(int M) { /// M：最高次 f[0]=f[1]=g[0]=g[1]=1; rep(i,2,M+5) f[i]=f[i-1]*i%mod; g[M+4]=powmod(f[M+4],mod-2); per(i,1,M+4) g[i]=g[i+1]*(i+1)%mod; } ll polysum(ll n,ll *arr,ll m) { // a[0].. a[m] \sum_{i=0}^{n-1} a[i] for(int i = 0; i <= m; i++) a[i] = arr[i]; a[m+1]=calcn(m,a,m+1); rep(i,1,m+2) a[i]=(a[i-1]+a[i])%mod; return calcn(m+1,a,n); } ll qpolysum(ll R,ll n,ll *a,ll m) { // a[0].. a[m] \sum_{i=0}^{n-1} a[i]*R^i if (R==1) return polysum(n,a,m); a[m+1]=calcn(m,a,m+1); ll r=powmod(R,mod-2),p3=0,p4=0,c,ans; h[0][0]=0;h[0][1]=1; rep(i,1,m+2) { h[i][0]=(h[i-1][0]+a[i-1])*r%mod; h[i][1]=h[i-1][1]*r%mod; } rep(i,0,m+2) { ll t=g[i]*g[m+1-i]%mod; if (i&1) p3=((p3-h[i][0]*t)%mod+mod)%mod,p4=((p4-h[i][1]*t)%mod+mod)%mod; else p3=(p3+h[i][0]*t)%mod,p4=(p4+h[i][1]*t)%mod; } c=powmod(p4,mod-2)*(mod-p3)%mod; rep(i,0,m+2) h[i][0]=(h[i][0]+h[i][1]*c)%mod; rep(i,0,m+2) C[i]=h[i][0]; ans=(calcn(m,C,n)*powmod(R,n)-c)%mod; if (ans<0) ans+=mod; return ans; } } int main() { polysum::init(N); while(~scanf("%d", &n)) { for (int i = 1; i <= n; ++i) scanf("%d", &a[i]); sort(a+1, a+1+n); LL pre = 1, ans = 0; for (int i = 1; i <= n; ++i) { b[0] = 0; for (int j = 1; j <= n-i+1; ++j) b[j] = (j*(polysum::powmod(j, n-i+1) - polysum::powmod(j-1, n-i+1))) % MOD; ans = (ans + pre*(polysum::polysum(a[i], b, n-i+1) - polysum::polysum(a[i-1], b, n-i+1))%MOD )%MOD; pre = (pre * a[i]) % MOD; } printf("%lld\n", (ans + MOD) % MOD); } return 0; }
杜教BM：
namespace linear_seq { const int N=10010; #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 mp make_pair #define all(x) (x).begin(),(x).end() #define SZ(x) ((int)(x).size()) typedef vector<ll> VI; typedef pair<ll,ll> PII; 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 res[N],base[N],_c[N],_md[N]; vector<ll> Md; void mul(ll *a,ll *b,int k) { rep(i,0,k+k) _c[i]=0; rep(i,0,k) if (a[i]) rep(j,0,k) _c[i+j]=(_c[i+j]+a[i]*b[j])%mod; for (int i=k+k-1;i>=k;i--) if (_c[i]) rep(j,0,SZ(Md)) _c[i-k+Md[j]]=(_c[i-k+Md[j]]-_c[i]*_md[Md[j]])%mod; rep(i,0,k) a[i]=_c[i]; } ll solve(ll n,VI a,VI b) { // a 系数 b 初值 b[n+1]=a[0]*b[n]+... // printf("%d\n",SZ(b)); ll ans=0,pnt=0; int k=SZ(a); assert(SZ(a)==SZ(b)); rep(i,0,k) _md[k-1-i]=-a[i];_md[k]=1; Md.clear(); rep(i,0,k) if (_md[i]!=0) Md.push_back(i); rep(i,0,k) res[i]=base[i]=0; res[0]=1; while ((1ll<<pnt)<=n) pnt++; for (int p=pnt;p>=0;p--) { mul(res,res,k); if ((n>>p)&1) { for (int i=k-1;i>=0;i--) res[i+1]=res[i];res[0]=0; rep(j,0,SZ(Md)) res[Md[j]]=(res[Md[j]]-res[k]*_md[Md[j]])%mod; } } rep(i,0,k) ans=(ans+res[i]*b[i])%mod; if (ans<0) ans+=mod; return ans; } VI BM(VI s) { VI C(1,1),B(1,1); int L=0,m=1,b=1; rep(n,0,SZ(s)) { ll d=0; rep(i,0,L+1) d=(d+(ll)C[i]*s[n-i])%mod; if (d==0) ++m; else if (2*L<=n) { VI T=C; ll c=mod-d*powmod(b,mod-2)%mod; while (SZ(C)<SZ(B)+m) C.pb(0); rep(i,0,SZ(B)) C[i+m]=(C[i+m]+c*B[i])%mod; L=n+1-L; B=T; b=d; m=1; } else { ll c=mod-d*powmod(b,mod-2)%mod; while (SZ(C)<SZ(B)+m) C.pb(0); rep(i,0,SZ(B)) C[i+m]=(C[i+m]+c*B[i])%mod; ++m; } } return C; } ll gao(VI a,ll n) { VI c=BM(a); c.erase(c.begin()); rep(i,0,SZ(c)) c[i]=(mod-c[i])%mod; return solve(n,c,VI(a.begin(),a.begin()+SZ(c))); } };
有递推式直接输入初始值和系数，不然打表打出前几项丢进去跑
代码：
#pragma GCC optimize(2) #pragma GCC optimize(3) #pragma GCC optimize(4) #include<bits/stdc++.h> using namespace std; #define y1 y11 #define fi first #define se second #define pi acos(-1.0) #define LL long long #define ll long long #define ls rt<<1, l, m #define rs rt<<1|1, m+1, r //#define mp make_pair #define pb push_back #define ULL unsigned LL #define pll pair<LL, LL> #define pli pair<LL, int> #define plii pair<LL,pii> #define pii pair<int, int> #define piii pair<pii, pii> #define pdi pair<double, int> #define pdd pair<double, double> #define mem(a, b) memset(a, b, sizeof(a)) #define debug(x) cerr << #x << " = " << x << "\n"; #define fio ios::sync_with_stdio(false);cin.tie(0);cout.tie(0); //head const int mod = 1e9 + 7; namespace linear_seq { const int N=10010; #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 mp make_pair #define all(x) (x).begin(),(x).end() #define SZ(x) ((int)(x).size()) typedef vector<ll> VI; typedef pair<ll,ll> PII; 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 res[N],base[N],_c[N],_md[N]; vector<ll> Md; void mul(ll *a,ll *b,int k) { rep(i,0,k+k) _c[i]=0; rep(i,0,k) if (a[i]) rep(j,0,k) _c[i+j]=(_c[i+j]+a[i]*b[j])%mod; for (int i=k+k-1;i>=k;i--) if (_c[i]) rep(j,0,SZ(Md)) _c[i-k+Md[j]]=(_c[i-k+Md[j]]-_c[i]*_md[Md[j]])%mod; rep(i,0,k) a[i]=_c[i]; } ll solve(ll n,VI a,VI b) { /// a 系数 b 初值 b[n+1]=a[0]*b[n]+... /// printf("%d\n",SZ(b)); ll ans=0,pnt=0; int k=SZ(a); assert(SZ(a)==SZ(b)); rep(i,0,k) _md[k-1-i]=-a[i];_md[k]=1; Md.clear(); rep(i,0,k) if (_md[i]!=0) Md.push_back(i); rep(i,0,k) res[i]=base[i]=0; res[0]=1; while ((1ll<<pnt)<=n) pnt++; for (int p=pnt;p>=0;p--) { mul(res,res,k); if ((n>>p)&1) { for (int i=k-1;i>=0;i--) res[i+1]=res[i];res[0]=0; rep(j,0,SZ(Md)) res[Md[j]]=(res[Md[j]]-res[k]*_md[Md[j]])%mod; } } rep(i,0,k) ans=(ans+res[i]*b[i])%mod; if (ans<0) ans+=mod; return ans; } VI BM(VI s) { VI C(1,1),B(1,1); int L=0,m=1,b=1; rep(n,0,SZ(s)) { ll d=0; rep(i,0,L+1) d=(d+(ll)C[i]*s[n-i])%mod; if (d==0) ++m; else if (2*L<=n) { VI T=C; ll c=mod-d*powmod(b,mod-2)%mod; while (SZ(C)<SZ(B)+m) C.pb(0); rep(i,0,SZ(B)) C[i+m]=(C[i+m]+c*B[i])%mod; L=n+1-L; B=T; b=d; m=1; } else { ll c=mod-d*powmod(b,mod-2)%mod; while (SZ(C)<SZ(B)+m) C.pb(0); rep(i,0,SZ(B)) C[i+m]=(C[i+m]+c*B[i])%mod; ++m; } } return C; } ll gao(VI a,ll n) { VI c=BM(a); c.erase(c.begin()); rep(i,0,SZ(c)) c[i]=(mod-c[i])%mod; return solve(n,c,VI(a.begin(),a.begin()+SZ(c))); } }; LL n, m; int main() { scanf("%lld %d", &n, &m); vector<LL> a(m, 0); a[0] = a[m-1] = 1; vector<LL> b(m, 1); printf("%lld\n", linear_seq::solve(n, a, b)); return 0; }
双hash:
typedef pair<int,int> hashv; const LL mod1=1000000007; const LL mod2=1000000009; hashv operator + (hashv a,hashv b) { int c1=a.fi+b.fi,c2=a.se+b.se; if (c1>=mod1) c1-=mod1; if (c2>=mod2) c2-=mod2; return {c1,c2}; } hashv operator - (hashv a,hashv b) { int c1=a.fi-b.fi,c2=a.se-b.se; if (c1<0) c1+=mod1; if (c2<0) c2+=mod2; return {c1,c2}; } hashv operator * (hashv a,hashv b) { return {1LL*a.fi*b.fi%mod1,1LL*a.se*b.se%mod2}; } hashv h[N], p[N]; hashv get_hash(int l, int r) { return h[r]-h[l-1]*p[r-l+1]; }
代码：
#pragma GCC optimize(2) #pragma GCC optimize(3) #pragma GCC optimize(4) #include<bits/stdc++.h> using namespace std; #define y1 y11 #define fi first #define se second #define pi acos(-1.0) #define LL long long //#define mp make_pair #define pb push_back #define ls rt<<1, l, m #define rs rt<<1|1, m+1, r #define ULL unsigned LL #define pll pair<LL, LL> #define pli pair<LL, int> #define pii pair<int, int> #define piii pair<int, pii> #define puu pair<ULL, ULL> #define pdd pair<long double, long double> #define mem(a, b) memset(a, b, sizeof(a)) #define fio ios::sync_with_stdio(false);cin.tie(0);cout.tie(0); #define fopen freopen("in.txt", "r", stdin);freopen("out.txt", "w", stout); //head const int N = 1e6 + 5; typedef pair<int, int> hashv; const LL mod1=1000000007; const LL mod2=1000000009; hashv operator + (hashv a,hashv b) { int c1=a.fi+b.fi,c2=a.se+b.se; if (c1>=mod1) c1-=mod1; if (c2>=mod2) c2-=mod2; return {c1,c2}; } hashv operator - (hashv a,hashv b) { int c1=a.fi-b.fi,c2=a.se-b.se; if (c1<0) c1+=mod1; if (c2<0) c2+=mod2; return {c1,c2}; } hashv operator * (hashv a,hashv b) { return {1LL*a.fi*b.fi%mod1,1LL*a.se*b.se%mod2}; } hashv h[N], p[N], base = {31, 47}; hashv get_hash(int l, int r) { return h[r]-h[l-1]*p[r-l+1]; } char s[N], t[N]; int n, m, c0 = 0, c1 = 0; int main() { scanf("%s", s+1); scanf("%s", t+1); n = strlen(s+1); m = strlen(t+1); for (int i = 1; i <= n; ++i) if(s[i] == '0') c0++; else c1++; p[0] = {1, 1}; h[0] = {0, 0}; for (int i = 1; i <= m; ++i) h[i] = h[i-1]*base + hashv{t[i]-'a'+1, t[i]-'a'+1}, p[i] = p[i-1]*base; int ans = 0; if(c0 == 0) { if(m%c1 == 0) { hashv h = {0, 0}; bool f = true; for (int i = c1; i <= m; i += c1) { if(h.fi == 0) h = get_hash(i-c1+1, i); else { if(get_hash(i-c1+1, i) != h) { f = false; break; } } } if(f) ans++; } } else if(c1 == 0) { if(m%c0 == 0) { hashv h = {0, 0}; bool f = true; for (int i = c0; i <= m; i += c0) { if(h.fi == 0) h = get_hash(i-c0+1, i); else { if(get_hash(i-c0+1, i) != h) { f = false; break; } } } if(f) ans++; } } else { for (int i = 1; i*c1 < m; ++i) { if((m-i*c1)%c0 == 0) { bool f = true; int l1 = i, l0 = (m-i*c1)/c0; hashv h0 = {0, 0}, h1 = {0, 0}; int now = 1; for (int i = 1; i <= n; ++i) { if(s[i] == '0') { if(h0.fi == 0) h0 = get_hash(now, now+l0-1); else { if(get_hash(now, now+l0-1) != h0) { f = false; break; } } now += l0; } else { if(h1.fi == 0) h1 = get_hash(now, now+l1-1); else { if(get_hash(now, now+l1-1) != h1) { f = false; break; } } now += l1; } } if(f && h0 != h1) ans++; } } } printf("%d\n", ans); return 0; }
Montgomery Modular Multiplication:
#define rep(i,a,n) for (int i=a;i<n;i++) 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; 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;
参考资料min25博客： https://min-25.hatenablog.com/entry/2017/08/20/171214
例题：电音之王
代码：
#pragma GCC optimize(2) #pragma GCC optimize(3) #pragma GCC optimize(4) #include<bits/stdc++.h> using namespace std; #define y1 y11 #define fi first #define se second #define pi acos(-1.0) #define LL long long //#define mp make_pair #define pb push_back #define ls rt<<1, l, m #define rs rt<<1|1, m+1, r #define ULL unsigned LL #define pll pair<LL, LL> #define pli pair<LL, int> #define pii pair<int, int> #define piii pair<int, pii> #define puu pair<ULL, ULL> #define pdd pair<long double, long double> #define mem(a, b) memset(a, b, sizeof(a)) #define fio ios::sync_with_stdio(false);cin.tie(0);cout.tie(0); #define fopen freopen("in.txt", "r", stdin);freopen("out.txt", "w", stout); //head #define rep(i,a,n) for (int i=a;i<n;i++) 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; 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; int T, k; ULL A0, A1, M0, M1, C, M; int main() { 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), res(1), tmp(0); res = res*a0; res = res*a1; for (int i = 2; i <= k; ++i) { tmp = m0*a1+m1*a0+c; res = res*tmp; a0 = a1; a1 = tmp; } printf("%llu\n", res.get()); } return 0; }
565
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
