[BZOJ3434][Wc2014]时空穿梭(莫比乌斯反演)

Address

https://www.luogu.org/problemnew/show/P4152
https://www.lydsy.com/JudgeOnline/problem.php?id=3434
http://uoj.ac/problem/54

Solution

考虑枚举第一个点的坐标 (x1,1,x1,2,...,x1,n) ( x 1 , 1 , x 1 , 2 , . . . , x 1 , n )
最后一个点的坐标 (xc,1,xc,2,...,xc,n) ( x c , 1 , x c , 2 , . . . , x c , n )
那么容易得出,这两个点连成的线段上的整点(不包括端点)个数为:

gcd(xc,1x1,1,xc,2x1,2,...,xc,nx1,n)1 gcd ( x c , 1 − x 1 , 1 , x c , 2 − x 1 , 2 , . . . , x c , n − x 1 , n ) − 1

对答案的贡献:
Cc2gcd(xc,1x1,1,xc,2x1,2,...,xc,nx1,n)1 C gcd ( x c , 1 − x 1 , 1 , x c , 2 − x 1 , 2 , . . . , x c , n − x 1 , n ) − 1 c − 2

考虑枚举 n n 维向量 z 满足:
zi=xc,1x1,1 z i = x c , 1 − x 1 , 1

那么 x1,i x 1 , i 就有 mizi m i − z i 种合法取值。
得到答案为:
z1=1m1z2=1m2...zn=1mnCc2gcd(z1,z2,...,zn)1i=1n(mizi) ∑ z 1 = 1 m 1 ∑ z 2 = 1 m 2 . . . ∑ z n = 1 m n C gcd ( z 1 , z 2 , . . . , z n ) − 1 c − 2 ∏ i = 1 n ( m i − z i )

这是一个显然的 Mobius 反演。考虑把 gcd gcd 解决掉:
dz1=1m1z2=1m2...zn=1mnCc2d1i=1n(mizi)[gcd(z1,z2,...,zn)=d] ∑ d ∑ z 1 = 1 m 1 ∑ z 2 = 1 m 2 . . . ∑ z n = 1 m n C d − 1 c − 2 ∏ i = 1 n ( m i − z i ) [ gcd ( z 1 , z 2 , . . . , z n ) = d ]

=dz1=1m1dz2=1m2d...zn=1mndCc2d1i=1n(mizid)k|gcd(z1,z2,...,zn)μ(k) = ∑ d ∑ z 1 = 1 ⌊ m 1 d ⌋ ∑ z 2 = 1 ⌊ m 2 d ⌋ . . . ∑ z n = 1 ⌊ m n d ⌋ C d − 1 c − 2 ∏ i = 1 n ( m i − z i d ) ∑ k | gcd ( z 1 , z 2 , . . . , z n ) μ ( k )

=dkμ(k)z1=1,k|z1m1dz2=1,k|z2m2d...zn=1,k|znmndCc2d1i=1n(mizid) = ∑ d ∑ k μ ( k ) ∑ z 1 = 1 , k | z 1 ⌊ m 1 d ⌋ ∑ z 2 = 1 , k | z 2 ⌊ m 2 d ⌋ . . . ∑ z n = 1 , k | z n ⌊ m n d ⌋ C d − 1 c − 2 ∏ i = 1 n ( m i − z i d )

=dkμ(k)Cc2d1i=1nj=1midk(mijdk) = ∑ d ∑ k μ ( k ) C d − 1 c − 2 ∏ i = 1 n ∑ j = 1 ⌊ m i d k ⌋ ( m i − j d k )

=ui=1nj=1miu(miuj)d|uCc2d1μ(ud) = ∑ u ∏ i = 1 n ∑ j = 1 ⌊ m i u ⌋ ( m i − u j ) ∑ d | u C d − 1 c − 2 μ ( u d )

=ui=1n(mimiumiu(miu+1)2u)d|uCc2d1μ(ud) = ∑ u ∏ i = 1 n ( m i ⌊ m i u ⌋ − ⌊ m i u ⌋ ( ⌊ m i u ⌋ + 1 ) 2 u ) ∑ d | u C d − 1 c − 2 μ ( u d )

ni=1(mimiumiu(miu+1)2u) ∏ i = 1 n ( m i ⌊ m i u ⌋ − ⌊ m i u ⌋ ( ⌊ m i u + 1 ) 2 u ) 是一个 n n 次多项式,可以用 O(n2) 的时间求出,设 i i 次项为 ai
=ui=0naiuid|uCc2d1μ(ud) = ∑ u ∏ i = 0 n a i u i ∑ d | u C d − 1 c − 2 μ ( u d )

预处理 f(u,c,i)=uid|uCc2d1μ(ud) f ( u , c , i ) = u i ∑ d | u C d − 1 c − 2 μ ( u d ) 后:
=ui=0naif(u,c,i) = ∑ u ∏ i = 0 n a i f ( u , c , i )

可以得出,由于多项式 a a 的系数只有对 u 整除(下取整)的结果,故 a a 的取值不超过 O(nm) 种。
对下界分块之后使用 f f 的前缀和计算。
复杂度 O(mnc+Tn3m) (这能过?)
~~辣鸡 BZOJ 评测机卡常数~~

Code

#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#define For(i, a, b) for (i = a; i <= b; i++)
#define Mobius(i, n) for (i = 1; i <= n;)
using namespace std;

inline int read()
{
    int res = 0; bool bo = 0; char c;
    while (((c = getchar()) < '0' || c > '9') && c != '-');
    if (c == '-') bo = 1; else res = c - 48;
    while ((c = getchar()) >= '0' && c <= '9')
        res = (res << 3) + (res << 1) + (c - 48);
    return bo ? ~res + 1 : res;
}

const int MaxN = 1e5, E = 23, P = 14, N = MaxN + 5, MOD = 10007;
int tot, pri[N], miu[N], C[N][E], pw[N][P], f[N][E], sum[N][E][P],
n, c, m[P], a[P], b[P];
bool mark[N];

int Min(int a, int b) {return a < b ? a : b;}

void sieve()
{
    int i, j, k;
    For (i, 0, MaxN) C[i][0] = pw[i][0] = 1;
    For (i, 1, MaxN)
    {
        For (j, 1, Min(18, i))
            C[i][j] = (C[i - 1][j] + C[i - 1][j - 1]) % MOD;
        For (j, 1, 11) pw[i][j] = 1ll * pw[i][j - 1] * i % MOD;
    }
    mark[0] = mark[miu[1] = 1] = 1;
    For (i, 2, MaxN)
    {
        if (!mark[i]) miu[pri[++tot] = i] = -1;
        For (j, 1, tot)
        {
            if (1ll * i * pri[j] > MaxN) break;
            mark[i * pri[j]] = 1;
            if (i % pri[j] == 0) break;
            else miu[i * pri[j]] = -miu[i];
        }
    }
    For (k, 2, 20)
        For (i, 1, MaxN)
        {
            int dl = MaxN / i;
            For (j, 1, dl)
                f[i * j][k] = (f[i * j][k] + MOD +
                    C[i - 1][k - 2] * miu[j]) % MOD;
            For (j, 0, 11)
                sum[i][k][j] = (sum[i - 1][k][j] +
                    f[i][k] * pw[i][j]) % MOD;
        }
}

void prod(int x, int y, int n)
{
    int i;
    a[n + 1] = b[0] = 0;
    For (i, 1, n + 1) b[i] = a[i - 1] * y % MOD;
    For (i, 0, n + 1) a[i] = (a[i] * x - b[i] + MOD) % MOD;
}

void work()
{
    int i, j, _min, ans = 0;
    n = read(); c = read();
    For (i, 1, n) m[i] = _min = read();
    For (i, 1, n) _min = Min(_min, m[i]);
    Mobius (i, _min)
    {
        int nxt = m[1] / (m[1] / i);
        For (j, 2, n) nxt = Min(nxt, m[j] / (m[j] / i)), a[j] = 0;
        a[a[0] = 1] = 0;
        For (j, 1, n) prod(1ll * m[j] * (m[j] / i) % MOD,
            (1ll * (m[j] / i) * (m[j] / i + 1) >> 1) % MOD, j - 1);
        For (j, 0, n) ans = (ans + a[j] * (sum[nxt][c][j] -
            sum[i - 1][c][j] + MOD)) % MOD;
        i = nxt + 1;
    }
    printf("%d\n", ans);
}

int main()
{
    int T = read();
    sieve();
    while (T--) work();
    return 0;
}
  • 0
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值