YeXiaoRain Blog

發表於 更新於 分類於 Disqus： 文章字數： 2.9k 所需閱讀時間 ≈ 10 分鐘

问题$\pi (x)$

$[1,n]$ 计算之间素数的个数,其中$n$可能达到$10^{12}$

足够小的n(Sieve_of_Eratosthenes)

https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes

也就是 每次把 $i^2,i(i+1),i(i+2),\cdots$ 设为False

1
2
3
4
5
6
7
8
def primes_with_sieve(N):
    is_prime = [1 for i in range(N+1)]
    is_prime[0]=is_prime[1]=0
    for i in range(2,N+1):
        if is_prime[i]:
            for j in range(i*i,N+1,i):
                is_prime[j] = 0
    return is_prime

勒让德(Legendre’s Formula)

https://mathworld.wolfram.com/LegendresFormula.html

核心想法是容斥原理

$\phi(x,a) =$ 在$\le x$ 中,与前$a$个质数均互质的数, $x\in \mathbb{R},a\in \mathbb{Z}$

显然$\phi(x,a) = \phi(\lfloor x\rfloor,a)$

显然$\phi(x,\ge \pi(x))=\phi(x,\pi(x))=1$

$\displaystyle \phi(x,a)=\lfloor x\rfloor-\sum \lfloor \frac{x}{p_i}\rfloor+\sum \lfloor \frac{x}{p_ip_j}\rfloor - \sum \lfloor \frac{x}{p_ip_jp_k}\rfloor+\cdots$

而类似与上面的筛数法，也就是 $\phi(x,a=\pi(\sqrt{x}))$ 留下的是$1$以及 $(\sqrt{x},x]$之间的质数的个数，所以有

$\phi(x,a=\pi(\sqrt{x})) = \pi(x)-\pi(\sqrt{x})+1$

因此$\pi(x) = \phi(x,\pi(\sqrt{x}))+\pi(\sqrt{x})-1$

满足递归关系 $\phi(x,a)=\phi(x,a-1)-\phi(\frac{x}{p_a},a-1)$

proof:

根据定义 $\phi(x,a)$和$\phi(x,a-1)$ 中间的差值，其实就是$[1,x]$中 质因子含有$p_a$ 且不含有$p_{<a}$ 的数

也就是$v = p_a\cdot t \le x$且$p_{<a}\not{|}t$

$t\le \frac{x}{p_a}$

也就是 $\phi(\frac{x}{p_a},a-1)$

令$m_k =\prod_{i=1\cdots k} p_i$ 由前$k$个质数乘起来的

$\phi(m_k,k) = \phi(m_k)$ 也就是欧拉函数(Totient Function) 表示 不大于$x$且和$x$互质的数的个数

如果把$x$表示成 $x=sm_k+t$, $t\in[0,m_k]$

$\phi(sm_k+t,k) = s\phi(m_k)+\phi(t,k)$

这很显然，因为和前$k$个质数互质的数 一定形成 $[0,m_k)$的周期

对于$\displaystyle t\in (\frac{m_k}{2},m_k]$

$\phi(t,k)=\phi(m_k)-\phi(m_k-t-1,k)$

这个其实就是 $[1,m_k]$ 这一段切割+翻转了一下，$[1,t],[t+1,m_k]$,后面一段可以对应$[0,m_k-(t+1)]$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
def primes_with_sieve(N):
    is_prime = [1 for i in range(N+1)]
    is_prime[0]=is_prime[1]=0
    for i in range(2,N+1):
        if is_prime[i]:
            for j in range(i*i,N+1,i):
                is_prime[j] = 0
    return is_prime

def primes_from_prime_bool(arr):
    primes = [0]
    for i in range(2,len(arr)):
        if arr[i]:
            primes.append(i)
    return primes

def Sieve_of_Eratosthenes(N):
    return sum(primes_with_sieve(N))

print(f"{Sieve_of_Eratosthenes(10**2)=}")
print(f"{Sieve_of_Eratosthenes(10**4)=}")
print(f"{Sieve_of_Eratosthenes(10**6)=}")


def Legendre_s_Formula(N):
    import math
    sqrt_N=math.isqrt(N)
    # sieve for sqrt_N
    pi = primes_with_sieve(sqrt_N)
    primes = primes_from_prime_bool(pi)
    for i in range(2,sqrt_N+1):
        pi[i]+=pi[i-1]
    pi_sqrt_N=pi[sqrt_N]
    m_k = [1]
    for i in range(1,len(primes)):
        m_k.append(m_k[-1]*primes[i])
        if m_k[-1] > N:
            break

    cnt = pi_sqrt_N-1 # pi(N) = pi(sqrt(N))-1 +phi(n,pi(sqrt(n)))
    task = [(1,N,pi_sqrt_N)] # coef, x, a
    def euler_phi(val,k):
        for i in range(1,k+1):
            val*=primes[i]-1
            val//=primes[i]
        return val

    itr = 0
    while itr < len(task):
        coef,x,a = task[itr]
        # task.pop(0)
        if a==1:
            cnt += coef*(x-x//2)
        elif x<=1:
            cnt += coef*x;
        elif a>=x:
            assert x<=sqrt_N
            cnt += coef*1
        elif a < len(m_k) and x > m_k[a]:
            # phi(s m_k+t,k) = s phi(m_k,k) + phi(t,k)
            cnt += coef*(x//m_k[a])*euler_phi(m_k[a],a)
            task.append((coef            ,x%m_k[a],a))
        elif a < len(m_k) and x*2 > m_k[a] and x < m_k[a]:
            cnt += coef*euler_phi(m_k[a],a)
            task.append((coef*-1,m_k[a]-x-1,a))
        else:
            task.append((coef*-1,x//primes[a],a-1))
            task.append((coef   ,x               ,a-1))
        itr+=1

    return cnt


print(f"{Legendre_s_Formula(10**2)=}")
print(f"{Legendre_s_Formula(10**4)=}")
print(f"{Legendre_s_Formula(10**6)=}")
print(f"{Legendre_s_Formula(10**8)=}")

Meissel’s Formula

https://mathworld.wolfram.com/MeisselsFormula.html

$\lfloor x\rfloor=1+\sum_{1\le i\le a}\lfloor\frac{x}{p_i}\rfloor-\sum_{1\le i < j\le a}\lfloor\frac{x}{p_ip_j}\rfloor+\cdots+(\pi(x)-a)+P_2(x,a)+P_3(x,a)+\cdots$
$P_i(x,a) =$ 恰有$i$个质因子的数的个数, 且范围$\le x$,质因子不属于前$a$个

前x = 1 + (前a个质数的倍数的数) + (剩余质数) + (剩余和数)

例如

$P_2(x,a) =\sum_{i\ge a+1}[\pi(\frac{x}{p_i})-(i-1)]$

其实就是 枚举 $i \le j$中的最小的$i$

那么值是 $p_ip_{j\ge i} \le x$

也就有了上面 $p_j \le \frac{x}{p_i}$ 且 $j\ge i$

$P_3(x,a)=\sum_{i\ge a+1} P_2(\frac{x}{p_i},i-1)$ 同样的 先枚举最小的,这里似乎mathworld写错了，写成a了?

$\displaystyle =\sum_{i=a+1}^c\sum_{j=i}^{\pi(\sqrt{\frac{x}{p_i}})}[\pi(\frac{x}{p_ip_j})-(j-1)]$

那么这里其实 注意到 如果再向$P_4$ 走那么 $p_i \le x^{\frac{1}{4}}$

同时注意到$p_a < p_i \le x^{\frac{1}{4}}$

所以其实通过控制$a$的大小 可以控制$P_?$的零项

例如$a=\pi(\sqrt{x})$ 那么就是上面的勒让德，因为$P_{\ge 2}=0$

而Meissel’s的 表达式，就是 调整到让$P_{\ge 3}=0$

也就是$a=\pi(x^{\frac{1}{3}})$

回到上面 也就是

前x = 1 + (前a个质数的倍数的数) + (剩余质数) + (剩余和数仅含2个质因子的和数)
$\lfloor x\rfloor=1+\sum_{1\le i\le a}\lfloor\frac{x}{p_i}\rfloor-\sum_{1\le i < j\le a}\lfloor\frac{x}{p_ip_j}\rfloor+\cdots+(\pi(x)-a)+P_2(x,a)$
$\phi(x,a)=1+\pi(x)-a+P_2(x,a)$

$\pi(x) = \phi(x,a)-1+a-P_2(x,a)$

$\displaystyle P_2(x,a)=\sum_{i\ge a+1}^{b=\pi(x^{\frac{1}{2}})}[\pi(\frac{x}{p_i})-(i-1)]$

$\displaystyle =\sum_{i\ge a+1}^{b=\pi(x^{\frac{1}{2}})}\pi(\frac{x}{p_i})-\sum_{i\ge a+1}^{b=\pi(x^{\frac{1}{2}})}(i-1)$

$\displaystyle =\sum_{i\ge a+1}^{b=\pi(x^{\frac{1}{2}})}\pi(\frac{x}{p_i})-\frac{1}{2}(a+b-1)(b-a)$ 这里似乎mathworld也没写对?

$\displaystyle \pi(x) = \phi(x,a)-1+a -\sum_{i\ge a+1}^{b}\pi(\frac{x}{p_i})+\frac{1}{2}(a+b-1)(b-a)$,其中$a=\pi(x^{\frac{1}{3}}),b=\pi(x^{\frac{1}{2}})$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
def primes_with_sieve(N):
    is_prime = [1 for i in range(N+1)]
    is_prime[0]=is_prime[1]=0
    for i in range(2,N+1):
        if is_prime[i]:
            for j in range(i*i,N+1,i):
                is_prime[j] = 0
    return is_prime

def primes_from_prime_bool(arr):
    primes = [0]
    for i in range(2,len(arr)):
        if arr[i]:
            primes.append(i)
    return primes

def meissel(N):
    # $\pi(x) = \phi(x,a)-1+a -\sum_{i\ge a+1}^{b}\pi(\frac{x}{p_i})+\frac{1}{2}(a+b-1)(b-a)$,
    # 其中$a=\pi(x^{\frac{1}{3}}),b=\pi(x^{\frac{1}{2}})$
    import math
    sqrt_N=math.isqrt(N)
    # sieve for sqrt_N
    pi = primes_with_sieve(sqrt_N)
    primes = primes_from_prime_bool(pi) # 1-index 即是个数也是下标
    for i in range(2,sqrt_N+1):
        pi[i]+=pi[i-1]

    def meissel_recursive(x):
        # a = pi(x^1/3)
        a = 1
        for i in range(1,len(primes)):
            if primes[i]**3 > x:
                assert i > 1
                a = i-1
                break
        # assert a == pi[math.floor(x**(1/3))] 浮点数问题

        # b = pi(x^1/2)
        b = pi[math.isqrt(x)]

        # (a,b]
        cnt = a-1+(a+b-1)*(b-a)//2
        # -sum pi(x/p_i),i=a+1..b
        for p in primes[a+1:b+1]:
            if x//p < len(pi):
                cnt -= pi[x//p]
            else:
                cnt -= meissel_recursive(x//p)

        # +phi(x,a)
        task = [(1,x,a)] # coef, x, a
        def euler_phi(val,k):
            for i in range(1,k+1):
                val*=primes[i]-1
                val//=primes[i]
            return val

        itr = 0
        while itr < len(task):
            coef,x,a = task[itr]
            # task.pop(0)
            if a==1:
                cnt += coef*(x-x//2)
            elif x<=1:
                cnt += coef*x;
            elif a>=x:
                assert x<=sqrt_N
                cnt += coef*1
            else:
                task.append((coef*-1,x//primes[a],a-1))
                task.append((coef   ,x               ,a-1))
            itr+=1
        return cnt

    return meissel_recursive(N)

print(f"{meissel(10**2)=}")
print(f"{meissel(10**4)=}")
print(f"{meissel(10**6)=}")
print(f"{meissel(10**8)=}")

Lehmer’s Formula

https://mathworld.wolfram.com/LehmersFormula.html

目录分级，是Meissei的变化

其实就是Lehmer是控制$a$的范围让$P_{\ge 3}=0$

而这里实际就是控制$a,b$的范围，让$P_{\ge 4}=0$

$\displaystyle \pi(x) = \phi(x,a)-1+a-\sum_{a<i\le b} (\pi(\frac{x}{p_i})-(i-1))-\sum_{i=a+1}^{c}\sum_{j=i}^{\pi(\sqrt{\frac{x}{p_i}})}[\pi(\frac{x}{p_ip_j})-(j-1)]$

$a=\pi(x^{\frac{1}{4}})$

$b=\pi(x^{\frac{1}{2}})$

$c=\pi(x^{\frac{1}{3}})$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
import time
def primes_with_sieve(N):
    is_prime = [1 for i in range(N+1)]
    is_prime[0]=is_prime[1]=0
    for i in range(2,N+1):
        if is_prime[i]:
            for j in range(i*i,N+1,i):
                is_prime[j] = 0
    return is_prime

def primes_from_prime_bool(arr):
    primes = [0]
    for i in range(2,len(arr)):
        if arr[i]:
            primes.append(i)
    return primes

def Sieve_of_Eratosthenes(N):
    return sum(primes_with_sieve(N))

def n_1_div_pwr(n,pwr): # python3 只有isqrt没有1/3?
    # return int(n**(1/pwr)):
    l = 0
    r = n+1
    while r-l>1:
        mid = (l+r)//2
        if mid**pwr > n:
            r = mid
        else:
            l = mid
    return l


def meissel_lehmer(N):
    # $\pi(x) = \phi(x,a)-1+a
    #   - \sum_{i=[a+1,b]}(\pi(\frac{x}{p_i})-(i-1))
    #   - \sum_{i=[a+1,c]} \sum_{j=[i,\pi(\sqrt(x/p_i))(\pi(\frac{x}{p_ip_j})-(j-1))
    # 其中 $a=\pi(x^{\frac{1}{4}}),
    #       b=\pi(x^{\frac{1}{2}})$
    #       c=\pi(x^{\frac{1}{3}})$
    import math
    sqrt_N=math.isqrt(N)
    pi = primes_with_sieve(sqrt_N)
    primes = primes_from_prime_bool(pi) # 1-index 即是个数也是下标
    for i in range(2,sqrt_N+1):
        pi[i]+=pi[i-1]

    def inner_recursive(x):
        if x < len(pi):
            return pi[x]
        a = max(pi[n_1_div_pwr(x,4)],1)
        c = max(pi[n_1_div_pwr(x,3)],a)
        b = max(pi[n_1_div_pwr(x,2)],c)

        cnt = a-1
        # -sum pi(x/p_i)-(i-1),i=a+1..b
        for i in range(a+1,b+1):
            cnt -= inner_recursive(x//primes[i]) - (i-1)

        #   - \sum_{i=[a+1,c]} \sum_{j=[i,\pi(\sqrt(x/p_i))(\pi(\frac{x}{p_ip_j})-(j-1))
        for i in range(a+1,c+1):
            mxj = pi[math.isqrt(x//primes[i])]
            for j in range(i,mxj+1):
                cnt -= inner_recursive(x//primes[i]//primes[j])-(j-1)

        # +phi(x,a)
        task = [(1,x,a)] # coef, x, a
        itr = 0
        while itr < len(task):
            coef,x,a = task[itr]
            # task.pop(0)
            if a==1:
                cnt += coef*(x-x//2)
            elif x<=1:
                cnt += coef*x;
            elif a>=x:
                assert x<=sqrt_N
                cnt += coef*1
            else:
                task.append((coef*-1,x//primes[a],a-1))
                task[itr] = (coef   ,x           ,a-1)
                continue
            itr+=1
        return cnt

    return inner_recursive(N)


for i in range(1,11):
    start_time = time.time()
    print(f"meissel_lehmer\t{i=},{meissel_lehmer(10**i):15}\t{time.time()-start_time}s")

我尝试在这个里也加了m_k但更慢了

Lucy_Hedgehog

Lucy的代码本身解决的是一定范围内 质数的和，和这里其实差不多，个数无非是每个权重为1

$S(x,t)$ 为$[2,v]$ 所有整数中, 通过筛法循环筛完 前$t(\ge 0)\in \mathbb{Z}$个质数 后仍然剩余的数的个数

显然,这些数要么是$\le p_t$的素数, 要么这些数的最小素因子$> p_t$

$ans=S(n,\pi(\sqrt{n}))$

  1. $S(x,0) = \sum_{i=2}^x 1 = x-1$ 所有数字都还没被筛掉
  2. $p_t^2 > x$, $S(x,t)=S(x,t-1)$
  3. $p_t^2 \le x$ 这种情况相对复杂

$\displaystyle S(x,t) = S(x,t-1) - \sum_{j\in[2,x],p_t是j的最小素因子}1$, 就是考虑 多删除$p_t$的倍数会多删除什么

$\displaystyle S(x,t) = S(x,t-1) - \sum_{j\in[1,\lfloor \frac{x}{p_t}\rfloor],j的最小素因子 \ge p_t}1$ 考虑$j_{new}=j_{old}/p_t$

$\displaystyle S(x,t) = S(x,t-1) - (S(\lfloor \frac{x}{t} \rfloor,t-1) - (t-1))$

其实这里已经可以看出，lucy的代码数学上并不更优秀，而是lucy用了记忆化和递推 来完成了而大量的cache，因为注意到 x 的值都是x除以多个p得到的

那么其实只有 [0...sqrt(x)][x/sqrt(x)...x/1], $O(\sqrt{x})$个可能的取值，于是可以优化

而实际上，当我们单点求值时 还是像 勒让德 类似的方法只是增加了$[x/sqrt(x)…x/1]$的记忆化

而 例如pe501需要 多点求值 时，其实更多的还是 只记忆化 $x^{1/2}$或者$x^{2/3},x^{3/4}$ 来完成

多点求值 会更优吗？

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
#define rep(i,a,n) for (ll i=a;i<(ll)n;i++)
#define per(i,a,n) for (ll i=n;i-->(ll)a;)
template<class T>using minPQ=priority_queue<T,vector<T>,greater<T>>;
#define i64 long long

inline i64 v2idx(i64 v, i64 N, i64 Ndr, i64 nv) {
  return v >= Ndr ? (N/v - 1) : (nv - v);
}

i64 primesum(i64 N) {
  i64 r = (i64)sqrt(N);
  i64 Ndr = N/r;
  assert(r*r <= N and (r+1)*(r+1) > N);
  i64 nv = r + Ndr - 1;
  i64 *V = new i64[nv];
  i64 *S = new i64[nv];
  for (i64 i=0; i<r; i++) V[i] = N/(i+1); // [N/1...N/r]
  for (i64 i=r; i<nv; i++) V[i] = V[i-1] - 1; // [N/r-1...1]
  for (i64 i=0; i<nv; i++) S[i] = V[i]  - 1; // [N/1-1...N/r-1 N/r-2...0]

  for (i64 p=2; p<=r; p++) if (S[nv-p] > S[nv-p+1]) { // is prime
    i64 sp = S[nv-p+1]; // S(p,p) = pi(p)
    i64 p2 = p*p;
    for (i64 i=0; i<nv; i++) {
      if (V[i] >= p2) {
        S[i] -= (S[v2idx(V[i]/p, N, Ndr, nv)] - sp);
      } else {
        break;
      }
    }
  }

  return S[0];
}

int main() {
  i64 N = 100'0000'0000;
  printf("%lld\n", primesum(N));
  return 0;
}

完整 代码

练习题目

project euler 501

codeforces 665F

参考

wikipedia prime counting function

mathworld.wolfram prime counting function

求十亿内所有质数的和，怎么做最快? - 穆罕默德·李的回答 - 知乎

codeforces 665F tutorial

codeforces blog counting primes

oiwiki Min_25筛

Lucy_Hedgehog’s method

Meisell-Lehmer’s method

Youtube- The Pattern to Prime Numbers?

https://codeforces.com/blog/entry/44466

https://www.cnblogs.com/LzyRapx/p/8453170.html

github primecount

library-checker

發表於 更新於 分類於 ， Disqus： 文章字數： 1.6k 所需閱讀時間 ≈ 5 分鐘
E(二进制,最短路)F()G()

https://atcoder.jp/contests/abc355

E - Guess the Sum

交互题

有隐藏数组 a[0,2^n-1]

每次只能询问 sum a[2^i * j,2^i *(j+1)] mod 100

求 sum a[l..r] % 100

要求 询问次数最小

我的思路

首先 进行 二进制拆解有一个方案是

1
2
3
4
5
6
7
8
9
query(l,r,ql,qr,i,j):
  assert(l <= ql and qr <= r)
  if(l == ql and qr == r) return ask(i,j);
  m = (l+r)/2
  if(ql < m) ret += query(l,m,ql,min(m,qr),i-1,j*2) 
  if(qr > m) ret += query(m,r,max(m,ql),qr,i-1,j*2+1) 
  return ret

ans = query(0,1 << n,l,r,n,0)

这个问题是 不能保证最小次数

例如 求a[0..14],那么 只需要 a[0..15] - a[15..15] 两次 而不是3次

这个从二进制上看单侧

14 = 1110, 也就是连续的最多需要两次可以计算出

所以想的是 l 的二进制 变为 r的二进制过程中， 只能 加 或 减 2的幂次，且幂次不能超过当前值的最高位的1

所以 l,r的高位相同的地方忽略掉，处理低不同的位 例如1100 -> 1110，高位1100相同，只用考虑00->10

所以贪心吗，感觉”实现题”写着都好难受

閱讀全文 »

發表於 更新於 分類於 ， Disqus： 文章字數： 570 所需閱讀時間 ≈ 2 分鐘
G(Dilworth, 网络流, kmp)

https://atcoder.jp/contests/abc354

G - Select Strings

n个字符串si,每个字符串有价值ai,选一个字符串集合，使得集合中两两没有子串(连续)关系, 求最大价值和

n 100

sum |si| 5000

ai [1,1e9]

2s

1024mb

我的思路

这字符串可以 kmp si#s0#s1#s2.. 计算出 是否是子串关系

但之后怎么求不知道了

如果是依赖的选择 可以转化成 网络流

这里是互斥的选择

閱讀全文 »

發表於 更新於 分類於 ， Disqus： 文章字數： 4.8k 所需閱讀時間 ≈ 16 分鐘
G(top tree, static top tree)

https://atcoder.jp/contests/abc351

G - Hash on Tree

给定n点树有根数

值数组a[n]

f(u) = a[u] , 如果u是叶子

f(u) = a[u] + prod f(child of u)如果u非叶子

q次操作,每次 修改 a[pos]=val,每次修改后输出 f(1)

n,q 2e5

4s

1024mb

我的思路

首先这个是收到深度控制的,如果暴力的更新的话

那么一个想法是能否flat掉这个行为

[a,b,c]x[x,y,1] = [0,a*1+b*y,1]

leaf = [0,ai,1]

非叶子 = [ai,1,0]

然而没找到 矩阵方法或者满足“结合率”的东西

如果能有 矩阵或者结合律满足，就是个线段树维护，就好了

其中 上面的乘还可以拆成左右的包裹，总之如果有办法就好了

但是在尝试中发现要得到 a1+by 这种总是伴随秩的下降，没啥办法

第二个想的是上生成方程, 也没有弄出东西

那么最后就是 能否树的结构通过 轻重链 来完成快速计算(q log n)? 也没想到具体维护方法

閱讀全文 »

發表於 更新於 分類於 ， Disqus： 文章字數： 3k 所需閱讀時間 ≈ 10 分鐘
G(概率论)H(二分,树)I(TODO)

https://codeforces.com/contest/1951

G Clacking Balls(2s,256mb)

1~m篮子 顺时针 环

n个球，第i个在 篮子ai中，每个篮子至多1个球

alice 可以进行1s操作

  • [1,n]等概率 选i
  • 如果 球已经不存在，则跳过（依然耗时）
  • 如果 球存在，则沿着把球顺时针移动一个篮子，（如果目标篮子有球，把目标篮子中原来的球抛弃）

直到只剩下一个球停止，求操作时间的期望值 mod 1e9+7

n 3e5

n,m 1e9

我的思路

不能说没有想法，只能说一点想法也没有

如果一个球消失，那么说明有前面的球超过了它，换个角度 如果 球i移动 导致球j消失，也就是i走到j的位置，那么换成i消失是等价的后续，不影响概率

所以转换了题意的第三条（谁移动谁消失）以后

如果说 最后剩下的是 第i个球，那么这个球没有触碰下一个球，而前面的所有球触碰了后一个球

所以是个min/max 容斥+期望 吗？

计算每个球的 消失(触碰后面球)的期望次数，再min-max 容斥掉?

感觉也不对

再换个角度， 如果i留下

1
[i-1 ................i] .....

没啥想法，总的来说，显然直接状态描述肯定不行的，不论是点坐标还是距离

那么要么就是孤立每个 间隔的期望什么的，要么就是孤立的什么点期望，总之需要一个办法把其中 关联的东西拆掉，

而且不只如此，还有环

然后就是考虑简化的问题， 同样是环，但是只有两个点

那么 Exp(state) = Exp(点的距离)

然而转化构成了 一个转化矩阵

$A\lambda=\lambda$的形式，已知$A$要求$\lambda$

另一个是展开环到无限的数轴，指定目标点，那么得到每个点的移动方案数？

閱讀全文 »
0%