数学学习笔记4

数学学习笔记4

二阶递推数列

若一个递推式 $f_i=af_{i-1}+bf_{i-2}(0<a,b)$,那么这个递推式就是一个二阶递推数列。

我们以斐波那契数列为例来探讨一下某个特殊类型的二阶递推数列的性质。

$$
f_i=f_{i-1}+f_{i-2}
$$

通项公式(可由特征方程推出,此处不详细介绍):
$$f_n=\dfrac{1}{\sqrt{5}}(\dfrac{(1+\sqrt{5})^n}{2^n}-\dfrac{(1-\sqrt{5})^n}{2^n})$$

性质

  1. 卡西尼性质:$f_{i-1}f_{i+1}-f_i^2=(-1)^i$。
  2. 附加性:$f_{n+k}=f_nf_{k-1}+f_{n+1}f_k$。
  3. $\forall k \in \mathbf{N},f_n \mid f_{nk}$。
  4. 若 $f_a \mid f_b$ 则 $a \mid b$。
  5. 最大公因数:$\gcd(f_a,f_b)=\gcd(a,b)$。
  6. 以斐波那契数列相邻两项作为输入会使欧几里德算法达到最坏复杂度。

求出第 $i$ 个值

我们可以用矩阵乘法,即:
$$\begin{bmatrix}f_i & f_{i+1}\end{bmatrix} \times \begin{bmatrix}0 & 1 \\ 1 & 1\end{bmatrix}=\begin{bmatrix}f_{i+1} & f_{i+2}\end{bmatrix}$$

时间复杂度可以看做是 $O(\log n)$。

特别的如果 $5$ 在 $\bmod$ 下有二次剩余的话,我们可以直接通过通项公式解决,时间复杂度是小常数的 $O(\log n)$。

求出值为 $x$ 的最前下标

即求出最小的 $i$ 使得 $f_i=x$。

首先这个情况必须满足 $5$ 有二次剩余,否则会很麻烦,此处假设 $5$ 有二次剩余。

下面令 $x \gets x \times \sqrt{5}$,通过通项公式计算,设 $A=\dfrac{(1+\sqrt{5})^n}{2^n},B=\dfrac{(1-\sqrt{5})^n}{2^n}$,那么有 $AB=(-1)^n$ 并且 $(A-B)^2=x^2$。

于是知二求一即可,分讨 $4$ 种情况,求一下二次剩余和 bsgs 即可。(可能还需要用到原根,详见笔记 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
#include<bits/stdc++.h>
#define ll long long
#define mod 1000000009
#define phi (mod-1)
#define inv 383008016
#pragma GCC optimize("Ofast")
#include<ext/pb_ds/assoc_container.hpp>
#include<ext/pb_ds/hash_policy.hpp>
#define g 13
using namespace std;
using namespace __gnu_pbds;
ll n,a,b,ans=LLONG_MAX,cnt,cnt1,cnt2,cntt;
__gnu_pbds::gp_hash_table<ll,ll> op;
ll qmi(ll a,ll b,ll p){
ll res = 1%p,t=a;
while(b){
if(b&1) res=res*t%p;
t=t*t%p;
b>>=1;
}
return res;
}
ll bsgs(ll a,ll b,ll p){
op.clear();
ll t = ceil(sqrt(p)),now = 1;
for(ll j=0;j<t;j++) op[b*now%p] = j,now = now*a%p;
ll cnt = 1;
for(ll i=0;i<=t;i++){
ll ans;
if(op.find(cnt)!=op.end()) ans=op[cnt];
else ans=-1;
if(ans>=0&&i*t-ans>0) return i*t-ans;
cnt = cnt*now%p;
}
return p+1;
}
int main(){
cin>>n;
n=n*inv%mod;
a = (1+inv)*qmi(2,mod-2,mod)%mod,b = (1-inv+mod)*qmi(2,mod-2,mod)%mod;
cout<<(qmi(a,2,mod)+qmi(b,2,mod))%mod<<endl;
//n is odd
cnt=n*n%mod,cnt=(cnt-4+mod)%mod,cntt=bsgs(g,cnt,mod);
if(cntt%2==0){
cnt1=(qmi(g,cntt/2,mod)+n)*qmi(2,mod-2,mod)%mod,cnt2=(qmi(g,cntt/2+phi/2,mod)+n)*qmi(2,mod-2,mod)%mod;
ans=min(ans,bsgs(a*a%mod,cnt1*a%mod,mod)*2-1);
ans=min(ans,bsgs(a*a%mod,cnt2*a%mod,mod)*2-1);
}
//n is oven
cnt=n*n%mod,cnt=(cnt+4+mod)%mod,cntt=bsgs(g,cnt,mod);
if(cntt%2==0){
cnt1=(qmi(g,cntt/2,mod)+n)*qmi(2,mod-2,mod)%mod,cnt2=(qmi(g,cntt/2+phi/2,mod)+n)*qmi(2,mod-2,mod)%mod;
ans=min(ans,bsgs(a*a%mod,cnt1,mod)*2);
ans=min(ans,bsgs(a*a%mod,cnt2,mod)*2);
}
if(ans<=mod) cout<<ans<<endl;
else cout<<"-1\n";
return 0;
}

模意义下周期性

先说结论,如果 $\sqrt{5}$ 在模 $p$ 下有可以表示的数,那么周期一定整除 $p-1$,否则周期一定整除 $2p+2$。($p$ 是大于 $5$ 的质数)

并且周期是一个环,不是 $\rho$ 形状的。

证明

设 $A=\dfrac{(1+\sqrt{5})^n}{2^n},B=\dfrac{(1-\sqrt{5})^n}{2^n}$,如果 $\sqrt{5} \in Z_p$,那么 $A \in Z_p,B \in Z_p$,根据费马小定理,那么 $A^p \equiv A,B^p \equiv B$,立即得 $d(p) \mid p-1$。

如果 $\sqrt{5} \not\in Z_p$,那么 $A,B \in Z_p[\sqrt{5}]$,此处涉及到“扩域”的知识,大概就是设 $z=ax+by$ 其中 $y=\sqrt{5}$,$z \in Z_p[\sqrt{5}]$,我们可以对于它定义加法乘法运算,并且满足交换律,结合律,分配律,故我们可以在这个域里面执行四则运算。

但是 $AB$ 仍然等于 $-1$,所以:

$$
\begin{aligned}
A^p &= (-BA)A^p\\
&=-A^{p+1}B\\
&=-B((\dfrac12)^2-(\dfrac{\sqrt{5}}{2})^2)\\
&=B
\end{aligned}
$$

解释一下为什么 $A^{p+1}=((\dfrac12)^2-(\dfrac{\sqrt{5}}{2})^2)$,首先我们知道 $A^{p-1}$ 因为 $5$ 在模 $p$ 意义下没有二次剩余,所以 $A^{p-1}$ 也没有可以表示出来的数,但是因为 $p$ 是奇素数,所以 $p-1$ 是偶数,那么 $\sqrt{5}^{p-1}=5^{\frac{p-1}{2}}$ 是可以表示出来的,它的平方就是 $5^{p-1}$,根据费马小定理可得 $5^{p-1} \equiv 1$,所以 $5^{\frac{p-1}{2}}=1$ 或者 $-1$,如果它等于 $1$,那么 $\sqrt{5}$ 就可以表示出来,所以它等于 $-1$。(需要二项式展开,此处略过)

欧拉准则

这里就可以推出来,如果 $p$ 是奇质数,那么当 $m^{\frac{p-1}{2}} \equiv 1$ 时,$m$ 有二次剩余,当它等于 $-1$ 时,$m$ 没有二次剩余。

所以 $A^{p-1}= -1$,所以 $A^{p+1} = -A^2=AB=((\dfrac12)^2-(\dfrac{\sqrt{5}}{2})^2)=-1$。

同理 $B^p=A$,所以 $A^{2p+3}=A^{2p+2}A=(-1)^2A=A$,于是找到一个循环节 $2p+2$。

证毕。

如果要找最小周期的话,像阶一样求解即可,但是判断的时候用矩阵乘法即可。

这个周期称作皮萨诺周期,如果我们要求其它数的周期的话,我们需要先分解,然后直接用 $\operatorname{lcm}$ 合并。

有性质 $d(p^k)=p \times d(p^{k-1})(2 \le k)$ 且 $p$ 是质数,特别的当 $p \le 5$ 的时候需要手动找循环节,$d(2)=3,d(3)=8,d(5)=20$。

那么合数的循环节也可以求了,并且可以证明 $d(m) \le 6m$,且只有在满足 $m=2\times 5^k$ 的形式时才取到等号。

特别的,如果是其它二阶递推我们把特征方程写出来之后也含有根号,如果这个根号能够在模 $p$ 意义下表示出来,也需要计算最小周期,不能认为周期一定是 $p-1$ 或者 $2p+2$。

寻找周期的代码(需要特判 $1 \le p \le 5$):

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
92
93
94
95
96
97
98
99
100
101
102
103
104
#include<bits/stdc++.h>
#define ll long long
#pragma GCC optimize("Ofast")
using namespace std;
ll T,mod,n,i,ans,tempp,d[1000005],tot,a[2][2],maps[2][2],temp[2][2];
inline ll gcd(ll x,ll y){
if(y==0) return x;
return gcd(y,x%y);
}
inline ll qmi(ll a,ll b,ll p){
ll res = 1%p,t=a;
while(b){
if(b&1) res=res*t%p;
t=t*t%p;
b>>=1;
}
return res;
}
inline char nc(){
static char buf[1000000],*p=buf,*q=buf;
return p==q&&(q=(p=buf)+fread(buf,1,1000000,stdin),p==q)?EOF:*p++;
}
inline ll read(){
ll res = 0;
char c = nc();
while(c<'0'||c>'9')c=nc();
while(c<='9'&&c>='0')res=res*10+c-'0',c=nc();
return res;
}
char obuf[1<<21],*p3=obuf;
inline void pc(char c){
p3-obuf<=(1<<20)?(*p3++=c):(fwrite(obuf,p3-obuf,1,stdout),p3=obuf,*p3++=c);
}
inline void write(ll x){
if(x<0) pc('-'),x=-x;
if(x>9) write(x/10);
pc(x%10+'0');
}
inline void times(ll a[2][2],ll b[2][2]){
temp[0][0] = (a[0][0]*b[0][0]+a[0][1]*b[1][0])%mod;
temp[0][1] = (a[0][0]*b[0][1]+a[0][1]*b[1][1])%mod;
temp[1][0] = (a[1][0]*b[0][0]+a[1][1]*b[1][0])%mod;
temp[1][1] = (a[1][0]*b[0][1]+a[1][1]*b[1][1])%mod;
a[0][0] = temp[0][0],a[0][1] = temp[0][1],a[1][0] = temp[1][0],a[1][1] = temp[1][1];
}
inline ll getf(ll p){
a[0][0] = 0,a[0][1] = 1,a[1][0] = a[1][1] = 0;
maps[0][0] = 0,maps[0][1] = 1,maps[1][1] = 1,maps[1][0] = 1;
while(p){
if(p&1) times(a,maps);
times(maps,maps);
p>>=1;
}
if(a[0][0]==0&&a[0][1]==1) return 1;
return 0;
}
inline ll solve(ll p){
ll now = p;
tot = 0;
for(ll i=2;i*i<=p;i++) while(p%i==0) d[++tot] = i,p/=i;
if(p>1) d[++tot] = p;
for(ll i=1;i<=tot;i++) if(getf(now/d[i])) now/=d[i];
return now;
}
inline bool check(ll x,ll y){
if(qmi(x,(y-1)/2,y)==1) return 0;
return 1;
}
inline ll sol(ll x){
mod = x;
if(x==1) return 1;
if(x==2) return 3;
if(x==3) return 8;
if(x==4) return 6;
if(x==5) return 20;
return check(5,x)?solve(2*x+2):solve(x-1);
}
mt19937 rnd(141414);
int main(){
T=read();
while(T--){
ans=1;
n=read();
if(n==1) ans=1;
if(n==2) ans=3;
if(n==3) ans=8;
if(n==4) ans=6;
if(n==5) ans=20;
if(n>5){
for(i=2;i*i<=n;i++){
if(n%i==0){
tempp = sol(i);
n /= i;
while(n%i==0) n/=i,tempp*=i;
ans = ans*tempp/gcd(ans,tempp);
}
}
if(n>1) tempp=sol(n),ans = ans*tempp/gcd(ans,tempp);
}
write(ans),pc('\n');
}
fwrite(obuf,p3-obuf,1,stdout);
return 0;
}