21 Days Training // Day 16

T3. gcd

My idea : decompose

Notice that f(n) is multiplicative(?) .

And,

1. $f(1)=1$ (Useless…)
2. $f(p)=2$ Here $p$ is a prime.
3. $f(p)=2\sum_{i=0}^{\lfloor\frac{q}{2}\rfloor} p^i$ Here, $p$ is a prime and $q$ is an odd.
4. $f(p)=2\sum_{i=0}^{\frac{q}{2}} p^i – p^{\frac{q}{2}}$ Here, $p$ is a prime and $q$ is an even.

Use these properties, we can get f(n) quickly.

(Useless) Trick :

1. You only need to calculate all primes from 1 to 10^6 , since if a number cannot divide by these primes, it is definitely a prime.
2. For N \le 10^6 , you can directly use the algorithm with a linear-like complexity to calculate it; For N>10^6 , you only need to calculate from 10^6 to N since the answer of N \le 10^6 can be preprocessed before.

Code :

#include <bits/stdc++.h>
using namespace std;
const long long LIM=1000000;
bitset<LIM+10> isprime, vis;
long long lprime[LIM+11], plist[80000], pcount;

void preprocess() {
    lprime[1]=0;
    isprime[2]=1, vis[2]=1;
    lprime[2]=2;
    plist[pcount++]=2;
    for (long long j=4;j<=LIM+10;j+=2) vis[j]=1, lprime[j]=2;
    for (long long i=3;i<=LIM+10;i++) {
        if (vis[i]) continue;
        plist[pcount++]=i;
        lprime[i]=i, isprime[i]=1, vis[i]=1;
        for (long long j=2*i;j<=LIM+10;j+=i) vis[j]=1, lprime[j]=(lprime[j]?lprime[j]:i);
    }
}

int main() {
    freopen("gcd.in", "r", stdin);
    freopen("gcd.out", "w", stdout);
    preprocess();
    long long N;cin >> N;
    long long ans=0;
    if (N<=LIM) {
        for (long long i=1;i<=N;i++) {
            if (i==1) {
                ans++;
                continue;
            }
            if (isprime[i]) {
                ans+=2;
                continue;
            }
            long long bk=i;
            map<long long, int> pfac;
            while (bk!=1) {
                if (pfac.find(lprime[bk])==pfac.end()) pfac[lprime[bk]]=0;
                pfac[lprime[bk]]++;
                bk /= lprime[bk];
            }
            long long f=1;
            for (map<long long, int>::iterator it=pfac.begin();it!=pfac.end();it++) {
                long long expsum=0, cexp=1;
                for (int k=0;k<=it->second/2;k++) {
                    expsum += cexp;
                    cexp*=it->first;
                }
                cexp/=it->first;
                expsum *= 2;
                if (!(it->second&1)) expsum -= cexp;
                f *= expsum;
            }
            ans += f;
        }
    }
    else {
        long long pfx=39656708;
        for (long long i=LIM;i<=N;i++) {
            if (i==1) {
                ans++;
                continue;
            }
            long long bk=i;
            map<long long, int> pfac;
            for (int j=0;j<pcount&&bk>1;j++) {
                if (plist[j]>sqrt(N)) break;
                if (!(bk%plist[j])) {
                    pfac[plist[j]]=0;
                    while (!(bk%plist[j])) {
                        pfac[plist[j]]++;
                        bk /= plist[j];
                    }
                }
            }
            if (bk!=1) {
                pfac[bk]=1;
            }
            long long f=1;
            for (map<long long, int>::iterator it=pfac.begin();it!=pfac.end();it++) {
                long long expsum=0, cexp=1;
                for (int k=0;k<=it->second/2;k++) {
                    expsum += cexp;
                    cexp*=it->first;
                }
                cexp/=it->first;
                expsum *= 2;
                if (!(it->second&1)) expsum -= cexp;
                f *= expsum;
            }
            ans += f;
        }
        ans += pfx;
    }
    cout << ans;
    return 0;
}

However, it still can only get 20…