Loner


  • 首页

  • 归档

  • 搜索

  • 隐藏背景动画

    显示背景动画

单位根反演

发表于 2018-08-13 | 分类于 学习笔记 |

引子

假设现在我们要求一个式子

不妨构造二项式

我们现在将$x=1$,$x=-1$ 带入,则有

两式相加,则有

问题已经解决,$Ans=2^{2*n-1}$。

那如果我们现在想计算

该怎么办?

简介

单位根反演就是做以上工作的东西。

简而言之,如果能快速求出某个多项式的点值,那我们就可能能求出$x$的特定倍数幂的系数和。

形式化的,我们可以求

知识储备

单位根

根据$FFT$的前置知识,我们可以知道,对于一个$n$,若有$x^n-1=0$,则称$x$是$n$次单位根,记为$\omega _ n$。

单位根在复平面上的分布将单位圆均分为$n$份,在第一象限的最靠近实轴的单位根记为$\omega n^1$,一般也可记为$\omega {n}$。其余的按逆时针排序,分别为$\omega _ n^{2 \cdots n-1}$

所以我们不难发现如下几个性质:

  1. $ \omega n^i +\omega n^{i+n/2} (2|n)$

  2. $\omega n^i=\omega n^{i \% n}$

  3. $\omega {kn}^{ki}=\omega {n}^{i}$

由于性质3的存在,我们可以定义本原单位根。

本原单位根的定义是:若$\omega,n$满足$\omega^n-1=0$,且$n$是最小满足条件的,则称$\omega$是$n$的本原单位根。

我们不难发现,对于单位根$\omega {n}^{i}$,若$\gcd(n,i)==1$,则$\omega n^i$一定是n的本原单位根。

否则,设$d=\gcd(i,n)$,则$\omega {n}^{i}=\omega {n’d}^{i’d}=\omega {n’}^{i’}$,$(\omega {n’}^{i’})^{n’}-1=0$。并不是$n$的本原单位根。

性质:

证明:

若$n \mid k$,则有

否则,这玩意是个等比数列求和,根据高中知识,我们有

得证。

应用

现在我们考虑一开始提出的那个问题,求

这个就相当于在求(设$N=nk$)

的次数为$k$的倍数的项的系数和。

我们不妨设

那我们只需要求

即可。

同时注意到

例题

[LOJ 6485] LJJ 学二项式定理

Link

题意

给定$n,s,a 0,a 1,a 2,a 3$,求

题解

我们不妨不关注$a$数组部分,先构造二项式

现在只要分别求出$f(x)$除$4$余$0 \cdots 3$的系数之和即可,直接套用之前所讲的方法即可。

需要注意,在模$P$意义下,设$g$为$P$的原根,则$\omega _ {n}^{1}=g^{\frac{P-1}{n}}$。

第二种构造方法:

注意,这里使用了$n-i$替换了$i$,所以要处理一下$a$数组的下标。

常系数齐次线性递推学习笔记

发表于 2018-06-28 | 分类于 学习笔记 |

简介

什么是常系数齐次线性递推?

设有数列${a 0,a 1,a 2 \cdots }$满足递推关系$a n=\sum\limits {i=1}^{k}a {n-i}*f _ i$,则称该数列满足k阶齐次线性递推关系。

矩阵乘法

这个很普及的东西还是提一下吧。

假设我们有一个满足k阶齐次线性递推关系的数列${a 0,a 1,a 2 \cdots}​$,它满足的齐次线性递推关系为$a n=\sum\limits {i=1}^{k}a {n-i}*f i​$,现在我们要求$a n​$。

假设我们现在维护着一个列矩阵,它的行数为k。

我们考虑如何让这个矩阵的所有元素都向前递推一格,不难设计出k行k列的转移矩阵。

我们可以得到

现在我们要求$a _ n​$,根据矩阵乘法的结合律,易得到:

所以我们只需要算出$M^n \times A​$,然后取最后一个数就可以了。

可以使用矩阵快速幂,复杂度$\mathcal{O}(k^3log _ 2n)$

特征多项式

特征值,特征向量:

若有常数$\lambda$,向量$\vec v$,满足$\lambda\vec v=A\vec v$ ,则称向量$\vec v$为矩阵$A$的一组特征向量,$\lambda$为矩阵$A$的一组特征值。

秩为k的矩阵有k组线性不相关的特征向量。

特征多项式

对关系式进行变换:$(\lambda I-A)\vec v=0$。

这个等式有解的充要条件是$det(\lambda I-A)=0$,大致可以看做向量被拍扁之类的。

可以得到一个$k$次多项式,这个多项式的值等于$0$时把这个方程称为矩阵$A$的特征方程。这个多项式称为矩阵$A$的特征多项式。

特征多项式记为$f(x)=det(\lambda I-A)$。

其中,$det()$为行列式函数,$I$为单位矩阵。

$k$个解自然就是$k$个特征值。(并不需要解出来,怎么处理后面会说)

跟据算数基本定理,最终的多项式有$k$个解,可以写作$f(x)=\prod i(\lambda i-x)$。

Cayley-Hamilton定理

根据Cayley-Hamilton定理,可知$f(A)=O$,$O$为0矩阵。

下面给出一个简单证明:

$f(A)=\prod\limits {i}(\lambda {i} I - A)$

考虑将$f(A)$得到的矩阵分别乘特征向量,如果最后都得到了$0$矩阵,因为这几个特征向量线性不相关,那么可证明$f(A)$乘以任意向量都会得到$0$矩阵,从而$f(A)$为$0$矩阵。

现在问题转换为证明$f(A)​$乘任意特征向量都会得到$0​$矩阵。

先证明:$(\lambda i I-A) \times (\lambda j I -A) = (\lambda j I -A) \times (\lambda i I -A)$

展开即可,不再赘述。

现在考虑任意一个特征向量 $\vec v _ i$

证毕。

优化递推

设转移矩阵为$M$。

根据矩阵快速幂那套理论,我们只要求出来$M^n$就可以了。

我们考虑$M$的特征多项式$f(x)$,这是一个$k$次多项式。

我们将$M$带入,就会得到一个关于$M$的$k$次多项式。

我们可以将$M^n$分解为$M^n=f(M) \times g(M) +R(M)$。

由于$f(M)=O$,所以$M^n=R(M)$,$R(M)$是一个次数不超过$k-1$的多项式。

也就是说,我们只要求出来$M^n \% f(M)$ 就万事大吉了。

但是我们怎么求呢?

我们考虑快速幂的过程(就是倍增)。

假设我们现在已知 $g(M)=M^{2^i} \% f(M)$,现在要求$h(M)= M^{2^{i+1}} \% f(M)$。

一个直接的想法就是令$H(M)=g(M) \times g(M)$。

但是这样做,$H(M)$ 的次数是$2k-2$的。

那我们考虑原本的递推关系,$a n=\sum\limits {i=1}^{k}a {n-i}*f i$。

不难得到$M^n=\sum\limits {i=1} ^{k} M^{n-i} \times f {i}$。

所以我们可以用这个式子将多余出来的系数都向前压一位。

这样我们就得到了一个$\mathcal{O}(k^2\log _ 2 n)$的做法。

那有没有优化的余地呢?
我们从倍增的过程入手,可以发现$H(M)=g(M) \times g(M)$的过程可以由FFT加速至$\mathcal{O}(klog _ 2k)$。

现在只要解决压系数的过程就行了。

我们只要让$h(M) =H(M) \% f(M)$就行了。

等等,$f(M)$怎么求?

根据定义,$f(x)=det(xI - M)$,得到

我们对其进行第一列展开,得到

$M _ {i,j}$代表$M$的算术余子式。

只要直接上多项式取模就行了。

最后的总复杂度$\mathcal{O}(k \log 2 k \log 2 n)$

我们手动模拟一下多项式取模的过程,其实可以发现我们手动向前压的过程就是在暴力取模。

附上BZOJ4161的$\mathcal{O}(k^2\log _ 2 n)$代码

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
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cctype>
#include<cmath>
#include<iostream>
#include<algorithm>
#include<vector>
#include<set>
#include<map>
#include<queue>
#include<stack>
#include<cassert>
typedef long long ll;
typedef unsigned long long ull;
using namespace std;

const int P=1000000007;
const int MAXN=4010;

int n,k,ans;
int f[MAXN],h[MAXN];

struct Matrix{
int a[MAXN];
Matrix (){memset(a,0,sizeof a);}
int& operator [] (const int &i) {return a[i];}
int operator [] (const int &i) const {return a[i];}
inline Matrix operator * (const Matrix &rhs) const
{
Matrix ret;
for(int i=0;i<k;i++)
for(int j=0;j<k;j++)
(ret[i+j]+=1ll*a[i]*rhs[j]%P)%=P;
for(int i=2*k-2;i>=k;ret[i--]=0) // notice ret[i--]=0
for(int j=1;j<=k;j++)
(ret[i-j]+=1ll*ret[i]*f[j]%P)%=P;
return ret;
}
}res;

Matrix ksm(Matrix a,int b) // 严格来说,不是快速幂
{
Matrix ret;
ret[0]=1;
for(;b;a=a*a,b>>=1) if(b&1) ret=ret*a;
return ret;
}

int main()
{
scanf("%d%d",&n,&k);
for(int i=1;i<=k;i++) scanf("%d",&f[i]),f[i]=f[i]>0?f[i]:f[i]+P;
for(int i=0;i<k;i++) scanf("%d",&h[i]),h[i]=h[i]>0?h[i]:h[i]+P;
if(n<k) printf("%d\n",h[n]);
res[1]=1;ans=0;
res=ksm(res,n);
for(int i=0;i<k;i++) ans=(ans+1ll*res[i]*h[i]%P)%P;
printf("%d\n",ans);
#ifdef LOCAL
system("pause");
#endif
}

杜教筛初步

发表于 2018-06-23 | 分类于 学习笔记 |

之前学杜教筛的时候感觉没有学的很明白,现在重新学习一遍。

杜教筛简介

何为杜教筛?
说直白一点就是在低于线性的时间内求出某个积性函数的前缀和。

假设现在我们有一个积性函数$f(i)$,设$S(i)=\sum\limits_{j=1}^{i} f(i)$。

我们要求$S(n)$,假设我们找到了另一个任意数论函数$g(i)$ 。

那么可以很轻松地得到如下的推导:

然后我们做进一步推导:

也就是说,如果我们能够找到一个积性函数$g$,并且我们能够快速计算$(f*g)$的前缀和,$g$的前缀和,那么我们的问题就可以用我们最后得到的那个式子进行计算,具体来说就是可以用数论分块,然后最后一项进行递归处理。

例题

BZOJ 3944

Link

题意

求$\sum\limits{i=1}^{n} \varphi(i)$和$\sum\limits{i=1}^{n} \mu(i)$。

$n\leqslant 2^{31}-1$

题解

我们首先考虑$\sum\limits_{i=1}^{n} \varphi(i)$。

按照刚才那套理论,我们需要找到一个函数$g$,并且$(g*\varphi)$和$g$的前缀和均能够容易求出。

我们不难想到,$\varphi *I=id$,并且$id$和$1$的函数均可以很容易地求出。

所以不妨令$g$函数就是恒等函数$I$。

那么我们就有:

然后就可以前半部分用公式,后半部分数论分块+递归。

我们再来考虑$\sum\limits_{i=1}^{n} \mu(i)$。

按照刚才那套理论,我们需要找到一个函数$g$,并且$(g*\mu)$和$g$的前缀和均能够容易求出。

我们不难想到$\mu *1=\varepsilon$。

那么我们不妨令$g$为恒等函数$I$。

于是我们带入得到:

直接数论分块+递归即可。

洛谷P3768 简单的数学题

Link

题意

求$\sum\limits{i=1} ^{n}\sum\limits\limits {j=1} ^{n} ij*gcd(i,j)$

$n \leqslant 10^{10}$

题解

我们首先对原式进行反演。

现在的问题就是快速求出$\varphi(d)d^2$的前缀和。

我们设$f(i)=\varphi(i)*i^2$,$S(i)=\sum\limits_{j=1}^{i} f(j)$。

现在我们要找到函数$g$,使得我们能够快速计算$(f*g)$的前缀和,$g$的前缀和。

也就是说,我们需要能够快速求以下这个式子:

所以我们只需要使$g(i)=i^2$,然后继续使用刚才的刚才的方法就可以完成此题。

[51NOD 1594] Gcd and Phi

发表于 2018-05-17 | 分类于 题解 |

题目大意:给定n,求$F(n) = \sum{i=1}^{n} \sum{j=1}^{n} \varphi( \gcd(\varphi(i),\varphi(j)))​$ ,多组数据。

$n \leqslant 10^5 ,T \leqslant 5$

不妨设$p[x]=\sum\limits_{i=1}^{n} [\varphi(i)==x]$

不妨枚举$\varphi$的取值,易作如下变形

不妨设$G(x)=\sum\limits{i=1}^{\lfloor \frac{n}{x} \rfloor } \sum\limits{j=1}^{\lfloor \frac{n}{x} \rfloor} p[ix] \times p[jx]​$

显然,

这个东西可以$\mathcal{O}(n \ log_2 n)$求。

我们看看$F(n)$变成了什么

这个东西也可以$\mathcal{O}(n \log_2 n)$求。

好,做完了。

代码如下

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
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cctype>
#include<cmath>
#include<iostream>
#include<algorithm>
#include<vector>
#include<set>
#include<map>
#include<queue>
#include<stack>
typedef long long ll;
typedef unsigned long long ull;
using namespace std;
const int MAXN=2E6+10;
int phi[MAXN],mu[MAXN],n,pri[MAXN],tot,p[MAXN],T;
ll f[MAXN],ans;
bool vis[MAXN];
void solve()
{
scanf("%d",&n);tot=0;
memset(vis,0,sizeof vis);
memset(f,0,sizeof f);
memset(p,0,sizeof p);
phi[1]=mu[1]=1;
for(int i=2;i<=n;i++)
{
if(!vis[i]) pri[++tot]=i,mu[i]=-1,phi[i]=i-1;
for(int j=1;j<=tot&&pri[j]*i<=n;j++)
{
vis[pri[j]*i]=1;
if(i%pri[j]==0) {mu[i*pri[j]]=0;phi[i*pri[j]]=phi[i]*pri[j];break;}
else mu[i*pri[j]]=-mu[i],phi[i*pri[j]]=phi[i]*(pri[j]-1);
}
}
for(int i=1;i<=n;i++) ++p[phi[i]];
ans=0;
for(int i=1;i<=n;i++)
for(int j=i;j<=n;j+=i)
f[i]+=p[j];
for(int i=1;i<=n;i++) f[i]=f[i]*f[i];
for(int i=1;i<=n;i++)
if(mu[i])
for(int j=1;i*j<=n;j++)
ans+=mu[i]*f[i*j]*phi[j];
printf("%lld\n",ans);
}
int main()
{
scanf("%d",&T);
while(T--) solve();
}

APIO2018 游记

发表于 2018-05-13 | 分类于 游记 |

这篇游记是回石家庄的火车上写的,所以可能有口胡

Day0 2018.5.10

从石家庄坐火车去北京,还是一样的火车站,还是一样的地铁站。

话说前一天竟然将电源落在机房了,不爽

到宾馆的时候大概是11点左右,但是那时候房间都没有腾出来,不爽。

所以闭(包)队先带我们去吃了个饭(话说还真有酱爆鸡丁这种操作)。

吃完,回酒店,竟然还没好,差评。

没办法,只好先去试机。

学校好远。。。

先敲了个后缀自动机,略微调了一会,还好。

又写了写后缀数组和$NTT$,基本上一遍过,感觉$OK$。

话说ztb和czy都调了不短的时间,23333。

临走时顺便白嫖了点水

房间终于好了,不愧是五星级酒店,爽到

Day1 2018.5.11

上午一开始先来了张刺激人心的图片,如下

哇,竟然有秦岳,害怕

首先先讲了讲折纸与信息技术,感觉很懵逼。

接下来是秦岳的图片拼接算法和(游戏)人工智能。

感觉挺有意思的。

话说首师大附中的午饭还不错

下午首先讲的是二分,竟然有HEOI2016原题,感觉听懂了不少。

讲课的时候竟然还能开弹幕,学到了。(性感AP,在线IO)

然后是hht的玄学一般图匹配。

晚上在宾馆摸着。

明天考试好虚啊。

Day2 2018.5.12

今天考试。

考成mengbier。

没了。

Day 3 2018.5.13

今天是北大的人来讲课啊。。。

上午讲了些什么杜教筛,洲阁筛和$Min_25$筛和比较高端的数论知识。

下午去80中,讲了些后缀自动机。(话说你画我猜还真是好玩)

晚上坐火车回石家庄。

2018.6.22 @ 高铁

123…6
Zhang_RQ

Zhang_RQ

In solitude, where we are least alone.

27 日志
8 分类
32 标签
友链
  • KingSann
© 2021 Zhang_RQ
由 Hexo 强力驱动
|
主题 — NexT.Pisces v6.0.1