【cogs2593】幂，暴搜+容斥

传送门
思路：
这个题搞了好几天，终于卡过去了，但我的做法并不是特别正规，如果大家想了解更加科学的做法，可以去TA爷的51nod题解
说一下部分分好了：
10分随便做
30分我们可以把每一个$x(x\in[2,n])$拆成$e^k$的形式，其中$e$我称之为单位元（或者说是最简形式），$e$所有的质因子次数的$\gcd=1$，显然单位元不同的数字不可能相同，那么对于每一个$x^y=(e^k)^{y}$，那么就是考虑其能否拆成$(e^{k'})^{y'}$的形式，其中$k'<k,y'\leq m$，如果可以拆，那么就是重复的数字
其实大多数的暴力，就算能跑$n=500000000,m=1000000000$，也是30分做法
先说一下之前写的一个不科学做法，基本思路是把所有数都计算进答案($nm$)，然后减去那些能被更小的数+更大的幂表示的数，容易发现质因子$\gcd=1$，那么一定不会被重复计算，所以只要考虑那些$\gcd>1$的情况就可以了，而且$10^9$内这些数的数量级大约是$\sqrt n$,所以我决定写一发，但调了很久才发现，这玩意好像只能做到$O(m\log^2n)$，显然不能满足要求
聪哥看到这道题以后，告诉我这和51nod上的”矩阵中不重复的元素“类似，而且最难的$V3$被TA爷用暴搜艹过去了= =，顺便聪爷跟我说科学的思路——考虑单位元是没有问题的，但要转换成”乘法矩阵”的形式
这里所谓的”乘法矩阵”，每个单位元$e$都有一个宽为$m$，长为$⌊\log_en⌋$的矩阵，第$i$行$j$列的数就是$e^{ij}$，所有矩阵的大小之和正好是$(n-1)m$，也就是说，我们只要能统计出这些矩阵中的不重复元素，就可以很简单地统计出答案了。容易发现很多单位元的乘法矩阵大小是一致的，所以，而且大于$\sqrt n$的单位元矩阵只有一行，元素显然不重复，所以只用考虑小于等于$\sqrt n$的单位元就可以了；另外，矩阵宽一定，长的范围是$[1,29]$，所以方法就是预处理这些矩阵中的不重复元素，然后暴搜$\sqrt n$以内的单位元个数，统计答案，最后再加上剩下的那些大于$\sqrt n$的单位元的矩阵大小之和就可以了
后面那一步很好搞，但怎么预处理不重复元素呢
聪爷给出了一个非常简单的容斥，但这是错误的233
考虑前缀累加的情况下只用算当前列i关于之前列（包含本身）的不重复元素就可以了，自己想了一下，我们要考虑$i$与$[1,i-1]$的$lcm$来容斥，而且容斥时还要考虑选到的最小列，比如我们选了$a_1,a_2..a_k,i$，$a_1$是最小的列数，那么它的贡献就是$\displaystyle (-1)^k\frac {a_1m}{lcm(a_1,a_2..a_k,i)}$，这样裸做最坏复杂度是求$e=2$的情况，为$O(2^{\log_2n})≈O(n)$
然后我就拼命想优化，代码中$ok$数组就是当时的一个优化，不过好像只能让程序快$0.1s～0.2s$(?)
比较重要的一个优化是枚举最小列数j，如果$[j+1,i-1]$中有$lcm(i,j)$的因子，那么我们就不用考虑$j$的贡献了，因为容斥贡献的绝对值大小只与最小列数和lcm有关，如果有数字可选且为$lcm(i,j)$的因子，那么选它和不选它构成两个集合，贡献的符号相反，大小相同，所以一定是0
但是还是只能跑$2^{28}$，如果$n>2^{28}$就GG了
不过，我突然惊奇地想到，$29$是个质数？！
那我们直接在之前的$[1,28]$的计算中加上对$29$的计算，只是符号相反，多除一个29而已，最后还要考虑一下选不选1的问题，就是$m-\lfloor\frac m {29}\rfloor$了
最后用奇怪的姿势过了这道题，啧啧，感觉有些……赛艇？

代码：

#include<cstdio>
#include<iostream>
#include<cmath>
#define min(x,y) ((x)<(y)?(x):(y))
#define LL long long
using namespace std;
int n,m,lim,XX;
int prime[32000],gcd[65][65];
bool vis[32000];
LL ans,tot,sum[65],ok[33];
inline LL G(LL x,LL y){return y?G(y,x%y):x;}
void init()
{
    lim=sqrt(n);
    for (int i=2;i<=lim;++i)
    {
        if (!vis[i]) prime[++prime[0]]=i;
        for (int j=1;i*prime[j]<=lim&&j<=prime[0];++j)
        {
            vis[i*prime[j]]=1;
            if (i%prime[j]==0) break;
        }
    }
    for (int i=0;i<=60;++i)
        for (int j=0;j<=60;++j)
            gcd[i][j]=G(i,j);
}
inline void dfs1(int x,LL lcm,int mi,int val)
{
    if (1LL*x*m/lcm<1||x<=mi)
    {
        ans+=(1LL*mi*m)/lcm*val;
        sum[29]-=(1LL*mi*m)/lcm/29*val;
        return;
    }
    ok[x]=lcm;
    dfs1(x-1,lcm,mi,val);
    ok[x]=0;
    if (lcm%x)
        ok[x]=lcm,
        dfs1(x-1,lcm/G(x,lcm)*x,mi,-val),
        ok[x]=0;
    else
    {
        bool flag=0;
        for (int i=XX-1;i>x;--i)
            if (lcm%i==0)
            {
                if (lcm==ok[i]) flag=1;
                break;
            }
        if (!flag)
            ok[x]=lcm,
            dfs1(x-1,lcm,mi,-val);
            ok[x]=0;
    }
}
inline void dfs2(int x,LL cnt,int g)
{
    if (x>prime[0]||cnt*prime[x]>lim)
    {
        if (g!=1) return;
        int t=(log(n)+1e-6)/log(cnt);
        ans+=sum[t];
        tot-=1LL*t*m;
        return;
    }
    LL t=1;
    for (int i=1;t*cnt<=lim;++i)
    {
        dfs2(x+1,cnt*t,gcd[g][i-1]);
        t*=prime[x];
    }
}
main()
{
    scanf("%d%d",&n,&m);
    init();
    sum[1]=m;
    for (int i=2;1<<i<=n&&i<=28;++i)
    {
        ans=sum[i-1];
        XX=i;
        dfs1(i-1,i,i,1);
        for (int j=i-1;j>=1;--j)
        {
            bool flag=0;
            for (int k=j+1;k<i;++k)
                if (i/gcd[i][j]*j%k==0)
                {
                    flag=1;
                    break;
                }
            if (!flag) dfs1(i-1,i/gcd[i][j]*j,j,-1);
        }
        sum[i]=ans;
    }
    sum[29]+=sum[28]+m-m/29;
    tot=1LL*n*m-m+1;
    ans=0;
    dfs2(1,1,0);
    printf("%lld\n",ans+tot);
}