「51nod 1847」奇怪的数学题

令 $\textrm{sgcd} (i, j) $ 为 $i,j$ 的次大公约数,计算:

$$\sum\limits_{i = 1}^n\sum\limits_{j = 1} ^n \textrm{sgcd} (i, j)^k$$

$n \leq 10^9,k \leq 50。$

令 $f(d)$ 为 $d$ 除以其最小质因子的值,特别的,$f(1) =0$。

简单变形得到:

$$\textrm {ans} = \sum\limits_{d = 1}^n f(d)^k\sum\limits_{i = 1}^{\lfloor \frac{n}{d} \rfloor} \mu (i) \lfloor \frac{n}{di} \rfloor ^2$$

利用杜老师筛,第二个 $\sum$ 后面的部分用 $n^{\frac{3}{4}}$ 的时间计算,考虑怎么快速计算前半部分。

容易发现这个部分可以方便的用埃筛计算,但是单次计算时间复杂度太高。

我们考虑埃筛计算的预处理部分,我们直接利用预处理中容斥的过程预处理出答案即可,可以理解为每次枚举最小的质因子,具体细节参考代码。

另外一个细节就是因为模数不是质数,无法利用插值得到自然数幂和,考虑到 $k$ 不大,利用第二类斯特林数求自然数幂和即可。

因为不难理解,直接给出推导:

$$\sum\limits_{i = 1}^n i^m = \sum\limits_{i = 1}^n \sum\limits_{j = 0}^m S(m, j) j!C(i, j) $$

$$ = \sum\limits_{j = 0}^m S(m, j)j!\sum\limits_{i = 0}^nC(i,j)$$

$$= \sum\limits_{j = 0}^m S(m, j)j!C(n + 1, j + 1)$$

即可在单次 $\mathcal{O} (m^2)$ 的时间内不需要求逆元计算前缀和。

总复杂度 $\mathcal{O} (\sqrt n k^2 + n^{\frac{3}{4}} + n^{0.7})$。

#include <bits/stdc++.h>

using namespace std;

#define ll             long long
#define db            double
#define uint         unsigned int
#define up(i,j,n)        for (int i = j; i <= n; i++)
#define down(i,j,n)    for (int i = j; i >= n; i--)
#define cadd(a,b)        a = add (a, b)
#define cpop(a,b)        a = pop (a, b)
#define cmul(a,b)        a = mul (a, b)
#define pr            pair<int, int>
#define fi            first
#define se            second
#define SZ(x)        (int)x.size()
#define bin(i)        (1 << (i))
#define Auto(i,node)    for (int i = LINK[node]; i; i = e[i].next)

template<typename T> inline bool cmax(T & x, T y){return y > x ? x = y, true : false;}
template<typename T> inline bool cmin(T & x, T y){return y < x ? x = y, true : false;}
template<typename T> inline T dmax(T x, T y){return x > y ? x : y;}
template<typename T> inline T dmin(T x, T y){return x < y ? x : y;}

const int MAXN = 1e6 + 5;
const int MAXM = 4e4 + 5;
const int L = 1e6;

int N, K, prime[MAXN], cnt;

uint pk[MAXN];

bool vis[MAXN];

namespace calc_miu {

    uint f[MAXN], g[MAXN];
    bool iscalc[MAXN];

    uint get(int n){
        if (n <= L) return f[n];
        if (iscalc[N / n]) return g[N / n];
        iscalc[N / n] = 1;
        uint sum = 1;
        for (int cl = 2, cr; cl <= n; cl = cr + 1){
            cr = n / (n / cl);
            sum -= get(n / cl) * (cr - cl + 1);
        }
        return g[N / n] = sum; 
    }

    void main(){
        f[1] = 1; 
        up (i, 2, L) {
            if (!vis[i]) {
                prime[++cnt] = i;
                pk[cnt] = 1;
                up (j, 1, K) pk[cnt] *= i;
                f[i] = -1;
            }
            up (j, 1, cnt) {
                int x = i * prime[j];
                if (x > L) break;
                vis[x] = 1;
                if (i % prime[j] == 0) break;
                f[x] = -f[i]; 
            }
        }
        up (i, 1, L) f[i] += f[i - 1];
    }
}

namespace calc_f {

    uint S[105][105];

    uint f0[MAXM], f1[MAXM], g0[MAXM], g1[MAXM], ff[MAXM], gg[MAXM];

    int m;

    uint calc(int n){
        uint sum = 0;
        up (i, 0, K) {
            int x = (n + 1) / (i + 1) * (i + 1);
            uint cur = x / (i + 1);
            up (j, 0, i) cur *= (n + 1 - j) == x ? 1 : (n + 1 - j);
            sum += cur * S[K][i];
        }
        return sum;
    }

    inline uint get(int n) { return n <= m ? ff[n] : gg[N / n]; }

    void main(){
        S[0][0] = 1;
        up (i, 0, K) up (j, 1, i) S[i][j] = S[i - 1][j - 1] + S[i - 1][j] * j;
        m = sqrt(N + 0.5);
        up (i, 1, m) {
            f0[i] = i - 1; f1[i] = calc(i) - 1;
            g0[i] = N / i - 1; g1[i] = calc(N / i) - 1;
        }
        up (o, 1, cnt) {
            if (prime[o] > m) break;
            uint x = pk[o], d0 = f0[prime[o] - 1], d1 = f1[prime[o] - 1];
            int mi = m / prime[o];
            up (i, 1, mi) {
                g0[i] -= g0[i * prime[o]] - d0;
                g1[i] -= (g1[i * prime[o]] - d1) * x;
                gg[i] += g1[i * prime[o]] - d1;
            }
            int tmp = N / prime[o], mx = dmin(N / prime[o] / prime[o], m);
            up (i, mi + 1, mx) {
                g0[i] -= f0[tmp / i] - d0;
                g1[i] -= (f1[tmp / i] - d1) * x;
                gg[i] += f1[tmp / i] - d1;
            }
            if ((ll)prime[o] * prime[o] <= m) {
                mx = prime[o] * prime[o];
                down (i, m, mx) {
                    f0[i] -= f0[i / prime[o]] - d0;
                    f1[i] -= (f1[i / prime[o]] - d1) * x;
                    ff[i] += f1[i / prime[o]] - d1;
                }
            }
        }
        up (i, 1, m) ff[i] += f0[i], gg[i] += g0[i];
    }    

}

uint calc(int n){
    uint sum = 0, cur = 0, last = 0;
    for (int cl = 1, cr; cl <= n; cl = cr + 1) {
        cr = n / (n / cl);
        cur = calc_miu::get(cr);
        sum += (uint)(n / cl) * (n / cl) * (cur - last);
        last = cur;
    }
    return sum;
}

int main(){
#ifdef cydiater
    freopen("input.in", "r", stdin);
#endif
    scanf("%d%d", &N, &K);
    calc_miu::main();
    calc_f::main();
    uint ans = 0, last = 0, cur = 0;
    for (int cl = 1, cr; cl <= N; cl = cr + 1) {
        cr = N / (N / cl);
        cur = calc_f::get(cr); 
        ans += (cur - last) * calc(N / cl);
        last = cur;
    }
    cout << ans << endl;
    return 0;
}