Atcoder abc300

發表於 更新於 分類於 ， Disqus： 文章字數： 2.4k 所需閱讀時間 ≈ 8 分鐘
Ex(meet-in-middle,常系数齐次线性递推,Botsan-Mori)

https://atcoder.jp/contests/abc300/tasks

G P-smooth number

<= n 的，且所有质因子 <= p的数的个数

n 1e16

p 100

4s

1024mb

我的思路

dfs(上限, 质数列表idx) = sum dfs(上限/ 质数[idx]^pwr, idx-1)

然后对于底部做cache

然而并不理想，有cache和没cache 1e15都需要 8~12s,

跟别说1e18

另一个角度是，既然题目都告诉了最大的范围的答案大约是2e10, 所以它大也不算大,

所以考虑meet-in-middle

似乎就过了

代码

https://atcoder.jp/contests/abc300/submissions/me

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
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
#define rep(i,a,n) for (ll i=a;i<(ll)(n);i++)
ll read(){ll r;scanf("%lld",&r);return r;}

int main(){
  vector<int> p100;
  rep(i,2,100) { // 暴力
    bool ok=true;
    rep(j,2,i) if(i%j==0) ok =false;
    if(ok) p100.push_back(i);
  }
  map<pair<ll,int>,ll> ans; // ans[upbound,idx]

  ll n=read();
  auto g=[&](const vector<int> &pp){
    vector<ll> res = {1};
    for(auto p:pp){
      int sz = res.size();
      for(ll t=p;t<=n;t*=p){
        rep(j,0,sz) {
          if(res[j] > n/t) break;
          res.push_back(res[j]*t);
        }
      }
      sort(begin(res),end(res));
    }
    return res;
  };

  auto calc=[&](const vector<ll>&v0,const vector<ll>&v1) -> ll{
    ll ans = 0;
    int j=v1.size()-1;
    rep(i,0,size(v0)){
      while(j>=0 and (v0[i] > n/v1[j] or v0[i]*v1[j] > n) ) j--; // overflow
      ans += j+1;
    }
    return ans;
  };

  int p=read();
  vector<int> pp[2] = {{},{}};
  rep(i,0,size(p100)) {
    pp[i%2].push_back(p100[i]);
    if(p100[i] == p) break;
  }
  printf("%lld\n",calc(g(pp[0]),g(pp[1])));
  return 0;
}

Ex - Fibonacci: Revisited

$a_n = \begin{cases} 1 & (0 \leq n \lt K) \\ \displaystyle{\sum_{i=1}^K} a_{n-i} & (K \leq n). \\ \end{cases}$

给定$N$, 求 $\displaystyle \sum_{m, m \mathrm{ AND } N = m} a_m$

$K \in [1,5\cdot 10^4]$

$N\in [0,10^{18}]$

6s

1024mb

我的思路

完全没思路啊

只能说m 是N二进制表示的子集

如果 k = 2, 那就是经典的fib, 矩阵快速幂 能算一个，但是不能 同时类和

而这里注意到 m 有 N最高位 和 m 没有N最高位的一一对应

f(v) = (v,v-1,v-2,…,v-k+1) 的行向量

b 是长度k的列向量=(1,0,…,0)^T

$a_v = f(v) \cdot b$

$t_i=2^i$, 如果$v$的$i$位为$0$, 那么

$A$是转移矩阵

$a_{v | t_i} = f(v | t_i) \cdot b = (f(v) \cdot A^{2^i}) \cdot b$

因此,$I$为单位矩阵

$\displaystyle \sum_{m,m\mathrm{AND} N}a_m = f(0) \prod_{i, N \mathrm{AND} 2^i = 2^i} (I+A^{2^i}) \cdot b$

$f(0)=b^T$

问题变成如何计算 转移矩阵的 $2^i$ 幂次

似乎缺少线性代数的知识

因为这里K的确有点大 5e4了都, 别说做乘法，光是表示就很大，虽然本身是稀疏的，但是幂次以后应该就不是稀疏的了

题解

Digit DP with carries

有 $a_1,a_2,…,a_k$面额的硬币, 多少中方案组成N元(mod 998244353)

$k \le 100, S=\sum_{i} a_i \le 1000,N\le 10^{18}$

可以变成digit dp, 大概橙色难度

对于硬币a: newdp[x] = dp[x]+dp[x-a]+dp[x-2a]+dp[x-3a]+...

但对于使用次数可以拆分成 2的幂次, 于是变成 for i=0..infty, for a, newdp[x] = dp[x]+dp[x-a*2^i], 注意这样拆分后外层是循环幂次，而内层是循环a数组

而在执行到i时，注意到对于给定的N mod 2^i 不再会变

而如果 dp[N] 的贡献来源有dp[v],那么dp[S = v + 2^i(1?+2?+4?+8?+...)(a_1?+a_2?+...+a_k)]

及v mod 2^i == N mod 2^i

所以其它v mod 2^i != N的直接忽略掉

所以规模也可以下降了

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
// carry[i] = carry[i*2^t + (原始N)%(2^t)]
carry = vector<int>(2S+1,0) // 2S是因为每轮 新增长<=S,而结果每次除以2, l/2+S <= l
carry[0] = 1
while(N){ // 每轮处理完 需要处理掉 mod 2 != N mod 2的, 题解这里写成M了，看了半天
	for i = 1 to K { for k = 2S - a_i to 0 { // 倒续为了不重复
		// 实际对应的是 += [k - 2^(t=轮次) * a[i]], 而轮次在每轮规模压缩下被除以2了
		carry[k + a_i] += carry[k]
	} }
	for k = 0 to 2S + 1 {
		p = 2k + (N mod 2) // mod 2 != N mod 2的, 不会对答案有贡献 直接忽略
		carry[k] = (p < 2S+1) ? carry[p] : 0
	}
	N = floor(N / 2)
}
return carry[0]

解释digit DP

显然可以用生成函数表示上面的问题, $a_i$的选择方案就是

$\displaystyle 1+x^{a_i}+x^{2a_i}+x^{3a_i}+\cdots = \frac{1}{1-x^{a_i}}$

答案显然就是 $\displaystyle [x^{N}] \prod_{i\in[1,k]} \frac{1}{1-x^{a_i}}$

又$\displaystyle \frac{1}{1-x}=\frac{1+x}{1-x^2}=\frac{(1+x)(1+x^2)}{1-x^4}=\frac{(1+x)(1+x^2)(1+x^4)}{1-x^8}=\prod_{i\ge 0}(1+x^{2^i})$, 神奇的变换

$\displaystyle F(x)=\prod_{i\in[1,K]}(1+x^{a_1})$, 这里其实就对应于上面的每一轮处理了,

那么答案$\displaystyle [x^N] \prod_{i\ge 0} F(x^{2^{i}})$

引理 B,Q,S_i

B能展开成 $B(x) = \prod_{k \geq 0} Q_{k}(x^{2^k})$

S_i=0/1 为取奇偶列$\mathcal{S}i \left(\sum{n \geq 0} a_n x^n \right) = \sum_{m \geq 0} a_{2m + i} x^m$

令$B_n(x) = \prod_{k \geq 0} Q_{n + k}(x^{2^k})$ 这对应的就是每轮的剩余部分 去除以2的动作

有$\lbrack x^N \rbrack A(x) B_n(x)$

$= \lbrack x^{\lfloor N/2 \rfloor} \rbrack \mathcal{S}_{N \bmod 2}(A(x) B_n(x))$

$=\lbrack x^{\lfloor N/2 \rfloor} \rbrack \mathcal{S}{N \bmod 2} (A(x) Q_n(x)) \prod{k \geq 1} Q_{n+k}(x^{2^{k-1}})$

$= \lbrack x^{\lfloor N/2 \rfloor} \rbrack \mathcal{S}{N \bmod 2}(A(x) Q_n(x)) B{n+1}(x).$

回到digit DP 也就是令$Q_n=F(x)$, 更精简的伪代码

1
2
3
4
5
6
7
8
A = 1 // 生成函数
k = 0 // int
while N > 0 {
	A = S(N mod 2, convolution(A, Q_k))
	k += 1
	N = floor(N / 2)
}
return [x^0] A

Bostan-Mori algorithm and Application to linear recurrence relation

解决问题 $[x^{n}]\frac{P(x)}{Q(x)}$,其中 $P(x),Q(x)$的幂次比较小,而$n$十分大的情况

对于, linear recurrence relation. 同样适用

因为上面的转换相当于 $\displaystyle a_n = \sum_{i=1}^d c_ia_{n-i}, n\ge d$

所以对于$n \ge d$有 $\displaystyle [x^n](1-\sum_{i=1}^d c_ix^i)A(x) = 0$

令$Q(x)= 1-\sum_{i=1}^d c_ix^i$,注意到$A(x)Q(x)$的最大非零幂次$\le d-1$

$P(x) = A(x)Q(x) = (A(x) \mod x^d)Q(x) \mod x^d$

$A(x) = \frac{(A(x) \mod x^d)Q(x) \mod x^d}{Q(x)}$

Botsan-Mori algorithm

把$a_0,a_1,…$ 按照奇偶划分开

偶$a_0,a_2,a_4,\cdots$

奇$a_1,a_3,a_5,\cdots$

首先显然有 $\lbrack x^N \rbrack A(x) = \lbrack x^{\lfloor N/2 \rfloor} \rbrack \mathcal{S}_{N \bmod 2} (A(x))$

$\displaystyle \mathcal{S}_i\left(\frac{P(x)}{Q(x)}\right)$

$\displaystyle =\mathcal{S}_i\left(\frac{P(x)Q(-x)}{Q(x)Q(-x)}\right)$

注意到$Q(x)Q(-x)$ 如果是偶项和奇项相乘的 一正一负抵消，所以结果只有偶次项, 所以分母不改变奇偶

$\displaystyle =\frac{\mathcal{S}_i (P(x)Q(-x))}{\mathcal{S}_0 (Q(x)Q(-x))}$

$\displaystyle =\frac{\mathcal{S}_i (P(x)Q(-x))}{Q(x^{\frac{1}{2}})Q(-x^{\frac{1}{2}})}$

分子 就是 d log d 的卷积

1
2
3
4
5
6
7
8
9
input : P(x), Q(x), N`
output : [x^N] P(x) / Q(x)`

while N > 0 {
	P(x) = S(N mod 2, P(x) Q(-x))`
	Q(x) = S(0, Q(x) Q(-x))`
	N = floor(N / 2)`
}
return P(0) / Q(0)`

另一方面

$\displaystyle \frac{P(x)}{Q(x)}$

$\displaystyle = \frac{P(x)Q_0(-x)}{Q_0(x)Q_0(-x)}$

$\displaystyle = P(x)Q_0(-x) \times \frac{1}{Q_1(x^2)}$

$\displaystyle = P(x)Q_0(-x) \times \frac{Q_1(-x^2)}{Q_1(x^2)Q_1(-x^2)}$

$\displaystyle = P(x)Q_0(-x) \times Q_1(-x^2) \times \frac{1}{Q_2(x^4)}$

$\displaystyle = P(x) Q_0(-x) \times \prod_{n \geq 1} Q_n(-x^{2^n})$

至此 题目中 如果是求单点， 都可以 $O(K\log K \log N)$算出来

回到题目

上面的 Botsan-Mori 中

在循环里, $A(x) = S_{N\mod 2}(A(x))$

而这里要 $m \mathrm{AND} N = m$

所以也就是改成

$A(x) = S_0(A(x))+S_1(A(x)), N\mod 2=1$

$A(x) = S_0(A(x)), N\mod 2=0$

即可

代码

https://atcoder.jp/contests/abc300/submissions/42091386

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
#include <bits/stdc++.h>
#include <atcoder/modint>
#include <atcoder/convolution>
const int MOD=998244353;
using mint = atcoder::static_modint<MOD>;
using namespace std;
typedef long long ll;
#define rep(i,a,n) for (ll i=a;i<(ll)(n);i++)
ll read(){ll r;scanf("%lld",&r);return r;}
vector<mint> S(int i,const vector<mint>&A){
  vector<mint> res;
  rep(j,0,size(A)) if((j&1)==i) res.push_back(A[j]);
  return res;
}
vector<mint> operator+(const vector<mint>&A,const vector<mint>&B){
  vector<mint> res(max(size(A),size(B)),0);
  rep(i,0,size(res)) res[i] = (i<(int)size(A)?A[i]:0)+(i<(int)size(B)?B[i]:0);
  return res;
}
vector<mint> negX(const vector<mint>&A){ // A(-x)
  vector<mint> res=A;
  rep(i,0,size(A)) if(i&1) res[i] = -res[i];
  return res;
}

int main(){
  ll k=read();
  ll n=read();
  // Q(x) = 1 - sum_{i=1..k} ci x^i
  vector<mint> Q(k+1,-1);
  Q[0] = 1;
  // P(x) = ((A mod x^k) * Q) mod x^k
  auto P = atcoder::convolution(Q, vector<mint>(k, 1));
  P.resize(k);
  // 基于 Botsan_Mori 修改
  while(n){ // m & n == n, sum [x^m] P(x)/Q(x)
    auto negxq = negX(Q);
    auto p_negxq = atcoder::convolution(P,negxq); // P(x)Q(-x)
    // 原始Botsan_Mori, P=S(n&1,p_nexq);
    P = S(0,p_negxq) + ((n&1)?S(1,p_negxq):vector<mint>{});
    Q = S(0,atcoder::convolution(Q,negxq)); // Q(x)=S(0,Q(x)Q(-x))
    n/=2;
  }
  printf("%d\n",(P[0]/Q[0]).val());
  return 0;
}

总结

G

卡了一下，还是自己想出来, 不过atcoder的机子真的快，我这边2s+,atcoder的只有955ms，不过也可能是我有几个编译参数的原因

???什么 有人用和我想的一样的 记忆化 下降 卡过了??? 而且似乎它记忆的内容更多

Ex

题解各种 书写错误 看自闭了

第一次见这个 Digit DP with carries, 常系数齐次线性递推

线性递推 = 生成方程 的卷积 感觉现在看很显然，自己却没想到

bostan-mori 也是第一次见

似乎这个用线性代数也有办法优化，但是线性代数忘完了

参考

官方题解

luogu 常系数齐次线性递推 P4723

Ryuhei Mori’s article (Japanese)

Codeforces: Bostan-Mori, an elegant algorithm to compute k-th term of linear recurrence in O(dlgdlgk)

bostan-mori, arxiv 2008: A Simple and Fast Algorithm for Computing the N-th Term of a Linearly Recurrent Sequence