NKOJ 3961 zMy的lcm （Pollard_Rho+欧拉函数+高精度）

P3961zMy的lcm

问题描述

---

首先注意到，题目中的异或是不存在的，因为出题人懒得写高精度异或。
先推导一番

$$Ans=\sum_{i=1}^{n}lcm(n,i)=n\sum_{i=1}^{n}\frac{i}{gcd(n,i)}=n\sum_{d|n}\sum_{i=1}^{n}\frac{i}{d}[gcd(n,i)==d]=n\sum_{d|n}\sum_{i=1}^{\frac{n}{d}}i[gcd(\frac{n}{d},i)==1]$$

这里用到一个公式，$\sum_{i=1}^{n}i[gcd(n,i)==1]=\frac{n\varphi(n)}{2}$
关于证明，考虑到如果$gcd(i,n)=1$，那么$gcd(n-i,n)=1$，由此得证。

因此

$$Ans=\frac{n}{2}\sum_{d|n}\frac{n}{d}\varphi(\frac{n}{d})=\frac{n}{2}\sum_{d|n}d\varphi(d)，令f(n)=\sum_{d|n}d\varphi(d)$$
考虑如何求$f(n)$，观察发现，$f(n)$可能是积性函数，尝试证明一下。
$$f(a)=\sum_{i|a}i\varphi(i)，f(b)=\sum_{j|b}j\varphi(j)，gcd(a,b)=1$$
那么显然 $$f(ab)=\sum_{i|a}\sum_{j|b}ij\varphi(ij)=\sum_{i|a}\sum_{j|b}ij\varphi(i)\varphi(j)=\sum_{i|a}i\varphi(i)\sum_{j|b}j\varphi(j)=f(a)f(b)$$
因此可以分解质因数后求$f(n)$，只需要考虑求$f(p^k)$，推导一下容易发现$f(p^k)=\frac{p^{2k}-1}{p+1}p+1$

因此只需要用$Pollard\ Rho$求质因数分解，然后高精度算一算就行了。
总复杂度$O(Tn^{\frac{1}{4}}+高精度复杂度)$

---

代码：

#include<stdio.h>
#include<algorithm>
#include<cstring>
#include<cstdlib>
#include<ctime>
#define ll long long
using namespace std;
const ll mod=1e9;
struct Big
{
    ll cur,S[500];
    void Output()
    {
        printf("%lld",S[cur]);
        ll i,k;
        for(i=cur-1;i>=0;i--)
        {
            k=mod/10;
            while(k>S[i])putchar('0'),k/=10;
            if(k)printf("%lld",S[i]);
        }
    }
    void add(ll k)
    {
        S[0]+=k;ll i=0;
        while(S[i]>=mod)S[i+1]+=S[i]/mod,S[i++]%=mod,cur=max(i,cur);
        while(S[cur+1])cur++;
    }
    void Add(const Big& o)
    {
        ll i,r=max(o.cur,cur);
        for(i=0;i<=r;i++)
        {
            S[i]+=o.S[i];
            if(S[i]>=mod)S[i+1]+=S[i]/mod,S[i]%=mod;
        }
        cur=r+5;while(cur&&S[cur]==0)cur--;
    }
    void multiply(ll k)
    {
        for(ll i=0;i<=cur;i++)S[i]*=k;
        for(ll i=0;i<cur;i++)if(S[i]>=mod)S[i+1]+=S[i]/mod,S[i]%=mod;
        while(S[cur]>=mod)S[cur+1]+=S[cur]/mod,S[cur]%=mod,cur++;
    }
    void Multiply(const Big o,Big& E)
    {
        ll i,j;memset(&E,0,sizeof(E));
        for(i=0;i<=cur;i++)
        for(j=0;j<=o.cur;j++)
        {
            E.S[i+j]+=S[i]*o.S[j];
            if(E.S[i+j]>=mod)E.S[i+j+1]+=E.S[i+j]/mod,E.S[i+j]%=mod;
        }
        E.cur=cur+o.cur+5;while(E.cur>0&&E.S[E.cur]==0)E.cur--;
    }
    void divide(ll k)
    {
        for(ll i=cur;i>0;i--)
        {
            S[i-1]+=S[i]%k*mod;S[i]/=k;
            if(S[cur]==0)cur--;
        }
        S[0]/=k;
    }
    void operator=(const Big& B)
    {
        memset(S,0,sizeof(S));cur=B.cur;
        for(ll i=cur;i>=0;i--)S[i]=B.S[i];
    }
}Ans,A,B,D,Zero;
ll T,P[5]={2,3,5,7,11},C[233],tot;
ll gcd(ll a,ll b)
{return b==0?a:gcd(b,a%b);}
ll QC(ll a,ll b,ll c)
{
    ll o=0;a%=c;
    while(b)
    {
        if(b&1)o+=a,o-=o>=c?c:0;
        b>>=1;a+=a;a-=a>=c?c:0;
    }
    return o;
}
ll QM(ll a,ll b,ll c)
{
    ll o=1;
    while(b)
    {
        if(b&1)o=QC(o,a,c);
        b>>=1;a=QC(a,a,c);
    }
    return o;
}
ll QQ(ll a,ll b)
{
    ll o=1;
    while(b)
    {
        if(b&1)o*=a;
        b>>=1;a*=a;
    }
    return o;
}
bool MR(ll x)
{
    if(!x&1)return x==2;
    if(x==1)return 0;
    ll i,j,k,p,t=0,m,pre;
    m=x-1;while(!m&1)m>>=1,t++;
    for(i=0;i<5;i++)
    {
        p=P[i];if(x%p==0)return x==p;
        p=QM(p,m,x);
        if(p==1)continue;
        for(j=1;j<=t;j++)
        {
            pre=p;p=QC(p,p,x);
            if(p==1)
            {
                if(pre!=1&&pre!=x-1)return 0;
                break;
            }
        }
        if(p!=1)return 0;
    }
    return 1;
}
void RHO(ll x)
{
    if(MR(x)){C[++tot]=x;return;}
    ll i,k,a,b,c,t;
    while(1)
    {
        i=1;k=2;
        a=b=1ll*rand()*rand()%x;
        c=1ll*rand()*rand()%x;
        while(i++)
        {
            b=(QC(b,b,x)+c)%x;
            if(a==b)break;
            t=gcd(a>b?a-b:b-a,x);
            if(t!=1){RHO(t);RHO(x/t);return;}
            if(i==k)k<<=1,a=b;
        }
    }
}
int main()
{
    srand(time(NULL));
    ll i,j,tmp;ll n;
    scanf("%lld",&T);
    while(T--)
    {
        scanf("%lld",&n);
        if(n==1){printf("1\n");continue;}
        Ans=Zero;Ans.add(1ll);
        tot=0;RHO(n);
        sort(C+1,C+tot+1);
        for(i=1;i<=tot;i=j)
        {
            j=i+1;while(j<=tot&&C[j]==C[i])j++;
            tmp=QQ(C[i],j-i);
            if(j-i>1)
            {
                A=Zero;B=Zero;
                A.add(tmp+1);
                B.add(tmp-1);
                A.Multiply(B,D);
                D.multiply(C[i]);
                D.divide(C[i]+1);
                D.add(1);
                Ans.Multiply(D,B);
                Ans=B;
            }
            else
            {
                A=Zero;B=Zero;
                A.add(tmp);
                B.add(tmp-1);
                A.Multiply(B,D);
                D.add(1);
                Ans.Multiply(D,B);
                Ans=B;
            }
        }
        A=Zero;
        A.add(n);
        Ans.divide(2);
        Ans.Multiply(A,B);
        B.add(n);
        B.Output();putchar('\n');
    }
}