【线性筛法求解积性函数】Archer

计算sigma（a=1~n）sigma（b=1~n）lcm(a,b).

oimaster的好题。

考虑a与b枚举是同域的，对于一个a只需统计比他小的的b，然后乘2即可，

则构造函数f[a]=-a+2*lcm(a,b)【b=1~a】；

设gcd（a,b）=d

f[a]=-a+2a*sigma(d|a)sigma(b<=a and gcd(a,b)=d) b/d

令b/d=k,

f[a]=-a+2a*sigma(d|a)sigma(k<=a/d and gcd(a/d,k)=1) k

因为（a/d）|a

则f[a]=-a+2a*sigma(d|a)sigma(k<=d and gcd(d,k)=1) k

做了这么多，其实只是想用一个公式，当n>1,1~n中与n互质的数的和为：（n*phi（n））/2

简略证明：对于一个与n互质的数d，n-d与n也互质，所以我们两两相加，一共有phi（n）/2对，而每对的和为n，得证。

f[a]=-a+2a*((sigma(d|a and d>1) d*phi(d)/2 )+ 1)

因为phi(1)=1

化解为f[a]=a*sigma(d|a) d*phi(d);

n为积性，phi为积性，所以f[a]为积性，

f[p^k]=(p^(3k+1)+p^k)/(p+1) 这一步我推了两次，每次都推错，明明是体力活的说。

得出递推式f[p^k]=p^3*f[p^(k-1)]-p^k*(p-1)

剩下就是线性筛法的事，论文上没说，我自己推了一下，却做不到线性。

如果是质数，就直接算。

如果正好包含p^1,则乘以f[p]

如果包含p^k,则设A*p^k=n,

现在已知f[A*p^(k-1)],设A*p^(k-1)=i;

f[n]=f[A]*f[p^k]=f[A]*（p^3*f[p^(k-1)]-p^k*(p-1)）

=p^3*f[A]*f[p^(k-1)]-f[A]*p^k*(p-1)

=p^3*f[i]-f[A]*p^k*(p-1)

推不下去了，所以log的分解算了。

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
const int maxn=1000000;
long long f[maxn+1],s[maxn+1],phi[maxn+1],prim[maxn+1],tot,check[maxn+1];
long long sqt(long long x) {return x*x*x;}
void make(long long x,long long k,long long &ai,long long &bi)
{
  ai=bi=0;
  for (;x%k==0;) {
    if (!bi) bi=(k-1);else bi*=k;
    x/=k;
  }
  ai=x;
}
void origin()
{
  long long ne,i,ai,bi,j;
  phi[1]=1,f[1]=1;
  tot=0;
  for (i=2;i<=maxn;i++) {
    if (!check[i]) {
      prim[++tot]=i;
      phi[i]=i-1;
      f[i]=sqt(i)*f[1]-i*phi[i];
    }
    for (j=1;j<=tot;j++) {
      ne=i*prim[j];
      if (ne>maxn) break;
      check[ne]=1;
      if (i%prim[j]==0) {
        phi[ne]=phi[i]*prim[j];
        make(ne,prim[j],ai,bi);
        f[ne]=sqt(prim[j])*f[i]-prim[j]*f[ai]*bi;
//        f[ne]=sqt(prim[j])*f[i]-(prim[j]-1)*ai*bi;
        break;
      }
      else {
        phi[i*prim[j]]=phi[i]*(prim[j]-1);
        f[ne]=f[i]*f[prim[j]];
      }
    }
  }
  for (i=1;i<=maxn;i++) s[i]=s[i-1]+f[i];
}
int main()
{
  freopen("Archer.in","r",stdin);
  freopen("Archer.out","w",stdout);
    long long t,i,x;
    origin();
    scanf("%I64d\n",&t);
    for (i=1;i<=t;i++) {
      scanf("%I64d\n",&x);
      printf("%I64d\n",s[x]);
    }
  return 0;
}