【51Nod 1363】最小公倍数之和

Description

给出一个n，求1-n这n个数，同n的最小公倍数的和。
例如：n = 6，1,2,3,4,5,6 同6的最小公倍数分别为6,6,6,12,30,6，加在一起 = 66。
由于结果很大，输出Mod 1000000007的结果。

Solution

做这道题，这是历经波澜！
$ans=\sum_{i=1}^n {i*n\over gcd(i,n)}$
设$f(d)=\sum_{i=1}^n i(gcd(i,n)=d)$
$ans=\sum_{d=1}^n {f(d)*n\over d}$
$f(d)=\sum_{i=1}^n i(gcd(i,n)=d)=\sum_{i=1}^n i(gcd(i/d,n/d)=1)={\varphi(n/d)*(n/d)\over 2}$
所以$ans=n+{n\over 2}\sum_{d|n,d≠n} \varphi(d)*d$
但是直接这样打只会对一半的点，时间复杂度不够大。
现在可以考虑，顺着质因数去线性算$\varphi(i)*i$
但是还可以转化一下公式。
设$n=\prod_{i=1}^B p[i]^{a[i]}$（p[i]是互异的质因数）
$\sum_{d|n,d≠n} \varphi(d)*d=\prod_{i=1}^B\sum_{j=0}^{a[i]} \varphi(p[i]^j)* p[i]^j$（把所有的质数，然后枚举其次方，所有的情况乘起来组成所有的因数）
$=\prod_{i=1}^B 1+\sum_{j=1}^{a[i]}(p[i]-1)p[i]^{j-1}* p[i]^j$
$=\prod_{i=1}^B 1+\sum_{j=1}^{a[i]}(p[i]-1)p[i]^{2j-1}$
$=\prod_{i=1}^B 1+(p[i]-1){p[i]^{2*a[i]+1}-p[i]\over {p[i]^2-1}}$
$=\prod_{i=1}^B 1+{p[i]^{2*a[i]+1}-p[i]\over {p[i]+1}}$
然后可以快速分解，打了一下不过由于不优美烂了。
然后想了一下。
预处理出根号范围内的质数，然后分解质因数，本来是根号的速度的，但是根号的那么范围$s\sqrt x$因为被分解后是不断的减小的，所以大部分的数据是可以过得。

Code

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define fo(i,a,b) for(i=a;i<=b;i++)
using namespace std;
const int mo=1000000007,ni2=500000004;
typedef long long ll;
ll i,j,k,l,t,n,m,cas;
ll ans,p[50000],ph[50007],pp[50007];
ll zhi[50007];
bool bz[5000001];
ll qsm(ll x,ll y){
    ll z=1;
    while(y){
        if(y&1)z=z*x%mo;
        x=x*x%mo;
        y/=2;
    }
    return z;
}
ll doing(ll x,ll y){
    ll o=(qsm(x,2*y+1)-x+mo)%mo;
    o=o*qsm(x%mo+1,mo-2)%mo;
    o=(o+mo)%mo;
    return o;
}
int main(){
    fo(i,2,40000){
        if(!bz[i])zhi[++zhi[0]]=i;
        fo(j,1,zhi[0]){
            t=zhi[j]*i;if(t>40000)break;bz[t]=1;
            if(!(i%zhi[j])){;break;}
        }
    }
    for(scanf("%lld",&cas);cas;cas--){
        scanf("%lld",&n);
        ll x=n;pp[0]=0;
        fo(i,1,zhi[0]){
            if(zhi[i]*zhi[i]>x)break;
            if(x%zhi[i]==0)pp[++pp[0]]=zhi[i],p[pp[0]]=0;
            while(x%zhi[i]==0)x/=zhi[i],p[pp[0]]++;
        }
        if(x>1)pp[++pp[0]]=x,p[pp[0]]=1;
        ans=1;
        fo(i,1,pp[0]){                 
            ans=ans*(doing(pp[i],p[i])+1)%mo;
        }
        ans=ans*n%mo*ni2%mo;
        ans=(ans+n%mo*ni2%mo)%mo;
        printf("%lld\n",ans);
    }
}