Codeforces Round 530

發表於 2019-01-08 更新於 2024-10-23 分類於 Codeforces ， Div1 Disqus：文章字數： 572 所需閱讀時間 ≈ 2 分鐘

B 数学题简单反证

题目HERE

求一个n x m <= 300 000的二维字符数组

满足任何2x2的字符都包含A T G C

并且和给你的数组不同的字符尽量少

我的思路

我有发现如果已经有

1
2

AG
TC

那么右侧一列,一定是AT,顺序不定

下面两个一定是AG顺序不定,

又发现如果对于任意2*3

1
2

XYZ
???

如果XYZ出现了3个字母

那么???下方的一定是XYZ且顺序不变

因为如果Y不在中间,那么会变成

1
2
3

XYZ
???
Y?Y

第3行填X 或 Z都是不行的

又发现2x3如果XYZ 3个字母不同,具有性质,

1
2

XYZ
???

下面的??? 是确定的,因为X列和Z列是相等的,所以一定是,其中O表示第4个字母

1
2

XYZ
ZOX

然后我就现在想如何构造去找,因为找到3个不同,那么整列都确定了

然后就自闭了

其实可以把上面的事实展开

XYZ
ZOX
XYZ
ZOX
XYZ
ZOX
XYZ
ZOX

有结论,如果某行中出现连续3个不同字符,那么这3个字符对应的列,每列包含且仅包含两个字符,且任意行包含连续这3个不同字符

交换行列有,如果某列中出现连续3个不同字符,那么这3个字符对应的行,每行包含且仅包含两个字符,且任意列包含连续这3个不同字符

注意到上面两条是事实,但是是互斥的

那么目前有3种情况

所有行列,的每一个行列,仅仅包含两个字符[且两两交替.废话]
存在一列包含3个不同字符,那么所有列都包含这3个不同字符,所有行的每一行仅仅包含两个字符,
和2把行列换一下

综上↓

以下是一句话题解

如果你发现了事实:

要么每一行只有两种字母

或

要么每一列只有两种字母

没了然后暴力枚举就完了这个题评分果然有2300,感觉这种题又是需要多枚举几个比如能枚举到5*5的能看出规律就好了

Möbius inversion formula

發表於 2019-01-06 更新於 2024-10-23 分類於 algorithm Disqus：文章字數： 877 所需閱讀時間 ≈ 3 分鐘

Möbius function

${\displaystyle \mu (n)=\delta _ {\omega (n)}^{\Omega (n)}\lambda (n)}$

${\displaystyle \delta }$ 是 Kronecker delta, λ(n) 是 Liouville 函数, ω(n) 是n的不同的质因数个数，Ω(n) 是质因子个数

定义

μ(n) = 0 如果n的质因子有幂次大于1的

μ(n) = 1 如果n由偶数个不同质数相乘

μ(n) = −1 如果n由奇数个不同质数相乘

性质

${\displaystyle \sum _ {d\mid n}\mu (d)={\begin{cases}1&{\text{if }}n=1,\\0&{\text{if }}n>1.\end{cases}}}$

有的地方写作

$\mu * 1 = \epsilon$ 其中星号表示狄利克雷卷积,正好对应的是所求和的d都是n的因子

2023年新的证明

现在看了一些数论基础，认识到之前的证明方式不够系统，像是没学九九乘法表硬算乘法一样

系统一点学习(https://yexiaorain.github.io/Math/number_theory_2/),应该按照以下步骤:

数论函数定义与常见数论函数(I,u,mu)
Dirichlet卷积与广义Dirichlet卷积(结合率,交换率等)
可乘函数与完全可乘函数定义与性质

这样的化学完以后,其实就有了

$\mu * u = I$

$f * I = f$

$g = f * u$也就是所谓的莫比乌斯变换

$g * \mu = (f * u) * \mu = f * (\mu * u) = f * I = f$ 也就有了莫比乌斯反演

性质的证明

首先要μ(x) != 0

需要x的各个质因子次数恰好为1

假设n的所有质因子有t个,既 $n = p_1^{a_1} * p_2^{a_2}…p_t^{a_t}$

那么所有是n的因子的x 且$\mu(x) \neq 0$ 的x,则为这t个质因子的组合

注意到不包含平方数的 Möbius function也可以写成

$\mu(n) = (-1)^{t}$, 其中n不包含平方数,t为其质因子个数

${\displaystyle \sum _{d\mid n}\mu (d) = \sum _ {k=0}^t {t \choose k}(-1)^{k} = (1-1)^t = 0 }$

证毕

积性函数

Möbius function 是积性函数!! 根据积性函数定义如果gcd(a,b) == 1 有 f(ab)=f(a)*f(b)

$\mu(ab) = \mu(a) \mu(b)$ ,当 a和b互质

wikipedia上,还有写些其它性质

不过和本文的关系不大,就没有 copy paste过来

Möbius inversion formula

如果 f和g都是算数函数,且

$g(n)=\sum_{d,\mid ,n}f(d)\quad\text{对所有整数 }n\ge 1$

g(n)表示它所有因子对应的f的和,所以一旦有题目F(x) = sum 所有f(y),y是x的因子，就可以联想到反演

那么有

${\displaystyle f(n)=\sum _{d,\mid ,n}\mu (d)g\left({\frac {n}{d}}\right)\quad {\text{对所有整数 }}n\geq 1}$

证明

${\displaystyle \sum _{n\mid x}\mu (n)g\left({\frac {x}{n}}\right)}$

${\displaystyle =\sum _{n\mid x}\mu (n)\sum _{m\mid {\frac {x}{n}}}f\left(m\right)}$

${\displaystyle=\sum _{m\mid x}f\left(m\right)\sum _{n\mid \frac{x}{m}}\mu (n)}$ (n能取到所有x的因子,m也能取到,且满足n,m其中一个确定时,另一个取值使得n*m为x的因子)

${\displaystyle=f(x)}$

见上面Möbius function的性质,也就是仅在m=x时右侧的sum 才不为0,且为1

实现

线性筛

const int maxn = 100000000;
int prime[maxn], tot, mu[maxn];
bool check[maxn];
void calmu(){
  mu[1] = 1;
  rep(i,2,maxn){
    if (!check[i]){
      prime[tot++] = i;
      mu[i] = -1;
    }
    rep(j,0,tot){
      if (i * prime[j] >= maxn) break;
      check[i * prime[j]] = true;
      if (i % prime[j] == 0){
        mu[i * prime[j]] = 0;
        break;
      }else
        mu[i * prime[j]] = -mu[i];
    }
  }
}

练习题目

CF Hello 2019 F

CF 548 Div2 D

参考

Möbius inversion formula

Möbius function

OEIS A008683

https://mathlesstraveled.com/2016/11/29/the-mobius-function-proof-part-1/

https://mathlesstraveled.com/2016/09/08/post-without-words-11/

http://2000clicks.com/MathHelp/NumberTh06MobiusFunction.aspx

https://www.youtube.com/watch?v=XKjQcPNWMo0

https://www.youtube.com/watch?v=-blqpqbgu0Q

FFT大数乘法(DFT,IDFT)

發表於 2019-01-02 更新於 2024-10-23 分類於 algorithm Disqus：文章字數： 3k 所需閱讀時間 ≈ 10 分鐘

to
FFT,IFT,DFT,IDFT

FFT 大数乘法

不论FFT还是FFT大数乘法,网上都有”大量”的资料

但很多都看得云里雾里,不太能真正的理解

有几篇快要讲清了,但实际靠的举例来描述其可行性,

本篇希望从个人理解以及偏数学的角度来解释

FFT 大数乘法要解决的问题

小学生都知道,乘法就是列竖式,然后进行一个数字一个数字的乘法再相加

那么对于一个长度a和长度b的两个数字相乘,进行了$O(ab)$次的运算,为了描述简便,假设两个数字的长度都为$n=max(a,b)$,不足的补零,也就是$O(n^2)$的时间复杂度

那么如果题目是长$50000$乘长$50000$的运算,那么用原始的乘法是无法在1s内解出的,

FFT可以做到$O(n \cdot log(n))$时间复杂度

计算逻辑总览

要计算 $a \cdot b$

且存在函数$f:x \to X$,能使

$f:x \to X$ 的时间复杂度 $\le O(n \cdot log (n))$
$f^{-1}:X \to x$ 的时间复杂度 $\le O(n \cdot log (n))$, 注:是逆运算不是$-1$次方
$f^{-1}(g(f(a),f(b))) = a\cdot b$, 其中$g$的时间复杂度 $\le O (n \cdot log(n))$

用文字描述,大自为

如果$a,b$可以通过$f$函数变为$A,B$
且$A,B$通过$g$的运算变为$C$
$C$再通过$f$的逆函数变为$c$
且$c == a \cdot b$的
那么原本$O(n^2)$的乘法,就可以用这样的方法替代,变为$O(n \cdot log(n))$

观察上面的式子,可以发现,如果有以下两个映射也可以实现

$f:x \to X$ 的时间复杂度 $\le O(n \cdot log (n))$
$g(f(a),f(b)) = a\cdot b$, 其中$g$的时间复杂度 $\le O(n \cdot log (n))$

之所以写成3条是因为我们具体用得方法, 和3条对应

所以简单的理解,为如果把$a,b$进行低时间复杂度的变形,再进行某种低时间复杂度的运算能变为$a\cdot b$,那么即可在低时间复杂度内算出$a\cdot b$

先抛开时间复杂度谈方法

$ f * g= \mathcal{F}^{-1} \{ \mathcal{F} \{ f \} \cdot \mathcal{F} \{ g \} \}$

以上就是整个FFT大数乘法的所有数学公式,实际问题是离散傅立叶变换

分别解释符号

$*$ 卷积不是乘号
$\mathcal{F}$傅里叶变换
$\cdot$ 乘
$f,g$ 这里看作$n$元向量

卷积(连续的)

定义

${\displaystyle (f * g)(t){\stackrel {\mathrm {def} }{=}}\ \int_{-\infty }^{\infty }f(\tau )g(t-\tau ),d\tau }$

简单讲,卷积的结果$(f * g)(t)$ 等于 $f(r)\cdot g(t-r)$的积分

这里就和乘法$a \cdot b = c$ 对应上了,在乘法过程中,不考虑进位, 每一位的取值都是所有整数

结果的第$t$位的原始值 = 所有$a$的$i$位乘上 $b$的$t-i$位的和

$c_t = \sum_{i=0}^t a_i\cdot b_{t-i}$

傅里叶变换(连续的)

也就是上面的$\mathcal{F}$

${\displaystyle {\hat {f}}(\xi )=\int_{-\infty }^{\infty }f(x)\ e^{-2\pi ix\xi },dx,}$ (定义)

逆变换

${\displaystyle f(x)=\int_{-\infty }^{\infty }{\hat {f}}(\xi )\ e^{2\pi ix\xi },d\xi ,}$ (由正向定义推导)

证明卷积的等式

要证明 $f * g= \mathcal{F}^{-1} \{ \mathcal{F} \{ f \} \cdot \mathcal{F} \{ g \} \}$

抄写自wiki

$\mathcal{F}\{f * g\}(\nu ) $

$= \int_{\mathbb{R}^n} \{f * g\} (z) , e^{-2 \pi i z\cdot \nu}, dz$ (傅立叶定义)

$= \int_{\mathbb{R}^n} \int_{\mathbb{R}^n} f(x) g(z-x), dx, e^{-2 \pi i z\cdot \nu}, dz.$ (卷积定义)

$= \int_{\mathbb{R}^n} f(x)\left(\int_{\mathbb{R}^n} g(z-x)e^{-2 \pi i z\cdot \nu},dz\right),dx.$(积分顺序)

$= \int_{\mathbb{R}^n} f(x) \left( \int_{\mathbb{R}^n} g(y) e^{-2 \pi i (y+x)\cdot\nu},dy \right) ,dx$,(换元$z=x+y$)

$= \int_{\mathbb{R}^n} f(x)e^{-2\pi i x\cdot \nu} \left( \int_{\mathbb{R}^n} g(y) e^{-2 \pi i y\cdot\nu},dy \right) ,dx$ (提取与$y$无关乘法因式)

$= \int_{\mathbb{R}^n} f(x)e^{-2\pi i x\cdot \nu},dx \int_{\mathbb{R}^n} g(y) e^{-2 \pi i y\cdot\nu},dy.$

$= {\mathcal {F}}\{f\}(\nu ) \cdot {\mathcal {F}}\{f\}(\nu)$

得证

小总结

至此,我们有了

乘法中按位表示和卷积对应
上面要找的f,和傅里叶变换对应,$f^{-1}$和傅里叶逆变换对应
有 $ f * g $ 等式

所以

乘法按位表示 $\to$ 卷积 $\to$ ((函数的傅里叶变换)的积)的傅里叶逆变换

当然这里证明的是连续的傅立叶和卷积,而我们下面用的是离散傅立叶(讲道理还是要单独证明的,这里并没有证明

接下来谈一谈,如何实现,并且能在时间复杂度内

先处理傅里叶变换的部分,首先离散傅里叶变换DFT的定义为 $x \to X$

$X_{k}=\sum _ {n=0}^{N-1}x _ {n}e^{-\frac{2 \pi i}{N}kn}\qquad k=0,\dots ,N-1.$

IDFT

$X_{k}=\frac{1}{N} \sum _ {n=0}^{N-1}x _ {n}e^{\frac{2 \pi i}{N}kn}\qquad k=0,\dots ,N-1.$

注: 有的地方系数不是$1$和$\frac{1}{N}$,是两个 $\frac{1}{\sqrt{N}}$, 本质上没有区别,也不影响后面的推导的核心内容

其中$i,\pi,e$为意义,单位虚数,元周率,自然对数, 这里主要是复平面的单位向量的知识点

$N$取$2^{\lceil log_2n \rceil}$ 也就是长度取到2的幂次, 不足的补零

把这个式子用矩阵表示$X = Wx$,发现其实是一个”系数矩阵$W$”

${\displaystyle W=\left(\omega^{jk}\right)_ {j,k=0,\ldots ,N-1}}$

${\displaystyle W={\begin{bmatrix}
1&1&1&1&\cdots &1\
1&\omega &\omega ^{2}&\omega ^{3}&\cdots &\omega ^{N-1}\
1&\omega ^{2}&\omega ^{4}&\omega ^{6}&\cdots &\omega ^{2(N-1)}\
1&\omega ^{3}&\omega ^{6}&\omega ^{9}&\cdots &\omega ^{3(N-1)}\
\vdots &\vdots &\vdots &\vdots &\ddots &\vdots \
1&\omega ^{N-1}&\omega ^{2(N-1)}&\omega ^{3(N-1)}&\cdots &\omega ^{(N-1)(N-1)}
\end{bmatrix}},{\displaystyle \omega =e^{-\frac{2\pi i}{N}}}}$

目标是求向量$X$,虽然求这个矩阵需要$O(n^2)$的时间,但可以不求矩阵,直接求向量$X$

希望能$O(n\cdot log(n))$求出列向量$X$

行列按照$0$-index作为下标, 那么有

矩阵的 $上一列_ i \cdot \omega^i = 下一列_ i$
在偶数列向量中,列上$i$和$i+\frac{n}{2}$互为相反数, 因为$\omega^{\frac{N}{2}} = -1$ (或见youtube

因此有递推式

$f(x,N,idx)$ = 初始向量为$x$, 大小为$N$(即要乘上$N$阶矩阵), 结果的第idx 个的值

对于 $idx < \frac{N}{2}$

$odd(x)$ 表示 $x$的奇数位置上的值构成的长度$\frac{N}{2}$的向量

$even(x)$ 表示 $x$的偶数位置上的值构成的长度$\frac{N}{2}$的向量

$f(x,N,idx) = f(even(x),\frac{N}{2},idx) + w^k \cdot f(odd(x),\frac{N}{2},idx)$

$f(x,N,idx+\frac{N}{2}) = f(even(x),\frac{N}{2},idx) - w^k \cdot f(odd(x),\frac{N}{2},idx)$

函数是折半的,复杂度显然

至此傅里叶变换可以在$O(n\cdot log(n))$时间复杂度内实现

为了帮助理解上面这个递归式, 这里举个例子, 比如$N=8, idx = 2$ 求$f(x,8,2)$

$X_2 = 1\cdot x_0 + w^2\cdot x_1 + w^4 \cdot x_2 + w^6 \cdot x_3 + w^8 \cdot x_4 + w^{10}\cdot x_5 + w^{12} \cdot x_6 + w^{14} \cdot x_7 $

$X_{2+\frac{8}{2}} = X_{6} = 1\cdot x_0 + w^6\cdot x_1 + w^{12} \cdot x_2 + w^{18} \cdot x_3 + w^{24} \cdot x_4 + w^{30}\cdot x_5 + w^{36} \cdot x_6 + w^{42} \cdot x_7 $

注意到在$N=8$时,$w^4 = -1$, 所以

$X_{2+4} = 1\cdot x_0 - w^2\cdot x_1 + w^4 \cdot x_2 - w^6 \cdot x_3 + w^8 \cdot x_4 - w^{10}\cdot x_5 + w^{12} \cdot x_6 - w^{14} \cdot x_7 $

这不就是

$X_2 = (w^0 \cdot x_0 + w^4 \cdot x_2 + w^8 \cdot x_4 + w^{12} \cdot x_6) + w^2 (w^0 \cdot x_1 + w^4 \cdot x_3 + w^8\cdot x_5 + w^{12} \cdot x_7 ) $

$X_{2+4} = (w^0 \cdot x_0 + w^4 \cdot x_2 + w^8 \cdot x_4 + w^{12} \cdot x_6) - w^2 (w^0 \cdot x_1 + w^4 \cdot x_3 + w^8\cdot x_5 + w^{12} \cdot x_7 ) $

注意到$w_{2N}^{2k} = w_{N}^k$

而其中的 $w^0 \cdot x_0 + w^4 \cdot x_2 + w^8 \cdot x_4 + w^{12} \cdot x_6$ 正是$x$的偶数项在 $N=4$ 时的$X_2$, 即$f(even(x), 4, 2)$

$w^0 \cdot x_1 + w^4 \cdot x_3 + w^8\cdot x_5 + w^{12} \cdot x_7 $ 正是$x$的奇数项在 $N=4$ 时的$X_2$,即$f(odd(x), 4, 2)$

那么$g$也就是其中的按位乘也就是简单的$O(n)$

最后傅里叶逆变换, 可以注意到上面的公式只在系数上多了一个负号,所以同理可以在$O(n \cdot log (n))$时间复杂度内实现

递归+原地

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

const double pi = acos(-1.0);

void FFT(complex<double> * A, int n /* size */, int flag /* 1:FFT,-1:IFFT*/) {
    if (n == 1) return;
    // assert((n & (n-1)) == 0); 一定要是2的幂次
    // e^(i*x) = cos(x)+i*sin(x)
    complex<double> Wm(cos(2 * pi / n), -sin(2 * pi / n) * flag);
    vector<complex<double> > tmp(n);
    int u = 0;
    for (int k = 1; k < n; k += 2, u++) tmp[u] = A[k]; // 奇数列
    for (int k = 0; k < n; k += 2) A[k / 2] = A[k]; // 偶数列
    for (int k = u, i = 0; k < n && i < u; k++, i++) A[k] = tmp[i];
    FFT(A        , n / 2, flag);
    FFT(A + n / 2, n / 2, flag);
    complex<double> W(1, 0); // 运算中是 Wm 的i 次方
    for (int i = 0; i < n / 2 ; i++){
        int j = n / 2 + i;
        auto U = A[i];
        auto T = W * A[j];
        A[i] = U + T;
        A[j] = U - T;
        W *= Wm;
    }
}

int main(){
  // 123*456 = 56088
  const int SIZE = 8; // 一定要是2的幂次
  complex<double> * a = new complex<double>[SIZE]{3,2,1,0,0,0,0,0};
  complex<double> * b = new complex<double>[SIZE]{6,5,4,0,0,0,0,0};
  FFT(a, SIZE, 1);
  FFT(b, SIZE, 1);
  complex<double> * c = new complex<double>[SIZE]{0,0,0,0,0,0,0,0};
  for(int i = 0;i < SIZE;i++) c[i] = a[i] * b[i];
  FFT(c, SIZE, -1);
  for(int i = 0;i < SIZE;i++) c[i] /= SIZE;
  // print
  for(int i = 0;i < SIZE;i++) printf("(%.3lf,%.3lf)", c[i].real(), c[i].imag());
  printf("\n");
  for(int i = 0;i < SIZE-1;i++) { // 进位
    c[i+1] += int(c[i].real() / 10);
    c[i] -= int(c[i].real() / 10) * 10;
  }
  for(int i = 0;i < SIZE;i++) printf("(%d)", int(c[i].real()));
  printf("\n");
  return 0;
}

递推+原地

注意,以下代码其中 for m 的部分应该是从2 开始直到n,仅仅是看上去简便,这里处理的时候把所有对应位置都乘上了2,所以对应cos,sin等等与m相关的都变了

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

#define rep(i,a,n) for (int i=a;i < n;i++)
const double pi = acos(-1.0);

// bit 翻转
int rev(int x, int len) {
  int ans = 0;
  while(len -- ){
    ans <<= 1;
    ans |= x & 1;
    x >>= 1;
  }
  return ans;
}

int getlog2(int n){
  return 31 - __builtin_clz(n);
}

void FFT(vector<complex<double> > &A, int flag /* 1:FFT,-1:IFFT */) {
  int n = A.size();
  if(n == 1) return ;
  assert((n & (n-1)) == 0); // 2 的幂次
  int lgn = getlog2(n);
  // 相当于 上面 多次递归每个位置放到的最终位置, 直接计算位置, 而神奇的是刚好 位置 = 原位置的按照长度lgn的bit翻转
  rep(i, 0, n) {
    int j = rev(i, lgn);
    if (j > i) swap(A[i], A[j]);
  }
  // 逻辑和递归里一样了, 区别是 这里不像递归里下标连续, 需要计算下标, 好在还是顺序的
  rep(pwr, 0, lgn){
    int m = 1 << pwr;
    complex<double> Wm(cos(pi / m), -sin(pi / m) * flag);
    for (int k = 0; k < n; k += (m<<1)) {
      complex<double> W(1, 0);
      rep(j, 0, m) {
        auto U = A[k + j];
        auto T = W * A[k + j + m];
        A[k + j] = U + T;
        A[k + j + m] = U - T;
        W *= Wm;
      }
    }
  }
}

int main(){
  // 123*456 = 56088
  const int SIZE = 8; // 一定要是2的幂次
  auto a = vector<complex<double> >{3,2,1,0,0,0,0,0};
  auto b = vector<complex<double> >{6,5,4,0,0,0,0,0};
  FFT(a, 1);
  FFT(b, 1);
  auto c = vector<complex<double> >(8,0);
  for(int i = 0;i < SIZE;i++) c[i] = a[i] * b[i];
  FFT(c,-1);
  for(int i = 0;i < SIZE;i++) c[i] /= SIZE;
  // print
  for(int i = 0;i < SIZE;i++) printf("(%.3lf,%.3lf)", c[i].real(), c[i].imag());
  printf("\n");
  for(int i = 0;i < SIZE-1;i++) { // 进位
    c[i+1] += int(c[i].real() / 10);
    c[i] -= int(c[i].real() / 10) * 10;
  }
  for(int i = 0;i < SIZE;i++) printf("(%d)", int(c[i].real()));
  printf("\n");
  return 0;
}

实例代码

CF R488 Div1 E ACCEPTED

练习题目

CF EDU57 G

CF R488 Div1 E

参考

YouTube - The Fast Fourier Transform Algorithm

模逆元/乘法逆元(CF 入门算法)

發表於 2018-12-30 更新於 2024-10-23 分類於 algorithm Disqus：文章字數： 854 所需閱讀時間 ≈ 3 分鐘

模逆元/乘法逆元

众所周知,如同线段树/FFT一样,乘法逆元是CF的入门必会算法,然而我只能凭空写线段树

每次遇到需要用乘法逆元,我都是快速搜一个函数代码,复制进来用

之前也反复学懂了三四遍,始终没记住(毕竟用的频率不算高)

以下大自整理总结一下

什么是乘法逆元

已知a和p求b, 其中gcd(a,p)=1

$a^{ -1 } \equiv b{\pmod {n}}$

求法

费马小定理

当p是质数时

$a^{p-2} \mod p $

扩展欧几里得

既解$\gcd=1$的扩展欧几里得方程, 已知$(a,p)$求一个$(b,k)$

$a\cdot b+p\cdot k = 1$, 其中$k$为某未知整数

把正常计算$\gcd$的过程的细节写出来int gcd(int a,int b){return b==0?a:gcd(b,a%b);}

$a \cdot b+p\cdot k = 1$,

初始为$(a,p)$

把$(d_i,d_{i+1})$ 变为 $(d_{i+1},d_{i+2})$

$d_1 = a + p\cdot c_1$ 分解式子1 ($c_1 = - a / p, d_1 = a \mod p$)

$d_2 = p + d_1\cdot c_2$ 分解式子2($c_2 = - p / d_1, d_2 = p \mod d_1$)

$d_3 = d_1 + d_2\cdot c_3$ 分解式子3

直到$d_i == 0$, 也就是$d_{i+1} = d_{i-1} + d_2 \cdot c_{i+1}$

$gcd(a,p) = gcd(a,p) \cdot 1 + 0 cdot 0$

在算法实现上也就是再分解, 然后反过来把式子i带入式子i+1

$d_n = d_{n-2} + d_{n-1} \cdot c_n$

$d_n = d_{n-2} + (d_{n-3}+d_{n-2} \cdot c_{n-1}) \cdot c_n$

然后我们只需要记录的是 $d_{n-2}, d_{n-1}$的系数 $gcd(a,p) = d_{n-2} \cdot k_0 + d_{n-1} \cdot k_1$

$gcd(a,p) = d_{n-2} \cdot k_0 + (d_{n-3} \cdot 1 + d_{n-2} \cdot c_{n-1}) \cdot k_1$ , 这个 $c_i = - a_i/b_i$ 在回溯过程中可以得到

$gcd(a,p) = d_{n-3} \cdot k_1 + d_{n-2} \cdot (k_0 - k_1 \cdot a_i / b_i )$

如此展开到最初的a和p即可

// 返回的 [x,y,z] 满足 a*x+b*y=z
tuple<int,int,int> exgcd(int a,int b){
  if(b==0)
    return {1,0,a};
  auto [x,y,z] = exgcd(b,a%b);
  return  {y,x-(a/b)*y,z};
}

以上代码需要C++17

要注意的是如果原来的数据有负数,零,需要在外部处理

扩展欧几里得的矩阵理解

${\begin{pmatrix} a & p \\ 1 & 0 \\ 0 & 1 \end{pmatrix} } {\begin{pmatrix}x \\ y\end{pmatrix}} \Rightarrow {\begin{pmatrix}{a \cdot x + b \cdot y = gcd(x,y) } \\ x \\ y \end{pmatrix}}$

所以可以建立 2*3矩阵,然后对它做列变换,直到得到gcd(a,p), 这时的列向量上的第二第三个值就是(x,y)

批量求 $< p $ 的模逆元

方法1

妇孺皆知,质数的欧拉筛法,又模逆元有可乘性,

所以如果 a=p1*p2,那么 a的模逆元 = p1的模逆元*p2的模逆元 mod p

只需要 O(p+ exgcd平均时间 * ~n/ln n )

方法2

线性递推

$x = p/a$ (整除)

$y = p \mod a < a$

有 $a \cdot x + y = p$

所以 $y = - a \cdot x (\bmod p)$

两边同时乘逆元 $inv_a \cdot inv_y$

$inv_a \cdot inv_y \cdot y = - a \cdot inv_a \cdot inv_y \cdot x (\bmod p)$

即 $inv_a = - x \cdot inv_y = - (p/a) \cdot inv_{p \bmod a} (\bmod p)$

注意到$ y = p \bmod a < a$,所以在计算$inv_a$时,$inv_y$ 已经计算好了

又注意处理负号,在mod p的意义下

$- (p/a) \cdot inv_{p \bmod a} = (p - (p/a)) \cdot inv_{p \bmod a} (\bmod p)$

综上

1
2
3

inv[1] = 1;
rep(i,2,p)
  inv[i] = ((p-p/i)*inv[p%i])%p;

例题

CF EDU57 F

参考

wikipedia

baidu百科

USACO 5.1.1 凸包模板题

發表於 2018-12-14 更新於 2024-10-23 分類於 USACO Disqus：文章字數： 1.1k 所需閱讀時間 ≈ 4 分鐘

USACO 5.1.1

简述

给你n(0<=n<=10000)个平面上点,求凸包的周长

坐标范围(-1'000'000<=x,y<=1'000'000)

解法

用官网介绍的方法

关于此题

实现和介绍方法的差异

首先USACO介绍的方法选取的一个中间的点,然后对周围的点排序,

关于这个,我这里实现是选取y最小(同理可以选取任何一个斜率上的边界点(因为它必定属于凸包(或位于凸包边上共线(也是属于啦))))

这就会导致所用的排序不太相同[我这样也可以用角度直接比]

因为是某斜率边界点,那么剩余的其它点必定在0至180度内,可以直接用叉乘比较,因此又需要注意,0和180度的比较上为了避免对称点的不定偏序关系,需要选取所选斜率上边界点上的边界点

有点拗口,具体举例,如果我选择y最小的点,假设有3个

A(-4,0),B(0,0),C(4,0),如果我选取B(0,0),作为比较的点的话,那么根据叉乘来排序,A(-4,0)和C(4,0)之间的左右关系是可颠倒的BA×BC=0=BC×BA

所以这种情况,建议选择y最小时,x最小或x最大的.[不过幸运的是USACO数据向来很弱,没有加处理仅仅取边界点也能AC]

而对于USACO介绍这种方法选中心点的时候,注意不应该用叉乘来作为比较函数,因为可能存在石头剪刀布的情况,(上面的0和180度也会发生)

举例有4个点A(-1,0),B(1,1),C(0,-1),O(0,0),当我们采取USACO介绍的方法进行排序时,假设选取O为中心的比较点,那么会有剩下3个点都在某另一个点左侧,也就没有偏序关系,无法排序,这种情况该用,计算角度进行排序

在选比较的点/中心点时,USACO还需要解决,绕了一圈以后,处理最初的点是否属于凸包的问题,而选取某一斜率边界点的话,不会有这个问题,因为该点一定在凸包边界上

实现

选取某斜率下的边界点(简单起见就y最小)

以该点为基础,排序(下面代码用逆时针排序),如果有多点共线,选取距离近的

for排序后的数组

逐个尝试push加入点
如果当前点已经有>1个点,且栈顶的两个点A=p[i-1],B=p[i],与新点Q,满足AB 不在 AQ右侧,(依赖前面的排序和距离从近到远),那么去掉栈顶也就是B,重新尝试push该Q

最后求周长就没啥说的了..

一点点正确性怀疑[都是被cf逼的],因为USACO要求的是精度到小数点后2位,所以我有点想知道,这么多浮点数运算损失了多少23333

代码

(没有处理上面说的0和180可能存在的问题,没有HACK的感觉真轻松,向CF日常200+测试数据低头)

/*
ID: 10157591
PROG: fc
LANG: C++
*/
#include <bits/stdc++.h>
#define rep(i,a,n) for (int i=a;i<n;i++)
#define per(i,a,n) for (int i=n-1;i>=a;i--)

using namespace std;

const string filename = "fc";

void usefile(){
  freopen((filename+".in").c_str(),"r",stdin);
  freopen((filename+".out").c_str(),"w",stdout);
}

int n;

pair<double,double> base = {0,0};
pair<double,double>p[10010];

// 右手定则, v1在v2右侧 为正, v1在v2左侧为负
double cross(const pair<double,double>& v1,const pair<double,double>& v2){
  return v1.first * v2.second - v1.second * v2.first;
}

pair<double,double> operator - (const pair<double,double>& p1,const pair<double,double>& p2){
  return {p1.first-p2.first,p1.second-p2.second};
}

double abs(const pair<double,double> & p){
  return sqrt(p.first*p.first + p.second*p.second);
}

// 逆时针 ,共线从近到远
bool cmp(const pair<double,double>& p1,const pair<double,double>& p2){
  double c = cross(p1-base,p2-base);
  if(c != 0) return c>0;
  return abs(p1-base) < abs(p2-base);
}

pair<double,double> ans[10010];

int main(){
  usefile();
  cin>>n;
  if(n < 2){
    cout<<0<<endl;
    return 0;
  }
  rep(i,0,n){
    scanf("%lf %lf",&p[i].first,&p[i].second);
  }
  base = p[0];
  rep(i,0,n){
    if(base.second > p[i].second){
      base = p[i];
    }
  }
  sort(p,p+n,cmp);
  ans[0] = p[0];
  int cnt = 1;
  rep(i,1,n){
    while(cnt>1 && cross(ans[cnt-1]-ans[cnt-2],p[i]-ans[cnt-1]) <= 0){ // 新加的点应该在最后一条射线左侧
      cnt--;
    }
    ans[cnt++] = p[i];
  }
  double output = 0;
  rep(i,0,cnt){
    output+=abs(ans[(i+1)%cnt]-ans[(i)%cnt]);
  }
  printf("%.2lf\n",output);

  return 0;
}

题目HERE

我的思路

以下是一句话题解

Möbius function

定义

性质

2023年新的证明

性质的证明

积性函数

Möbius inversion formula

实现

练习题目

参考

FFT 大数乘法

FFT 大数乘法要解决的问题

计算逻辑总览

先抛开时间复杂度谈方法

卷积(连续的)

傅里叶变换(连续的)

证明 卷积的等式

小总结

接下来谈一谈,如何实现,并且能在时间复杂度内

递归+原地

递推+原地

实例代码

相关延伸

练习题目

参考

模逆元/乘法逆元

什么是乘法逆元

求法

费马小定理

扩展欧几里得

扩展欧几里得 的矩阵理解

批量求 $< p $ 的 模逆元

方法1

方法2

例题

参考

USACO 5.1.1

简述

关于此题

实现

代码

证明卷积的等式

扩展欧几里得的矩阵理解

批量求 $< p $ 的模逆元