【51Nod 1238】最小公倍数之和 V3

http://www.51nod.com/onlineJudge/questionCode.html#!problemId=1238
\(A(n)=\sum\limits_{i=1}^n\frac{in}{(i,n)}\),则\(ans=\sum\limits_{i=1}^n\left(2A(i)-i\right)\)
\[ \begin{aligned} A(n)=&n\sum_{d|n}\sum_{i=1}^{\frac nd}i\left[\left(i,\frac nd\right)=1\right]\\ =&\frac n2\sum_{d|n}\left(\varphi(d)d+[d=1]\right)\\ =&\frac n2\sum_{d|n}\varphi(d)d+\frac n2\\ ans=&\sum_{i=1}^n\left(2A(i)-i\right)\\ =&\sum_{i=1}^ni\sum_{d|i}\varphi(d)d\\ =&\sum_{i=1}^ni\sum_{d=1}^{\left\lfloor\frac ni\right\rfloor}\varphi(d)d^2\\ \end{aligned} \]
\(S(n)=\sum\limits_{i=1}^n\varphi(i)i^2\),求出\(O\left(\sqrt n\right)\)不同下取整取值的\(S\)就可以算出答案了,所以现在重点是杜教筛\(S(n)\)
先让\(f(n)=\varphi(n)n^2\)卷上\(f(n)=n^2\)
\[\sum_{i=1}^{n}\sum_{d|i}\varphi(d)d^2\left(\frac id\right)^2=\sum_{i=1}^ni^2\sum_{d|i}\varphi(d)=\sum_{i=1}^ni^3\]
\[ \begin{aligned} \sum_{i=1}^ni^3=&\sum_{i=1}^{n}\sum_{d|i}\varphi(d)d^2\left(\frac id\right)^2\\ =&\sum_{i=1}^ni^2\sum_{d=1}^{\left\lfloor\frac ni\right\rfloor}\varphi(d)d^2\\ =&\sum_{i=1}^ni^2S\left(\left\lfloor\frac ni\right\rfloor\right)\\ \end{aligned} \]
如果要求\(S(n)\)\(S(n)=\sum\limits_{i=1}^ni^3-\sum\limits_{i=2}^ni^2S\left(\left\lfloor\frac ni\right\rfloor\right)\),时间复杂度\(O\left(n^{\frac 23}\right)\)
PS:\(\sum\limits_{i=1}^ni^3=\left(\frac{n(n+1)}{2}\right)^2\)

#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;

const int N = 100003;
const int B = 4641589;
const int p = 1000000007;
const int ni2 = 500000004;
const int ni6 = 166666668;

ll n;
bool notp[B + 1];
int phi[B + 1], sum[B + 1], num = 0, prime[B + 1];

void Euler_shai() {
    sum[1] = phi[1] = 1;
    for (int i = 2; i <= B; ++i) {
        if (!notp[i]) prime[++num] = i, phi[i] = i - 1;
        for (int j = 1; j <= num && prime[j] * i <= B; ++j) {
            notp[prime[j] * i] = true;
            if (i % prime[j] == 0) {
                phi[prime[j] * i] = prime[j] * phi[i];
                break;
            } else
                phi[prime[j] * i] = (prime[j] - 1) * phi[i];
        }
        sum[i] = (1ll * phi[i] * i % p * i % p + sum[i - 1]) % p;
    }
}

int ps[B];
ll ndx;
int S(ll x) {return (ndx = n / x) <= B ? sum[ndx] : ps[x];}

void DJ_shai() {
    for (ll i = n, y; i >= 1; i = n / (y + 1)) {
        y = n / i;
        if (y <= B) continue;
        int &s = ps[i];
        s = y % p * (y % p) % p * ((y + 1) % p) % p * ((y + 1) % p) % p * ni2 % p * ni2 % p;
        for (ll j = 2, pre = 1, spre = 1, sj; j <= y; spre = sj, pre = j, ++j) {
            j = y / (y / j);
            sj = j % p * ((j + 1) % p) % p * ((j * 2 + 1) % p) % p * ni6 % p;
            ((s -= (sj - spre + p) % p * S(i * j) % p) += p) %= p;
        }
    }
}

main() {
    scanf("%lld", &n);
    Euler_shai();
    
    int ans = 0;
    DJ_shai();
    
    for (ll i = 1, pre = 0; i <= n; pre = i, ++i) {
        i = n / (n / i);
        (ans += (i + pre + 1) % p * ((i - pre) % p) % p * ni2 % p * S(i) % p) %= p;
    }
    
    printf("%d\n", ans);
    return 0;
}

转载于:https://www.cnblogs.com/abclzr/p/6770714.html

  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值