【JZOJ4732】【NOIP2016提高A组模拟8.23】函数（线筛求欧拉函数+Pollard_Rho大数分解）

Problem

• 定义一个函数f(x)，定义域和值域均为$N^*$，且满足$\sum_{d|n} f(d)=n$。
• 现给定N个幸运数字$A_1,A_2,...,A_n$，求$\sum_{i=1}^n f(A_i)$。

Hint

Solution

解谜f(x)——欧拉函数！！！

• 首先，f(x)实际上就是欧拉函数$\varphi(x)$，这是欧拉函数的性质。
• 证明起来也很简单。在$\sum_{d|n} f(d)=n$这个式子中，我们枚举d，不妨就看做枚举某个数x与n的最大公因数的逆元$\div n$，即$\frac nd=gcd(x,n)$。
• 既然这样，那么$f(d)=\varphi(d)$即为所有≤d且与d互质的数的个数。不妨设数a在其中之一，即a满足$a≤d\ ∧\ gcd(a,d)=1$，则$gcd(a\times \frac nd,n)=\frac nd$。
• 因此，满足$gcd(x,n)=\frac nd$的x的个数即为$\varphi(d)$。那么我们枚举每一种gcd（实际上也就是枚举n的所有约数），方案数总和即为n。

算法I：线筛求欧拉函数

• 我们知道$\varphi$有一种极其简单的线性求法。不知道的戳这里。
• 我们可以预处理出区间$[1,10^7]$中所有数的欧拉函数值。这样，就破掉第0、1、2、4、5、6、7个点。

---

• 时间复杂度：$O(10^7+N)$。
• 期望得分：70points。

第3个点

• Ai全是7，因此，答案为$3*10^7*\varphi(7)=3*10^7*6=1.8*10^8$。
• 注意N太大了，不能读入所有Ai；而观察到只有这个点的N为$3*10^7$，因此可以特判+输出+return。

---

• 时间复杂度：$O(1)$。
• 期望得分：10points。

算法II：Pollard_Rho大数分解

• 观察到第8、9个点的Ai都大的一匹，但是N却小得可怜。这启示我们应将所有Ai分解质因数，然后用最傻逼的那种方法求$\varphi$。
• 发现Ai有18位数的，因此，$O(\sqrt{A_i})$的分解吃不消。
• 所幸，我们有独特的Pollard_Rho大数分解算法，不懂的戳这里。

---

• 时间复杂度：$O(\sum_{i=1}^n A_i^{\frac 14})$。
• 直接用的话，估计只有20points；结合之前的算法，期望得分：100points。

Code

#include <bits/stdc++.h>
#define fo(i,a,b) for(i=a;i<=b;i++)
using namespace std;
typedef long long ll;
typedef long double ld;

const ll A=1e7;
int i,j,n,cnt,p[A],phi[A],x;
bool com[A+1];
ll a,ans,pf[100],d[5]={2,3,5,7,10007};
ld cur,tmp;

ll ti(ll x,ll y,ll m)
{
    ld d=1;
    d=d*x/m*y;
    return ((x*y-((ll)d)*m)%m+m)%m;
}
ll fpow(ll x,ll y,ll m)
{
    ll ans=1;
    for(;y;y>>=1)
    {
        if(y&1) ans=ti(ans,x,m);
        x=ti(x,x,m);
    }
    return ans;
}
bool miller_rabin(ll k)
{
    int i; ll t,n;
    fo(i,0,4)
    {
        if(k==d[i]) return 1;
        if(k%d[i]==0) return 0;
        n=k-1;
        while((n&1)==0) n>>=1;
        for(t=fpow(d[i],n,k); n<k-1&&t!=1&&t!=k-1; n<<=1) t=ti(t,t,k);
        if(~n&1&&t<k-1) return 0;
    }
    return 1;
}
ll gcd(ll x,ll y) 
{
    return y?gcd(y,x%y):x;
}
ll get(ll k)
{
    ll a=rand()%(k+1), b=a, c=rand()%k,r,i=1,n=2;
    while(1)
    {
        i++;
        a=(ti(a,a,k)+c)%k;
        if(a==b) return k;
        r=gcd(abs(a-b),k);
        if(1<r&&r<k) return r;
        if(i==n) b=a, n<<=1;
    }
}
void pollard_rho(ll k)
{
    if(k==1) return;
    if(miller_rabin(k)) {pf[++pf[0]]=k; return;}
    ll p=k;
    while(p==k) p=get(k);
    pollard_rho(k/p),pollard_rho(p);
}

int main()
{
    srand(unsigned(time(0)));
    scanf("%d",&n);
    if(n==3e7) {printf("%d",6*n); return 0;}
    phi[1]=1;
    fo(i,2,A)
    {
        if(!com[i]) phi[p[++cnt]=i]=i-1;
        for(j=1; j<=cnt&&1ll*i*(x=p[j])<=A; j++)
        {
            com[i*x]=1;
            if(i%x==0) {phi[i*x]=x*phi[i]; break;}
            phi[i*x]=phi[i]*(x-1);
        }
    }
    fo(i,1,n)
    {
        scanf("%lld",&a);
        if(a<=A) ans+=1ll*phi[a];   
        else
        {
            pf[0]=0;
            pollard_rho(a);
            sort(pf+1,pf+pf[0]+1);
            pf[0]=unique(pf+1,pf+pf[0]+1)-pf-1;
            cur=a;
            while(pf[0]) 
            {
                tmp=(ld)1/(ld)pf[pf[0]--];
                cur*=((ld)1-tmp);
            }
            ans+=cur;
        }
    } 
    printf("%lld",ans);
}