一些记号：
- 记 \(\mathbb{P}\) 为质数集，\(p_i\) 表示第 \(i\) 个质数。
- 记 \(\operatorname{lpf}(x)\) 表示 \(x\) 的最小质因数为第几个质数。
- 特别地，令 \(\operatorname{lpf}(1)=\infty\)。
设 \(c_i\) 表示 \(p_i\) 在 \(m!\) 中的次数，注意到当 \(p_i>\sqrt m\) 时，\(c_i=\lfloor\dfrac{m}{p_i}\rfloor\)，所以 \(c_i\) 的取值只有 \(O(\sqrt m)\) 种，而且 \(c_i\) 单调不降。
设 \(F(i,n)\) 表示不考虑前 \(i\) 个质数时的答案，即：
\[F(i,n)=\sum_{x=1}^n \sum_{d|m!}[\operatorname{lpf}(xd)>i]\sigma_0(xd) \]注意 \(xd=1\) 的情况也会被算入。
考虑递推到 \(F(i-1,n)\)，考虑 \(p_i\) 分别在 \(x\) 和 \(d\) 中出现的次数：
\[\begin{aligned} F(i-1,n)&=\sum_{k\geq 0,p_i^k\leq n}\sum_{0\leq j\leq c_i}(k+j+1) F(i,\lfloor\frac{n}{p_i^k}\rfloor)\\ &=\sum_{k\geq 0,p_i^k\leq n}\frac{1}{2}(c_i+1)(c_i+2k+2)F(i,\lfloor\frac{n}{p_i^k}\rfloor) \end{aligned} \]如果我们能处理出 \(p_{i+1}^2>n\) 时的边界情况，这个递推的时间复杂度就和 min25 筛是一样的了。
当 \(p_{i+1}^2>n\) 时，直接考虑定义式，\(F(i,n)\) 即为：
\[\begin{aligned} F(i,n)&=\sum_{d|m!}[\operatorname{lpf}(d)>i]\sigma_0(d)+\sum_{j>i,p_j\leq n}\sum_{d|m!}[\operatorname{lpf}(d)>i]\sigma_0(p_jd) \end{aligned} \]我们先考虑一个子问题 \(\sum\limits_{d|T}\sigma_0(d)\)，其中设 \(T=\prod\limits_i p_i^{t_i}\)。
相当于对每一个质因数 \(p_i\) 选一个 \(t'_i\leq t_i\) 以确定 \(d\)，然后这个 \(d\) 的贡献是 \(\prod\limits_i (1+t_i')\)。于是：
\[\begin{aligned} \sum\limits_{d|T}\sigma_0(d)&=\prod_{i}\left(\sum_{t_i'=0}^{t_i}(1+t_i')\right)\\ &=\prod_i \frac{1}{2}(t_i+1)(t_i+2) \end{aligned} \]然后再考虑 \(\sum\limits_{d|T}\sigma_0(p_jd)\)：
\[\begin{aligned} \sum_{d|T}\sigma_0(p_jd)&=\sum_{d|(p_jT)}\sigma_0(d)-\sum_{d|(T/p_j^{t_j})}\sigma_0(d)\\ &=\sum_{d|T}\sigma_0(d)\frac{\frac{1}{2}(t_j+2)(t_j+3)-1}{\frac{1}{2}(t_j+1)(t_j+2)}\\ \end{aligned} \]再看回原式，\(\sum\limits_{d|m!}[\operatorname{lpf}(d)>i]\sigma_0(p_jd)\) 其实可以转化成 \(\sum\limits_{d|(m!/\prod_{k\leq i}p_k^{c_k})}\sigma_0(p_jd)\)，于是有：
\[\begin{aligned} F(i,n)&=\sum_{d|m!}[\operatorname{lpf}(d)>i]\sigma_0(d)+\sum_{j>i,p_j\leq n}\sum_{d|m!}[\operatorname{lpf}(d)>i]\sigma_0(p_jd)\\ &=\sum\limits_{d|(m!/\prod_{k\leq i}p_k^{c_k})}\sigma_0(d)+\sum_{j>i,p_j\leq n}\sum\limits_{d|(m!/\prod_{k\leq i}p_k^{c_k})}\sigma_0(p_jd)\\ &=\sum\limits_{d|(m!/\prod_{k\leq i}p_k^{c_k})}\sigma_0(d)+\sum_{j>i,p_j\leq n}\sum_{d|(m!/\prod_{k\leq i}p_k^{c_k})}\sigma_0(d)\frac{\frac{1}{2}(c_j+2)(c_j+3)-1}{\frac{1}{2}(c_j+1)(c_j+2)}\\ &=\left(\prod_{j>i}\frac{1}{2}(c_i+1)(c_i+2)\right)\left(1+\sum_{j>i,p_j\leq n}\frac{\frac{1}{2}(c_j+2)(c_j+3)-1}{\frac{1}{2}(c_j+1)(c_j+2)}\right) \end{aligned} \]相当于我在数轴上有 \(O(\sqrt m)\) 个关键点，关键点之间的质数的 \(c\) 都相等，我在数轴上还有 \(O(\sqrt n)\) 个关键点用于询问，于是我们把这 \(O(\sqrt n)+O(\sqrt m)\) 个关键点混在一起用 min25 筛出来它们的素数个数即可。
有关 min25 筛复杂度的证明：
考虑 \(F(i,j)\) 的状态数，第二维的取值为 \(n,\lfloor\frac{n}{2}\rfloor,\cdots,\lfloor\frac{n}{\sqrt n}\rfloor,\sqrt n,\cdots,1\) 共 \(O(\sqrt n)\) 种（整除分块），第一维的取值范围是 \(O(\pi(\sqrt j))\)，于是共有状态数：
\[\begin{aligned} &\sum_{i=1}^{\sqrt n}O(\pi(\sqrt i))+\sum_{i=1}^{\sqrt n}O(\pi(\sqrt{\frac{n}{i}}))\\ \approx&\sum_{i=1}^{\sqrt n}O\!\left(\frac{\sqrt i}{\ln \sqrt i}\right)+\sum_{i=1}^{\sqrt n}O\!\left(\frac{\sqrt{\frac{n}{i}}}{\ln \sqrt{\frac{n}{i}}}\right)\\ \approx&O\left(\int_1^{\sqrt n}\frac{\sqrt x}{\ln \sqrt x}\text{ d}x+\int_{1}^{\sqrt n}\frac{\sqrt{\frac{n}{x}}}{\ln \sqrt{\frac{n}{x}}}\text{ d}x\right)\\ \approx&O\left(\frac{n^{\frac{3}{4}}}{\ln n}\right) \end{aligned} \]转移时的复杂度忽略，否则太难分析了（
所以你可以把复杂度当成 \(O(n^{\frac{3}{4}})\)。
写代码时写 id 那一段时竟然绕着绕着给自己绕晕了，后来发现只需要记住一件事：id 只是一种给关键点的编码方式，所以你只需要找到一种编码方式保证不会存在两个关键点编号相同即可。
数论题代码果然及其难写，所以代码也写得及其丑，跑得也及其慢。看啥时候有空规范一下自己数论题的码风（
#include<bits/stdc++.h>
#define N 100010
#define ll long long
#define LNF 0x7fffffffffffffff
using namespace std;
namespace modular
{
const int mod=323232323,inv2=(mod+1)/2;
inline int add(int x,int y){return x+y>=mod?x+y-mod:x+y;}
inline int dec(int x,int y){return x-y<0?x-y+mod:x-y;}
inline int mul(int x,int y){return 1ll*x*y%mod;}
inline void Add(int &x,int y){x=x+y>=mod?x+y-mod:x+y;}
inline void Dec(int &x,int y){x=x-y<0?x-y+mod:x-y;}
inline void Mul(int &x,int y){x=1ll*x*y%mod;}
}using namespace modular;
inline int poww(int a,int b)
{
int ans=1;
while(b)
{
if(b&1) ans=mul(ans,a);
a=mul(a,a);
b>>=1;
}
return ans;
}
inline ll read()
{
ll x=0;
int f=1;
char ch=getchar();
while(ch<'0'||ch>'9')
{
if(ch=='-') f=-1;
ch=getchar();
}
while(ch>='0'&&ch<='9')
{
x=(x<<1)+(x<<3)+(ch^'0');
ch=getchar();
}
return x*f;
}
int tot;
int cnt,prime[N];
bool notprime[N];
ll n,m,sn,sm;
ll kp[N*5];
namespace ID
{
int idn1[N],idn2[N],idm1[N],idm2[N],idmm[N];
int &id(ll x)
{
if(x<=n&&x==n/(n/x)) return x<=sn?idn1[x]:idn2[n/x];
if(x<=m&&x==m/(m/x)) return x<=sm?idm1[x]:idm2[m/x];
return idmm[x];
}
}using ID::id;
int is1[N],is2[N];
ll get(ll p)
{
ll x=m,ans=0;
while(x)
{
x/=p;
ans+=x;
}
return ans;
}
ll c[N];
void init()
{
const int nn=100000;
is1[0]=1;
for(int i=2;i<=nn;i++)
{
if(!notprime[i])
{
prime[++cnt]=i;
c[cnt]=get(i);
is1[cnt]=mul(is1[cnt-1],poww(mul(inv2,mul((c[cnt]+1)%mod,(c[cnt]+2)%mod)),mod-2));
int a=dec(mul(inv2,mul((c[cnt]+2)%mod,(c[cnt]+3)%mod)),1);
int b=mul(inv2,mul((c[cnt]+1)%mod,(c[cnt]+2)%mod));
is2[cnt]=add(is2[cnt-1],mul(a,poww(b,mod-2)));
}
for(int j=1;j<=cnt&&i*prime[j]<=nn;j++)
{
notprime[i*prime[j]]=1;
if(!(i%prime[j])) break;
}
}
}
ll g[N*5];
void work1()
{
for(int i=1;i<=tot;i++) g[i]=kp[i]-1;
for(int i=1;i<=cnt;i++)
{
ll p=prime[i],p2=p*p;
for(int j=tot;p2<=kp[j];j--)
g[j]-=g[id(kp[j]/p)]-(i-1);
}
}
int f[N*5];
int sum1[N*5],sum2[N*5];
bool vis[N*5];
int bound(int i,ll n)
{
int a=mul(sum1[tot],is1[i]);
int b=prime[i]<=n?dec(sum2[id(n)],is2[i]):0;
return mul(a,add(b,1));
}
void work2()
{
for(int i=cnt;i>=1;i--)
{
ll p=prime[i],p2=p*p,p22=(i==cnt?LNF:1ll*prime[i+1]*prime[i+1]);
for(int j=tot;p2<=kp[j];j--)
{
if(!vis[j]) f[j]=bound(i,kp[j]),vis[j]=1;
Mul(f[j],mul(inv2,mul((c[i]+1)%mod,(c[i]+2)%mod)));
ll pk=p;
for(int k=1;pk<=kp[j];k++,pk*=p)
{
int coef=mul(inv2,mul((c[i]+1)%mod,(c[i]+2*k+2)%mod));
ll v=kp[j]/pk;
if(p22>v) Add(f[j],mul(coef,bound(i,v)));
else Add(f[j],mul(coef,f[id(v)]));
}
}
}
}
int main()
{
n=read(),m=read();
sn=sqrt(n),sm=sqrt(m);
init();
for(ll l=1,r;l<=n;l=r+1)
kp[++tot]=r=(n/(n/l));
for(ll l=1,r;l<=m;l=r+1)
kp[++tot]=r=(m/(m/l));
for(ll i=1;i*i<=m;i++)
kp[++tot]=i;
sort(kp+1,kp+tot+1);
tot=unique(kp+1,kp+tot+1)-kp-1;
for(int i=1;i<=tot;i++) id(kp[i])=i;
work1();
sum1[1]=1;
//kp[i-1]+1 ~ kp[i]
for(int i=2;i<=tot;i++)
{
ll np=g[i]-g[i-1];
ll c=get(kp[i]);
sum1[i]=mul(sum1[i-1],poww(mul(inv2,mul((c+1)%mod,(c+2)%mod)),np));
int a=dec(mul(inv2,mul((c+2)%mod,(c+3)%mod)),1);
int b=mul(inv2,mul((c+1)%mod,(c+2)%mod));
sum2[i]=add(sum2[i-1],mul(np,mul(a,poww(b,mod-2))));
}
work2();
printf("%d\n",f[id(n)]);
return 0;
}
/*
10 3
*/
/*
323232323 23
*/