类欧几里得算法
注意下面有(hen)的(duo)地方i和x打混了，懒得改了，应该挺容易看懂（吧）
以下的除法均为整除，$\lambda$表示（可以直接求出的）常量。
我们欲求的是$f(k_1,k_2,a,b,c,n)=\sum_{i=0}^n x^{k_1} ((ax+b)/c)^{k_2}$
如果$a \geq c$或$b \geq c$，那么
$f(k_1,k_2,a,b,c,n)=\sum_{i=0}^n x^{k_1} ((a/c)x+(b/c)+(a \bmod c \times x+b \bmod c)/c)^{k2}$
把后面那坨玩意儿二项式展开一下就是形如若干项$\lambda\sum_{i=0}^n x^{k_1} ((a\bmod c\times x+b \bmod c)/c)^{k_2}$之和的形式，求出$f(k_1,k_2,a \bmod c,b \bmod c,c,n)$即可。
如果a=0或k2=0或an+b<c，那么就是要求$\lambda\sum_{i=0}^n x^{k_1}$，插值一下就好了。
否则的话，a<c且b<c，$0^{k_2}=0$。
考虑进行精彩变换，把$x^{k_2}$变换成$\sum_{i=0}^{x-1} ((i+1)^{k_2}-i^{k_2})$。
令$m=(an+b)/c>0$，那么
$f(k_1,k_2,a,b,c,n)=\sum_{t=0}^{m-1} ((t+1)^{k_2}-t^{k_2}) \sum_{x=0}^n x^{k_1} [(ax+b)/c \geq t+1]$
$=\sum_{t=0}^{m-1} ((t+1)^{k_2}-t^{k_2}) \sum_{x=0}^n x^{k_1} [x> (tc+c-b-1)/a ]$
$=\sum_{t=0}^{m-1} ((t+1)^{k_2}-t^{k_2}) \sum_{x=0}^n x^{k_1}-\sum_{t=0}^{m-1} ((t+1)^{k_2}-t^{k_2})\sum_{g=0}^{(tc+c-b-1)/a} g^{k1}$
注意到$(t+1)^{k_2}-t^{k_2}$和$\sum_{g=0}^x g^{k_1}$均为多项式，假设分别为A和B。
$=(\sum_{t=0}^{m-1} ((t+1)^{k_2}-t^{k_2}) )(\sum_{x=0}^n x^{k_1})-\sum_{t=0}^{m-1} (\sum_{p=0}^{k_2} A_pt^p)(\sum_{q=0}^{k_1} B_q ((tc+c-b-1)/a)^q)$
前一部分$\sum_{t=0}^{m-1}((t+1)^{k_2}-t^{k_2})$和$\sum_{x=0}^n x^{k_1}$可以插值，后面一坨把每一项展开，就是形如$\lambda\sum_{t=0}^{m-1} t^p ((tc+c-b-1)/a)^q$这样，求$f(p,q,c,c-b-1,a,m-1)$即可。
实现的时候可以直接实现一个函数$g(a,b,c,n)$，返回所有$f(k_1,k_2,a,b,c,n)$的值。
观察到$(a,c) \rightarrow (a\bmod c,c) \rightarrow (c,a \bmod c)$，所以递归次数是$O(\log(c))$的。
#include <iostream> #include <stdio.h> #include <math.h> #include <string.h> #include <time.h> #include <stdlib.h> #include <string> #include <bitset> #include <vector> #include <set> #include <map> #include <queue> #include <algorithm> #include <sstream> #include <stack> #include <iomanip> using namespace std; #define pb push_back #define mp make_pair typedef pair<int,int> pii; typedef long long ll; typedef double ld; typedef vector<int> vi; #define fi first #define se second #define fe first #define FO(x) {freopen(#x".in","r",stdin);freopen(#x".out","w",stdout);} #define Edg int M=0,fst[SZ],vb[SZ],nxt[SZ];void ad_de(int a,int b){++M;nxt[M]=fst[a];fst[a]=M;vb[M]=b;}void adde(int a,int b){ad_de(a,b);ad_de(b,a);} #define Edgc int M=0,fst[SZ],vb[SZ],nxt[SZ],vc[SZ];void ad_de(int a,int b,int c){++M;nxt[M]=fst[a];fst[a]=M;vb[M]=b;vc[M]=c;}void adde(int a,int b,int c){ad_de(a,b,c);ad_de(b,a,c);} #define es(x,e) (int e=fst[x];e;e=nxt[e]) #define esb(x,e,b) (int e=fst[x],b=vb[e];e;e=nxt[e],b=vb[e]) typedef long long ll; const int MOD=1e9+7; inline ll qp(ll a,ll b) { ll x=1; a%=MOD; while(b) { if(b&1) x=x*a%MOD; a=a*a%MOD; b>>=1; } return x; } namespace Lagrange { ll x[23333],y[23333],a[23333],g[23333],h[23333],p[23333]; int N; void work() { for(int i=0;i<N;++i) a[i]=0; g[0]=1; for(int i=0;i<N;++i) { for(int _=0;_<=i;++_) h[_+1]=g[_]; h[0]=0; for(int _=0;_<=i;++_) h[_]=(h[_]-g[_]*(ll)x[i])%MOD; for(int _=0;_<=i+1;++_) g[_]=h[_]; } for(int i=0;i<N;++i) { for(int j=0;j<=N;++j) p[j]=g[j]; for(int j=N;j;--j) p[j-1]=(p[j-1]+p[j]*(ll)x[i])%MOD; ll s=1; for(int j=0;j<N;++j) if(i!=j) s=s*(x[i]-x[j])%MOD; s=y[i]*qp(s,MOD-2)%MOD; for(int _=0;_<N;++_) a[_]=(a[_]+p[_+1]*s)%MOD; } } vector<int> feed(vector<int> v) { N=v.size(); for(int i=0;i<N;++i) x[i]=i,y[i]=v[i]; work(); v.clear(); for(int i=0;i<N;++i) v.pb(a[i]); while(v.size()&&!v.back()) v.pop_back(); return v; } ll calc(vector<int>&v,ll xx) { ll s=0,gg=1; xx%=MOD; for(int i=0;i<N;++i) s=(s+gg*v[i])%MOD,gg=gg*xx%MOD; return s; } } using Lagrange::feed; using Lagrange::calc; //ps[k]=\sum_{i=0}^x i^k vector<int> ps[2333]; //rs[k]=\sum_{i=0}^x ((i+1)^k-i^k) vector<int> rs[2333]; struct arr{ll p[11][11];}; ll C[233][233]; arr calc(ll a,ll b,ll c,ll n) { arr w; if(n==0) a=0; if(a==0||a*n+b<c) { for(int i=0;i<=10;++i) { ll t=calc(ps[i],n),s=b/c; for(int j=0;i+j<=10;++j) w.p[i][j]=t,t=t*s%MOD; } return w; } for(int i=0;i<=10;++i) w.p[i][0]=calc(ps[i],n); if(a>=c||b>=c) { arr t=calc(a%c,b%c,c,n); ll p=a/c,q=b/c; for(int i=0;i<=10;++i) for(int j=1;i+j<=10;++j) { ll s=0,px=1; for(int x=0;x<=j;++x,px=px*p%MOD) { ll qy=1; for(int y=0;x+y<=j;++y,qy=qy*q%MOD) { //x^(i) (px)^x q^y ??^(j-x-y) s+=px*qy%MOD*C[j][x]%MOD*C[j-x][y] %MOD*t.p[i+x][j-x-y]; s%=MOD; } } w.p[i][j]=s; } return w; } ll m=(a*n+b)/c; arr t=calc(c,c-b-1,a,m-1); for(int i=0;i<=10;++i) for(int j=1;i+j<=10;++j) { ll s=calc(rs[j],m-1)*calc(ps[i],n)%MOD; for(int p=0;p<j;++p) { for(unsigned q=0;q<ps[i].size();++q) { ll v=C[j][p]*ps[i][q]%MOD; //v*t^p*((tc+c-b-1)/a)^q s-=t.p[p][q]*v; s%=MOD; } } w.p[i][j]=s%MOD; } return w; } int T,n,a,b,c,k1,k2; int main() { for(int i=0;i<=230;++i) { C[i][0]=1; for(int j=1;j<=i;++j) C[i][j]=(C[i-1][j-1]+C[i-1][j])%MOD; } for(int i=0;i<=10;++i) { ll sp=0,sr=0; vector<int> p,r; for(int j=0;j<=20;++j) sp+=qp(j,i),sr+=qp(j+1,i)-qp(j,i), sp%=MOD,sr%=MOD,p.pb(sp),r.pb(sr); ps[i]=feed(p); rs[i]=feed(r); } scanf("%d",&T); while(T--) { scanf("%d%d%d%d%d%d", &n,&a,&b,&c,&k1,&k2); arr s=calc(a,b,c,n); int p=s.p[k1][k2]; p=(p%MOD+MOD)%MOD; printf("%d\n",p); } }
