2019 ICPC上海网络赛部分题解

C. Triple

考虑分情况不同做法。当 n \le 1000 时,枚举 A_i, B_j,则根据不等式解出满足条件的 C_k \in [|A_i-B_j|,A_i+B_j]。当 n > 1000,考虑枚举不合法的方案数,即形如 A_i + B_j < C_k 的方案数,不等式左边可以通过FFT取得,然后前缀和再枚举 C_k 容斥掉即可。

大力施加常数优化,AC仅需1256ms啦啦啦~

#include <bits/stdc++.h>
using namespace std;
typedef long long lld;
typedef complex<double> cplx;
int revs[1<<19], A[300010]; cplx wmks[1<<19];
const double PI = acos(-1.0);

void dft(cplx a[], int DFT, int N) {
    int *rev = revs + N;
    for (int i = 0; i < N; i++)
        if (i < revs[N+i]) swap(a[i], a[rev[i]]);
    for (int m = 2, m2 = 1; m <= N; m <<= 1, m2 <<= 1) {
        cplx *wmk = wmks + m + (DFT==-1 ? m2 : 0), u, t;
        for (int k = 0; k < N; k += m)
            for (int j = 0; j < m2; j++)
                t = wmk[j] * a[k+j+m2], u = a[k+j],
                a[k+j] = u + t, a[k+j+m2] = u - t;
    }
    if (DFT == -1) for (int i = 0; i < N; i++) a[i] /= N;
}

int a[100010], b[100010], c[100010];

lld solveWithFFT(int n)
{
    static cplx f[262144], g[262144], h[262144];
    static lld bc[262144], ac[262144], ab[262144];
    int ma = 0, mb = 0, mc = 0, md, len = 1; cplx tmp, tmp2;
    for (int i = 0; i < n; i++) ma = max(ma, a[i]), mb = max(mb, b[i]), mc = max(mc, c[i]);
    md = max(ma+mb, max(mb+mc, mc+ma)) + 2; while (len < md) len <<= 1;
    for (int i = 0; i < len; i++) f[i] = g[i] = h[i] = 0;
    for (int i = 0; i < n; i++) f[a[i]] += 1, g[b[i]] += 1, h[c[i]] += 1;
    dft(f, 1, len), dft(g, 1, len), dft(h, 1, len);
    for (int i = 0; i < len; i++) tmp = f[i] * g[i], tmp2 = g[i] * h[i], g[i] = f[i] * h[i], f[i] = tmp2, h[i] = tmp;
    dft(f, -1, len), dft(g, -1, len), dft(h, -1, len);
    for (int i = 0; i < len; i++) bc[i] = f[i].real()+0.2, ac[i] = g[i].real()+0.2, ab[i] = h[i].real()+0.2;
    for (int i = 1; i < len; i++) bc[i] += bc[i-1], ac[i] += ac[i-1], ab[i] += ab[i-1];
    lld ans = 1LL * n * n * n;

    for (int i = 0; i < n; i++)
    {
        if (c[i] <= len) ans -= ab[c[i]-1];
        if (b[i] <= len) ans -= ac[b[i]-1];
        if (a[i] <= len) ans -= bc[a[i]-1];
    }

    return ans;
}

inline int abs(int x) { return x < 0 ? -x : x; }

lld solveWithBruteForce(int n)
{
    for (int i = 0; i <= 200001; i++) A[i] = 0;
    for (int i = 0; i < n; i++) A[c[i]]++;
    for (int i = 1; i <= 200001; i++) A[i] += A[i-1];
    lld ans = 0;
    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
            ans += A[a[i]+b[j]] - A[max(abs(a[i]-b[j])-1,0)];
    return ans;
}

void solve(int cas)
{
    int n; scanf("%d", &n);
    for (int i = 0; i < n; i++) scanf("%d", &a[i]);
    for (int i = 0; i < n; i++) scanf("%d", &b[i]);
    for (int i = 0; i < n; i++) scanf("%d", &c[i]);
    printf("Case #%d: %lld\n", cas, n <= 1000 ? solveWithBruteForce(n) : solveWithFFT(n));
}

int main()
{
    for (int N = 2, N2 = 1; N <= (1<<18); N <<= 1, N2 <<= 1)
    {
        int *rev = revs + N; cplx *wmk = wmks + N;
        for (int i = 0; i < N; i++)
            rev[i] = (rev[i>>1]>>1)|((i&1)?N2:0);
        for (int i = 0; i < N2; i++)
            wmk[i] = cplx(cos(PI*i/N2), sin(PI*i/N2));
        for (int i = 0; i < N2; i++)
            wmk[i+N2] = cplx(wmk[i].real(), -wmk[i].imag());
    }

    int T; scanf("%d", &T);
    for (int i = 1; i <= T; i++) solve(i);
    return 0;
}

D. Counting Sequences I

首先发现是几个本质不同的方案的排列,于是想到打表。需要一个快速出结果的剪枝。

#include <bits/stdc++.h>
using namespace std;
typedef long long lld;
const int MOD = 1e9+7;
int n, jspx[3010], val[3010];
lld fac[3010], inv[3010], invs[3010], ans;

void dfs(int cur, int last, lld mul, int sum)
{
    if (cur == n)
    {
        if (mul == sum)
        {
            for (int i = 1; i <= n; i++)
                if (jspx[i])
                    printf("%dx%d ", jspx[i], i);
            printf("\n");
            lld tot = fac[n];
            for (int i = 1; i <= n; i++)
                tot = tot * invs[jspx[i]] % MOD;
            ans = (ans + tot) % MOD;
        }
    }
    else if (mul * pow(last, n-cur) <= sum + last * (n - cur))
    {
        for (int i = last; i <= n; i++)
        {
            val[cur] = i;
            jspx[i]++;
            dfs(cur + 1, i, mul * i, sum + i);
            jspx[i]--;
        }
    }
}

int main()
{
    fac[0] = fac[1] = invs[0] = invs[1] = inv[1] = 1;
    for (int i = 2; i < 3010; i++)
        fac[i] = i * fac[i-1] % MOD,
        inv[i] = (MOD-MOD/i) * inv[MOD%i] % MOD,
        invs[i] = invs[i-1] * inv[i] % MOD;
    scanf("%d", &n);
    dfs(0, 1, 1, 0);
    printf("ANS = %lld\n", ans);
    return 0;
}

E. Counting Sequences II

不知道题解在写什么鬼东西……

考虑 [1,m] 中的每个数字 t 的排列型生成函数。如果 t 为奇数,则其生成函数为 f(x) = 1 + \frac{x}{1!} + \frac{x^2}{2!} + \cdots = e^x;否则只能出现偶数次,则为 g(x) = 1 + \frac{x^2}{2!} + \frac{x^4}{4!} + \cdots = \frac{e^x + e^{-x}}{2} = \ch(x)

那么答案相当于求

\left(\frac{e^x+e^{-x}}{2}\right)^{\left\lfloor \frac{m}{2} \right\rfloor} \left(e^x\right)^{\left\lceil \frac{m}{2} \right\rceil}

\frac{x^n}{n!} 项系数。注意到 \left\lceil \frac{m}{2} \right\rceil = \left\lfloor \frac{m}{2} \right\rfloor + (m \mod 2),令 k = \left\lfloor \frac{m}{2} \right\rfloor,则可以化简为

\frac{1}{2^k} \sum_{s=0}^k C_{k}^{s} e^{(2s+(m\mod 2))x}

即最后答案为

\frac{1}{2^k} \sum_{s=0}^k \frac{k! (2s+(m\mod 2))^n}{s! (k-s)!}

#include <bits/stdc++.h>
using namespace std;
typedef long long lld;
const int MOD = 1e9+7, MAXN = 2e5+5;
int fac[MAXN], inv[MAXN], invs[MAXN];

int fpow(lld a, lld k)
{
    lld b = 1; k %= MOD-1;
    for (; k; k >>= 1, a = a * a % MOD)
        if (k & 1) b = b * a % MOD;
    return b;
}

void solve()
{
    lld ans = 0, n; int m, k;
    scanf("%lld %d", &n, &m); k = m / 2;
    for (int s = 0; s <= k; s++)
        ans += 1ll * fac[k] * invs[s] % MOD * invs[k-s] % MOD * fpow(2*s+(m&1), n) % MOD;
    ans = ans % MOD * fpow((MOD+1)/2, k) % MOD;
    printf("%lld\n", ans);
}

int main()
{
    fac[0] = fac[1] = inv[1] = 1;
    invs[0] = invs[1] = 1;
    for (int i = 2; i < MAXN; i++)
        fac[i] = 1ll * i * fac[i-1] % MOD,
        inv[i] = 1ll * (MOD-MOD/i) * inv[MOD%i] % MOD,
        invs[i] = 1ll * invs[i-1] * inv[i] % MOD;
    int T; scanf("%d", &T);
    while (T--) solve();
    return 0;
}

H. Luhhy’s Matrix

考虑让每个矩阵只被运算一次,那么可以使用类似于分块的想法设立分界线。

#include <bits/stdc++.h>
using namespace std;
typedef long long lld;
const int MAXN = 2e5+5;
unsigned lastans, seed, pow17[16], pow19[16];
char cost[MAXN];

struct matrix
{
    char item[16][16];
    int col[16], row[16];
    // let row[k] = item[i][k], col[k] = item[k][i]

    void norm()
    {
        for (int k = 0; k < 16; k++)
        {
            row[k] = col[k] = 0;
            for (int i = 0; i < 16; i++)
                row[k] = (row[k] << 1) | (item[i][k] & 1),
                col[k] = (col[k] << 1) | (item[k][i] & 1);
        }
    }

    matrix operator*(const matrix &b) const
    {
        matrix ans;
        for (int i = 0; i < 16; i++)
            for (int j = 0; j < 16; j++)
                ans.item[i][j] = cost[col[i]&b.row[j]];
        ans.norm();
        return ans;
    }

    static matrix getOne()
    {
        matrix ans; memset(ans.item, 0, sizeof(ans.item));
        for (int i = 0; i < 16; i++) ans.item[i][i] = 1;
        ans.norm(); return ans;
    }

    unsigned getAns() const
    {
        unsigned ans = 0;
        for (int i = 0; i < 16; i++)
            for (int j = 0; j < 16; j++)
                ans += item[i][j] * pow17[i] * pow19[j];
        return ans;
    }

    void output() const
    {
        for (int i = 0; i < 16; i++)
            for (int j = 0; j < 16; j++)
                printf("%d%c", int(item[i][j]), " \n"[j==15]);
    }

    static matrix generate()
    {
        matrix ans;
        seed ^= lastans;
        for (int i = 0; i < 16; i++)
        {
            seed ^= seed * seed + 15;
            for (int j = 0; j < 16; j++)
                ans.item[i][j] = (seed >> j) & 1;
        }
        ans.norm();
        return ans;
    }
} Q[MAXN], tail;

void solve()
{
    int front = 0, rear = 0, cut = 0, n; lastans = 0;
    scanf("%d", &n); tail = matrix::getOne();

    while (n--)
    {
        int t; scanf("%d %u", &t, &seed);

        if (t == 1)
        {
            Q[rear++] = matrix::generate();
            if (rear == front+1) tail = matrix::getOne(), cut = rear;
            else tail = Q[rear-1] * tail;
        }
        else
        {
            front++;

            if (front == cut)
            {
                cut = rear;
                for (int i = rear-2; i >= front; i--) Q[i] = Q[i+1] * Q[i];
                tail = matrix::getOne();
            }
        }

        lastans = (front == rear) ? 0 : (tail * Q[front]).getAns();
        printf("%u\n", lastans);
    }
}

int main()
{
    for (int sz2 = 1; sz2 < 65536; sz2 <<= 1)
        for (int i = 0; i < sz2; i++)
            cost[sz2+i] = cost[i] ^ 1;
    pow17[0] = pow19[0] = 1;
    for (int i = 1; i < 16; i++)
        pow17[i] = pow17[i-1] * 17,
        pow19[i] = pow19[i-1] * 19;
    int T; scanf("%d", &T);
    while (T--) solve();
    return 0;
}

K. Peekaboo

圆上整点。直接分解 a,b 然后暴力枚举点对。看题解以为卡了64倍常数,结果交了就过了(大雾)

#include <bits/stdc++.h>
using namespace std;
typedef long long lld;
typedef pair<int,int> pear;
#define F first
#define S second
inline lld sqr(int a) { return 1ll * a * a; }
lld gcd(lld a, lld b) { return b ? gcd(b, a%b) : a; }

void gaussInteger(lld R, vector<pear> &tot)
{
    auto check = [&tot] (lld r, lld rr) -> void
    {
        if (r == 1 || r % 4 != 1) return;
        for (int n = 1; n*n*2 < r; n++)
        {
            int m = int(sqrt(r-n*n)+1e-5);
            if (m * m + n * n != r) continue;
            if (gcd(m,n) != 1 || m <= n) continue;
            tot.emplace_back(rr*(m*m-n*n), 2*n*m*rr);
        }
    };

    for (int i = 1; i * i <= R; i++)
    {
        if (R % i) continue;
        check(R/i, i);
        if (i * i != R) check(i, R/i);
    }
}

void fillOver(vector<pear> &proc, int r)
{
    int cnt = proc.size();
    proc.resize(cnt*8+4);
    for (int i = 0; i < cnt; i++)
        proc[i+cnt] = pear(proc[i].S, proc[i].F);
    for (int i = 0; i < 2*cnt; i++)
        proc[i+cnt*2] = pear(proc[i].F, -proc[i].S),
        proc[i+cnt*4] = pear(-proc[i].F, proc[i].S),
        proc[i+cnt*6] = pear(-proc[i].F, -proc[i].S);
    proc[cnt*8] = pear(0, r);
    proc[cnt*8+1] = pear(r, 0);
    proc[cnt*8+2] = pear(0, -r);
    proc[cnt*8+3] = pear(-r, 0);
}

pair<pear,pear> ans[100050];

void solve()
{
    vector<pear> possA, possB;
    int a, b, c, tot = 0;
    scanf("%d %d %d", &a, &b, &c);
    gaussInteger(a, possA), gaussInteger(b, possB);
    fillOver(possA, a), fillOver(possB, b);

    for (auto A : possA) for (auto B : possB)
    {
        lld r2 = sqr(c), r3 = sqr(A.F-B.F)+sqr(A.S-B.S);
        if (r2 == r3) ans[tot++] = make_pair(A, B);
    }

    sort(ans, ans+tot);
    printf("%d\n", tot);
    for (int i = 0; i < tot; i++) printf("%d %d %d %d\n", ans[i].F.F, ans[i].F.S, ans[i].S.F, ans[i].S.S);
}

int main()
{
    int T; scanf("%d", &T);
    while (T--) solve();
    return 0;
}

HDU6414 带劲的多项式

考虑形如

f(x) = \prod_{i=1}^r (x-\lambda_i)^{l_i}

形式的多项式,它的一阶导函数为

f'(x) = \left( \sum_{i=1}^{r}l_i(x-\lambda_i) \right) \left(\prod_{i=1}^r (x-\lambda_i)^{l_i -1}\right)

则可以发现它们之间有

\gcd(f(x), f'(x)) = \prod_{i=1}^r (x-\lambda_i)^{l_i-1}

则有

\frac{f(x)}{\gcd(f(x), f'(x))} = \prod_{i=1}^r (x-\lambda_i)

由于根重数不同,每个根的因子会在不同时刻从上式消失。每次将 \gcd(f, f’) 当作新的 f 进行迭代即可。

另外,STL好慢啊233。

#include <bits/stdc++.h>
using namespace std;
typedef long long lld;
const int MOD = 998244353;
typedef vector<int> Poly;

inline lld fpow(lld a, int k)
{
    lld b = 1;
    for (; k; k>>=1, a=a*a%MOD)
        if (k&1) b=b*a%MOD;
    return b;
}

void norm(Poly &a)
{
    lld qwq = fpow(a.back(), MOD-2);
    for (auto &x : a) x = x * qwq % MOD;
}

Poly diff(const Poly &a)
{
    Poly _a = Poly(a.size() - 1);
    for (int i = 0; i < _a.size(); i++)
        _a[i] = a[i+1] * (i+1LL) % MOD;
    return _a;
}

void divide(const Poly &a, const Poly &p, Poly &q, Poly &r)
{
    assert(p.size() >= 1 && p.back() != 0);
    q.clear(), r = a; int P = p.size();
    lld pq = fpow(p.back(), MOD-2);

    for (int i = a.size()-1; i >= P-1; i--)
    {
        q.push_back(r[i] * pq % MOD);
        for (int j = 0; j < P; j++)
            r[i-P+1+j] = (r[i-P+1+j] + MOD - 1LL * r[i] * pq % MOD * p[j] % MOD) % MOD;
    }

    reverse(q.begin(), q.end());
    while (q.back() == 0) q.pop_back();
    while (r.back() == 0) r.pop_back();
}

Poly gcd(Poly a, Poly b)
{
    Poly c, d;
    while (b.size()) divide(a, b, c, d), a = b, b = d;
    norm(a);
    return a;
}

void solve()
{
    int n; scanf("%d", &n); vector<pair<int,int>> ans;
    Poly orig(n+1, 0), dif, cbs, tmp, nil, nul;
    for (int i = 0; i <= n; i++) scanf("%d", &orig[i]);
    norm(orig); dif = diff(orig), tmp = gcd(orig, dif);
    divide(orig, tmp, cbs, nil); orig = tmp, norm(orig);

    for (int l = 1; dif.size(); l++)
    {
        dif = diff(orig), tmp = gcd(orig, dif);
        divide(orig, tmp, nul, nil);

        if (nul != cbs)
        {
            // an root appears!
            Poly a, b;
            divide(cbs, nul, b, a);
            assert(b.size() == 2);
            ans.emplace_back((MOD-b[0])%MOD, l);
            cbs = nul;
        }

        orig = tmp, norm(orig);
    }

    printf("%d\n", int(ans.size()));
    for (auto p : ans) printf("%d %d\n", p.first, p.second);
}

int main()
{
    int T; scanf("%d", &T);
    while (T--) solve();
    return 0;
}

2019 ICPC南昌网络赛部分题解

A. Enju with Math Problem

给定欧拉函数表中连续100项,求是否存在或从哪个点开始。

考虑和去年EC筛 \mu^2(n) 一样的思路,检查 5, 7, 9, 11, 13, 17, 23 的倍数所在位置,利用中国剩余定理求解可能位置,然后利用区间筛筛出该区间的欧拉函数,对比。

#include <bits/stdc++.h>
using namespace std;
typedef long long lld;
const int MAXN = 2e4+5;
int pri[MAXN], pcn; bool isnp[MAXN];
int test[100], qaq[100], phi[100];

void init_prime()
{
    for (int i = 2; i < MAXN; i++)
    {
        if (!isnp[i]) pri[pcn++] = i;
        for (int j = 0; j < pcn; j++)
        {
            if (i * pri[j] >= MAXN) break;
            isnp[i * pri[j]] = true;
            if (i % pri[j] == 0) break;
        }
    }
}

int check(int y)
{
    phi[0] = qaq[0] = y;
    for (int i = 1; i < 100; i++)
        phi[i] = qaq[i] = qaq[i-1] + 1;

    int stt = qaq[0];
    if (stt <= 0 || qaq[99] > 150000000) return -1;

    for (int i = 0; i < pcn; i++)
    {
        int ft = ((stt - 1) / pri[i] + 1) * pri[i] - stt;
        for (int j = ft; j < 100; j += pri[i])
        {
            phi[j] = phi[j] / pri[i] * (pri[i] - 1);
            while (qaq[j] % pri[i] == 0) qaq[j] /= pri[i];
        }
    }

    for (int i = 0; i < 100; i++)
        if (qaq[i] > 1)
            phi[i] = phi[i] / qaq[i] * (qaq[i] - 1);
    for (int i = 0; i < 100; i++)
        if (test[i] != phi[i])
            return -1;
    return stt;
}

lld extgcd(lld a, lld b, lld &x, lld &y)
{
    if (!b) return x = 1, y = 0, a;
    lld d = extgcd(b, a % b, y, x);
    return y -= (a / b) * x, d;
}

pair<lld,lld> modeqs(lld b[], const int w[], int k)
{
    lld bj = b[0], wj = w[0], x, y, d;
    for (int i = 1; i < k; i++)
    {
        b[i] %= w[i]; d = extgcd(wj, w[i], x, y);
        if ((bj - b[i]) % d != 0) return make_pair(0LL, -1LL);
        x = x * (b[i] - bj) / d % (w[i] / d), y = wj / d * w[i];
        bj = ((bj + x * wj) % y + y) % y, wj = y;
    }
    return make_pair(bj, wj);
}

int bs[8][25]; lld b[8];
const int m[8] = { 6, 4, 6, 10, 12, 16, 18, 22 };
const int W[8] = { 9, 5, 7, 11, 13, 17, 19, 23 };
vector<int> poss[8];

int dfs(int i)
{
    if (i == 8)
    {
        auto getAns = modeqs(b, W, 8);
        if (getAns.first == 0) return -1;
        return check(getAns.first);
    }
    else
    {
        for (auto v : poss[i])
        {
            b[i] = W[i] - v;
            int t = dfs(i+1);
            if (t > 0) return t;
        }

        return -1;
    }
}

int solve()
{
    memset(bs, 0, sizeof(bs));
    for (int i = 0; i < 8; i++)
        poss[i].clear();

    for (int i = 0; i < 100; i++)
    {
        scanf("%d", &test[i]);
        for (int j = 0; j < 8; j++)
            if (test[i] % m[j] == 0)
                bs[j][i % W[j]]++;
    }

    for (int i = 0; i < 8; i++)
        for (int j = 0; j < W[i]; j++)
            if (bs[i][j] >= 100 / W[i])
                poss[i].push_back(j);
    return dfs(0);
}

int main()
{
    init_prime();
    int T; scanf("%d", &T);

    while (T--)
    {
        int ans = solve();
        if (ans == -1) printf("NO\n");
        else printf("YES\n%d\n", ans);
    }

    return 0;
}

D. Interesting Series

考虑 F_n = 1 + a + a^2 + \cdots + a^{n-1} = \frac{a^n – 1}{a-1}

那么所求相当于是

Answer(K) = \frac{\sum_{s \in K} a^s – C_{N}^{|K|}}{a-1}

考虑构造多项式

f(n) = \prod_{i=1}^n (1 + a^s x)

那么所求 \sum_{s \in K} a^s 就是 f(n)x^{|K|} 项系数之和。

分治FFT计算即可。注意到模数比较小,所以考虑用double或者long double跑FFT,但是要注意预处理精度(感谢Claris姐姐呜呜呜),不能多次乘。

#include <bits/stdc++.h>
typedef long long lld;
const double PI = acos(-1.0);
const int MOD = 100003;
typedef std::complex<double> cplx;
int fac[MOD], inv[MOD], invs[MOD];

void dft(cplx a[], int DFT, int N) {
    static int rev[1<<20]; static cplx w[1<<20], ww[1<<20];
    for (int i = 0; i < N; i++) {
        rev[i] = (rev[i>>1]>>1)|((i&1)?N>>1:0);
        if (i < rev[i]) swap(a[i], a[rev[i]]);
        w[i] = cplx(cos(PI*i/N), sin(PI*i/N));
        ww[i] = cplx(w[i].real(), -w[i].imag());
    }
    for (int d = 0, ctz = __builtin_ctz(N); (1<<d) < N; d++) {
        cplx u, t; int m2 = 1<<d, m = m2<<1;
        for (int k = 0; k < N; k += m)
            for (int j = 0; j < m2; j++)
                t = (DFT==1?w:ww)[j<<(ctz-d)] * a[k+j+m2], u = a[k+j],
                a[k+j] = u + t, a[k+j+m2] = u - t;
    }
    if (DFT == -1) for (int i = 0; i < N; i++) a[i] /= N;
}

lld fpow(lld a, lld k) {
    lld b = 1;
    for (; k; k>>=1, a=a*a%MOD)
        if (k&1) b=b*a%MOD;
    return b;
}

int t[131072], ss[MOD], a;

void solve(int l, int r, int L, int R)
{
    if (R - L <= 31)
    {
        int len = R-L+1, sz = 1;
        for (int i = 0; i < len; i++) t[L+i] = !i;
        for (int w = l; w <= r; w++, sz++)
            for (int i = sz; i > 0; i--)
                (t[L+i] += t[L+i-1] * fpow(a, ss[w]) % MOD) %= MOD;
    }
    else
    {
        int mid = (L+R)>>1, m = (l+r)>>1, len=R-L+1;
        solve(l, m, L, mid), solve(m+1, r, mid+1, R);
        static cplx a[131072], b[131072];
        for (int i = 0; i < len/2; i++)
            a[i] = t[L+i], b[i] = t[mid+1+i],
            a[len/2+i] = b[len/2+i] = 0;
        dft(a, 1, len), dft(b, 1, len);
        for (int i = 0; i < len; i++) a[i] *= b[i];
        dft(a, -1, len);
        for (int i = 0; i < len; i++)
            t[L+i] = lld(a[i].real()+0.5) % MOD;
    }
}

int main()
{
    inv[0] = fac[0] = invs[0] = invs[1] = inv[1] = fac[1] = 1;
    for (int i = 2; i < MOD; i++)
        fac[i] = 1ll * fac[i-1] * i % MOD,
        inv[i] = 1ll * (MOD-MOD/i) * inv[MOD%i] % MOD,
        invs[i] = 1ll * invs[i-1] * inv[i] % MOD;

    int N, Q, s;
    scanf("%d %d %d", &N, &a, &Q);
    for (int i = 1; i <= N; i++)
        scanf("%d", &ss[i]);
    int sz = 30, st = 32;
    while (sz < N) sz <<= 1, st <<= 1;
    solve(1, N, 0, st-1);

    while (Q--)
    {
        scanf("%d", &s);
        lld ans = t[s] - 1ll * fac[N] * invs[s] * invs[N-s] % MOD;
        ans = (ans + MOD) * inv[a-1] % MOD;
        printf("%lld\n", ans);
    }

    return 0;
}

话说回来补题的时候一直没过竟然是因为……输出答案时那个组合数写错了呜呜呜,变量意义改变不改名,调试两行泪。

H. The Nth Item

考虑作出 f(n+x) = a \times f(n) + b \times f(n+1) 的表达式,然后分段打表。

#include <cstdio>
using namespace std;
typedef long long lld;
const int MAXN = 5e5+5;
const int MOD = 998244353;
int trans[1002][2], f[MAXN][2];

int main()
{
    trans[0][0] = trans[1][1] = f[0][1] = 1;

    for (int i = 2; i <= 1001; i++)
    {
        trans[i][0] = 2ll * trans[i-1][1] % MOD;
        trans[i][1] = (trans[i-1][0] + 3ll * trans[i-1][1]) % MOD;
    }

    for (int i = 1; i < MAXN; i++)
    {
        f[i][0] = (1ll * trans[1000][0] * f[i-1][0] + 1ll * trans[1000][1] * f[i-1][1]) % MOD;
        f[i][1] = (1ll * trans[1001][0] * f[i-1][0] + 1ll * trans[1001][1] * f[i-1][1]) % MOD;
    }

    int Q; lld N; scanf("%d %lld", &Q, &N);
    lld lastAns = N, tot = 0; N = 0;

    for (int i = 0; i < Q; i++)
    {
        lld now = lastAns ^ N; N = now; now %= (MOD-1)/2;
        lld ans = (1ll * trans[now%1000][0] * f[now/1000][0] + 1ll * trans[now%1000][1] * f[now/1000][1]) % MOD;
        tot ^= ans; lastAns = ans * ans;
    }

    printf("%lld\n", tot);
    return 0;
}

2019 ICPC徐州网络赛部分题解

E. XKC’s basketball team

就在区间 [i+1,R] 找最右边大于 a_i+m 的数的位置。然后考虑到线段树去查询,先查右子树,再查左子树,然后写着写着就成了二叉树了……

#include <bits/stdc++.h>
using namespace std;
const int MAXN = 5e5+5;
int mmax[MAXN*4], a[MAXN];
#define lson i<<1
#define rson i<<1|1
int n, m;

int query(int x, int y, int L, int R, int i)
{
    if (mmax[i] < x || R <= y) return -1; else if (L == R) return L;
    int mid = (L+R)>>1; int q = query(x, y, mid+1, R, rson);
    if (~q) return q; else return query(x, y, L, mid, lson);
}

void build(int L, int R, int i)
{
    if (L == R) { scanf("%d", &mmax[i]); a[L] = mmax[i]; return; }
    int mid = (L+R)>>1; build(L,mid,lson), build(mid+1,R,rson);
    mmax[i] = max(mmax[lson], mmax[rson]);
}

int main()
{
    scanf("%d %d", &n, &m);
    build(1, n, 1);
    for (int i = 1; i <= n; i++)
    {
        int loc = query(a[i]+m, i, 1, n, 1);
        if (loc > 0) loc = loc - i - 1;
        printf("%d%c", loc, " \n"[i==n]);
    }
    return 0;
}

H. function

想不到,也只能抄抄题解过日子了。

#include <bits/stdc++.h>
using namespace std;
typedef long long lld;
const int MAXN = 2e5+5;
const int MOD = 998244353;
const int INV2 = (MOD+1)/2;
int pri[MAXN], pcn, sp[MAXN];
int m, Sqrt, id1[MAXN], id2[MAXN];
bool isnp[MAXN]; lld g[MAXN], h[MAXN], w[MAXN], n;

int main()
{
    scanf("%lld", &n);

    for (int i = 2; i < MAXN; i++)
    {
        if (!isnp[i])
        {
            pcn++; pri[pcn] = i;
            sp[pcn] = (sp[pcn-1] + i) % MOD;
        }

        for (int j = 1; j <= pcn; j++)
        {
            if (i * pri[j] >= MAXN) break;
            isnp[i * pri[j]] = true;
            if (i % pri[j] == 0) break;
        }
    }

    m = 0, Sqrt = sqrt(n);

    for (lld L = 1, R; L <= n; L = R+1)
    {
        R = n / (n / L), w[++m] = n / L;
        (w[m] <= Sqrt ? id1[w[m]] : id2[R]) = m;
        g[m] = (w[m] - 1) % MOD;
        h[m] = (w[m] % MOD + 2) * (w[m] % MOD - 1) % MOD * INV2 % MOD;
    }

    for (int j = 1; j <= pcn; j++)
    {
        for (int i = 1, d; i <= m && pri[j] <= w[i]/pri[j]; i++)
        {
            d = (w[i]/pri[j]<=Sqrt) ? id1[w[i]/pri[j]] : id2[n/(w[i]/pri[j])];
            g[i] -= g[d] - j + 1;
            h[i] -= pri[j] * (h[d] - sp[j-1]) % MOD;
        }
    }

    lld ans = 0;

    for (int i = 1; i <= pcn; i++)
    {
        if (n / pri[i] >= pri[i])
        {
            for (lld p = pri[i]; p <= n; p *= pri[i])
            {
                lld fk = n/p%MOD;
                ans += (n+1)%MOD*fk%MOD;
                ans -= fk*(fk+1)%MOD*INV2%MOD*p%MOD;
            }
        }
    }

    for (int i = 1; i < Sqrt; i++)
    {
        lld sp2 = h[i] - h[i+1], sp1 = g[i] - g[i+1];
        ans += (n+1)%MOD*i%MOD*sp1%MOD;
        ans -= i*(i+1ll)/2%MOD*sp2%MOD;
    }

    ans %= MOD; if (ans < 0) ans += MOD;
    printf("%lld\n", ans);
    return 0;
}

J. Random Access Iterator

假设深度为树高的叶节点集合为 R,则相当于从起点按某种移动方式转移到 R 中的概率。考虑 dp[u] 为从 u 转移不到 R 中的概率,那么有

dp[u] = \frac{1}{d(u)} \sum_{(u,v)\in T} dp[v]

同时注意深度为树高的叶节的概率为 0,其他叶节点概率为 1,即DP初值。然后两遍DFS解决。

#include <vector>
#include <cstdio>
using namespace std;
typedef long long lld;
const int MAXN = 1e6+5;
const int MOD = 1e9+7;
vector<int> G[MAXN];
int maxdep;

lld fpow(lld a, lld k)
{
    lld b = 1;
    for (; k; k>>=1, a=a*a%MOD)
        if (k&1) b=b*a%MOD;
    return b;
}

void dfs1(int u, int p, int d)
{
    maxdep = max(maxdep, d);
    for (auto v : G[u]) if (v != p) dfs1(v, u, d+1);
}

lld dfs2(int u, int p, int d)
{
    int k = 0; lld dp = 0;
    for (auto v : G[u]) if (v != p) dp += dfs2(v, u, d+1), k++;
    if (k) return fpow(dp * fpow(k, MOD-2) % MOD, k);
    else return d != maxdep ? 1 : 0;
}

int main()
{
    int n; scanf("%d", &n);

    for (int i = 1, u, v; i < n; i++)
    {
        scanf("%d %d", &u, &v);
        G[u].push_back(v);
        G[v].push_back(u);
    }

    dfs1(1,0,1);
    printf("%lld\n", (1+MOD-dfs2(1,0,1))%MOD);
    return 0;
}

K. Center

#include <unordered_map>
#include <unordered_set>
using namespace std;
typedef long long lld;
typedef pair<int,int> pii;

const lld INF = 1e8;
const lld STD = 1e7;
#define X first
#define Y second

lld trans(int x, int y) { return x * INF + y; }
lld trans(pii n) { return n.X * INF + n.Y; }
pii inv(lld t) { return pii(t/INF, t%INF); }
pii operator*(pii a, int n) { return pii(a.X * n, a.Y * n); }
pii operator+(pii a, pii b) { return pii(a.X + b.X, a.Y + b.Y); }
pii operator-(pii a, pii b) { return pii(a.X - b.X, a.Y - b.Y); }

pii p[1010];
unordered_set<lld> existed;
unordered_map<lld, int> newCenter;

int main()
{
    int n, x, y, tot;
    scanf("%d", &n), tot = n;

    for (int i = 0; i < n; i++)
    {
        scanf("%d %d", &x, &y);
        x += STD, y += STD;
        p[i] = make_pair(x, y);
        existed.insert(trans(x, y));
    }

    for (int i = 0; i < n; i++)
    {
        int ans = 0;
        for (int j = 0; j < n; j++)
            if (!existed.count(trans(p[i] * 2 - p[j])))
                ans++;
        tot = min(tot, ans);
    }

    for (int i = 0; i < n; i++)
        for (int j = i+1; j < n; j++)
            newCenter[trans(p[i] + p[j])]++;

    int ptp = 0; pii nc;
    for (auto cnt : newCenter)
        if (cnt.second > ptp)
            nc = inv(cnt.first), ptp = cnt.second;
    int maybe = 0;
    for (int j = 0; j < n; j++)
        if (!existed.count(trans(nc - p[j])))
            maybe++;
    tot = min(tot, maybe);

    printf("%d\n", tot);
    return 0;
}

2019 ICPC南京网络赛部分题解

B. super_log

就考察一下你是不是真的理解了拓展欧拉定理。

#include <bits/stdc++.h>
using namespace std;
typedef long long lld;
const int MAXN = 1e6+5;
int pri[MAXN], pcn, phi[MAXN];
bool isnp[MAXN];

lld fpow(lld a, lld k, lld m)
{
    lld b = 1;
    for (; k; k>>=1, a=a*a%m)
        if (k&1) b=b*a%m;
    return b;
}

bool judge(lld a, lld k, lld m)
{
    lld b = 1;
    while (k)
    {
        b = b * a;
        if (b >= m) return true;
        k--;
    }
    return false;
}

lld getAns(int a, int b, int m)
{
    if (m == 1) return a < m ? 0 : m;
    if (b == 0) return 1;
    lld t = getAns(a, b-1, phi[m]);
    lld ans = fpow(a, t, m);
    bool lg = t < phi[m] ? judge(a, t, m) : true;
    if (lg) ans += m;
    //printf("f(%d,%d,%d)=%lld\n", a, b, m, ans);
    return ans;
}

void solve()
{
    int a, b, m;
    scanf("%d %d %d", &a, &b, &m);
    printf("%lld\n", getAns(a,b,m)%m);
}

int main()
{
    phi[1] = 1;
    for (int i = 2; i < MAXN; i++)
    {
        if (!isnp[i]) pri[pcn++] = i, phi[i] = i-1;
        for (int j = 0; j < pcn; j++)
        {
            int k = i * pri[j];
            if (k >= MAXN) break;
            isnp[k] = true;
            if (i % pri[j] == 0) { phi[k] = phi[i] * pri[j]; break; }
            else phi[k] = phi[i] * (pri[j] - 1);
        }
    }

    int T; scanf("%d", &T);
    while (T--) solve();
    return 0;
}

C. Tsy’s number 5

抄题解。令 f_k = \sum_{i=1}^n [\varphi(i)=k],则有

\begin{aligned} ans & = \sum_{i=1}^n \sum_{j=1}^n \varphi(i) \varphi(j) 2^{\varphi(i) \varphi(j)} \\ & = \sum_{i=1}^n \sum_{j=1}^n f_i f_j ij 2^{ij} \\ & = 2 \sum_{i=1}^n \sum_{j=1}^i f_i f_j ij 2^{ij} – \sum_{i=1}^n {f_i}^2 i^2 2^{j^2} \\ & = 2 \sum_{i=1}^n f_ii \sum_{j=1}^i f_jj {\sqrt2}^{2ij} – \sum_{i=1}^n {f_i}^2 i^2 2^{j^2} \\ & = 2 \sum_{i=1}^n f_ii {\sqrt2}^{i^2} \sum_{j=1}^i f_jj {\sqrt2}^{j^2} {\sqrt2}^{-(i-j)^2} – \sum_{i=1}^n {f_i}^2 i^2 2^{j^2} \end{aligned}

#include <bits/stdc++.h>
using namespace std;
typedef long long lld;
const int MAXN = 1e5+5;
const int MOD = 998244353, G = 3;
const int SQRT2 = 116195171;
int phi[MAXN], sp[MAXN], qwq[MAXN], A[262144], B[262144];

inline lld fpow(lld a, lld k)
{
    lld b = 1;
    for (; k; k>>=1, a=a*a%MOD)
        if (k & 1) b=b*a%MOD;
    return b;
}

inline void init_prime()
{
    static int pri[MAXN], pcn;
    static bool isnp[MAXN];

    phi[1] = sp[0] = 1; sp[1] = fpow(SQRT2, MOD-2);
    for (int i = 2; i < MAXN; i++)
    {
        if (!isnp[i]) pri[pcn++] = i, phi[i] = i-1;
        sp[i] = fpow(SQRT2, MOD-1-1ll*i*i%(MOD-1));
        for (int j = 0; j < pcn; j++)
        {
            if (i * pri[j] >= MAXN) break;
            isnp[i * pri[j]] = true;
            if (i % pri[j] == 0) { phi[i * pri[j]] = phi[i] * pri[j]; break; }
            else phi[i * pri[j]] = phi[i] * (pri[j] - 1);
        }
    }
}

inline void dft(int a[], int DFT, int N)
{
    static int rev[262144], wmk[262144];
    for (int i = 0; i < N; i++)
        rev[i] = (rev[i>>1]>>1)|((i&1)?(N>>1):0);
    for (int i = 0; i < N; i++)
        if (i < rev[i]) swap(a[i], a[rev[i]]);

    for (int m = 2, m2 = 1; m <= N; m <<= 1, m2 <<= 1) {
        wmk[0] = 1, wmk[1] = fpow(G, (MOD-1)/m*DFT+(DFT==-1?MOD-1:0));
        for (int j = 2; j < m2; j++) wmk[j] = 1ll * wmk[j-1] * wmk[1] % MOD;
        for (int k = 0; k < N; k += m)
            for (int j = 0, u, t; j < m2; j++)
                t = 1ll * wmk[j] * a[k+j+m2] % MOD, u = a[k+j],
                a[k+j] = (u+t)%MOD, a[k+j+m2] = (u-t+MOD)%MOD;
    }

    if (DFT == -1)
        for (int i = 0, invN = fpow(N, MOD-2); i < N; i++)
            a[i] = 1ll * a[i] * invN % MOD;
}

inline void solve()
{
    int n; scanf("%d", &n);
    int len = 1; while (len < n+1) len <<= 1; len <<= 1;
    for (int i = 0; i < len; i++) A[i] = B[i] = 0;
    for (int i = 0; i <= n; i++) B[i] = sp[i], A[phi[i]]++;
    for (int i = 0; i <= n; i++)
        qwq[i] = A[i] = 1ll * i * A[i] % MOD * fpow(SQRT2, 1ll*i*i%(MOD-1)) % MOD;
    dft(A, 1, len); dft(B, 1, len);
    for (int i = 0; i < len; i++) A[i] = 1ll * A[i] * B[i] % MOD;
    dft(A, -1, len); lld ans = 0;
    for (int i = 1; i <= n; i++)
        ans += (2ll * A[i] - qwq[i]) * qwq[i] % MOD;
    ans %= MOD; if (ans < 0) ans += MOD;
    printf("%lld\n", ans);
}

int main()
{
    init_prime();
    int T; scanf("%d", &T);
    while (T--) solve();
    return 0;
}

D. Robots

考虑从 u 出发,到达 n 的期望步数,逆向建图拓扑排序搞定。

dp[u] = \frac{1}{d(u) + 1} \left( dp[u] + \sum_{(u,v)\in G} dp[v] \right) + 1

再考虑从 u 出发,到达 n 的期望消耗天数累计值,则有

dp2[u] = \frac{1}{d(u) + 1} \left( dp2[u] + \sum_{(u,v)\in G} dp2[v] \right) + dp[u]

其实是乱搞猜结论出来的,但是后面那个 dp[u] 的意义其实可以理解为,我这一步到终点的天数,是由路径上的每条边叠加得到的。

这个期望的期望好像题解没给证明,先这样吧QAQ。

#include <bits/stdc++.h>
using namespace std;
typedef long long lld;
const int MAXN = 1e5+5;
vector<int> G[MAXN];
int d[MAXN], d2[MAXN], Q[MAXN];
double one[MAXN], dp[MAXN], dp2[MAXN];

void run(double dp[], double dist[], int n)
{
    int front = 0, rear = 0; Q[rear++] = n;
    memcpy(d, d2, sizeof(d));

    while (front < rear)
    {
        int u = Q[front++];
        if (u != n) dp[u] = (d2[u]+1.0)/d2[u]*dist[u] + dp[u]/d2[u];

        for (auto v : G[u])
        {
            d[v]--; if (d[v] == 0) Q[rear++] = v;
            dp[v] += dp[u];
        }
    }
}

void solve()
{
    int n, m; scanf("%d %d", &n, &m);
    for (int i = 1; i <= n; i++)
        dp2[i] = dp[i] = d2[i] = 0, G[i].clear(), one[i]=1;
    for (int i = 0, u, v; i < m; i++)
        scanf("%d %d", &u, &v), G[v].push_back(u), d2[u]++;
    run(dp, one, n);
    run(dp2, dp, n);
    printf("%.2lf\n", dp2[1]);
}

int main()
{
    int T; scanf("%d", &T);
    while (T--) solve();
    return 0;
}

E. K Sum

套路变换。

\begin{aligned} f_n(k) & = \sum_{d=1}^n d^2 \left({\sum_{l_i=1}^n}\right)^k [\gcd(l_i) = 1] \\ & = \sum_{d=1}^n d^2 \sum_{i=1}^{\lfloor n/d \rfloor} \mu(i) \left\lfloor \frac{n}{id} \right\rfloor^k \\ & = \sum_{s=1}^n \left(\sum_{id=s}d^2 \mu(i)\right)\left\lfloor\frac{n}{s}\right\rfloor^k \end{aligned}

合并,后面那个等比数列求和,注意特判 1,所以需要得到 k \mod 1e9+7k \mod 1e9+6

\sum_{i=2}^k f_n(i) = \sum_{s=1}^n g(s) \sum_{i=2}^k\left\lfloor\frac{n}{s}\right\rfloor^i

杜教筛部分看起来像是 id^2 * \mu,那么配上 1

S(n) = \sum_{i=1}^n g(i) = \frac{n(n+1)(2n+1)}{6} – \sum_{i=2}^n S\left(\left\lfloor\frac{n}{i}\right\rfloor\right)

#include <bits/stdc++.h>
using namespace std;
typedef long long lld;
const int MOD = 1e9+7;
const int MAXN = 3e6+5;
int _G1[MAXN];

void init_prime()
{
    static int pri[MAXN], pcn;
    static bool isnp[MAXN];

    _G1[1] = 1;
    for (int i = 2; i < MAXN; i++)
    {
        if (!isnp[i]) pri[pcn++] = i, _G1[i] = (1ll * i * i-1) % MOD;
        for (int j = 0; j < pcn; j++)
        {
            if (i * pri[j] >= MAXN) break;
            isnp[i * pri[j]] = true;
            if (i % pri[j] == 0)
            {
                _G1[i * pri[j]] = _G1[i] * 1ll * pri[j] % MOD * pri[j] % MOD;
                break;
            }
            else
            {
                _G1[i * pri[j]] = _G1[i] * (1ll * pri[j] % MOD * pri[j] - 1)% MOD;
            }
        }
    }

    for (int i = 1; i < MAXN; i++)
        _G1[i] = (_G1[i-1] + _G1[i]) % MOD;
}

lld sumG(int n)
{
    static map<int,lld> _G2;
    if (n < MAXN) return _G1[n];
    if (_G2.count(n)) return _G2[n];

    lld ans = n * (n+1ll) % MOD * (2ll*n+1) % MOD * ((MOD+1)/6) % MOD;
    for (int L = 2, R; L <= n; L = R+1)
    {
        R = n / (n / L);
        ans -= (R-L+1) * sumG(n/L) % MOD;
    }

    ans %= MOD; if (ans < 0) ans += MOD;
    return _G2[n] = ans;
}

inline lld fpow(lld a, lld k)
{
    lld b = 1;
    for (; k; k>>=1, a=a*a%MOD)
        if (k&1) b=a*b%MOD;
    return b;
}

inline lld powers(int ns, int k, int k2)
{
    if (ns == 1) return k2-1;
    lld ans = fpow(ns, k+1) - 1ll * ns * ns % MOD;
    ans = ans * fpow(ns-1, MOD-2) % MOD;
    return ans;
}

inline void solve()
{
    int n, k = 0, k2 = 0; lld ans = 0;
    static char kk[MAXN];
    scanf("%d %s", &n, kk);
    for (int i = 0; kk[i]; i++)
        k = (k * 10ll + kk[i] - '0') % (MOD-1),
        k2 = (k2 * 10ll + kk[i] - '0') % MOD;

    for (int L = 1, R; L <= n; L = R+1)
    {
        R = n / (n / L);
        ans += (sumG(R) - sumG(L-1)) * powers(n/L, k, k2) % MOD;
    }

    ans %= MOD; if (ans < 0) ans += MOD;
    printf("%lld\n", ans);
}

int main()
{
    init_prime();
    int T; scanf("%d", &T);
    while (T--) solve();
    return 0;
}

G. Quadrilateral

参考题解。我怎么写了一坨屎QAQ。

#include <bits/stdc++.h>
using namespace std;
typedef long long lld;
inline lld sum1(lld c) { return c*(c+1)/2; }
inline lld sum2(lld c) { return c*(c+1)*(2*c+1)/6; }

// count the pairs for [1,a] + [1,b] <= [1,c]
inline lld func1(lld a, lld b, lld c)
{
    if (a > b) swap(a, b); c--;
    if (a < 1 || b < 1) return 0;
    if (c <= 0) return 0;
    else if (c < a) return sum1(c)*(c+1) - sum2(c);
    else if (c <= b) return (c+1-a)*(c+1-a)*a - sum1(c-a)*a + sum1(a-1)*(c+1) - sum2(a-1);
    else if (c < a+b) return (c-b)*(a+b)*(c+1) - (a+b+c+1)*(sum1(c)-sum1(b)) + (sum2(c)-sum2(b)) + (b+1-a)*(c+1-a)*a - sum1(b-a)*a + sum1(a-1)*(c+1) - sum2(a-1);
    else return a*(a+b)*(c+1) - (a+b+c+1)*(sum1(a+b)-sum1(b)) + (sum2(a+b)-sum2(b)) + (b+1-a)*(c+1-a)*a - sum1(b-a)*a + sum1(a-1)*(c+1) - sum2(a-1);
}

// calc the valid pairs of ([1,a], [1,b], [1,c], d)
inline lld func2(lld a, lld b, lld c, lld d)
{
    return a * b * c
         - func1(b, c, a-d)
         - func1(a, c, b-d)
         - func1(a, b, c-d)
         - func1(a, b, d-1)
         + func1(a, b, d-c-1);
}

// calc the valid pairs of ([la,ra], [lb,rb], [lc,rc], d)
inline lld func3(pair<int,int> a, pair<int,int> b, pair<int,int> c, int d)
{
    return func2(a.second, b.second, c.second, d)
         - func2(a.second, b.first-1, c.second, d)
         - func2(a.first-1, b.second, c.second, d)
         - func2(a.second, b.second, c.first-1, d)
         + func2(a.first-1, b.first-1, c.second, d)
         + func2(a.first-1, b.second, c.first-1, d)
         + func2(a.second, b.first-1, c.first-1, d)
         - func2(a.first-1, b.first-1, c.first-1, d);
}

inline void solve()
{
    pair<int,int> iv[4]; lld ans = 0;
    for (int i = 0; i < 4; i++)
        scanf("%d %d", &iv[i].first, &iv[i].second);
    for (int i = iv[3].first; i <= iv[3].second; i++)
        ans += func3(iv[0], iv[1], iv[2], i);
    printf("%lld\n", ans);
}

int main()
{
    int T; scanf("%d", &T);
    while (T--) solve();
    return 0;
}

2018CCPC吉林区域赛 G. High Priestess

题目链接:HDU-6561

给不超过 10^4 个阻值为 1 \Omega 的电阻,通过串联和并联使得最后形成的电路有效阻值为 r

这道题可以考虑有理逼近中的连分数。

a_0 + \dfrac{1}{a_1 + \dfrac{1}{a_2 + \dfrac{1}{\ddots + \dfrac{1}{a_n}}}}

我们通常将它记为 [a_0, a_1, a_2, \cdots, a_n]

\frac{p_k}{q_k}[a_0, a_1, \cdots, a_k],我们称 \frac{p_k}{q_k}[a_0, a_1, a_2, \cdots, a_n] 的第 k 个渐进分数。

可以证明,所有有理数都可以表示为这样的一个有限连分数。

考虑斐波那契的 \phi = \frac{1 + \sqrt5}{2},我们可以令连分数序列为无限的 [1, 1, 1, \cdots, 1, \cdots]。而且我们可以发现,\frac{p_k}{q_k} 正好等于 \frac{F_{k+1}}{F_k}。数论中还会研究,任何无理数也都可以用连分数来逼近。

可以看出来,将这样的一个有理数摆进来,首先将它的整数部分去掉,然后再求倒数,就变成了 [a_1, a_2, \cdots, a_n],问题规模得以缩小。当我们迭代到没有小数部分的时候,这个有理数已经被精确表示了。类似于一个辗转相除法的例子,对不对?时间复杂度为 O(\log q)

我们需要(cai)知(chu)道(lai)的是,渐进级数越精确,误差越小。如果我们用某个连分数来表示这道题的答案呢?

一定有 a_0 = 0。那么我们不妨设定一个数列 \lbrace b_i \rbrace,其中

b_i = \begin{cases} \frac{1}{a_i + \frac{1}{b_{i+1}}}, & i = 2k+1 \\ a_i + b_{i+1}, & i=2k \end{cases}

那么可以认为 b_i 是由某一些电阻构成的复合电阻。当 i 为奇数时,我们将 a_i1 \Omega 的电阻与 b_{i+1} 并联起来;当 i 为偶数时,我们将 a_i1 \Omega 的电阻与 b_{i+1} 串联起来。注意当超过连分数长度 n 时,b_{i+1}0\infty,取决于奇偶性。

首先题目给定的 r 输入进来以后一定是一个有理数。那我们不妨先求出 [ r – \epsilon, r + \epsilon ] 中分母最小的分数 \frac{p_i}{q_i}

我们将 \frac{p_i}{q_i} 的连分数表示出来。但是由于题目要求所用电阻数量不能超过 10^4,我们可以将超过的部分略去,留下一堆多余的电阻,这样比略去的逼近效果略好。

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll EPS = 10000000;
const int MAXN = 10000;
int kase;

ll gcd(ll a, ll b) { return b ? gcd(b, a % b) : a; }

void euclid(ll a, ll b, ll x, ll y, ll &p, ll &q)
{
    ll K = a / b; a -= K * b, x -= K * y;
    if (x > y) p = q = 1; else euclid(y, x, b, a, q, p);
    p += K * q;
}

void solve()
{
    static char buf[1000];
    static ll lfs[1000];
    scanf("%s", buf);
    ll fz = 0, fm = 1, fzm, fmm, fzM, fmM;

    for (int i = 2; buf[i]; i++)
    {
        fz = fz * 10 + buf[i] - '0';
        fm *= 10;
    }

    if (fm < EPS) fz *= EPS / fm, fm = EPS;
    fzm = fz - fm / EPS, fmm = gcd(fzm, fm), fzm /= fmm, fmm = fm / fmm;
    fzM = fz + fm / EPS, fmM = gcd(fzM, fm), fzM /= fmM, fmM = fm / fmM;
    euclid(fzm, fmm, fzM, fmM, fz, fm);
    if (fmM < fm) fz = fzM, fm = fmM;
    if (fmm < fm) fz = fzm, fm = fmm;
    int tot = 1, sum = 0, idx = 0, idg, id1, id2, idq = -1, op;

    while (fz && sum < MAXN)
    {
        swap(fz, fm);
        sum += (lfs[tot++] = min(fz / fm, MAXN - sum + 0ll));
        fz %= fm;
    }

    // printf("[ 0"); for (int i = 1; i < tot; i++) printf(", %lld", lfs[i]); printf(" ]\n");
    printf("Case %d:\n%d %d\n", ++kase, sum, sum-1);

    idg = sum-1;
    for (int j = tot-1; j > 0; j--)
    {
        op = j & 1;
        id1 = idx++;

        for (int i = 1; i < lfs[j]; i++)
        {
            id2 = idx++;
            printf("%d %d %d\n", op, id1, id2);
            id1 = ++idg;
        }

        if (~idq)
        {
            id2 = idq;
            printf("%d %d %d\n", op, id1, id2);
            ++idg;
        }

        idq = idg;
    }
}

int main()
{
    int T; scanf("%d", &T);
    while (T--) solve();
    return 0;
}

这道题总而言之感觉还是有点怪怪的,例如 0.000001 一定是无解的,但是竟然也能 AC。就这样吧,应该是良心出题人没有出卡人数据。

2019杭电多校第五场 I. discrete logarithm problem

题目链接:HDU 6632

已知数字 a,b,p,并且 p=2^a3^b+1。想求解离散对数 a^x \equiv b \mod p

需要一个东西叫做 Pohlig-Hellman 算法。以下是来自维基百科内容的翻译,将使用群论语言。

p^e 阶循环群

这个算法用于计算 p^e 阶循环群 G=\langle g \rangle 的离散对数 g^x = h

  1. 初始化 x_0 := 0

  2. 计算 \gamma := g^{p^{e-1}},由拉格朗日配集分解计数定理知,ord~\gamma = p

  3. 对于 k \in \lbrace 0,\cdots,e-1 \rbrace,分别做

    1. 计算 h_k := (g^{-x_k}h)^{p^{e-1-k}}。由此构造可知,此元素的阶一定整除 p,即阶为 1p,因此 h_k \in \langle \gamma \rangle

    2. 使用 BSGS 算法,计算 d_k \in \lbrace 0, \cdots, p-1 \rbrace,其中 \gamma^{d_k} = h_k。这个操作的复杂度为 O(\sqrt{p})

    3. x_{k+1} := x_k + p^k d_k

  4. 返回 x_e

e \ll p 时,这个算法优于 BSGS,复杂度为 O(e\sqrt{p})

一般合数阶循环群

首先我们将 G = \langle g \rangle 的阶数刻画为 n = \prod_{i=1}^r {p_i}^{e_i}。现在我们解方程 g^x = h

  1. 对循环群阶数做因数分解

  2. 对于 i \in \lbrace 1, \cdots, r \rbrace,分别做

    1. 计算 g_i := g^{n / {p_i}^{e_i}},由拉格朗日定理知此元素的阶为 {p_i}^{e_i}

    2. 计算 h_i := h^{n/{p_i}^{e_i}},由此构造知 h_i \in \langle g_i \rangle

    3. 利用上面算法计算 x_i 使得 {g_i}^{x_i} = h_i

  3. 利用中国剩余定理解方程组 x \equiv x_i \mod {p_i}^{e_i}

  4. 返回 x

因此总的时间复杂度为

T(n) = O\left( \sum_{i=1}^r e_i (\log n + \sqrt{p_i}) \right)

回到本题

显然因为 p_i = 2, 3,我们不用 BSGS……好的。

#include <bits/stdc++.h>
using namespace std;
typedef long long lld;

lld fpow(lld _a, lld k, lld MOD) {
    __int128 b = 1, a = _a % MOD;
    while (k) {
        if (k & 1) b = b * a % MOD;
        a = a * a % MOD;
        k >>= 1;
    }
    return b;
}

__int128 extgcd(__int128 a, __int128 b, __int128 &x, __int128 &y) {
    if (!b) return x = 1, y = 0, a;
    __int128 d = extgcd(b, a % b, y, x);
    return y -= (a / b) * x, d;
}

int T, e2, e3, G;
lld MOD, n, p2, p3;

lld solvePE(lld g, lld h, lld p, lld pe, int e)
{
    lld x = 0, pp = 1;
    lld y = fpow(g, pe/p, MOD);

    for (int k = 0; k < e; k++)
    {
        __int128 hk = fpow(g, pe-x, MOD);
        hk = hk * h % MOD;
        hk = fpow(hk, pe/p/pp, MOD);
        __int128 d, yy;
        for (d = 0, yy = 1; d < p; d++, yy = yy * y % MOD)
            if (yy == hk) break;
        x = (x + pp * d) % MOD, pp *= p;
    }

    return x;
}

__int128 solveN(lld h)
{
    if (e3 == 0) return solvePE(G, h, 2, p2, e2);
    else if (e2 == 0) return solvePE(G, h, 3, p3, e3);
    __int128 x2 = solvePE(fpow(G, p3, MOD), fpow(h, p3, MOD), 2, p2, e2);
    __int128 x3 = solvePE(fpow(G, p2, MOD), fpow(h, p2, MOD), 3, p3, e3);
    __int128 invp2, invp3; assert(extgcd(p2, p3, invp2, invp3) == 1);
    return (x2 * p3 * invp3 + x3 * p2 * invp2) % (MOD-1);
}

int main()
{
    cin >> T;

    while (T--)
    {
        long long a, b;
        cin >> MOD >> a >> b;
        n = MOD-1, e2 = e3 = 0, p2 = p3 = 1;
        while (n % 2 == 0) n /= 2, p2 *= 2, e2++;
        while (n % 3 == 0) n /= 3, p3 *= 3, e3++;
        for (G = 2; (e2 > 0 && fpow(G, (MOD-1)/2, MOD) == 1)
                 || (e3 > 0 && fpow(G, (MOD-1)/3, MOD) == 1); G++);
        __int128 inda = solveN(a), indb = solveN(b), x, y;
        __int128 d = extgcd(inda, MOD-1, x, y);
        lld ans = indb % d ? -1
            : lld((indb / d * x % (MOD-1) + (MOD-1)) % ((MOD-1)/d));
        cout << ans << endl;
    }

    return 0;
}

我竟然都忘了 extgcd 是怎么求多解的,我好菜啊QAQ

min25筛

这也是个用于求积性函数前缀和的筛法。

对于积性函数 f(n),满足这样几个要求时可以用min25来做。

  • n \in prime 时, f(n) 可以使用关于 n 的多项式表示。

  • n \in prime 时,f(n^k) 要能够快速计算答案。

预处理 Σf(p)

在min25筛的过程中,我们需要预处理多项式的 x^k 相关的这个东西

\sum_{i=2}^n i^k [i \in prime]

那么记所有质数的集合为 P,第 j 个素数为 P_jx 的最小质因子为 minp(x),我们需要预处理这样一个东西

g(n,j) = \sum_{i=2}^n i^k [i \in P ~ or ~ minp(i) > P_j]

举个例子吧

\begin{aligned} g(n,0) &= 2^k + 3^k + 4^k + 5^k + 6^k + 7^k + \cdots + n^k \\ g(2n,1) &= 2^k + 3^k + 5^k + 7^k + 9^k + 11^k + \cdots + (2n-1)^k \\ g(6n,2) &= 2^k + 3^k + 5^k + 7^k + 11^k + \cdots + (6n-5)^k + (6n-1)^k \end{aligned}

相当于在这个求和中,g(n,j) 将前 j 个素数的倍数全部去掉了,只留下了那些素数的 k 次方和。

\begin{aligned} g(2n,0) – g(2n,1) &= 4^k + 6^k + 8^k + 10^k + 12^k + \cdots + (2n)^k \\ &= 2^k \left( 2^k + 3^k + 4^k + 5^k + 6^k + \cdots + n^k \right) \\ \\ g(6n,1) – g(6n,2) &= 9^k + 15^k + 21^k + 27^k + 33^k + \cdots + (6n-3)^k \\ &= 3^k \left( 3^k + 5^k + 7^k + 9^k + 11^k + \cdots + (2n-1)^k \right) \end{aligned}

那么从 g(n,j-1) 转移到 g(n,j) 的过程中,我们考虑去掉 p_j^k 的倍数 \sum (p_j \prod p_{j+i}^{e_i})^k

差不多就是有这样一个转移

g(n,j) = g(n,j-1) – p_j^k \left( g\left(\left\lfloor\frac{n}{p_j}\right\rfloor, j-1\right) – \sum_{x=1}^{j-1} p_x^k \right)

为什么要减去前 j-1 个素数 k 次方和的转移呢?因为 (p_x p_j)^k 在这之前都被筛掉了啊。

如果 p_j^2 > n,那肯定没有东西能被删掉了。所以有这样一个转移

g(n,j) = \begin{cases} g(n,j-1), &p_j^2 \gt n \\ g(n,j-1) – p_j^k \left( g(\lfloor n / p_j\rfloor, j-1) – \sum_{x=1}^{j-1} p_x^k \right), &p_j^2 \le n \end{cases}

我们需要一些转换,才能把上面那个东西从递归变成迭代。

我们会发现,我们可能需要经常性的计算 g(\lfloor n/x\rfloor, j) 的值,那我们就一起维护了吧。

这里是预处理的代码。

// euler sieve, spk for sum p^k from p_1 ~ p_pcn
int pri[MAXN], pcn, spk[MAXN];
int m, Sqrt, id1[MAXN], id2[MAXN];
long long gk[MAXN], w[MAXN];

inline long long sum_ik(lld n); // for i=2..n, sum i^k

void init()
{
    m = 0, Sqrt = sqrt(n);

    for (long long L = 1, R; L <= n; L = R+1)
    {
        R = n / (n / L), w[++m] = n / L;
        (w[m] <= Sqrt ? id1[w[m]] : id2[R]) = m;
        gk[m] = sum_ik(w[m]);
    }

    for (int j = 1; j <= pcn; j++) // p_pcn^2 should be greater than n
    {
        for (int i = 1; i <= m && pri[j] <= w[i]/pri[j]; i++)
        {
            int d = (w[i]/pri[j]<=Sqrt) ? id1[w[i]/pri[j]] : id2[n/(w[i]/pri[j])];
            gk[i] -= (spk[j] - spk[j-1]) * (gk[d] - spk[j-1]);
        }
    }
}

其中,id1[i] 表示数论分块中整除值 \lfloor n/l \rfloor = i 的块编号,id2[n/i] 表示整除值 \lfloor n / l \rfloor = i 的块编号,w[s] 表示数论分块第 s 个块的整除值。这样比较节省空间。

然后考虑下面那个迭代。首先,w_i 就是 g 函数的第一个参数。

p_j^2 > w_i,也就是第 i 个块整除值的时候,没有必要继续处理了。内层的 i=1 \cdots m 相当于从大往小枚举块,也就不用利用滚动数组什么的记录上一轮 j-1 时的具体值了。

最后,init执行完时,gk[i] 就是第 i 个整除值 w_i 时,g(w_i, pcn) 的值了,也就是我们想要的东西了。

由于那个多项式可能会有好几个 p^x,所以都维护出来吧。

求 Σf(i)

我们刚才去掉了合数处的点值,那么我们还要根据积性函数的性质把那些点值加回来。

那么我们现在构造这样一个函数

S(n,j) = \sum_{x=2}^n f(x)~[ minp(x) \ge p_j ]

例如

\begin{aligned} S(20,3) &= f(5) + f(7) + f(11) + f(13) + f(17) + f(19) \\ S(20,2) – S(20,3) &= f(3) + f(9) + f(15) \\ S(20,1) – S(20,2) &= f(2) + f(4) + f(6) + \cdots + f(18) + f(20) \end{aligned}

我们需要从最大的质数开始,一个个往回补上点值。

S(n,j) = \sum_{k \ge j} \sum_{e} f(p_k^e)\left(1 + S\left(\left\lfloor n/p_k^e\right\rfloor, k+1\right)\right)

考虑 p_k^e 的倍数,显然当 e=0 的时候不能算上贡献;当 e>0 的时候,\sum p_k^e \prod p_{k+i}^{e_i} 的贡献由 f(p_k^e) 乘上 \sum f(\prod p_{k+i}^{e_i}) 得到。为什么还有个 1+ 呢?S(x,k) 里面可没有 p_j 的东西,不加就算不上 \sum f(p_k^e) 的贡献了呀。

那我们整理一下这个鬼一样的式子。

\begin{aligned} S(n, j) &= \sum_{k \ge j} f(p_k) + \sum_{k \ge j} \sum_{e} \left( f(p_k^e) S(\lfloor n/p_k^e\rfloor, k+1) + f(p_k^{e+1}) \right) \\ &= g(n, |P|) – \sum_{i=1}^{j-1} f(p_i) + \sum_{k \ge j} \sum_{e} \left( f(p_k^e) S(\lfloor n/p_k^e\rfloor) + f(p_k^{e+1}) \right) \end{aligned}

所以这一部分的代码可以写成

long long S(long long x, int y)
{
    if (x <= 1 || pri[y] > x) return 0;
    int d = x <= Sqrt ? id1[x] : id2[n/x];
    long long ans = g(x) - g(pri[y]-1);

    for (int i = y; i <= pcn && pri[i] <= x / pri[i]; i++)
    {
        long long p1 = pri[i], p2 = p1 * p1;
        for (int e = 1; p1 <= x / pri[i]; ++e, p1 = p2, p2 *= pri[i])
            ans += S(x/p1, i+1) * f(i,e) + f(i,e+1);
    }

    return ans;
}

其中 g(n)\sum_{i=1}^n f(i) [i \in prime],也就是上面式子中的 g(n,|P|) 可以从上面预处理的式子中得到。而 f(i,e) 就是 f(p_i^e)

最后你要求的 \sum_{i=1}^n f(i) = S(n,1) + 1

举例来说,我已知积性函数 f(p^e) = C_{e+A-1}^{A-1},那么 f(p)=A,我可以预处理上面预处理部分 sp0, g0 等内容,那么 g(x) - g(pri[y]-1) 应该写为 A * (g0[d]-(y-1)),而 f(i,e) 则写为 C(e+A-1, A-1)。对应例题为 HDU6537。

举例来说,我已知积性函数 f(p^e) = Ae+1,那么 f(p)=A+1,我可以预处理 sp0, g0 等内容,那么 g(x) - g(pri[y]-1) 应该写为 (A+1) * (g0[d]-(y-1)),而 f(i,e) 则写为 A*e+1。对应具体例题为 SPOJ DIVCNT3。

再举例来说,我已知积性函数 f(p^e) = p \oplus e,那么可以在 p=2f(2)=3,然后其他时候有 f(p)=p-1,预处理 sp0, sp1, g0, g1 等内容,那么 g(x) - g(pri[y]-1) 应该写为 g1[d]-g0[d]-sp1[y-1]+y-1,并特判 if (y == 1) ans += 2;,而 f(i,e) 则写为 pri[i]^e。对应具体例题为 LOJ6053。

后记

所以说这是不是很像埃式筛啊,所以它又被称为拓展埃式筛。

据说时间复杂度是

O\left( \frac{n^{\frac{3}{4}}}{\log n} \right)

咱也不会证,咱也不敢证。

2019杭电多校第三场 E. Easy Math Problem

题目想要对于给定 nk

\sum_{i=1}^n \sum_{j=1}^n \gcd(i,j)^k ~ lcm(i,j) ~[\gcd(i,j) \in prime]

那我们考虑一下枚举某个素数 d = \gcd(i,j) 的情况吧。

\sum_{d \in P} d^{k-1} \sum_{i=1}^n \sum_{j=1}^n ij ~ [\gcd(i,j) = d]

我们先来考虑一下

f(N) = \sum_{i=1}^N \sum_{j=1}^N ij [\gcd(i,j)=1]

先根据对称性,可以得到

\sum_{i=1}^N \sum_{j=1}^i ij [\gcd(i,j)=1] = \frac{f(N)+1}{2}

再根据非常显然的 \gcd(i,j)=1 \Leftrightarrow \gcd(i,j-i)=1,可以配对

\sum_{j=1}^i j [\gcd(i,j)=1] = \frac{i \varphi(i) + [i=1]}{2}

那么就有

\begin{aligned} \frac{f(N)+1}{2} &= \sum_{i=1}^N \sum_{j=1}^i ij [\gcd(i,j)=1] \\ &= \frac{\sum_{i=1}^N i^2 \varphi(i) +1}{2} \end{aligned}

所以说刚才求的 f(N) 就相当于对 i^2 \varphi(i) 求了个前缀和。根据杜教筛一般套路,给配上个 g(n) = n^2 就可以筛筛筛了。

f(n) = \frac{n(n+1)(2n+1)}{6} – \sum_{i=2}^n i^2 f\left(\left\lfloor\frac{n}{i}\right\rfloor\right)

再回到原式,其实相当于求

\sum_{d \in P} d^{k+1} \sum_{i=1}^{\lfloor n/d\rfloor} i^2 \varphi(i)

那么我们对 n 进行数论分块,就可以获得答案了。考虑块 [L,R] 内前面一个东西贡献的值,其实相当于

\sum_{i=L}^R i^{k+1} [i \in prime]

那么把它变成一个前缀和的差,构造一个积性函数

h(p) = p^{k+1}; h(p^e) = 0, e>1

然后考虑用min25筛,就可以做了。

由于min25筛需要快速求出 p^k 的前缀和,那么就得上拉格朗日插值一下。

推荐看 这篇文章 学习min25。我也要好好花时间学习一下啦 QAQ

#include <bits/stdc++.h>
using namespace std;
typedef long long lld;
const int MAXN = 7e6+5;
const int MOD = 1e9+7;
lld N, K, ans;

lld fpow(lld a, lld k) {
    lld b = 1;
    while (k) {
        if (k & 1) b = b * a % MOD;
        a = a * a % MOD;
        k >>= 1;
    }
    return b;
}

namespace interpolation
{
    typedef vector<lld> Poly;

    Poly& operator*=(Poly &a, const Poly &b) {
        Poly c(a.size()+b.size()-1, 0);
        for (int i = 0; i < a.size(); i++)
            for (int j = 0; j < b.size(); j++)
                (c[i+j] += a[i] * b[j] % MOD) %= MOD;
        return a = c;
    }

    Poly& operator*=(Poly &a, lld qwq) {
        for (lld &ai : a)
            ai = ai * qwq % MOD;
        return a;
    }

    Poly& operator+=(Poly &a, const Poly &b) {
        a.resize(max(a.size(), b.size()));
        for (int i = 0; i < b.size(); i++)
            (a[i] += b[i]) %= MOD;
        return a;
    }

    Poly Lagrange(const Poly &X, const Poly &FX) {
        Poly c;
        int n = X.size();
        for (int i = 0; i < n; i++) {
            Poly x({FX[i]});
            for (int j = 0; j < n; j++)
                if (j != i)
                    x *= Poly({-X[j],1}),
                    x *= fpow(X[i]-X[j], MOD-2);
            c += x;
        }
        return c;
    }

    lld evaluate(const Poly &a, lld x) {
        lld xx = 1, ans = 0;
        for (int i = 0; i < a.size(); i++)
            (ans += a[i] * xx % MOD) %= MOD,
            xx = xx * x % MOD;
        return ans;
    }

    Poly init_lag() {
        Poly a, b;
        lld last = 0;
        for (int i = 0; i <= K+2; i++) {
            a.push_back(i);
            (last += fpow(i, K+1)) %= MOD;
            b.push_back(last);
        }
        return Lagrange(a, b);
    }
}

namespace euler_seive
{
    int pri[MAXN], _B1[MAXN], pcn;
    bool isnp[MAXN];
    unordered_map<lld, int> _B2;
    const lld inv6 = fpow(6, MOD-2);

    void init_prime()
    {
        _B1[1] = 1;

        for (int i = 2; i < MAXN; i++)
        {
            if (!isnp[i])
            {
                pri[++pcn] = i;
                _B1[i] = i-1;
            }

            for (int j = 1; j <= pcn; j++)
            {
                if (1ll * i * pri[j] >= MAXN) break;
                isnp[i * pri[j]] = true;
                if (i % pri[j] == 0) {
                    _B1[i * pri[j]] = _B1[i] * pri[j];
                    break;
                } else {
                    _B1[i * pri[j]] = _B1[i] * (pri[j]-1);
                }
            }
        }

        for (int i = 2; i < MAXN; i++)
            _B1[i] = (_B1[i-1] + 1ll * _B1[i] * i % MOD * i % MOD) % MOD;
    }

    inline lld sum2(lld a) { a %= MOD; return a * (a+1) % MOD * (2*a+1) % MOD * inv6 % MOD; }

    int sumId2Phi(lld n)
    {
        if (n < MAXN) return _B1[n];
        else if (_B2.count(n)) return _B2[n];
        lld ans = n % MOD * (n%MOD+1) % MOD * ((MOD+1)/2) % MOD;
        ans = ans * ans % MOD; // \sum n^3
        for (lld l = 2, r; l <= n; l = r+1)
            r = n / (n / l),
            ans -= (sum2(r) - sum2(l-1)) * sumId2Phi(n/l) % MOD;
        ans = (ans % MOD + MOD) % MOD;
        return _B2[n] = ans;
    }
}

namespace min25_seive
{
    using euler_seive::pri;
    using euler_seive::pcn;
    interpolation::Poly _sumik_f;
    int _sumpk[MAXN], m, Sqrt;
    int id1[MAXN], id2[MAXN], tot;
    lld g[MAXN], w[MAXN];

    inline lld sumik(lld n) {
        return interpolation::evaluate(_sumik_f, n % MOD) - 1;
    }

    inline void init()
    {
        for (int i = 1; i <= pcn; i++)
            _sumpk[i] = (_sumpk[i-1] + fpow(pri[i], K+1)) % MOD;
        _sumik_f = interpolation::init_lag();

        m = 0;
        Sqrt = sqrt(N);

        for (lld l = 1, r; l <= N; l = r+1)
        {
            r = N / (N / l), w[++m] = N / l;
            (w[m] <= Sqrt ? id1[w[m]] : id2[N/w[m]]) = m;
            g[m] = sumik(w[m]);
        }

        tot = lower_bound(pri+1, pri+1+pcn, Sqrt) - pri;
        for (int j = 1; j <= tot; j++)
        {
            for (int i = 1; i <= m && pri[j] <= w[i]/pri[j]; i++)
            {
                int k = (w[i]/pri[j]<=Sqrt) ? id1[w[i]/pri[j]] : id2[N/(w[i]/pri[j])];
                (g[i] -= (_sumpk[j] - _sumpk[j-1]) * (g[k] - _sumpk[j-1]) % MOD) %= MOD;
            }
        }
    }

    inline lld sumpk(lld n) {
        return g[n <= Sqrt ? id1[n] : id2[N/n]];
    }
}

void solve()
{
    scanf("%lld %lld", &N, &K);
    min25_seive::init();

    ans = 0;
    for (lld l = 1, r, s = 0; l <= N; l = r+1)
    {
        r = N / (N / l);
        lld ss = min25_seive::sumpk(r);
        ans += euler_seive::sumId2Phi(N/l) * (ss-s) % MOD;
        s = ss;
    }

    printf("%lld\n", (ans % MOD + MOD) % MOD);
}

int main()
{
    euler_seive::init_prime();
    int T;
    scanf("%d", &T);
    while (T--) solve();
    return 0;
}

我好菜啊,反演不会推,min25也不会写,一道题写了四个小时到处是bug,到处TLE和MLE。

std好快啊。qwq

2019杭电多校第一场 K. Function

题目想要求

\sum_{i=1}^n \gcd\left(\left\lfloor\sqrt[3]i\right\rfloor,i\right) mod~998244353

首先可以想到,利用 \lfloor\sqrt[3]{i}\rfloor 的不同值进行分块。我们先记 m=\lfloor\sqrt[3]{n}\rfloor, 那么原式变为

\sum_{c=1}^{m} \sum_{i=c^3}^{\min((c+1)^3-1, n)} \gcd(c,i)

看着那个min比较碍眼,不如考虑当 c 的情况,那么有

H(c)=\sum_{i=c^3}^{(c+1)^3-1} \gcd(c, i) = \sum_{d|c} d \sum_{i=c^3}^{(c+1)^3-1} \left[\gcd(c,i)=d\right]

看起来有个什么类似于前缀和的东西?再设

f(n,d) = \sum_{i=1}^{n} \left[ \gcd(d,i) = 1 \right]

根据莫比乌斯函数的容斥意义容易得到

f(n,d) = \sum_{g|d} \mu(g) \left\lfloor\frac{n}{g}\right\rfloor

回到上面那个式子

\begin{aligned} H(c) &= \sum_{d|c} d \times \left( f\left(\left\lfloor\frac{(c+1)^3-1}{d}\right\rfloor,\frac{c}{d}\right) – f\left(\left\lfloor\frac{c^3-1}{d}\right\rfloor,\frac{c}{d}\right) \right) \\ &= \sum_{d|c} d \sum_{g|\frac{c}{d}} \mu(g) \left( \left\lfloor\frac{(c+1)^3-1}{gd}\right\rfloor – \left\lfloor\frac{c^3-1}{gd}\right\rfloor \right) \\ &= \sum_{gkd=c} d ~\mu(g) \left( \left\lfloor\frac{(c+1)^3-1}{gd}\right\rfloor – \left\lfloor\frac{c^3-1}{gd}\right\rfloor \right) \\ &= \sum_{kt=c} \left( \sum_{gd=t} d ~\mu(g) \right) \left( \left\lfloor\frac{(c+1)^3-1}{gd}\right\rfloor – \left\lfloor\frac{c^3-1}{gd}\right\rfloor \right) \\ &= \sum_{kt=c} \varphi(t) \left( \left\lfloor\frac{(c+1)^3-1}{t}\right\rfloor – \left\lfloor\frac{c^3-1}{t}\right\rfloor \right) \end{aligned}

好像可以发现 c | (c+1)^3-1,那后面那个玩意又有点碍眼,看看好像只有 t=1 的时候不太一样?不如我们……

\begin{aligned} H(c) &= \sum_{kt=c} \varphi(t) \left( \frac{(c+1)^3-1}{t} – \left\lfloor\frac{c^3}{t} \right\rfloor + 1 \right) \\ &= \sum_{kt=c} \varphi(t)~k~(3c+3) + \sum_{kt=c} \varphi(t)k \\ &= c + 3(c+1) (\varphi * id)(c) \end{aligned}

那我们继续看后面那部分吧

\begin{aligned} H'(m) &= \sum_{d|m} d \sum_{i=m^3}^n [\gcd(m,i) = d] \\ &= \sum_{d|m} d \sum_{g|\frac{m}{d}} \mu(g) \left( \left\lfloor\frac{n}{gd}\right\rfloor – \left\lfloor\frac{m^3-1}{gd}\right\rfloor \right) \\ &= \sum_{t|m} \varphi(t) \left( \left\lfloor\frac{n}{t}\right\rfloor – \left\lfloor\frac{m^3-1}{t}\right\rfloor \right) \end{aligned}

再回去看看最开始要的东西

H'(m) + \sum_{i=1}^{m-1} H(i)

不过 H(i) 看起来可以线性筛再前缀和?

好像 O(n^{\frac{1}{6}}) 单组询问的复杂度可以接受?Let’s go~

#include <bits/stdc++.h>
typedef long long lld;
using namespace std;
const size_t MAXN = 1e7+5;
const int MOD = 998244353;
bool isnp[MAXN];
int pcn, pri[MAXN/10];
lld g[MAXN], s[MAXN];
int phi[MAXN], q[MAXN];
int mod;

void init_prime()
{
    g[1] = q[1] = phi[1] = 1;

    for (int i = 2; i < MAXN; i++)
    {
        if (!isnp[i])
        {
            g[i] = s[i] = 2*i-1;
            q[i] = i-1;
            phi[i] = i-1;
            pri[pcn++] = i;
        }

        for (int j = 0; j < pcn; j++)
        {
            int k = i * pri[j];
            if (k >= MAXN) break;
            isnp[k] = true;

            if (i % pri[j] == 0)
            {
                s[k] = (s[i] + q[i] + 0ll) * pri[j];
                q[k] = q[i] * pri[j];
                g[k] = g[i] / s[i] * s[k];
                phi[k] = phi[i] * pri[j];
                break;
            }
            else
            {
                phi[k] = phi[i] * (pri[j] - 1);
                q[k] = q[pri[j]];
                s[k] = s[pri[j]];
                g[k] = g[i] * s[k];
            }
        }
    }

    for (int i = 1; i < MAXN; i++)
    {
        g[i] = (g[i-1] + i + (3 * i + 3) * (g[i] % MOD)) % MOD;
    }
}

int two(__int128 n)
{
    if (n < 8) return 1;
    int l = 1, r = 1e9+7, m;
    while (r - l > 1)
    {
        m = (l + r) >> 1;
        if (m > n / m / m)
            r = m;
        else
            l = m;
    }
    return r-1;
}

void solve()
{
    __int128 n = 0;
    char buf[30];
    scanf("%s", buf);
    for (int i = 0; buf[i]; i++)
        n = n * 10 + (buf[i] - '0');
    int m = two(n);
    lld ans = g[m-1];
    __int128 qwq = m; qwq *= m; qwq *= m; qwq -= 1;
    for (int i = 1; i * i <= m; i++)
    {
        if (m % i == 0)
        {
            ans += phi[i] * ((n / i % MOD) - (qwq / i % MOD)) % MOD;
            if (m != i * i)
                ans += phi[m / i] * (n / (m / i) % MOD - (qwq / (m / i) % MOD)) % MOD;
        }
    }
    ans %= MOD;
    if (ans < 0) ans += MOD;
    printf("%lld\n", ans);
}

int main()
{
    init_prime();
    int T;
    scanf("%d", &T);
    while (T--) solve();
    return 0;
}