1 条题解

  • 0
    @ 2025-8-24 23:15:48

    自动搬运

    查看原文

    来自洛谷，原作者为

    avatar Naszt
    数论，我的挚爱 | 看个人介绍请删除 "display: none;" 元素

    搬运于2025-08-24 23:15:48，当前版本为作者最后更新于2025-03-06 20:44:23，作者可能在搬运后再次修改，您可在原文处查看最新版

    自动搬运只会搬运当前题目点赞数最高的题解，您可前往洛谷题解查看更多

    以下是正文

    互质与整除

    观前提醒

    原本这题是作为 D 题的，不过多元多项式优化其实优化并不大，代码却极其难写，有点故意套模板的嫌疑。

    所以本题就仅仅需要其思想，时间复杂度是 O(d3(n)+σ0(n)log⁡n)O(d_3(n)+\sigma_0(n)\log n) 的。
    加上 EI 的多元多项式优化可以做到 $\Theta(s\sigma_0(n)\log\sigma_0(n)+\sigma_0(n)\log n)$。
    ※ σ0(n):=∑ij=n1\sigma_0(n):=\sum_{ij=n}1，为因子个数函数。
    ※ d3(n):=∑ijk=n1d_3(n):=\sum_{ijk=n}1，为因子的因子个数函数。

    但 d3(n)d_3(n) 其实跑不满，实际上期望可能是 d3(n)log⁡n\frac{d_3(n)}{\log n} 的，而且常数相对较低。
    当然，后面也是会讲一下多元多项式优化做法的。

    思路分析

    先表示出答案函数 a(n)a(n)

    定义 ψ(n)\psi(n)：方程 φ(x)=n\varphi(x)=n 的 xx 的解数。
    答案 a(n)a(n) 即为 ∑d∣nψ(d)\sum_{d\mid n}\psi(d)。

    定义 ψ(n)\psi(n) 的 Dirichlet 生成函数：

    Ψ~(x)=∑n≥1ψ(n)nx\tilde{\Psi}(x)=\sum_{n\ge1}\frac{\psi(n)}{n^x}

    由于每有一个 φ(x)=n\varphi(x)=n 就对 ψ(n)\psi(n) 产生 11 的贡献，
    由此我们可以化简 Ψ~(x)\tilde{\Psi}(x)：

    $$\begin{aligned} \tilde{\Psi}(x)&=\sum_{n\ge1}\frac{\psi(n)}{n^x}=\sum_{n\ge1}\frac{1}{\varphi^x(n)}\\ &=\prod_{p\in\mathcal{P}}\left(1+\sum_{i\ge1}\frac{1}{\varphi^x(p^i)}\right)\\ &=\prod_{p\in\mathcal{P}}(1+(p-1)^{-x}+(p^2-p)^{-x}+(p^3-p^2)^{-x}+\dots)\\ &=\prod_{p\in\mathcal{P}}(1+(p-1)^{-x}(1+p^{-x}+p^{-2x}+\dots))\\ &=\prod_{p\in\mathcal{P}}\left(1+\frac{(p-1)^{-x}}{1-p^{-x}}\right)\\ &=\zeta(x)\prod_{p\in\mathcal{P}}(1+(p-1)^{-x}-p^{-x})\\ \end{aligned} $$

    后面的乘积难以化简，于是我们将其单独提出。 即为定义 f(n)f(n) 的 Dirichlet 生成函数：

    $$\tilde{F}(x)=\sum_{n\ge1}\frac{f(n)}{n^x}=\prod_{p\in\mathcal{P}}(1+(p-1)^{-x}-p^{-x}) $$

    那么答案函数 a(n)a(n) 的 Dirichlet 生成函数为：

    $$\tilde{\Alpha}(x)=\zeta(x)\tilde{\Psi}(x)=\zeta^2(x)\tilde{F}(x) $$

    又因为 ζ2(x)\zeta^2(x) 是 σ0(n)\sigma_0(n) 的 Dirichlet 生成函数，
    所以答案函数：

    $$a(n)=\sum_{d\mid n}f(d)\sigma_0\left(\frac{n}{d}\right) $$

    现在考虑快速计算 f(d)f(d)

    $$f(d)=[d^{-x}]\prod_{p\in\mathcal{P}}(1+(p-1)^{-x}-p^{-x}) $$

    考虑直接计算，其值域为 σ0(n)\sigma_0(n)，共做 O(σ0(n))O(\sigma_0(n)) 次狄利克雷卷积。
    每次狄利克雷卷积可以直接将其看做背包问题 dp，时间复杂度为其值域。

    那么时间复杂度为 O(σ02(n))O(\sigma_0^2(n))，当然也是跑不满的。

    一种很自然的想法是将其转为多项式，但这并不是多项式。

    我们模仿贝尔级数的思想，定义 xi=pi−xx_i=p_i^{-x}，其中 pip_i 为第 ii 小的质数。
    对任意一项因式分解就可以将其转为单项式：

    $$n^{-x}=\left(\prod_{i\ge1} p_i^{t_{i,n}}\right)^{-x}=\prod_{i\ge1} x_i^{t_{i,n}} $$

    于是可以把 Dirichlet 生成函数转成多元普通生成函数：

    $$f(d)=\left[\prod_{i\ge1} x_i^{t_{i,d}}\right]\prod_{p\in\mathcal{P}}\left(1+\prod_{i\ge1} x_i^{t_{i,p-1}}-\prod_{i\ge1} x_i^{t_{i,p}}\right) $$

    考虑 σ0(n)\sigma_0(n) 次暴力卷积，不过其实每次只会更新指数比它高的项，
    那么每次卷积的时间复杂度达不到 σ0(n)\sigma_0(n)。

    时间复杂度实际上为 ∑d∣nσ0(d)=d3(n)\sum_{d|n}\sigma_0(d)=d_3(n)，换句话就是因子的因子个数。

    代码实现

    以下代码有个 ss 的常数，你可以通过状压等方法去掉，不过没什么必要。

    出题人代码

    #include<bits/stdc++.h>
#define il inline __attribute__((__always_inline__))

typedef unsigned int i4;
typedef unsigned long long i8;

namespace Miller_Rabin{
  const i4 Pcnt=7;
  const i8 p[Pcnt]={2,325,9375,28178,450775,9780504,1795265022};
  il i8 mul(i8 a,i8 b,i8 p){
    return __int128(a)*b%p;
  }
  il i8 pow(i8 a,i8 b,i8 p){
    i8 ans=1;
    for(;b;a=mul(a,a,p),b>>=1)if(b&1)ans=mul(ans,a,p);
    return ans;
  }
  il bool prime(i8 x){
    if(x<3||x%2==0)return x==2;
    i8 u=x-1,t=0;
    while(u%2==0) u/=2,++t;
    i8 ud[]={2,325,9375,28178,450775,9780504,1795265022};
    for(i8 a:ud){
      i8 v=pow(a,u,x);
      if(v==1||v==x-1||v==0) continue;
      for(i4 j=1;j<=t;j++){
        v=mul(v,v,x);
        if(v==x-1&&j!=t){v=1;break;}
        if(v==1) return 0;
      }
      if(v!=1) return 0;
    }
    return 1;
  }
}using Miller_Rabin::prime;

const i4 MAXS=16,MAXN=104000,MAXA=64;
const i8 P=998244353;

i8 p[MAXA];i4 T,s,a[MAXA];
struct Vec{
  char t[MAXS];
  i4 M=0;
  il Vec(){memset(t,0,sizeof(char)*MAXS);}
  il bool next(Vec&v){
    for(i4 i=0,w=1;i!=s;w*=(a[i++]+1)){
      if(t[i]==v[i]){M-=t[i]*w,t[i]=0;continue;}
      M+=w,t[i]++;return 1;
    }
    return 0;
  }
  il char&operator[](char x){return t[x];}
  il Vec operator+(Vec v){Vec x;for(char i=0;i!=s;i++)x[i]=t[i]+v[i];return x;}
  il Vec operator-(Vec v){Vec x;for(char i=0;i!=s;i++)x[i]=t[i]-v[i];return x;}
};

il i4 M(Vec v){i4 x=0;for(i4 i=0,w=1;i!=s;w*=(a[i++]+1))x+=v[i]*w;return x;}
il i8 D(Vec v){i8 x=1;for(i4 i=0;i!=s;i++)x*=(v[i]+1);return x;}
il i8 V(Vec v){i8 x=1;for(i4 i=0;i!=s;i++)x=x*powl(p[i],v[i])+.5;return x;}

il void conv(i4*f,Vec n_vd,Vec n_vp,const i4 Vd,const i4 Vp,const bool mod1,const bool mod2){
  static i4 F[MAXN];Vec i;
  if(mod1)do F[i.M+Vd]=f[i.M];while(i.next(n_vd));
  if(mod2)do F[i.M+Vp]+=P-f[i.M];while(i.next(n_vp));
  if(mod1)do f[i.M+Vd]=(f[i.M+Vd]+F[i.M+Vd])%P,F[i.M+Vd]=0;while(i.next(n_vd));
  if(mod2)do f[i.M+Vp]=(f[i.M+Vp]+F[i.M+Vp])%P,F[i.M+Vp]=0;while(i.next(n_vp));
}
il i4 slove(Vec n){
  static i4 f[MAXN];f[0]=1;i4 ans=0;
  std::map<i8,char>Pn;Vec d;
  for(i4 i=0;i!=s;i++){
    Vec d_1;i8 fact=p[i]-1;bool fg=1;
    for(auto [p_w,w]:Pn)while(fact%p_w==0)fact/=p_w,fg&=(++d_1[w]<=a[w]);
    if(fact!=1)fg=0;d[i]=1;conv(f,n-d_1,n-d,M(d_1),M(d),fg,1);d[i]=0;Pn[p[i]]=i;
  }
  do if(!Pn.count(V(d)+1)&&prime(V(d)+1))conv(f,n-d,d,d.M,0,1,0);while(d.next(n));
  do ans=(ans+D(n-d)*f[d.M])%P,f[d.M]=0;while(d.next(n));return ans;
}
int main(){
  std::ios::sync_with_stdio(0);
  std::cin.tie(0),std::cout.tie(0);
  std::cin>>T;
  while(T--){
    std::cin>>s;Vec n;
    for(i4 i=0;i!=s;i++)std::cin>>p[i]>>a[i],n[i]=a[i];
    if(p[0]!=2){std::cout<<"2\n";continue;}
    std::cout<<slove(n)<<"\n";
  }
  return 0;
}

    验题人代码

    #include <cstdio>
#include <iostream>
#include <cmath>
#define ll unsigned long long
using namespace std;

ll FastMul(ll x, ll y, ll mod){return (__int128)x * y % mod;}
ll FastPow(ll a, ll b, ll mod){
  ll res = 1;
  while(b){
    if(b & 1) res = FastMul(res, a, mod);
    b >>= 1, a = FastMul(a, a, mod);
  }
  return res;
}

ll FastPowNoMod(ll a, int b){
  ll res = 1;
  while(b){
    if(b & 1) res *= a;
    b >>= 1, a *= a;
  }
  return res;
}

const int ListNum = 12, List[ListNum] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};
bool Miller_Rabin(ll x){
  if(x == 1) return false;
  ll D = x - 1;
  while(!(D & 1)) D >>= 1;
  for(int i = 0; i < ListNum; ++i){
    if(x == List[i]) return true;
    ll k = FastPow(List[i], D, x), d = D;
    if(k != 1){
      while(d < x - 1){
        if(k == x - 1) break;
        d <<= 1, k = FastMul(k, k, x);
      }
      if(d == x - 1) return false;
    }
  }
  return true;
}

const int Mx = 103700, Nx = 16, Mod = 998244353;

ll N;
ll prime[Nx], ksm[Nx]; // power[i][j] = prime[i] ^ j
int alpha[Nx], cnt; // N = prod_{i \le cnt} prime[i] ^ alpha[i]
int sigma0, prod[Nx]; // Find position for vector bbeta[]

int F[Mx], G[Mx]; // ans = sum_{d|n} f(d) sigma0(n/d)

int bbeta[Nx], id; ll now; // Numerate divisor of n, now = prod_{i \le cnt} prime[i] ^ bbeta[i]
int gama[Nx], ids; // Numerate pres such that pre * now | n 
inline void getpos(ll x){
  id = 0;
  for(int i = 1; i <= cnt; ++i){
    int sum = 0;
    while(x % prime[i] == 0) x /= prime[i], ++sum;
    id += prod[i - 1] * sum, bbeta[i] = sum;
  }
}
inline bool pushup(){
  int pos = 1; ++id;
  while(bbeta[pos] == alpha[pos] && pos <= cnt) bbeta[pos] = 0, now /= ksm[pos], ++pos;
  return (pos == cnt + 1) ? 0 : (++bbeta[pos], now *= prime[pos], 1);
}
inline bool reduce(){
  int pos = 1;
  while(gama[pos] == 0 && pos <= cnt) gama[pos] = alpha[pos] - bbeta[pos], ids += prod[pos - 1] * gama[pos], ++pos;
  return (pos == cnt + 1) ? 0 : (--gama[pos], ids -= prod[pos - 1], 1);
}
inline bool add(){ // now = divisor(now)
  int pos = 1; ++id;
  while(bbeta[pos] == alpha[pos] && pos <= cnt) bbeta[pos] = 0, now /= (alpha[pos] + 1), ++pos;
  return (pos == cnt + 1) ? 0 : (++bbeta[pos], now += now / bbeta[pos], 1);
}

//int list[Mx];
//void print(){ cout << F[7] << endl; return;
//	for(int i = 0; i < sigma0; ++i) cout << list[i] << ' ' << F[i] << '\n';
//	cout << "===============\n";
//}

int Solve(){
  if(prime[1] != 2) return 2;
  prod[0] = N = 1;
  for(int i = 1; i <= cnt; ++i){
    prod[i] = prod[i - 1] * (alpha[i] + 1);
//		power[i][0] = 1;
    ksm[i] = FastPowNoMod(prime[i], alpha[i]);
//		for(int j = 1; j <= alpha[i]; ++j) power[i][j] = power[i][j - 1] * prime[i];
    N *= ksm[i];
  }
  sigma0 = prod[cnt];
  F[0] = 2, F[1] = -1;
  for(int i = 2; i < sigma0; ++i) F[i] = 0;
//	for(int i = 1; i <= cnt; ++i) bbeta[i] = 0;
//	now = 1, id = 0;
//	do list[id] = now; while(pushup());
//	print();
  for(int i = 2; i <= cnt; ++i){
    for(int j = 0; j < sigma0; ++j) G[j] = F[j];
    int idp = prod[i - 1];
    for(int j = 1; j <= cnt; ++j) gama[j] = alpha[j] - (i == j), bbeta[j] = (i == j);
    ids = sigma0 - 1 - idp;
    do F[ids + idp] = (F[ids + idp] - F[ids]) % Mod; while(reduce());
//		cout << "prime " << prime[i] << " mod " << N % (prime[i] - 1) << endl;
//		print();
    if(N % (prime[i] - 1) == 0){
      getpos(prime[i] - 1);
      for(int j = 1; j <= cnt; ++j) gama[j] = alpha[j] - bbeta[j];
      ids = sigma0 - 1 - id;
//			cout << "yesmod " << ids << ' ' << id << endl;
      do F[ids + id] = (F[ids + id] + G[ids]) % Mod; while(reduce());
    }
//		cout << "prime " << prime[i] << " don " << idp << endl;
//		print();
  }
  
  for(int i = 1; i <= cnt; ++i) bbeta[i] = 0; id = 0; now = 1;
  while(pushup()) if(N % (now + 1) != 0 && Miller_Rabin(now + 1)){
    for(int i = 1; i <= cnt; ++i) gama[i] = alpha[i] - bbeta[i];
    ids = sigma0 - 1 - id;
    do F[ids + id] = (F[ids + id] + F[ids]) % Mod; while(reduce());
//		cout << "divisor " << now << endl;
//		print();
  }
//	print();
  for(int i = 1; i <= cnt; ++i) bbeta[i] = 0; id = 0; now = 1;
  ll ans = 0;
  do ans += 1ll * now * F[sigma0 - 1 - id]; while(add());
  ans %= Mod;
  return ans < 0 ? ans + Mod : ans;
}

int main(){
  ios::sync_with_stdio(0);
  cin.tie(0), cout.tie(0);
  int T; cin >> T;
  while(T--){
    cin >> cnt;
    for(int i = 1; i <= cnt; ++i) cin >> prime[i] >> alpha[i];
    cout << Solve() << endl;
  }
  return 0;
}

    多元多项式优化方法

    观察多元普通生成函数：

    $$f(d)=\left[\prod_{i\ge1} x_i^{t_{i,d}}\right]\prod_{p\in\mathcal{P}}\left(1+\prod_{i\ge1} x_i^{t_{i,p-1}}-\prod_{i\ge1} x_i^{t_{i,p}}\right) $$

    仿照 P4389 付公主的背包 优化 0/1 背包：

    $$f(d)=\left[\prod_{i\ge1} x_i^{t_{i,d}}\right]\exp\left(\sum_{p\in\mathcal{P}}\ln\left(1+\prod_{i\ge1} x_i^{t_{i,p-1}}-\prod_{i\ge1} x_i^{t_{i,p}}\right)\right) $$

    对其展开成形式幂级数：

    $$\begin{aligned} f(d)&=\left[\prod_{i\ge1} x_i^{t_{i,d}}\right]\exp\left(\sum_{p\in\mathcal{P}}\sum_{k\ge1}\frac{(-1)^{k+1}}{k} \left(\prod_{i\ge1} x_i^{t_{i,p-1}}-\prod_{i\ge1} x_i^{t_{i,p}}\right)^k\right)\\ &=\left[\prod_{i\ge1} x_i^{t_{i,d}}\right]\exp\left(\sum_{p\in\mathcal{P}}\sum_{\substack{j,k\ge0\\[3pt]j+k\ne0}}\frac{(-1)^{j+1}}{j+k}\binom{j+k}{j}\left(\prod_{i\ge1} x_i^{t_{i,p-1}}\right)^j\left(\prod_{i\ge1} x_i^{t_{i,p}}\right)^k\right) \end{aligned} $$

    但是这样会有一个小问题：
    当 p=2p=2 时要计算 ln(2−x1)ln(2-x_1)
    但是很可惜，不能按上方一般的形式展开，所以我们可以将其单独提出：

    $$f(d)=\left[\prod_{i\ge1} x_i^{t_{i,d}}\right]\exp\left(\sum_{\substack{p\in\mathcal{P}\\[3pt]p\ne2}}\sum_{\substack{j,k\ge0\\[3pt]j+k\ne0}}\frac{(-1)^{j+1}}{j+k}\binom{j+k}{j}\left(\prod_{i\ge1} x_i^{t_{i,p-1}}\right)^j\left(\prod_{i\ge1} x_i^{t_{i,p}}\right)^k\right)(2-x_1) $$

    由于对于所有的 d∣nd\mid n，都有 ti,d<ti,nt_{i,d}<t_{i,n}
    且已经给出了 nn 的标准质因子分解形式 n=∏i=1spiαin=\prod_{i=1}^s p_i^{\alpha_i}。

    所以我们仅需计算 ss 维，每一维都分别的次数都分别低于 αi\alpha_i 的多项式的 exp。

    这里如果暴力递推计算 exp 仍然是 d3(n)d_3(n) 的时间复杂度。

    多元多项式

    接下来最后一步就是多元多项式 exp 了，
    我这里主要是复读 > EI 大神的 blog < 和有关 Uoj#596 的 blog，
    在此之前我们先浅谈一下多元多项式。

    这里约定多元多项式为 ss 维，每一维都分别的次数都分别 ≤αi\le\alpha_i 的多项式。
    那么它的值域 N=∏i=1s(αi+1)=σ0(n)N=\prod_{i=1}^s(\alpha_i+1)=\sigma_0(n)

    多元多项式 DFT/IDFT

    多元多项式 DFT 相当于按顺序依次对其每一维 DFT。
    多元多项式 IDFT 相当于反向按顺序的依次对其每一维 IDFT。

    那么为什么多元多项式 DFT 就是依次对其每一维 DFT 呢？
    因为做 DFT 相当于把系数表示法转为了点值表示法，
    所以依次对其每一维带入单位根即可计算点值。

    由于多元多项式 DFT 是线性算法。
    有「求逆原理」：$(A_1A_2\dots A_s)^{-1}=A_s^{-1}A_{s-1}^{-1}\dots A_1^{-1}$
    这样就得到了多元函数的 IDFT 算法。

    多元多项式乘法

    我们计算卷积 F(x1,x2,…,xs)⋅G(x1,x2,…,xs)F(x_1,x_2,\dots,x_s)\cdot G(x_1,x_2,\dots,x_s)，
    可以发现卷积后的值域为 ∏i=1s(2αi+1)=O(2sN)\prod_{i=1}^s(2\alpha_i+1)=O(2^sN)

    然而实际上大多数情况下我们仅需计算 $\bmod(x_1^{\alpha_1},x_2^{\alpha_2},\dots,x_s^{\alpha_s})$ 意义下的值。
    若 ss 较大，值域将会急剧膨胀，这被称为维度爆炸。
    而 EI 以一种构造性的算法解决了这个问题：

    显然，高维多项式是要避免的。
    于是我们将下标 i⃗=(i1,i2,…,is)\vec{i}=(i_1,i_2,\dots,i_s) 编码为一个 (α1+1,α2+1,...,αs+1)(\alpha_1+1,\alpha_2+1,...,\alpha_s+1) 进制数。
    即令 $\lvert\vec{i}\rvert=\sum\limits_{t=1}^si_t\prod_{n=1}^{t-1}(\alpha_n-1)=i_1+i_2*(\alpha_1+1)+...+i_s*(\alpha_1+1)*(\alpha_2+1)*...*(\alpha_{s-1}+1)$。
    这样就把多维映射到了一维上，下标 i⃗\vec{i} 和编码 ∣i⃗∣\lvert\vec{i}\rvert 是一一对应的。
    而且可以发现大多数时候有 $\lvert\vec{i}+\vec{j}\rvert=\lvert\vec{i}\rvert+\lvert\vec{j}\rvert$。
    所以多元多项式可以被这样的一元多项式来代替：

    $$F(x_1,x_2,\dots,x_s)=\sum_{\substack{\vec{i}=(i_1,\dots,i_s)\\i_1\le \alpha_1,\dots,i_s\le \alpha_s}}^{}f_{\vec{i}}x_1^{i_1}\dots x_s^{i_s}\textcolor{red}\to F(x)=\sum_{i=0}^{N-1}f_ix^i $$

    其中 fi⃗=f∣i⃗∣f_{\vec{i}}=f_{\lvert\vec{i}\rvert}。

    只有 i⃗+j⃗\vec{i}+\vec{j} 超过了模意义的范围时 $\lvert\vec{i}+\vec{j}\rvert=\lvert\vec{i}\rvert+\lvert\vec{j}\rvert$ 才不成立，
    所以我们要将超出范围（产生进位）的贡献去除（我们要计算模意义下的卷积）。
    考虑占位函数 χ(n)\chi(n)，使得 $\lvert\vec{i}+\vec{j}\rvert=\lvert\vec{i}\rvert+\lvert\vec{j}\rvert$ 成立当且仅当 χ(i)+χ(j)=χ(i+j)\chi(i)+\chi(j)=\chi(i+j)。
    这样，我们计算二元幂级数 ∑ifixiyχ(i)\sum_i f_ix^iy^{\chi(i)} 的卷积，然后利用第二维去除无用贡献即可。

    接下来就是 EI 的一个精妙构造了：

    由于 $[i_t+j_t\le \alpha_t]=\left\lfloor\frac{i+j}{\prod_{n=1}^{t-1}(\alpha_n+1)}\right\rfloor-\left\lfloor\frac{i}{\prod_{n=1}^{t-1}(\alpha_n+1)}\right\rfloor-\left\lfloor\frac{j}{\prod_{n=1}^{t-1}(\alpha_n+1)}\right\rfloor\in\{0,1\}$

    令 $\chi(i)=\sum\limits_{t=1}^{s-1}\left\lfloor \frac{i}{\prod_{n=1}^{t}(\alpha_n+1)}\right\rfloor=\left\lfloor \frac{i}{(\alpha_1+1)}\right\rfloor+...+\left\lfloor \frac{i}{(\alpha_1+1)(\alpha_2+1)...(\alpha_{s-1}+1)}\right\rfloor$

    这个占位函数是极好的，因为 χ(i+j)−χ(i)−χ(j)∈[0,s)\chi(i+j)-\chi(i)-\chi(j)\in[0,s)。
    所以 $\lvert\vec{i}+\vec{j}\rvert=\lvert\vec{i}\rvert+\lvert\vec{j}\rvert$ 成立当且仅当 χ(i)+χ(j)≡χ(i+j)(mods)\chi(i)+\chi(j)\equiv\chi(i+j)\pmod{s}。

    我们令 1≡ys1\equiv y^s，相当于模 ys−1y^s-1 的循环卷积。
    这样，值域也被这么压缩到了 sNsN。

    多元多项式全家桶

    多元多项式的其他运算可以直接套一元多项式的模板。

    对其求导可以考虑一个特殊的微分算子 D=xddx\mathfrak{D}=x\frac{d}{dx}，
    这里的 xx 是压缩成一维时的变量。

    可以发现它满足常见的导数性质，
    这样就解决了多项式 ln 和牛顿迭代。

    时间复杂度

    时间复杂度分为：

    1. 判 σ0(n)\sigma_0(n) 次质数（所有的 d+1d+1）：Θ(σ0(n)log⁡n)\Theta(\sigma_0(n)\log n)
    2. 计算 exp 之前的级数：O(σ0(n))O(\sigma_0(n))
    3. 多元多项式 exp：Θ(sσ0(n)log⁡σ0(n))\Theta(s\sigma_0(n)\log\sigma_0(n))

    第 22 处的时间复杂度证明（不严谨）：
    若所有质因子减 11 都还是因子：O(slog⁡2(n))O(s\log^2(n))
    若所有因子加 11 都是质数（除去上方情况）：A309891(n)+σ0(n)A309891(n)+\sigma_0(n)
    当 σ0(n)\sigma_0(n) 较大时这里的时间复杂度就为 O(σ0(n))O(\sigma_0(n))
    较小时就没必要考虑时间复杂度了

    第 33 处的时间复杂度以及实现方法，对于卷积部分：
    做 ss 个长为 2N2N 的 DFT，
    然后在 yy 维暴力进行卷积，相当于进行了 DFT+IDFT。
    再做 ss 个长为 2N2N 的 IDFT，
    卷积顺序正确性可以看上方的证明。
    复杂度为 $\Theta(sN\log N)+\Theta(s^2N)=\Theta(sN\log N)=\Theta(s\sigma_0(n)\log\sigma_0(n))$。
    而牛顿迭代并不会增加时间复杂度。

    总时间复杂度为：

    $$\Theta(s\sigma_0(n)\log\sigma_0(n)+\sigma_0(n)\log n) $$

    暴力 exp 版本代码

    是 Θ(sd3(n)+σ0(n)log⁡n)\Theta(sd_3(n)+\sigma_0(n)\log n) 的，过不了这题。

    #include<bits/stdc++.h>
namespace Miller_Rabin{
  const int Pcnt=7;
  const int64_t p[Pcnt]={2,325,9375,28178,450775,9780504,1795265022};
  int64_t mul(int64_t a,int64_t b,int64_t p){
    return __int128(a)*b%p;
  }
  int64_t pow(int64_t a,int64_t b,int64_t p){
    int64_t ans=1;
    for(;b;a=mul(a,a,p),b>>=1)if(b&1)ans=mul(ans,a,p);
    return ans;
  }
  bool prime(int64_t x){
    if(x<3||x%2==0)return x==2;
    int64_t u=x-1,t=0;
    while(u%2==0) u/=2,++t;
    int64_t ud[]={2,325,9375,28178,450775,9780504,1795265022};
    for(int64_t a:ud){
      int64_t v=pow(a,u,x);
      if(v==1||v==x-1||v==0) continue;
      for(int j=1;j<=t;j++){
        v=mul(v,v,x);
        if(v==x-1&&j!=t){v=1;break;}
        if(v==1) return 0;
      }
      if(v!=1) return 0;
    }
    return 1;
  }
}using Miller_Rabin::prime;

typedef std::complex<double> C;
typedef int32_t i32;
typedef int64_t i64;

const double PI2=6.2831853071795864769l;
const i32 MAXS=16,MAXN=104000,MAXA=64;
const i64 P=998244353;

i64 T,s,p[MAXA],a[MAXA];
i64 nCr[MAXA<<1][MAXA],Inv[MAXA<<1];
struct Vec{
  i32 t[MAXS];
  Vec(){memset(t,0,sizeof(i32)*MAXS);}
  bool next(Vec v,i32 w=-1){for(i32 i=0;i!=s;i++){if(i==w||t[i]==v[i]){t[i]=0;continue;}t[i]++;return 1;}return 0;}
  i32&operator[](i32 x){return t[x];}
  bool operator+=(Vec v){for(i32 i=0;i!=s;i++)if((t[i]+=v[i])>a[i])return 0;return 1;}
  Vec operator-(Vec v){Vec x;for(i32 i=0;i!=s;i++)x[i]=t[i]-v[i];return x;}
};
i64 fpow(i64 a,i64 n=P-2){i64 x=1;for(;n;(a*=a)%=P,n>>=1)if(n&1)(x*=a)%=P;return x;}

i64 M(Vec v){i64 x=0;for(i64 i=0,w=1;i!=s;w*=(a[i++]+1))x+=v[i]*w;return x;}
i64 D(Vec v){i64 x=1;for(i64 i=0;i!=s;i++)x*=(v[i]+1);return x;}
i64 S(Vec v){i64 x=0;for(i64 i=0;i!=s;i++)x+=v[i];return x;}
i64 V(Vec v){i64 x=1;for(i64 i=0;i!=s;i++)x=x*powl(p[i],v[i])+.5;return x;}

void Init(){
  for(i32 i=1;i!=MAXA<<1;i++)Inv[i]=fpow(i);
  for(i32 i=0,j=1;i!=MAXA<<1;i++)
    for(nCr[i][0]=1,j=1;j<=i&&j!=MAXA;j++)
      nCr[i][j]=(nCr[i-1][j-1]+nCr[i-1][j])%P;
}
void exp(int*f,Vec n){
  static int g[MAXN];Vec d,i;g[0]=1,d.next(n);
  do{do g[M(d)]=(g[M(d)]+S(i)*f[M(i)]%P*g[M(d-i)])%P;
    while(i.next(d));g[M(d)]=g[M(d)]*fpow(S(d))%P;
  }while(d.next(n));
  do f[M(d)]=g[M(d)],g[M(d)]=0;while(d.next(n));
}
int slove(Vec n){
  static int f[MAXN];
  std::map<i64,i32> Pn;Pn[2]=0;
  for(int i=1;i!=s;i++){
    Vec d_1,t;i64 fact=p[i]-1,j=0,k=0;
    for(auto [p_w,w]:Pn)
      while(fact%p_w==0)d_1[w]++,fact/=p_w;
    if(fact!=1)d_1[0]=MAXA;
    do{do f[M(t)]=(f[M(t)]+(j&1?1ll:P-1)*Inv[j+k]%P*nCr[j+k][j])%P;
      while((t[i]=++k)<=a[i]);t[i]=k=0;++j;
    }while(t+=d_1);Pn[p[i]]=i;
  }
  Vec d;d.next(n);do{
    if(Pn.count(V(d)+1)||!prime(V(d)+1))continue;
    Vec t;int j=1;
    for(;t+=d;j++)f[M(t)]=(f[M(t)]+(j&1?1ll:P-1)*Inv[j])%P;
  }while(d.next(n));
  exp(f,n);i64 ans=0;
  do ans=(ans+D(n-d)*(2ll*f[M(d)]+(M(d)%(a[0]+1)?P-f[M(d)-1]:0)))%P;
  while(d.next(n));memset(f,0,sizeof(i32)*D(n));return ans;
}
int main(){
  std::ios::sync_with_stdio(0);
  std::cin.tie(0),std::cout.tie(0);
  Init();std::cin>>T;
  while(T--){
    std::cin>>s;Vec n;
    for(int i=0;i!=s;i++)std::cin>>p[i]>>a[i],n[i]=a[i];
    if(p[0]!=2){std::cout<<"2\n";continue;}
    std::cout<<slove(n)<<"\n";
  }
  return 0;
}

    后话

    我其实还有几个小疑问：

    本题中 αi\alpha_i 都很小，是否有更快的方法进行卷积？
    总时间复杂度能不能做到 Θ(σ0(n)log⁡n)\Theta(\sigma_0(n)\log n)？

    本题中的几个数论函数如何快速求前缀和？各种筛法似乎都不可用。

    不过这就在我的能力之外了，欢迎各位大佬解决。

    • 1

    [Math×Girl] 互质与整除

    信息

    ID
    11389
    时间
    2000ms
    内存
    512MiB
    难度
    6
    标签
    递交数
    0
    已通过
    0
    上传者
    1. 关于
    2. 联系我们
    3. 隐私
    4. 服务条款
    5. Language
    6. 主题
    3. Worker 0, 41ms