bzoj 5104 Fib数列

Posted on | Views:

二次剩余 + $BSGS$ .

这个题可以暴力水过去.先找出模 $P=10^9+9​$ 意义下的循环节,发现是 $\frac {P-1} 3​$ .于是 $O(\frac {P-1} 3)​$ 暴力判.

  • 考虑正经一点的,比较优秀的做法.记 $\phi=\frac {\sqrt 5+1} 2​$ ,则数列第 $n​$ 项 $F_n=\frac 1 {\sqrt 5} \cdot (\phi^n-(-\frac 1 \phi)^n)​$ .
  • 注意 $5$ 在模 $P$ 意义下可以开根号,记开出来的根为 $k=383008016$ ,答案为 $x$ ,那么要解的方程化为,

$$
\phi^x-(-\frac 1 \phi)^x=k\cdot n
$$

  • 换元,令 $t=\phi ^x$ ,则方程化为,

$$
t^2-(-1)^x\cdot=(k\cdot n)t
$$

  • 得到了 $t$ 的一元二次方程,对 $x$ 为奇数/偶数分别求解.若 $\Delta$ 不为二次剩余,则无解.
  • 否则,求根公式解得 $t​$ ,再根据 $t=\phi ^x​$ 用 $BSGS​$ 解出最小 $x​$ .
  • 瓶颈在 $BSGS​$ 上,时间复杂度 $O(\sqrt P)​$ .

暴力

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
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
inline int read()
{
	int out=0,fh=1;
	char jp=getchar();
	while ((jp>'9'||jp<'0')&&jp!='-')
		jp=getchar();
	if (jp=='-')
		fh=-1,jp=getchar();
	while (jp>='0'&&jp<='9')
		out=out*10+jp-'0',jp=getchar();
	return out*fh;
}
const int P=1e9+9;
inline int add(int a,int b)
{
	return (a+b>=P)?(a+b-P):(a+b);
}
int solve(int n)
{
	if(n==1)
		return 1;
	int len=333333336;
	int a=0,b=1,c;
	for(int i=2;i<=len;++i)
	{
		c=add(a,b);
		if(c==n)
			return i;
		a=b;
		b=c;
	}
	return -1;
}
int main()
{
	int n=read();
	cout<<solve(n)<<endl;
	return 0;
}

正经做法

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
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
inline int read()
{
	int out=0,fh=1;
	char jp=getchar();
	while ((jp>'9'||jp<'0')&&jp!='-')
		jp=getchar();
	if (jp=='-')
		fh=-1,jp=getchar();
	while (jp>='0'&&jp<='9')
		out=out*10+jp-'0',jp=getchar();
	return out*fh;
}
const int P=1e9+9,inv2=(P+1)>>1;
const int k=383008016,phi=691504013;
const int Len=(P+1)/3;
inline int add(int a,int b)
{
	return ((a+b)%P+P)%P;
}
inline int mul(int a,int b)
{
	return 1LL * a * b % P;
}
int fpow(int a,int b)
{
	int res=1;
	while(b)
	{
		if(b&1)
			res=mul(res,a);
		a=mul(a,a);
		b>>=1;
	}
	return res;
}
int inv(int x)
{
	return fpow(x,P-2);
}
bool EulerJudge(int x)
{
	return fpow(x,(P-1)>>1)!=(P-1);
}
int wi;
struct Complex
{
	int x,y;//x+yw
	Complex(int x=0,int y=0):x(x),y(y) {}
	Complex operator * (const Complex &rhs) const
	{
		int tx=mul(mul(y,rhs.y),wi);
		tx=add(tx,mul(x,rhs.x));
		int ty=add(mul(x,rhs.y),mul(y,rhs.x));
		return Complex(tx,ty);
	}
	friend Complex operator ^ (Complex a,int b)
	{
		Complex res=Complex(1,0);
		while(b)
		{
			if(b&1)
				res=res*a;
			a=a*a;
			b>>=1;
		}
		return res;
	}
};
int Cipolla(int n)
{
	if(!EulerJudge(n))
		return -1;
	srand(time(NULL));
	int a;
	while("RLDAKIOI")
	{
		a=rand()%P;
		wi=add(mul(a,a),-n);
		if(!EulerJudge(wi))
			break;
	}
	Complex res=Complex(a,1);
	res=res^((P+1)>>1);
	res.x=add(res.x,P);
	return min(res.x,P-res.x);
}
map<int,int> mp;
int BSGS(int a,int b)//a^x=b
{
	mp.clear();
	int m=ceil(sqrt(P));
	for(int j=0; j<m; ++j)
		mp[mul(b,fpow(a,j))]=j;
	for(int i=1; i<=m; ++i)
	{
		int j,val=fpow(a,i*m);
		if(mp.find(val)!=mp.end())
		{
			j=mp[val];
			return i*m-j;
		}
	}
	return P;
}
int n;
int solve(int b,int c,int t)
{
	if(t==-1)
		return P;
	t=add(t,-b);
	t=mul(t,inv2);
	int x=BSGS(phi,t);
	if(c==1 && x%2==0)
		return P;
	if(c==-1 && x%2==1)
		return P;
	return x;
}
int calc(int b,int c)
{
	int Delta=add(mul(b,b),mul(-4,c));
	int t=Cipolla(Delta);
	if(t==-1)
		return P;
	return min(solve(b,c,t),solve(b,c,add(P,-t)));
}
int main()
{
	n=read();
	int ans=P;
	int b=mul(P-k,n);
	ans=min(ans,calc(b,1));//odd
	ans=min(ans,calc(b,-1));//even
	if(ans<P)
		printf("%d\n",ans);
	else
		puts("-1");
	return 0;
}