• 【UR #5】怎样跑得更快 神秘做法


    链接

    自己一直没有想到 \(\gcd\) 技巧,想了一种神秘的方法,记一下并且补一下证明。

    首先是个特殊矩阵高斯消元。

    \(A_{i, j} = \gcd(i, j)^{c-d} \times i^d \times j^d\)

    我的想法是利用性质快速高斯消元,因为相当于那个矩阵系数给你了,你可以手玩,对于每个 \(k\),要求出 \(x_k\),那么可以选若干行,使得他们的线性组合恰好第 \(k\) 列的系数不是 \(0\),剩下基本是 \(0\)

    先考虑 \(c = 1, d = 0\) 的情况,这时候系数矩阵就是 \(\gcd(i, j)\)

    然后我手玩了 \(k = n\) 的情况,如果 \(k\) 是质数,设 \(A_i\) 为第 \(i\) 行,那么 \(A_k - A_1\) 即可解决。

    我就猜是不是类似可以套一个莫比乌斯容斥一下,考虑设置这样的一个东西 \(\sum_{i|k} A_{k/i} \times \mu(i)\)

    猜的啊,发现这个玩意,只有 \(ik\) 项系数不是 \(0\),并且第 \(ik\) 项的系数恰好是第 \(k\) 项的 \(i\) 倍,那么就可以递推了!

    考虑 \(c = n, d = 0\) 的情况,我发现还是满足上面性质的!

    然后考虑 \(d \not =0\),我稍微试了一下 \(\sum_{i|k} A_{k/i} \times \mu(i) \times i^d\) 发现也是对的!也就找到了普使的规律!

    这为啥对呢?试着感性理解一下。

    看了 vfk 题解感觉非常深刻!你看成一个函数 \(A_{i, j} = f(\gcd(i, j)) \times i^d \times h(j)\)

    考虑对于一个 \(k\) 时,第 \(x\) 列的系数。

    这个 \(h(x)\) 就是常量了,不用管。

    考虑一个 \(t = \gcd(k/i, x)\) 的贡献,假设不是 \(0\),那么你考虑

    首先如果 \(t = k\),那么只能是 \(x = bk (b \in N^+), i = 1\)。即是倍数,贡献是 \(1\)

    否则 \(t < k\),考虑因为有贡献,所以对于每个质因子,如果没被 \(x\) 限制一定得是 \(k / i\) 限制,这样就得到了最小的 \(i\),肯定有贡献,然后你考虑肯定存在一个质因子现在这个 \(k / i_{min}\) 可以再将一格,否则 \(t\) 就比 \(k\) 大了,这样就利用 \(\mu\) 的贡献抵消变成 \(0\) 了。

    然后加入 \(i^d\) 的影响,你考虑那个 \(\sum\)\(\times i^d\),事实上是补齐了这个东西都是 \(k^d\)

    所以我这个东西只在 vfk 模型中的 \(g(i)\) 完全积性才可以,感觉非常的不实用。

    不过我能想出来这么奇怪的东西感觉很神秘,想了 1h。

    这个式子可以抽象等价为 \([k | x] = [\sum_{i | k} \gcd(k/i,x) \times \mu(i) \not= 0]\) 是啥子呢?

    复杂度 \(O(n \log n)\)

    // Skyqwq
    #include <bits/stdc++.h>
    
    #define pb push_back
    #define fi first
    #define se second
    #define mp make_pair
    
    using namespace std;
    
    typedef pair<int, int> PII;
    typedef long long LL;
    
    template <typename T> bool chkMax(T &x, T y) { return (y > x) ? x = y, 1 : 0; }
    template <typename T> bool chkMin(T &x, T y) { return (y < x) ? x = y, 1 : 0; }
    
    template <typename T> void inline read(T &x) {
        int f = 1; x = 0; char s = getchar();
        while (s < '0' || s > '9') { if (s == '-') f = -1; s = getchar(); }
        while (s <= '9' && s >= '0') x = x * 10 + (s ^ 48), s = getchar();
        x *= f;
    }
    
    const int N = 1e5 + 5, P = 998244353;
    
    int n, c, d, q, b[N];
     
    bool st[N];
    
    int mu[N], pr[N], t;
    
    vector<int> D[N];
    
    int gcd(int a, int b) {
        return b ? gcd(b, a % b) : a;
    }
    
    int inline power(int a, int b) {
        if (b < 0) return power(power(a, P - 2), -b);
        int res = 1;
        while (b) {
            if (b & 1) res = (LL)res * a % P;
            a = (LL)a * a % P;
            b >>= 1;
        }
        return res;
    }
    
    void inline prework(int n) {
        for (int i = 1; i <= n; i++) mu[i] = 1;
        for (int i = 2; i <= n; i++) {
            if (!st[i]) pr[++t] = i, mu[i] = -1;
            for (int j = 1; pr[j] * i <= n; j++) {
                st[i * pr[j]] = 1;
                if (i % pr[j] == 0) {
                    mu[i * pr[j]] = 0;
                    break;
                }
                mu[i * pr[j]] = -mu[i];
            }
        }
        for (int i = 1; i <= n; i++) {
            if (mu[i]) for (int j = i; j <= n; j += i) D[j].pb(i);
        }
    }
    
    void inline solve1() {
        for (int i = 2; i <= n; i++) {
            if (b[i] != b[1]) {
                puts("-1");
                return;
            }
        }
        printf("%d", b[1]);
        for (int i = 2; i <= n; i++) printf(" %d", 0);
        puts("");
    }
    
    void inline add(int &x, int y) {
        x += y;
        if (x >= P) x -= P;
    }
    
    void inline del(int &x, int y) {
        x -= y;
        if (x < 0) x += P;
    }
    
    int ans[N];
    
    int pc[N], pd[N];
    
    int inline A(int i, int j) {
        return 1ll * pc[i] * pd[j] % P;
    }
    
    void inline solve2() {
        for (int k = n; k; k--) {
            int u = 0, v = 0;
            for (int i: D[k]) {
                int xs = 1ll * mu[i] * pd[i] % P;
                if (xs < 0) xs += P;
                add(u, 1ll * xs * A(k / i, k) % P);
                add(v, 1ll * xs * b[k / i] % P);
            }
            for (int z = 2; z * k <= n; z++) {
                int zx = 1ll * u * pd[z] % P;
                del(v, 1ll * ans[z * k] * zx % P);
            }
            if (!u) {
                if (!v) continue;
                else {
                    puts("-1");
                    return;
                }
            }
            ans[k] = 1ll * v * power(u, P - 2) % P;
        }
        for (int i = 1; i <= n; i++) printf("%d ", ans[i]);
        puts("");
    }
    
    int main() {
        read(n), read(c), read(d), read(q);
        prework(n);
        for (int i = 1; i <= n; i++) pc[i] = power(i, c), pd[i] = power(i, d);
        while (q--) {
            for (int i = 1; i <= n; i++) read(b[i]);
            solve2();
        }
        return 0;
    }
    
  • 相关阅读:
    Window PHP 使用命令行模式
    LNMP ftp 可以登录无权限操作?
    linux 允许mysql用户远程访问
    解决报错:scandir() has been disabled for security reasons
    LNMP 配置二级域名
    MUI 图片上传剪切预览,可选(拍照+系统相册)
    MUI 单个图片上传预览(拍照+系统相册):先选择->预览->上传提交
    MUI 单图片压缩上传(拍照+系统相册): 选择立即上传
    循环递归的区别?
    如何让自己的广播只让指定的 app 接收?
  • 原文地址:https://www.cnblogs.com/dmoransky/p/15820481.html
Copyright © 2020-2023  润新知