bzoj3512 DZY Loves Math IV

题目链接：bzoj3512
题目大意：
给定n,m,求$\sum_{i=1}^{n}\sum_{j=1}^{m}\phi(i j)$模$10^9+7$。
1<=n<=10^5，1<=m<=10^9，本题共4组数据。

话说，样例可真良心。。直接给你个极限的
Input
100000 1000000000
Output
857275582

题解：
杜教筛、数论
设$S(n,m)=\sum_{i=1}^{m}\phi(ni)$。
那么有$Ans=\sum_{i=1}^{n}S(i,m)$。
膜了一发题解，大概有两种吧。有种用到了一个显然的但是我不会的一个结论。
所以我去学了另一种虽然有点慢但是是我不会上面那个还能做的(期间还膜了膜奥爷爷和男神orz)。
主要是考虑如何弄开$S(n,m)=\sum_{i=1}^{m}\phi(ni)$。
那么重点就是这个$n$，将其分解质因数，分类：
- 若$n$有质因子$p$的指数大于1，设为$p^{r},r>1$，那么$\phi(ni)=p^{r-1}*\phi(\frac{n}{p^{r-1}}i)$。设$P$为$n$的所有质因子指数超过1的那部分，即$\frac{n}{P}=p_1p_2…p_n$，就是剩下的质因子指数都为1了。所以就有$S(n,m)=P*S(\frac{n}{P},m)$。

<script type="math/tex; mode=display" id="MathJax-Element-87"></script>
- 对于满足 $|\mu(n)|=1$的 $n$，设其有一质因子$p$，则有 $$S(n,m)=\phi(p)*S(\frac{n}{p},m)+S(n,\lfloor\frac{m}{p}\rfloor)$$ 为什么呢？
- 从这个式子 $S(n,m)=\sum_{i=1}^{m}\phi(ni)$出发
- 对于 $p$不能整除$i$的， $\phi(ni)=\phi(p)*\phi(\frac{n}{p}i)$，即 $$\sum_{i=1}^{n}\phi(ni)=\phi(p)*S(\frac{n}{p},m)=(p-1)*S(\frac{n}{p},m)$$
- 那么对于 $p|i$， $\phi(ni)=p*\phi(\frac{n}{p}i)$，即 $$\sum_{i=1}^{n}\phi(ni)=p*S(\frac{n}{p},m)$$
- 即对所有 $p|i$，都少了 $S(\frac{n}{p},m)=\sum_{i=1}^{m}\phi(\frac{n}{p}i)$，换元把 $\frac{n}{p}->i$，即 $\sum_{i=1}^{\lfloor\frac{m}{p}\rfloor}\phi(ni)=S(n,\lfloor\frac{m}{p}\rfloor)$
- 综上，就有 $$S(n,m)=\phi(p)*S(\frac{n}{p},m)+S(n,\lfloor\frac{m}{p}\rfloor)$$

然后呢，就是杜教筛的帕了，用其处理求欧拉函数的前缀和。
同样记忆化搜索。
对于$S$函数，要考虑边界判断，n=1时m=1时等等等等。
也是记忆化递归。
不同的是，这个不能用$map$来记忆化了，会MLE，状态数太多了。
所以就开个$f[i]$只存$S(i,m)$的。
怎么说，这种做法真的贼慢。。十几秒。。

#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<map>
#include<cmath>
#include<iostream>
#include<algorithm>
using namespace std;
typedef long long LL;
#define N 1000100

const LL ny=500000004;
const LL mod=1000000007;
map<int,LL> Phi;bool ispri[N];
int cnt,M,pri[N/10],p[N],minp[N];LL phi[N],lim,sum[N];
LL mymax(LL x,LL y){return (x>y)?x:y;}
int mymin(int x,int y){return (x<y)?x:y;}
void pre()
{
    cnt=0;phi[1]=sum[1]=1;p[1]=1;minp[1]=1;
    for (int i=2;i<=lim;i++)
    {
        if (!ispri[i]) {pri[++cnt]=i;phi[i]=(i-1)%mod;minp[i]=i;p[i]=1;}
        for (LL j=1;j<=cnt && i*pri[j]<=lim;j++)
        {
            ispri[i*pri[j]]=true;
            minp[i*pri[j]]=pri[j];
            if (i%pri[j]==0)
            {
                p[i*pri[j]]=p[i]*pri[j]%mod;
                phi[i*pri[j]]=phi[i]*pri[j]%mod;
                break;
            }
            p[i*pri[j]]=p[i];
            phi[i*pri[j]]=phi[i]*(pri[j]-1)%mod;
        }
    }
    for (LL i=2;i<=lim;i++) sum[i]=(sum[i-1]+phi[i])%mod;
}
LL solve(int n)
{
    if (n<=lim) return sum[n];
    if (Phi.count(n)) return Phi[n];
    LL ans=(LL)n*(n+1)%mod*ny%mod;int i,r;
    for (int i=2;i<=n;i=r+1)
    {
        r=mymin(n,n/(n/i));
        ans-=(LL)(r-i+1)*solve(n/i)%mod;
        ans%=mod;if (ans<0) ans+=mod;
    }
    Phi[n]=ans;
    return ans;
}
LL f[N];
LL s(int n,int m)
{
    if (m==0) return 0;
    if (m==M && f[n]) return f[n];
    if (n==1) return solve(m);
    if (m==1) return phi[n];
    LL k=minp[n];
    LL ans=((LL)(k-1)*s(n/k,m)%mod+s(n,m/k)%mod)%mod;
    if (m==M) f[n]=ans;
    return ans;
}
int main()
{
    //freopen("a.in","r",stdin);
    //freopen("a.out","w",stdout);
    int n,m,i;LL ans=0;
    scanf("%d%d",&n,&m);M=m;
    lim=(LL)trunc(pow(m,2.0/3.0));
    lim=mymax(lim,(LL)n);pre();
    for (i=1;i<=n;i++)
     ans=(ans+(LL)p[i]*s(i/p[i],m)%mod)%mod;
    if (ans<0) ans+=mod;
    printf("%lld\n",ans);
    return 0;
}