类欧几里得算法

https://loj.ac/problem/138

注意下面有(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 mod c imes x+b mod c)/c)^{k2}$

把后面那坨玩意儿二项式展开一下就是形如若干项$lambdasum_{i=0}^n x^{k_1} ((amod c imes x+b mod c)/c)^{k_2}$之和的形式，求出$f(k_1,k_2,a mod c,b mod c,c,n)$即可。

如果a=0或k2=0或an+b<c，那么就是要求$lambdasum_{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}$可以插值，后面一坨把每一项展开，就是形如$lambdasum_{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)  ightarrow (amod c,c) ightarrow (c,a mod 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
",p);
    }
}