HDU 5528 Count a * b(线性筛+积性函数)

$去年长春赛区的B题，金牌数论题$
$我用了比较丑陋的方法过的，其实这题可以推导$
$但是看了人家推的，除了叉姐的我看得懂，其他人的我都看不懂$
$先打个表看下里面0和非0元素的个数把$
$很快就发现，如果一个数字不是全是一个因子的次方的话，拆成两个互质的数相乘即可$
$如果是x^n的话，这得好好观察，经验来说一般有公式$
$f[x^{n+1}]=x*f[x^n]+m$
$凑一下这个m就好了$
$当然这题还没那么简单，n是10^9，所以需要先筛因子，然后对于每个因子，求这个公式$
$考虑到因子很大，可以用递归，100W以下打表，以上的话就筛最小的素因子$
$跑得比较慢，但是考虑到比赛时候是应该能过的$

---

代码

#include <map>
#include <set>
#include <stack>
#include <queue>
#include <cmath>
#include <string>
#include <vector>
#include <cstdio>
#include <cctype>
#include <cstring>
#include <sstream>
#include <cstdlib>
#include <iostream>
#include <algorithm>
#pragma comment(linker,"/STACK:102400000,102400000")

using namespace std;
#define   MAX           1000005
#define   MAXN          1000005
#define   maxnode       105
#define   sigma_size    30
#define   lson          l,m,rt<<1
#define   rson          m+1,r,rt<<1|1
#define   lrt           rt<<1
#define   rrt           rt<<1|1
#define   middle        int m=(r+l)>>1
#define   LL            long long
#define   ull           unsigned long long
#define   mem(x,v)      memset(x,v,sizeof(x))
#define   lowbit(x)     (x&-x)
#define   pii           pair<int,int>
#define   bits(a)       __builtin_popcount(a)
#define   mk            make_pair
#define   limit         10000

//const int    prime = 999983;
const int    INF   = 0x3f3f3f3f;
const LL     INFF  = 0x3f3f;
const double pi    = acos(-1.0);
const double inf   = 1e18;
const double eps   = 1e-8;
const LL     mod   = 1e9+7;
const ull    mx    = 133333331;

/*****************************************************/
inline void RI(int &x) {
      char c;
      while((c=getchar())<'0' || c>'9');
      x=c-'0';
      while((c=getchar())>='0' && c<='9') x=(x<<3)+(x<<1)+c-'0';
 }
/*****************************************************/

bool prime[MAX];
int pr[MAX];
int num[MAX];
ull f[MAX];
int tot;

void init(){
    mem(prime,0);tot=0;
    f[1]=1;
    for(int i=2;i<MAX;i++){
        if(!prime[i]){
            pr[tot++]=i;
            f[i]=2*i-1;
            num[i]=i;
        }
        for(int j=0;j<tot&&pr[j]*i<MAX;j++){
            prime[pr[j]*i]=1;
            if(i%pr[j]==0){
                if(i==num[i]){
                    num[i*pr[j]]=num[i]*pr[j];
                    f[i*pr[j]]=f[i]*pr[j]+(ull)(pr[j]-1)*i;
                }
                else{
                    f[i*pr[j]]=f[i/num[i]]*f[pr[j]*num[i]];
                    num[i*pr[j]]=num[i]*pr[j];
                }
                break;
            }
            else{
                num[i*pr[j]]=pr[j];
                f[i*pr[j]]=f[i]*f[pr[j]];
            }
        }
    }    
}

vector<pii> v;
int divi[MAX];
int xx;
void dfs(int x,int y){
    if(x==v.size()){
        divi[xx++]=y;
        return;
    }
    for(int i=0;i<=v[x].second;i++){
        dfs(x+1,y);
        y*=v[x].first;
    }
}

ull dfs(int x){
    if(x<MAX) return f[x];
    int k=x;
    for(int i=0;i<tot&&pr[i]*pr[i]<=k;i++){
        if(k%pr[i]==0){
            int ret=1;
            while(k%pr[i]==0){
                k/=pr[i];
                ret*=pr[i];
            }
            if(k==1){
                return dfs(x/pr[i])*pr[i]+(ull)(pr[i]-1)*(x/pr[i]);
            }
            else return dfs(k)*dfs(x/k);
            break;
        }
    }
    if(k!=1) return 2*k-1;
}
int main(){
    init();
    int t;
    cin>>t;
    while(t--){
        int n;
        scanf("%d",&n);
        int x=n;
        v.clear();
        for(int i=0;i<tot&&pr[i]*pr[i]<=x;i++){
            if(x%pr[i]==0){
                int j=0;
                while(x%pr[i]==0){
                    x/=pr[i];
                    j++;
                }
                v.push_back(mk(pr[i],j));
            }
        }
        if(x!=1) v.push_back(mk(x,1));
        xx=0;
        dfs(0,1);
        ull ans=0;
        for(int i=0;i<xx;i++){
            ull tmp=dfs(divi[i]);
            ans+=(ull)divi[i]*divi[i]-tmp;
            //cout<<tmp<<endl;
        }
        printf("%I64u\n",ans);
    }
    return 0;
}