A simple Blog for wyx I've been down the bottle hoping.
真·杜教筛 51nod 1237 1238
发表于: | 分类: Oi | 评论:0 | 阅读:59

以前刷的都是裸题然后根本不用推导背式子就行的那种

这两天开始每天做两道数学题提高智商

首先是这个

$\sum_{i=1}^n\sum_{j=1}^n \gcd(i,j)$

先推柿子

$=\sum_{d=1}d\sum_{i=1}^n\sum_{j=1}^n [\gcd(i,j)==d]$

$=\sum_{d=1}d\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{n}{d}}[\gcd(i,j)==1]$

$=\sum_{d=1}d(2*\varphi(\frac{n}{d})-1)$

然后根号分块求个$\varphi$的前缀和就行了

#include <map>
#include <stdio.h>
#include <string.h>
#include <iostream>
#include <algorithm>
using namespace std;
const int mod = 1e9+7;
const int inv2 = 500000004;
const int N = 1e7+5;
const int Max = 1e7;
typedef long long LL;

inline int calc_1(LL x) {
    x %= mod;
    return x * (x+1) % mod * inv2 % mod;
}

int prime[N], cnt, phi[N];
bool F[N];

inline void init() {
    phi[1] = 1;
    for(int i=2;i<=Max;++i) {
        if(!F[i]) {
            prime[++cnt] = i;
            phi[i] = i - 1;
        }
        for(int j = 1; prime[j] * i <= Max; ++ j) {
            F[i*prime[j]] = 1;
            if(i % prime[j] == 0) {
                phi[i*prime[j]] = phi[i] * prime[j];
                break;
            } 
            phi[i*prime[j]] = phi[i] * (prime[j]-1);
        }
    }
    for(int i=2;i<=Max;++i) (phi[i] += phi[i-1]) %= mod;
}

map <LL,int> mp;

inline int calc_s(LL x) {
    if(x <= Max) return phi[x];
    if(mp.count(x)) return mp[x];
    LL temp1 = calc_1(x);
    LL temp2 = 0, last = 0;
    for(LL i = 2; i <= x; i = last + 1) {
        last = x / (x/i);
        (temp2 += 1ll * (last - i + 1) % mod * calc_s(x / i) % mod) %= mod;
    }
    temp1 = temp1 - temp2;
    temp1 += mod;
    temp1 %= mod;
    return mp[x] = temp1;
}

inline int calc(LL x) {
    LL temp1 = 0, temp2 = calc_1(x), last = 0;
    for(LL i = 1; i <= x; i = last + 1) {
        last = x / (x/i);
        (temp1 += 1ll * (calc_1(last) - calc_1(i-1) + mod) % mod * calc_s(x / i) % mod) %= mod;
    }
    temp1 = temp1 * 2 % mod;
    temp1 = temp1 - temp2 + mod;
    temp1 %= mod;
    return temp1;
}

int main() {
    init();
    LL x;
    cin >> x;
    cout << calc(x) << endl;
}

然后是这个

$\sum_{i=1}^n\sum_{j=1}^n {Lcm(i,j)}$

$=\sum_{i=1}^n\sum_{j=1}^n{\frac{ij}{\gcd(i,j)}}$

$=\sum_{d=1}^n\frac{1}{d}\sum_{i=1}^n\sum_{j=1}^n{ij[\gcd(i,j)==d]}$

$=\sum_{d=1}^n\frac{1}{d}\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{n}{d}}{d^2\times ij[\gcd(i,j)==d]}$

$=\sum_{d=1}^n d\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{n}{d}} ij[\gcd(i,j)==1]$

$=\sum_{d=1}^n d\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{n}{d}} j[\gcd(i,j)==1]$

$=\sum_{d=1}^n d\sum_{i=1}^{\frac{n}{d}}i^2\varphi(i)$

$=\sum_{d=1}^n\sum_{i=1}^{\frac{n}{d}}i^2 \varphi(i)$

然后就是套路的卷积之后杜教筛咯……

$f(x)=x*x,g(x)=x^2\varphi(x)$

$f*g(n)=\sum_{x=1}^n\sum_{d|x}f(d)g(\frac{x}{d})$

$=\sum_{x=1}\sum_{d|x}d^2(\frac{x}{d})^2\varphi(\frac{x}{d})$

$=\sum_{x=1}^n x*x \sum_{d|x}\varphi(\frac{x}{d})$

$=\sum_{x=1}^n x^3$

所以有$\sum_{x=1}^n x^3 = \sum_{i=1}^n S(\frac{n}{i})$

$S(n)=\sum_{x=1}^n x^3 - \sum_{i=2}^n S(\frac{n}{i})$

呼结束啦

#include <map>
#include <stdio.h>
#include <string.h>
#include <iostream>
#include <algorithm>
using namespace std;
typedef long long LL;
const LL N = 1e7+5;
const LL Max = 1e7;
const LL mod = 1000000007;
const LL inv4 = 250000002;
const LL inv6 = 166666668;
const LL inv2 = 500000004;
map <LL,LL> mp;

LL phi[N], prime[N], cnt;
bool F[N];

inline void init() {
    phi[1] = 1;
    for(LL i=2;i<=Max;++i) {
        if(!F[i]) {
            prime[++cnt] = i;
            phi[i] = i - 1;
        }
        for(LL j = 1; prime[j] * i <= Max; ++j) {
            F[prime[j]*i] = 1;
            if(i%prime[j] == 0) {
                phi[prime[j]*i] = prime[j] * phi[i];
                break;
            }
            phi[prime[j]*i] = phi[i] * (prime[j]-1);
        }
    }
    for(LL i=1;i<=Max;++i) phi[i] = (LL)phi[i] * i % mod * i % mod;
    for(LL i=1;i<=Max;++i) (phi[i] += phi[i-1])%=mod;
}

inline LL calc_1(LL x) {
    x %= mod;
    return x * (x+1) % mod * inv2 % mod;
}

inline LL calc_2(LL x) {
    x %= mod;
    return x * (x+1) % mod * (2*x+1) % mod * inv6 % mod;
}

inline LL calc_3(LL x) {
    x %= mod;
    return x * (x+1) % mod * x % mod * (x+1) % mod * inv4 % mod;
}

inline LL calcs(LL x) {
    if(x <= Max) return phi[x];
    if(mp.count(x)) return mp[x];
    LL temp1 = calc_3(x), temp2 = 0;
    LL last = 0;
    for(LL i = 2; i <= x; i = last + 1) {
        last = x / (x/i);
        (temp2 += 1ll*((calc_2(last)-calc_2(i-1)+mod)%mod)*calcs(x/i)%mod) %= mod;
    }
    temp1 = (temp1 - temp2 + mod) % mod;
    return mp[x]= temp1;
}

inline LL calc(LL x) {
    LL last, temp = 0;
    for(LL i = 1; i <= x; i = last + 1) {
        last = x / (x/i);
        (temp += 1ll*((calc_1(last) - calc_1(i-1) + mod) % mod) * calcs(x/i) % mod) %= mod;
    }
    return temp;
}

int main() {
    init();
    LL x;
    cin >> x;
    cout << calc(x) << endl;
}

Title - Artist
0:00

站点地图 网站地图
Copyright © 2015-2017 A simple Blog for wyx
Powered by Typecho自豪的采用Sgreen主题

TOP