登录
当前位置： 首页 > 工具软件 > misaka > 使用案例 >

「LibreOJ Round #11」Misaka Network 与求和（杜教筛 + Min_25）

晏经武
2023-12-01

#572. 「LibreOJ Round #11」Misaka Network 与求和

推式子

∑ i = 1 n ∑ j = 1 n f ( g c d ( i , j ) ) k ∑ d = 1 n f ( d ) k ∑ i = 1 n d ∑ j = 1 n d [ g c d ( i , j ) = 1 ] ∑ d = 1 n f ( d ) k ∑ K = 1 n d μ ( k ) ( n K d ) 2 t = K d ∑ t = 1 n ( n t ) 2 ∑ d ∣ t f ( d ) k μ ( t d ) 我 们 记 f ( x ) k = F ( x ) 上 面 式 子 后 半 部 分 是 一 个 迪 利 克 雷 卷 积 形 式 : F ∗ μ 所 以 我 们 卷 上 一 个 I ， 有 F ∗ μ ∗ I = F ∗ ϵ = F 得 到 后 半 部 分 的 前 缀 和 S ( n ) = ∑ i = 1 n F ( i ) − ∑ i = 2 n S ( n i ) \sum_{i = 1} ^{n} \sum_{j = 1} ^{n} f(gcd(i, j)) ^ k\\ \sum_{d = 1} ^{n} f(d) ^k \sum_{i = 1} ^{\frac{n}{d}} \sum_{j = 1} ^{\frac{n}{d}}[gcd(i, j) = 1]\\ \sum_{d = 1} ^{n} f(d) ^k \sum_{K = 1} ^{\frac{n}{d}} \mu(k) \left( \frac{n}{Kd} \right) ^2\\ t = Kd\\ \sum_{t = 1} ^{n} \left(\frac{n}{t} \right) ^ 2 \sum_{d \mid t} f(d) ^ k \mu(\frac{t}{d})\\ 我们记f(x) ^ k = F(x)\\ 上面式子后半部分是一个迪利克雷卷积形式:F * \mu\\ 所以我们卷上一个I，有F * \mu * I = F * \epsilon = F\\ 得到后半部分的前缀和S(n) = \sum_{i = 1} ^{n} F(i) - \sum_{i = 2} ^{n} S(\frac{n}{i})\\ i=1∑n​j=1∑n​f(gcd(i,j))kd=1∑n​f(d)ki=1∑dn​​j=1∑dn​​[gcd(i,j)=1]d=1∑n​f(d)kK=1∑dn​​μ(k)(Kdn​)2t=Kdt=1∑n​(tn​)2d∣t∑​f(d)kμ(dt​)我们记f(x)k=F(x)上面式子后半部分是一个迪利克雷卷积形式:F∗μ所以我们卷上一个I，有F∗μ∗I=F∗ϵ=F得到后半部分的前缀和S(n)=i=1∑n​F(i)−i=2∑n​S(in​)

化简到这里只需要跟上面一题类似用Min_25求 ∑ i = 1 n F ( i ) \sum\limits_{i = 1} ^{n} F(i) i=1∑n​F(i)，然后用杜教筛求 S ( n ) S(n) S(n)即可得到答案。

代码

/*
  Author : lifehappy
*/
#pragma GCC optimize(2)
#pragma GCC optimize(3)
#include <bits/stdc++.h>

#define mp make_pair
#define pb push_back
#define endl '\n'
#define mid (l + r >> 1)
#define lson rt << 1, l, mid
#define rson rt << 1 | 1, mid + 1, r
#define ls rt << 1
#define rs rt << 1 | 1

using namespace std;

typedef long long ll;
typedef unsigned long long ull;
typedef pair<int, int> pii;

const double pi = acos(-1.0);
const double eps = 1e-7;
const int inf = 0x3f3f3f3f;

inline ll read() {
    ll f = 1, x = 0;
    char c = getchar();
    while(c < '0' || c > '9') {
        if(c == '-')    f = -1;
        c = getchar();
    }
    while(c >= '0' && c <= '9') {
        x = (x << 1) + (x << 3) + (c ^ 48);
        c = getchar();
    }
    return f * x;
}

#define uint unsigned int

const int N = 1e6 + 10;

uint quick_pow(uint a, int n) {
    uint ans = 1;
    while(n) {
        if(n & 1) ans = ans * a;
        a = a * a;
        n >>= 1;
    }
    return ans;
}

namespace Min_25 {
    uint prime[N], g[N], sum[N], f[N], calc[N];

    int a[N], id1[N], id2[N], n, m, k, cnt, T;

    bool st[N];

    int ID(int x) {
        return x <= T ? id1[x] : id2[n / x];
    }

    void init() {
        cnt = m = 0;
        T = sqrt(n + 0.5);
        for(int i = 2; i <= T; i++) {
            if(!st[i]) {
                prime[++cnt] = i;
                f[cnt] = quick_pow(i, k);
                sum[cnt] = sum[cnt - 1] + 1;
            }
            for(int j = 1; j <= cnt && 1ll * i * prime[j] <= T; j++) {
                st[i * prime[j]] = 1;
                if(i % prime[j] == 0) {
                    break;
                }
            }
        }
        for(ll l = 1, r; l <= n; l = r + 1) {
            r = n / (n / l);
            a[++m] = n / l;
            if(a[m] <= T) id1[a[m]] = m;
            else id2[n / a[m]] = m;
            g[m] = a[m] - 1;
        }
        for(int j = 1; j <= cnt; j++) {
            for(int i = 1; i <= m && 1ll * prime[j] * prime[j] <= a[i]; i++) {
                g[i] -= g[ID(a[i] / prime[j])] - sum[j - 1];
            }
        }
        for(int i = 1; i <= T; i++) {
            st[i] = 0;
        }
        /*
            非递归版本
        */

        // for(int j = cnt; j >= 1; j--) {
        //     for(int i = 1; i <= m && 1ll * prime[j] * prime[j] <= a[i]; i++) {
        //         for(ll k = prime[j]; k * prime[j] <= a[i]; k *= prime[j]) {
        //             calc[i] += calc[ID(a[i] / k)] + (g[ID(a[i] / k)] - sum[j - 1]) * f[j];
        //         }
        //     }
        // }
    }

    // uint solve(int x) {
    //     if(x <= 1) return 0;
    //     return calc[ID(x)] + g[ID(x)];
    // }

    /*
    下面是递归版本
    */
    // uint solve(int n, int m) {
    //     if(n < prime[m] || n <= 1) return 0;
    //     uint ans = 0;
    //     for(int j = m; j <= cnt && 1ll * prime[j] * prime[j] <= n; j++) {
    //         for(ll i = prime[j]; i * prime[j] <= n; i *= prime[j]) {
    //             ans += solve(n / i, j + 1) + (g[ID(n / i)] - sum[j - 1]) * f[j];
    //         }
    //     }
    //     return ans;
    // }

    // uint solve(int n) {
    //     if(n <= 1) return 0;
    //     return solve(n, 1) + g[ID(n)];
    // }
}

unordered_map<int, uint> ans_s;

uint S(int n) {
    if(ans_s.count(n)) return ans_s[n];
    uint ans = Min_25::solve(n);
    for(uint l = 2, r; l <= n; l = r + 1) {
        r = n / (n / l);
        ans -= (r - l + 1) * S(n / l);
    }
    return ans_s[n] = ans;
}

int main() {
    // freopen("in.txt", "r", stdin);
    // freopen("out.txt", "w", stdout);
    // ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
    uint n = read(), k = read();
    Min_25::n = n, Min_25::k = k;
    Min_25::init();
    uint ans = 0;
    for(uint l = 1, r; l <= n; l = r + 1) {
        r = n / (n / l);
        ans += (n / l) * (n / l) * (S(r) - S(l - 1));
    }
    cout << ans << endl;
    return 0;
}

写法二

∑ i = 1 n ∑ j = 1 n f ( g c d ( i , j ) ) k ∑ d = 1 n f ( d ) k ∑ i = 1 n d ∑ j = 1 n d [ g c d ( i , j ) = = 1 ] ∑ d = 1 n f ( d ) k ( 2 ∑ i = 1 n d ϕ ( i ) − 1 ) \sum_{i = 1} ^{n} \sum_{j = 1} ^{n}f(gcd(i, j)) ^k\\ \sum_{d = 1} ^{n} f(d) ^ k \sum_{i = 1} ^{\frac{n}{d}} \sum_{j = 1} ^{\frac{n}{d}}[gcd(i, j) == 1]\\ \sum_{d = 1} ^{n} f(d) ^ k(2 \sum_{i = 1} ^ {\frac{n}{d}} \phi(i) - 1)\\ i=1∑n​j=1∑n​f(gcd(i,j))kd=1∑n​f(d)ki=1∑dn​​j=1∑dn​​[gcd(i,j)==1]d=1∑n​f(d)k(2i=1∑dn​​ϕ(i)−1)
充分利用上面递推Min_25，得到一个复杂度更优的算法。

/*
  Author : lifehappy
*/
#pragma GCC optimize(2)
#pragma GCC optimize(3)
#include <bits/stdc++.h>

#define mp make_pair
#define pb push_back
#define endl '\n'
#define mid (l + r >> 1)
#define lson rt << 1, l, mid
#define rson rt << 1 | 1, mid + 1, r
#define ls rt << 1
#define rs rt << 1 | 1

using namespace std;

typedef long long ll;
typedef unsigned long long ull;
typedef pair<int, int> pii;

const double pi = acos(-1.0);
const double eps = 1e-7;
const int inf = 0x3f3f3f3f;

inline ll read() {
    ll f = 1, x = 0;
    char c = getchar();
    while(c < '0' || c > '9') {
        if(c == '-')    f = -1;
        c = getchar();
    }
    while(c >= '0' && c <= '9') {
        x = (x << 1) + (x << 3) + (c ^ 48);
        c = getchar();
    }
    return f * x;
}

#define uint unsigned int

const int N = 1e6 + 10;

uint quick_pow(uint a, int n) {
    uint ans = 1;
    while(n) {
        if(n & 1) ans = ans * a;
        a = a * a;
        n >>= 1;
    }
    return ans;
}

namespace Min_25 {
    uint prime[N], g[N], sum[N], f[N], calc[N];

    int a[N], id1[N], id2[N], n, m, k, cnt, T;

    bool st[N];

    int ID(int x) {
        return x <= T ? id1[x] : id2[n / x];
    }

    void init() {
        cnt = m = 0;
        T = 1000000;
        for(int i = 2; i <= T; i++) {
            if(!st[i]) {
                prime[++cnt] = i;
                f[cnt] = quick_pow(i, k);
                sum[cnt] = sum[cnt - 1] + 1;
            }
            for(int j = 1; j <= cnt && 1ll * i * prime[j] <= T; j++) {
                st[i * prime[j]] = 1;
                if(i % prime[j] == 0) {
                    break;
                }
            }
        }
        for(ll l = 1, r; l <= n; l = r + 1) {
            r = n / (n / l);
            a[++m] = n / l;
            if(a[m] <= T) id1[a[m]] = m;
            else id2[n / a[m]] = m;
            g[m] = a[m] - 1;
        }
        for(int j = 1; j <= cnt; j++) {
            for(int i = 1; i <= m && 1ll * prime[j] * prime[j] <= a[i]; i++) {
                g[i] -= g[ID(a[i] / prime[j])] - sum[j - 1];
            }
        }
        for(int i = 1; i <= T; i++) {
            st[i] = 0;
        }
        for(int j = cnt; j >= 1; j--) {
            for(int i = 1; i <= m && 1ll * prime[j] * prime[j] <= a[i]; i++) {
                for(ll k = prime[j]; k * prime[j] <= a[i]; k *= prime[j]) {
                    calc[i] += calc[ID(a[i] / k)] + (g[ID(a[i] / k)] - sum[j - 1]) * f[j];
                }
            }
        }
    }

    uint solve(int x) {
        if(x <= 1) return 0;
        return calc[ID(x)] + g[ID(x)];
    }
}

namespace Djs {
    uint prime[N], phi[N], cnt;

    bool st[N];

    void init() {
        phi[1] = 1;
        for(int i = 2; i < N; i++) {
            if(!st[i]) {
                prime[++cnt] = i;
                phi[i] = i - 1;
            }
            for(int j = 1; j <= cnt && 1ll * i * prime[j] < N; j++) {
                st[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 = 1; i < N; i++) {
            phi[i] += phi[i - 1];
        }
    }

    unordered_map<int, uint> ans_s;

    uint S(int n) {
        if(n < N) return phi[n];
        if(ans_s.count(n)) return ans_s[n];
        uint ans = 1ll * n * (n + 1) / 2;
        for(uint l = 2, r; l <= n; l = r + 1) {
            r = n / (n / l);
            ans -= (r - l + 1) * S(n / l);
        }
        return ans_s[n] = ans;
    }
}

int main() {
    // freopen("in.txt", "r", stdin);
    // freopen("out.txt", "w", stdout);
    // ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
    uint n = read(), k = read();
    Min_25::n = n, Min_25::k = k;
    Min_25::init();
    Djs::init();
    uint ans = 0;
    for(uint l = 1, r; l <= n; l = r + 1) {
        r = n / (n / l);
        ans += (2 * Djs::S(n / l) - 1) * (Min_25::solve(r) - Min_25::solve(l - 1));
    }
    cout << ans << endl;
    return 0;
}
 类似资料：

相关阅读

设置杜松子酒找不到的路线numpy：有效地与索引数组求和如何从Java筛选器获取请求URL？@FXML批注和FXMLLoader类未解析为Java 11和JavaFX 11中的类型Angular将并行和链接请求与$ http.then（）和$ q.all（）组合

相关文章

PyCharm下载和安装教程C++11的std::function和std::bind用法详解DTD示例11Tableau条件筛选器Spring MVC请求转发与重定向

相关问答

python案例求教？IntelliJ无法识别JavaFX 11和OpenJDK 11JavaFX 11不能与JDK 11、intellij idea一起工作配置Jetty 10/11请求日志Firebase工具和Java 11

相关文档

Biopython 教程与手册 v1.61+Python 和 Pandas 数据分析教程C++ 11 FAQ 中文版PHP 异步编程：手把手教你实现 co 与 KoaSwift 教程
快捷导航： 新手教程 算法原理 架构设计 Java进阶 数据库进阶 大厂专栏 面试经验 编程笔记 编程问答 所有专题 文档资料 工具软件 电子书籍 小牛导航
在线工具： 房贷计算器 个税计算器 Linux命令查询 Json格式化 正则表达式 颜色转换 AES加解密 SHA1加密 MD5加密 毒鸡汤 字数统计 随机密码生成 进制转换 Base64编解码 励志句子
Copyright © 2019-2024 小牛知识库@xnip.cn. All Rights Reserved.