题解 | #Bash Plays with Functions#

提供一篇不使用dp的解法。利用贝尔级数不难证明$f_0=\frac{1+x}{1-x}={\mu }^2\ast 1f\_0=\backslash frac\{1+x\}\{1-x\}=\backslash mu^2*1$，则$f_r={\mu }^2\ast \underset{r+1\mathrm{个}}{\underbrace{1\ast 1\ast \cdots \ast 1}}=\frac{1+x}{(1-x{)}^{r+1}}f\_r=\backslash mu^2*\backslash underbrace\{1*1*\backslash cdots\; *1\}\_\{r+1个\}=\backslash frac\{1+x\}\{(1-x)^\{r+1\}\}$。设$g(x)=\frac{1}{(1-x{)}^{r+1}}g(x)=\backslash frac\{1\}\{(1-x)^\{r+1\}\}$，那么有$g(x)={\prod }_{i=1}^s\left(\genfrac{}{}{0px}{}{r+k_i}{k_i}\right),x=p_1^{k_1}{p}_2^{k_2}\cdots p_s^{k_s}g(x)=\backslash prod\backslash limits\_\{i=1\}^\{s\}\backslash tbinom\{r+k\_i\}\{k\_i\},x=p\_1^\{k\_1\}p\_2^\{k\_2\}\backslash cdots\; p\_s^\{k\_s\}$，具体证明过程可以在oiwiki的普通生成函数板块找到，剩下的事情就只有求积性函数一类的工作了。

#if __has_include(<bits/stdc++.h>)
#include<bits/stdc++.h>
#endif
#include<iostream>
#include<bitset>
#include<map>
#include<vector>
#include<unordered_map>
#include<algorithm>
#include<sstream>
#include<stack>
#include<queue>
#include<cmath>
using namespace std;
typedef long long ll;
typedef long double ld;
typedef pair<ll,ll> pll;
template<class T,class U>
ostream& operator<<(ostream& os,const pair<T,U>& p){os<<p.first<<" "<<p.second;return os;}
#define debug(x) cerr<<#x<<":"<<x<<endl
#define debug2(x,y) cerr<<#x<<":"<<x<<"  "<<#y<<":"<<y<<endl
#define debug3(x,y,z) cerr<<#x<<":"<<x<<"  "<<#y<<":"<<y<<"  "<<#z<<":"<<z<<endl
#define debuga(a) for(int i=0;i<a.size();++i) debug2(i,a[i]);
#define debugar(a,b) for(int i=0;i<b;++i) debug2(i,a[i]);
using namespace std;
typedef long long ll;
const ll N=1e6+1000,N_=2e6+5000,M=1e9+7;
int cnt=0,prime[N],part[N],prod[N_];
int pfm[N][10],pfmsize[N];
int INV[N_];
void pre(){
    ios::sync_with_stdio(0);
    cin.tie(0);
    cout.tie(0);
    INV[0]=INV[1]=part[1]=prod[0]=prod[1]=1;
    for(int i=2;i<N;++i){
        prod[i]=(ll)i*prod[i-1]%M;
        INV[i]=(M-M/i)*INV[M%i]%M;
        if(!part[i]){
            prime[cnt++]=i;
            part[i]=i;
        }
        for(int j=0;i*prime[j]<N;++j){
            part[i*prime[j]]=i%prime[j]?prime[j]:part[i];
            if(i%prime[j]==0) break;
        }
    }
    for(ll i=N;i<N_;++i){
        prod[i]=i*prod[i-1]%M;
        INV[i]=(M-M/i)*INV[M%i]%M;
    }
    for(ll i=2;i<N_;++i) INV[i]=(ll)INV[i]*INV[i-1]%M;
    for(int i=2;i<N;++i){
        ll t=part[i],p=i,c=0;
        while(t!=1){
            p/=part[p];
            ++c;
            if(part[p]!=t){
                pfm[i][pfmsize[i]++]=c;
                t=part[p];
                c=0;
            }
        }
    }
}
ll comb(ll n,ll m){return (ll)prod[m]*INV[n]%M*INV[m-n]%M;}
void foreachp(ll r,ll n,ll ind,ll& res,ll t){
    if(ind==pfmsize[n]){
        res=res+t;
        if(res>=M) res%=M;
        return;
    }
    foreachp(r,n,ind+1,res,t*(comb(pfm[n][ind],r+pfm[n][ind])+comb(pfm[n][ind]-1,r-1+pfm[n][ind]))%M);
}
signed main(){
    pre();
    ll q,r,n;
    cin>>q;
    while(q--){
        cin>>r>>n;
        ll res=0;
        foreachp(r,n,0,res,1);
        cout<<res<<'\n';
    }
}