现在定义函数(F_m(n))表示将(n)表示为若干(m)的非负整数次幂的和的方案数
定义(G_m^k(n))为(k)个(F_m(n))卷积起来的结果，现给定(n,m,k)，求(sumlimits_{i=0}^n G_m^k(i))，答案对(10^9+7)取模
数据范围：(0<=n<=10^{18},2<=m<=10^{18},1<=k<=20)
Solution
神仙题qwq
考虑将这个问题形象化一点，说白了就是用(1,m,m^2,m^3...m^w)去凑(n)，并且要求一个方案中前面的数不能比后面的数大，求方案数
然后卷积起来的话就是。。每个数都有(K)个（为了与下面的变量区分开来，题目中的(k)统一写成大写），但是每个数中的(k)个是不一样的，可以算作不同的方案
我们先考虑比较简单的(k=1)的情况：
总共有(log)个可以用的数字(1,m,m^2,...,m^w)，我们将这些数按照从小到大的顺序存到一个数组(a)里面去，我们令(f[i][j])表示用了(a)中的前(i)个数字，当前的和为(j)的方案数，这个状态显然不够优秀因为第二维太大了
考虑这样一个事情：假设我们现在只要求一个位置的值，比如说只用求(F(n))（(m)那个下标就不写了，默认就是(m)，后面关于(g)函数的表示同），那么很明显并不是所有的(F(x))都要用到的，我们可以只计算我们需要用到的位置的值
注意到在转移的时候，我们从小到大枚举(i)，(a_i)一定是(a_{i-1})的倍数，我们将第二维(j)写成(k*a_{i-1}+r)的形式（其中(r<a_i)，其实也就是商和余数这样的形式），我们将(j)放在模(a_{i-1})意义下看，会发现在转移的时候，因为加进来的(a_i)一定是(a_{i-1})的倍数，所以模(a_{i-1})的余数(r)是不会变的，也就是说如果说我们最后要求(F(n))，那么中间会用到的(j)一定是形如(k*a^{i-1}+(n\%a_{i-1}))
所以我们可以将dp数组的第二维设成(a_{i-1})前面的那个系数，记这个新的dp数组为(h[i][j])，表示可用前(i)个(a)数组中的值，当前的和为(j*a_{i-1}+(n\%a_{i-1}))的方案数，那么有转移：
然而现在的问题是。。这个第二维还是太大了
然后这个时候我们有一个结论（我不会证），当(i)固定的时候这个(h[i])是一个关于第二维的(i)次多项式
那就。。暴力维护前(i+1)个值，第二维比较大的情况就直接插值就好了
但是现在的问题是，这样我们求出来的是单点的值，答案要我们求前缀和
这里又有一个很神的操作，我们给(a)数组里面多加一个(1)，然后就相当于给这个多项式乘上了一个((1+x+x^2+x^3+...))，然后这个时候得到的就是前缀和了（具体实现上只要将(h[0])初始化为每一位都是(1)就行）
复杂度是(O(log_mn(nlog_mn+(log_mn)^2))的，其中((log_mn)^2)那项是拉格朗日插值的复杂度，还是有点爆炸，这个时候我们可以选择用线性插值但是我不会qwq
注意到这里插值的点((x,y))中的(x)是连续的（从(0)到(i)），所以我们预处理一下就可以做到(O(log_mn))插值查询了（其实也不能说是把多项式的系数插了出来。。只是能算值而已不过那样在这题里面就够了）
代码大概长这样
#include<iostream>
#include<cstdio>
#include<cstring>
#include<vector>
#define ll long long
using namespace std;
const int N=2010,MOD=1e9+7;
ll a[N],pw[N];
ll fac[N],invfac[N];
int K,ans,lg,cnt;
ll m,n;
int plu(int x,int y){return (1LL*x+y)%MOD;}
int mul(int x,int y){return 1LL*x*y%MOD;}
int ksm(int x,int y){
int ret=1,base=x;
for (;y;y>>=1,base=mul(base,base))
if (y&1) ret=mul(ret,base);
return ret;
}
struct T1{/*{{{*/
ll X[N],Y[N];
ll pre[N],suf[N],fm[N];
int n;
void calc(){
int tmp;
for (int i=1;i<=n;++i){
tmp=mul(invfac[i-1],invfac[n-i]);
tmp=mul(tmp,Y[i]);
if ((n-i)%2==1)
fm[i]=(MOD-tmp)%MOD;
else
fm[i]=tmp;
}
}
int get_val(ll x){
if (x+1<=n) return Y[x+1];
x%=MOD;
int ret=0,tmp=1;
pre[0]=1; suf[n+1]=1;
for (int i=1;i<=n;++i) pre[i]=mul(pre[i-1],plu(x,MOD-X[i]));
for (int i=n;i>=1;--i) suf[i]=mul(suf[i+1],plu(x,MOD-X[i]));
for (int i=1;i<=n;++i)
ret=plu(ret,mul(fm[i],mul(pre[i-1],suf[i+1])));
return (ret+MOD)%MOD;
}
}h[N];/*}}}*/
void prework(){
pw[0]=1; lg=0;
while (pw[lg]<=n/m) pw[lg+1]=pw[lg]*m,++lg;
cnt=0;
for (int i=0;i<=lg;++i)
for (int j=1;j<=K;++j)
a[++cnt]=pw[i];
a[0]=1;
fac[0]=1;
for (int i=1;i<=cnt;++i) fac[i]=mul(fac[i-1],i);
invfac[cnt]=ksm(fac[cnt],MOD-2);
for (int i=cnt-1;i>=0;--i) invfac[i]=mul(invfac[i+1],i+1);
}
void dp(){
h[0].X[1]=0; h[0].Y[1]=1; h[0].n=1; h[0].calc();
ll r;
for (int i=1;i<=cnt;++i){
h[i].n=i+1;
r=n%a[i]/a[i-1];
for (int j=0;j<=i;++j){
h[i].X[j+1]=j;
h[i].Y[j+1]=plu(h[i].Y[j],h[i-1].get_val(a[i]/a[i-1]*j+r));
}
h[i].calc();
}
printf("%d
",h[cnt].get_val(n/a[cnt]));
}
int main(){
#ifndef ONLINE_JUDGE
freopen("a.in","r",stdin);
#endif
scanf("%lld%lld%d",&n,&m,&K);
ans=0;
prework();
dp();
}