JZOJ 5025 Sum

感性理解后面部分是关于 i 的积性函数（设为 f(i)），然后考虑质数幂处的值。

f(pc)=pc∑j=1pc∑k=1lcm(gcd(pc,j),gcd(pc,k))=c∑i=0pipc∑j=1pc∑k=1[lcm(gcd(pc,j),gcd(pc,k))=pi]=c∑i=0pipc∑j=1pc∑k=1([lcm(gcd(pc,j),gcd(pc,k))∣pi]−[lcm(gcd(pc,j),gcd(pc,k))∣pi−1])=c∑i=0pi[(i∑d=0φ(pc−d))2−(i−1∑d=0φ(pc−d))2]=c∑i=0pi[(pc−[i<c]pc−i−1)2−(pc−pc−i)2]=(2c+1)(p2c−p2c−1)+pc−1

然后 Min_25 筛即可。

代码：

#include <cmath>
#include <cstdio>
using namespace std;
const int MX = 1e5;
long long n;
int lim;
int vis[MX + 5],cnt,prime[MX + 5];
unsigned Gprime[MX + 5][2];
int tot,le[MX + 5],ge[MX + 5];
long long lis[2 * MX + 5];
unsigned G[2 * MX + 5][3],Fprime[2 * MX + 5];
inline int &id(long long x)
{
    return x <= lim ? le[x] : ge[n / x];
}
unsigned F(int k,long long n)
{
    if(n < prime[k] || n <= 1)
        return 0;
    unsigned ret = Fprime[id(n)] - (3 * Gprime[k - 1][1] - 3 * Gprime[k - 1][0] + (k - 1));
    for(register int i = k;i <= cnt && (long long)prime[i] * prime[i] <= n;++i)
    {
        long long pw0 = 1,pw = prime[i],pw2 = (long long)prime[i] * prime[i];
        for(register int c = 1;pw2 <= n;pw0 = pw,pw = pw2,pw2 *= prime[i],++c)
            ret += ((2 * c + 1) * (pw * pw - pw * pw0) + pw0) * F(i + 1,n / pw) + ((2 * c + 3) * (pw2 * pw2 - pw2 * pw) + pw);
    }
    return ret;
}
int main()
{
    freopen("sum.in","r",stdin),freopen("sum.out","w",stdout);
    scanf("%lld",&n),lim = sqrt(n);
    for(register int i = 2;i <= lim;++i)
    {
        if(!vis[i])
            prime[++cnt] = i,Gprime[cnt][0] = Gprime[cnt - 1][0] + i,Gprime[cnt][1] = Gprime[cnt - 1][1] + i * i;
        for(register int j = 1;j <= cnt && i * prime[j] <= lim;++j)
        {
            vis[i * prime[j]] = 1;
            if(!(i % prime[j]))    
                break;
        }
    }
    for(register long long l = 1,r;l <= n;l = r + 1)
    {
        r = n / (n / l);
        lis[id(n / l) = ++tot] = n / l;
        G[tot][0] = n / l - 1;
        if(n / l & 1)
            G[tot][1] = (n / l + 1) / 2 * (n / l) - 1;
        else
            G[tot][1] = (n / l) / 2 * (n / l + 1) - 1;
        long long tmp[3];
        tmp[0] = n / l,tmp[1] = n / l + 1,tmp[2] = 2 * (n / l) + 1;
        for(register int i = 0;i < 3;++i)
            if(!(tmp[i] % 2))
            {
                tmp[i] /= 2;
                break;
            }
        for(register int i = 0;i < 3;++i)
            if(!(tmp[i] % 3))
            {
                tmp[i] /= 3;
                break;
            }
        G[tot][2] = tmp[0] * tmp[1] * tmp[2] - 1;
    }
    for(register int k = 1;k <= cnt;++k)
    {
        int p = prime[k];
        long long s = (long long)p * p;
        for(register int i = 1;lis[i] >= s;++i)
            G[i][0] -= G[id(lis[i] / p)][0] - (k - 1),
            G[i][1] -= (long long)p * (G[id(lis[i] / p)][1] - Gprime[k - 1][0]),
            G[i][2] -= s * (G[id(lis[i] / p)][2] - Gprime[k - 1][1]);
    }
    for(register int i = 1;i <= tot;++i)
        Fprime[i] = 3 * G[i][2] - 3 * G[i][1] + G[i][0];
    printf("%u\n",F(1,n) + 1U);
}