Atcoder abc222

發表於 2022-08-08 更新於 2024-04-11 分類於 Atcoder ， ABC Disqus：文章字數： 4.1k 所需閱讀時間 ≈ 14 分鐘

G(数学, 欧拉定理,phi,miller-robin,pollard-rho)H(题意转化,DP,FFT,生成函数,拉格朗日反演,P-recursive )

G - 222

https://atcoder.jp/contests/abc222/tasks/abc222_g

在数列2,22,222,2222,22222,….中

N个X, 首个是 Xi的倍数的下标是?, 或者不存在

范围

N 200

Xi [1,1e8]

我的思路

一眼看上去很数学, 很像Project Euler的题

$2222222 = 2 * 1111111 = 2 * \frac{(10^7 - 1)}{9}$

其实就是问对于x

是否 2 * (10^7 - 1) = 9 k x

首先x的2的幂次为0/1

好像有点绕

$kx = 1111111 = 10^0+10^1+10^2+\cdots$

右边虽然项数为合数时可以拆分, 例如$6 = 3 * 2$, $111111 = 111 \cdot 1001 = 11 \cdot 10101$

但不知道是否能拆出所有

另一个就是对于比较小的11111的部分,可以pollard-rho分解

考虑长除法?

每次除法取mod 乘10 加1

但1e8 不知道效率怎么样

PE 129 做过类似的, 但是问题是首个让是它倍数的最小$111\cdots 111$长度超过百万的是哪个因子

而有一些可用的结论

除了上面$2,5$因子外$kx = 111\cdots 111$始终有解, 且$111\cdots 111$ 的长度不超过$n$ (因为模数随着长度变化成环)

因此如果暴力的话, 期望值是在 $O(NAi)$ 的

想了下打表超过1e6的记录下来, 未超过的现场算, 但很多超过1e6的

int two(int v){
  int c = 0;
  ll m = 0;
  do{
    m*=10;
    m+=2;
    c++;
    m%=v;
  }while(m!=0);
  return c;
}

void calc(){
  rep(i,1000000, 100000000+1){
    if(i % 1000000 == 0) printf("progress %lld\n",i/1000000);
    if(i % 4 == 0 || i % 5 == 0)continue;
    int res = two(i);
    if(res > 1000000) printf("ans[%lld] = %d\n",i,res);
  }
}

另一个就是根据PE129的证明过程, 反正有$\phi(n)$ 或者$\phi(9n)$ 是一个解

那么可以找$\phi(n) , \phi(9n)$的因子尝试, 但这样是否能保证是最小的呢????? 根据倍数, 显然最小的是这个解的因子

$\phi(n) = n \cdot (1-1/p1) \cdot (1-1/p2) \cdots$

似乎可做?

题解

题解说不需要这么多, 就欧拉定理+暴力找phi就够了

看来我用高级算法乱搞了太多

代码

https://atcoder.jp/contests/abc222/submissions/33867609

16ms 还不是最快的, 有人10ms

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

typedef long long ll;
typedef __int128 lll;
#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;)
#define pb push_back

ll read(){ll r;scanf("%lld",&r);return r;} // read


ll quick_p(ll b, ll p,ll mod){
  ll r = 1;
  while(p){
    if(p%2)(r*=b)%=mod;
    (b*=b)%=mod;
    assert(r>0);
    assert(b>0);
    p/=2;
  }
  return r%mod;
}

bool is_prime_32(ll v){
  if(v == 2)return true;
  if(v < 2)return false;
  if(v%2 == 0)return false;
  ll test_g[] = {2, 325, 9375, 28178, 450775, 9780504, 1795265022};
  ll startp = v-1;
  while(startp%2 == 0)startp>>=1;
  rep(i,0,7){
    ll p = startp;
    ll base = test_g[i];
    // don't break may cause 4033 bug
    if(base % v == 0)continue;
    bool find = false;
    ll r = quick_p(base,p,v);
    while(p != v-1){
      if(r == v-1){
        find = true;
        break;
      }
      // -1 开始的序列, 或全1序列
      if(r == 1){
        if(p == startp){
          find = true;
          break;
        }
        return false;
      }
      p*=2;
      (r*=r)%=v;
    }
    if(!find){
      return false;
    }
  }
  return true;
}

ll my_sqrt(ll v){
  assert(v > 1);
  ll l = 1;
  ll r = v; // care overflow
  ll ret = 1;
  while(l < r){
    ll m = (l+r)/2;
    ll m2 = m*m;
    if(m2 == v) return m;
    if(m2 < v) {
      ret = m;
      l = m + 1;
    } else {
      r = m - 1;
    }
  }
  return ret;
}

ll randint(ll low,ll hi){
  return low + (rand() % static_cast<int>(hi - low + 1));
}

ll Pollard_Rho(ll N) { // 返回一个> 1的因数
  assert(N > 1);
  if (N == 4) return 2;
  ll ret = my_sqrt(N); // 质数平方 效率低 提前判断
  if(ret * ret == N) return ret;
  while(true) {
    ll c = randint(1, N - 1); // 生成随机的c
    auto f = [=](ll x) { return ((lll)x * x + c) % N; }; // ll 表示__int128，防溢出
    ll t = 0, r = 0; // 初始两个相同
    do{
      t = f(t); // 1倍速度
      r = f(f(r)); // 2倍速度
      ll d = gcd(abs(t - r), N);
      if (d > 1 && d < N) return d;
    }while (t != r);
  }
}

// 分解x为质因数, sorted, {prime,power}
vector<pair<ll,int> > fenjie(ll x) {
  vector<int> res = {};
  deque <ll> arr = {x};
  while(arr.size()){
    ll v = arr.front();
    arr.pop_front();
    if(v == 1) continue;
    if(is_prime_32(v)) {
      res.pb(v);
      continue;
    }
    ll divisor = Pollard_Rho(v);
    arr.push_back(divisor);
    arr.push_back(v/divisor);
  }
  sort(res.begin(),res.end());
  vector<pair<ll,int> > ret;
  rep(i,0,res.size()){
    if(i == 0|| res[i] != res[i-1]) ret.push_back({res[i], 1});
    else ret.back().second++;
  }
  return ret;
}

ll phi(ll n){
  auto primes = fenjie(n);
  // printf("%lld =",n);
  // for(auto [v,pwr]:primes) printf("[%lld %d]",v,pwr);
  // printf("\n");
  ll ret = n;
  for(auto [v,pwr]:primes) ret = (ret/v)*(v-1);
  return ret;
}

// -------------------------- lib -------------------------- 

int n ;

void dfs(int idx, vector<pair<ll,int>> primes, ll mul, ll & ans){
  if(idx == (int)primes.size()){
    // test pwr 10^p = k 9v + 1,  10^p % 9v == 1
    if(quick_p(10,mul,9*n) == 1) ans = min(ans,mul);
    return ;
  }
  rep(pwr,0,primes[idx].second+1){
    if(mul > ans) return;
    dfs(idx+1,primes,mul,ans);
    mul *= primes[idx].first;
  }
}

void w(){
  n = read();
  if(n % 4 == 0 || n % 5 == 0){
    printf("-1\n");
    return ;
  }
  if(n % 2 == 0) n /= 2;
  ll phin = n % 3 == 0 ? phi(9*n) : phi(n);
  if(phin == 1){
    printf("1\n");
    return ;
  }
  auto primes = fenjie(phin);
  ll ans = phin;
  dfs(0, primes, 1, ans);
  printf("%lld\n",ans);
}


int main(){
  int t = read();
  while(t--) w();
  return 0;
}

H - Beautiful Binary Tree

https://atcoder.jp/contests/abc222/tasks/abc222_h

给定N, 问多少个满足条件的有根二叉树, 每个点上有数字 0 或 1, 叶子点上都是1

至多n-1次操作, 让根上值为N, 其它所有点的值为0

每次操作, 把一个点的值全部加到它的父节点或,父节点的父节点上

答案mod 998244353

限制

N 1e7

1024 mb

我的思路

首先叶子上全是1, 且n-1次内全部移动完, 限制了最大的树的高度

而这里二叉树还可能有多个点只有一个子节点的

至于如何操作呢

那感觉上也是贪心从下向上,

而如果叶子向上看是 x-1-?的形式, 那么中间的一定会操作, 所以叠加上去, 0-(1+x)-?

而如果是 x-0-?, 那就直接跳过

这里要注意的是可能有长成这样的

因此顺序应该是从深度从大到小,而不是所有叶子做bfs

再看如果给定图做dfs的话,

dfs(i) 表示把低层的都收集到i的次数

那么对于一个节点 u-v-k

v 原来是1, 那么次数 = dfs(v) + 1

如果v 原来是0, 那么次数 = sum (dfs(k) + 1), k 是v的所有子节点

换句话说, dfs过程中一部分是在合并和, 还有一部分是在+1

所以本质上,能让所有的和 = N, 就要看所有+1的来源, 当然根上可以直接放1

又注意到上面写的 dfs转移方程式, 其实每次+1, 对应一个次移动

那么一共n-1次+1, 也就意味根上一定是1

然后感觉上, 可以考虑左右树拆分

左树贡献 i的话, 右侧贡献为 n-1 - i

两边独立, 似乎就可以fft/ntt 来搞了

然后如何变成和主问题等价的子问题呢?

考虑其中一个子节点让它对根贡献i, 记作$h(i)$

子节点为空, 则贡献i = 0,方案1

子节点为1, 则贡献为 i = 子树贡献(和原问题等价) + 1

子节点为0, 则考虑它的子节点, 因为它不能是叶子,它至少有一个子节点

那么1个的情况 i = 子树贡献() + 1, 方案数 x 2

那么2个的情况 i = 左子树贡献x + 1 + 右子树贡献y + 1, 方案数加和,

有点问题是这样下面贡献可能根是非1的, 因此f(x) 的意义改成产生的+1贡献, 根也是0的情况

$f(x) = \sum_{i=0}^{x} h(i)\cdot h(x-i)$,

$h(0) = 1$ // 对应无节点

$h(x) = f(x-1) + 2 * f(x-1) + \sum_{i=0}^{x-2} f(i)\cdot f(x-2-i), x > 0$

答案就是$f(n-1)$

这种自身相互依赖的用cdq二分好像能做?吗?

** 好像我的过程漏掉了总和, 只考虑操作步数……… 推了半天推了个锤子,白推了 **

官方题解

满足条件的树的充要

根和叶子都是1
所有1的和 = N
0不连续

证明

首先所有1 要汇总到根, 所有不在根上的1, 至少操作1次, 而最多n-1次,因此根有1,且所有其它1恰好被操作1次,而非1的地方不被操作

因此也不能有连续的0

问题变成统计上面的树的个数了

定义,对于$i > 0$

$a_i = i$个点是1,根也是1的满足要求的树的方案数

$b_i = i$个点是1,根是0的,满足剩余要求的树的方案

因为没有连续零,那么bi 要么有单个子树 $2 a_i$, 要么两个子树都不为空

所以 $b_i = 2a_i + \sum_{j=1}^{i-1} a_ja_{i-j}$

类似的,对于$a_i$

$a_1 = 1$

一个子节点时 $2(a_{i-1}+b_{i-1})$

所以 $a_i = 2(a_{i-1}+b_{i-1}) + \sum_{j=1}^{i-2} (a_j+b_j)(a_{i-1-j} + b_{i-1-j}), i > 1$

~~一点简化?(并不是) 是不是令$a_0 = 1,b_0 = 0$ 可以让上面变成完全的求和式子~~

这里也说分治类fft 可以做到 $N log^2 N$, 虽然没试过两个怎么做分治, 但会超时

生成方程

$a_0 = b_0 = 0$

分别把$a_i$和$b_i$作为系数做它们的生成方程$A(x),B(x)$

那么第一个表达式和$B=2A+A^2$等价

第二个和$A = x + 2x(A+B) + x(A+B)^2$

然后两个生成式带入一下

$A = x(1+A+B)^2 = x(1+3A+A^2)^2$

用 Newton’s algorithm 据说可以 $O(N log N)$, 也会超时

Lagrange inversion theorem 拉格朗日反演?

解决的问题, 给定F(x)

找G(x) 使得 G(F(x)) = x

这里有一点要用到,但是没有证明的是 $G(F(x)) = x$ 则 $F(G(x)) = x$, 只能说在$F(x)$的值域内$F(G(F(x)) = F(x)$即$F(G(x)) = x$可以证明, 但在值域外不知道如何证明…

如果 F,G满足, 且它们0次项系数均为0, 1次项系数均非0

那么有$\lbrack x^n \rbrack G(x) = \frac{1}{n} \lbrack x^{n-1} \rbrack \left(\frac{x}{F(x)}\right)^n.$

即是G的n次项 = $(\frac{x}{F(x)})^n$ 的$n-1$次项的$\frac{1}{n}$

证明

辅助lemma: 对于任何0次系数为0,存在非0次系数不为0的$F(x)$, 有对于整数$k$

$\lbrack x^{-1} \rbrack F’(x) F(x)^k = \lbrack k = -1\rbrack$, 即是$k=-1$时$-1$次系数为$1$,否则$-1$次系数为$0$

证明lemma:

对于$k\neq -1$, 显然求导法则$F’(x) F(x)^k = \frac{\left ( F(x)^{k+1} \right)’}{k+1}$

对于$k = -1$, $F(x) = \sum_{i>0} a_i x^i$

$\frac{F’(x)}{F(x)} = \frac{a_1+2a_2x+3a_3x^2+\cdots}{x(a_1+a_2x+a_3x^2+\cdots}= x^{-1} \frac{1 + 2\frac{a_2}{a_1}x + \cdots}{1 + \frac{a_2}{a_1}x + \cdots}$

也就是右侧这个分式除完以后是$1+k_1x+k_2x^2+\cdots$的样子, 因此 -1 次方的系数是 1, lemma 证毕.

因此$G(F(x)) = x$, 的$G$ 满足条件

$G’(F)\cdot F’ = 1$ ( 同时求导

展开$\sum_i i(\lbrack x^i\rbrack G(x) ) F^{i-1} F’ = 1$, (基本的求导法则 $(ax^i)’ = iax^{i-1}$, 注意到前两项都是系数而非生成函数

$\sum_{i} i (\lbrack x^i \rbrack G (x)) F^{i-1-n} F’ = F^{-n}.$ ( 同乘上$F^{-n}$

$\lbrack x^{-1}\rbrack \left(\sum_{i} i (\lbrack x^i \rbrack G (x)) F^{i-1-n} F’\right) = \lbrack x^{-1}\rbrack \left(F^{-n}\right).$ (提取$-1$次项目的系数

因为左侧$i (\lbrack x^i \rbrack G (x)) $ 整个都是系数,以及上面的lemma, 左侧只有$i=n$时生成系数的内容才为$1$,其它则是$0$

$n[x^n]G = [x^{-1}]F^{-n}$

变形一下,就有了最初要证明的$\lbrack x^n \rbrack G(x) = \frac{1}{n} \lbrack x^{n-1} \rbrack \left(\frac{x}{F(x)}\right)^n.$

这玩意百科上说建立了函数方程和幂级数之间的联系

使用示例: 卡特兰数Catalan number

$c_0 = 1$

$c_{n+1} = \sum_{i=0}^n c_{i} \cdot c_{n-i}$

令$C$为以卡特兰数 1,1,2,5,14为系数的生成方程

令$F(x) = C(x) - 1$, 保证0次项系数为0, 1次系数非0

那么$F(x) = x(F(x)+1)^2$ , 根据卡特兰数本身推导的定义的到的

令$G(x) = \frac{x}{(x+1)^2}$

那么有$G(F(x)) = \frac{F(x)}{(F(x)+1)^2} = x$

那么$c_n$ 就有了

因此

$
\begin{aligned}
[x^n]F(x) &= \frac{1}{n} \lbrack x^{n-1} \rbrack \left(\frac{x}{G(x)}\right)^n \\
&= \frac{1}{n} \lbrack x^{n-1} \rbrack (x + 1)^{2n} \\
&= \frac{1}{n} \binom{2n}{n-1} = \frac{1}{n+1} \binom{2n}{n},
\end{aligned}$

回到原问题

$A(x) = x(1+3A(x)+A(x)^2)^2$

令$G(x) = \frac{x}{(1+3x+x^2)^2}$

感觉到一点点套路了, 就是如果本身A(x) 的等式里是 $A(x) = x W(A)$的形式,直接$G(x) = \frac{x}{W(x)}$ 就行了,因为这样就有 $G(A) = \frac{A}{W(A)} = \frac{xA}{xW(A)} = \frac{xA}{A} = x$

同时$A(G(x)) = x$, 可以验证$x = A(G) = G(1+3A(G)+A(G)^2)^2 = \frac{x}{(1+3x+x^2)^2}(1+3x+x^2)^2$

$\lbrack x^n \rbrack A(x) = \frac{1}{n}\lbrack x^{n-1}\rbrack \left(\frac{x}{G(x)}\right)^{2n}= \frac{1}{n}\lbrack x^{n-1}\rbrack (1+3x+x^2)^{2n}$

现在只要找到$(1+3x+x^2)^{2N} $的$N-1$次的系数,再除以$N$就是要的答案了

$\left((1+3x+x^2)^{2N}\right)’ = 2N(1+3x+x^2)^{2N-1}(1+3x+x^2)’$

$(1+3x+x^2)\left((1+3x+x^2)^{2N}\right)’ = 2N(1+3x+x^2)^{2N}(1+3x+x^2)’$

把这个用生成方程表示$\sum u_k x^k = (1+3x+x^2)^{2N}$

$(1+3x+x^2)\left(\sum u_k x^k\right)’ = 2N(1+3x+x^2)’ (\sum u_k x^k)$

$(1+3x+x^2)(\sum ku_k x^{k-1}) = 2N(3+2x)(\sum u_k x^k)$

考虑两边$x^{k-1}$次项系数

$ku_k + 3(k-1)u_{k-1} + (k-2)u_{k-2} = 2N(3u_{k-1}+2u_{k-2})$

$u_k = \frac{(6N-3k+3)u_{k-1} + (4N - k + 2)u_{k-2}}{k}$

这样可以O(N) 递推, ??? 这样会触发很多乘法逆元的出现吗? 还是用分数做中间过程?

不做递推的直接求

$\begin{aligned}
u_k &= \lbrack x^k \rbrack (1+3x+x^2)^{2N} \\
&= \lbrack x^k \rbrack ( (1+3x)+x^2)^{2N} \\
&=\lbrack x^k \rbrack \sum_{0 \leq j \leq 2N} \binom{2N}{j}x^{2j} (1+3x)^{2N-j} \\
&=\sum_{0 \leq j \leq 2N} \binom{2N}{j}\lbrack x^{k-2j} \rbrack (1+3x)^{2N-j} \\
&=\sum_{0 \leq j \leq 2N} \binom{2N}{j} \binom{2N-j}{k-2j} 3^{k-2j},
\end{aligned}$

至此可以$O(N)$, 算出

再次注意,要算的是$N-1$次项系数再除以$N$

继续优化 P-recursive

TODO, orz

据说能做到$O(\sqrt{N} log N)$

代码

基于最后那个非递推求的

maspy 的只有46ms

https://atcoder.jp/contests/abc222/submissions/33870319

300+ms

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

typedef long long ll;
#define MOD 998244353
#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;)
#define pb push_back

ll read(){ll r;scanf("%lld",&r);return r;} // read

int n;

int qpow(ll v,int pwr){
  ll r = 1;
  while(pwr){
    if(pwr%2) (r*=v)%=MOD;
    (v*=v)%=MOD;
    pwr/=2;
  }
  return r;
}

int fac[20000010] = {1};
int ifac[20000010];

ll binom(ll n,ll m){
  if(n < 0 || m < 0 || m > n) return 0;
  return fac[n]*(ll)ifac[m]%MOD*(ll)ifac[n-m]%MOD;
}

int main(){
  int n = read();
  rep(i,1,2*n+1) fac[i] = fac[i-1]*(ll)i % MOD;
  ifac[2*n] = qpow(fac[2*n],MOD-2);
  per(i,0,2*n) ifac[i] = ifac[i+1] * (ll)(i+1) % MOD;
  ll ans = 0;
  ll p3 = qpow(3,n-1);
  ll inv3_sq = qpow(3*3, MOD-2);
  // sum_{i=[0..2n]}binom(2n,i)binom(2n-i,n-1-2i) 3^{n-1-2i}
  rep(i,0,2*n+1) {
    if(n-1-2*i < 0) break;
    (ans += binom(2*n,i)*binom(2*n-i,n-1-2*i)%MOD*p3) %= MOD;
    (p3 *= inv3_sq)%=MOD;
  }
  printf("%lld\n",ans*qpow(n,MOD-2) % MOD); // 1/n
  return 0;
}

使用了atcoder modint

#include <bits/stdc++.h>
#include <atcoder/modint>
using mint = atcoder::modint998244353;
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;)

ll read(){ll r;scanf("%lld",&r);return r;} // read

int n;

mint fac[20000010] = {1};
mint ifac[20000010];

mint binom(ll n,ll m){
  if(n < 0 || m < 0 || m > n) return 0;
  return fac[n] * ifac[m] * ifac[n-m];
}

int main(){
  int n = read();
  rep(i,1,2*n+1) fac[i] = fac[i-1] * i ;
  ifac[2*n] = fac[2*n].inv();
  per(i,0,2*n) ifac[i] = ifac[i+1] * (i+1);
  mint ans = 0;
  mint p3 = mint(3).pow(n-1);
  mint inv3_sq = mint(3*3).inv();
  // sum_{i=[0..2n]}binom(2n,i)binom(2n-i,n-1-2i) 3^{n-1-2i}
  rep(i,0,2*n+1) {
    if(n-1-2*i < 0) break;
    ans += binom(2*n,i)*binom(2*n-i,n-1-2*i)*p3;
    p3 *= inv3_sq;
  }
  printf("%d\n",(ans/n).val()); // 1/n
  return 0;
}