YZOJ P4587 斐波那契数列
时间限制:1234MS 内存限制:43210KB
难度:\(6.5\) (既然是自己搬的题还是正常一点吧w)
-
题目描述
定义模意义下的递推数列 \(\displaystyle f_n=\left\{ {\begin{array}{*{20}{c}} 1&{,n \le 2}\\ {{f_{n – 1}} + {f_{n – 2}}}&{,n > 2} \end{array}} \right.\),其中模数为 \(1000000009\) 。
给定整数 \(c\)(\(0 \leq c < 1000000009\)),求出它最早出现在数列的哪个位置,并输出下标。
若 \(c\) 永远不会出现在此数列的任一位置,则输出 \(-1\) 。
-
输入格式
多组数据。
第一行一个正整数 \(T\) (\(0 < T \leq 100\)) 表示 \(T\) 组数据。
接下来 \(T\) 行每行一个数表示每组数据的 \(c\) 。
-
输出格式
对于每组数据,输出一行一个数表示答案。
-
样例输入
1 2 3 |
2 987 321 |
-
样例输出
1 2 |
16 -1 |
Source: BZOJ 5104
在其他地方偷偷藏了点私货www,然而这里没有(不需要藏)
简约概述:有公式 \(\displaystyle f_n=\frac{1}{\sqrt5}\left(\frac{1+\sqrt5}{2}\right)^n – \frac{1}{\sqrt5}\left(\frac{1-\sqrt5}{2}\right)^n\) ,可以套用 BSGS 以及 Cipolla 解出 \(n\) 。
极致内容:
-
Section A
首先,\(f_n\) 的通项公式可以通过简单的数学构造得到。
因为一阶线性常系数齐次递推的通项公式为等比数列,所以同理设 \(t, q\),使得 \(\displaystyle f_n-tf_{n-1} = q\left(f_{n-1}-tf_{n-2}\right)\) 。
展开整理得 \(\displaystyle f_n = \left(q+t\right)f_{n-1}-qtf_{n-2}\) ,即 \(\displaystyle \left\{ {\begin{array}{*{20}{c}}{q + t = 1}\\{qt = – 1}\end{array}} \right.\) 。
构造二次函数 \(x^2-x-1=0\) ,由韦达定理可知 \(q, t\) 分别为此方程两根,即 \(\displaystyle \frac{1+\sqrt5}{2}\) 或 \(\displaystyle \frac{1-\sqrt5}{2}\) 。
然后套很多个等比数列求和公式,或者直接设 \(\displaystyle f_n = ax_1^n + bx_2^n\) 代入特殊值解方程组,都可以得出 \(\displaystyle f_n=\frac{1}{\sqrt5}\left(\frac{1+\sqrt5}{2}\right)^n – \frac{1}{\sqrt5}\left(\frac{1-\sqrt5}{2}\right)^n\) 。
-
Section B
接下来对这个公式进行一点处理。
设 \(c=f_n\),\(\displaystyle t=\frac{1+\sqrt5}{2}\) ,则可知 \(\displaystyle \frac{1-\sqrt5}{2} = -\frac{1}{t}\) ,所以有: \(\displaystyle \sqrt 5 c = {\left( t \right)^n} – {\left( { – {1 \over t}} \right)^n}\) 。
换元,再设 \(x=t^n\) ,则有: \(\displaystyle \sqrt 5 c = x – {\left( { – 1} \right)^n}\frac{1}{x}\) 。
再整理,可得以 \(x\) 为主元的一元二次方程 \(\displaystyle x^2 – \sqrt5cx – {\left( { – 1} \right)^n} = 0\) 。
此方程 \(\displaystyle \Delta = 5c^2 + 4\left( { – 1} \right)^n\) ,然后套用一元二次方程求根公式可得 \(\displaystyle x=\frac{{\sqrt 5 c \; \pm \sqrt \Delta }}{2}\) 。
公式中其中涉及到 \(\sqrt \Delta\) 的部分,可以采用 Cipolla 求解。
其中 \(\left(-1\right)^n\) 可以根据 \(n\) 的奇偶性分成 \(-1\) 和 \(1\) 两种取值,分类讨论。
又因为 \(t^n = x\) ,所以跑 BSGS 求出 \(n\) 的最小正整数解,并与其他情况求出的 \(n\) 取最小即可得出答案。
-
Section C
会 BSGS 的人可以去 AC 了。
BSGS (Baby-step giant-step) 是求解模意义下方程 \(\displaystyle {a^x} \equiv b\pmod p\) (这里假设 \(p\) 为质数)中 \(x\) 的最小正整数解的一种方法。
首先,最小正整数解 \(x\) 不会超过 \(p-1\) 。
若 \(x \geq p\) ,那么可以设 \(\displaystyle x=i\left(p-1\right)+j\) ,其中 \(i \geq 1, 0 \leq j < p-1\) 。
那么有 \(\displaystyle {a^x} = {a^{i\left(p-1\right) + j}} = {a^{i\left(p-1\right)}}{a^j} \equiv b\pmod p\) 。
因为 \(p\) 为质数,\(a < p\) ,所以由费马小定理得 \(\displaystyle a^{p-1} \equiv 1\pmod p\) 。
所以有 \(\displaystyle {a^x} = {a^{ip}}{a^j} \equiv {a^j} \equiv b\pmod p\) ,即 \(\displaystyle a^x \equiv a^{x \bmod \left(p-1\right)} \pmod p\) 。
所以只需要在 \((0, p-1]\) 内寻找 \(x\) 即可。
暴力枚举不现实,尝试将 \(p\) 分成两部分,即设 \(\displaystyle s=\lceil \sqrt p \rceil\),那么 \(x\) 可以分解为 \(is-j\) ,其中 \(0 < i \leq s, 0 \leq j < s\) 。
那么有 \(\displaystyle {a^x} = {a^{is – j}} \equiv b\pmod p\) ,即 \(\displaystyle {a^{is}} \equiv b{a^j}\pmod p\) 。
枚举 \(j\) ,把等号右边的值存到 hash_table(\(O(1)\)) 或 std::map (\(O(\log s)\))中,然后枚举 \(i\) 配对。
因为 \(i, j\) 最多只到 \(\displaystyle s=\lceil \sqrt p \rceil\) ,所以复杂度为 \(O(\sqrt p)\) 或 \(O(\sqrt p \log p)\)。
-
Section D
会 Cipolla 的人可以去 AC 了。
Cipolla 是一种求解二次剩余的高效算法。
虽然求解二次剩余也可以使用 BSGS,但是使用 Cipolla 时间复杂度会比之优秀(也不知道我能不能成功卡掉)。
对于奇质数 \(p\) 及整数 \(n\),若存在整数 \(x, 0 \leq x < p\) 且满足 \(\displaystyle x^2 \equiv n\pmod p\) ,那么 \(n\) 就称作在模 \(p\) 意义下的二次剩余。
定义数域 \(\displaystyle \mathbb {F}_p\) 包含 \([0, p-1]\) 中的整数,其中的运算均在模意义下进行。若无特殊说明,下列操作均在 \(\displaystyle \mathbb {F}_p\) 下进行。
先来证明几个引理。
\(Lemma.1\):对于方程 \(\displaystyle x^2 \equiv n \pmod p\) ,总共只有 \(\displaystyle \frac{p-1}{2}\) 个 \(n\)(\(0<n<p\)) 使得 \(x\) 有解。
\(Proof\):假设存在两个不同的数 \(u, v\) 使得 \(\displaystyle u^2 \equiv v^2 \pmod p\) ,那么有 \(\displaystyle p\; |\; \left(u^2-v^2\right)\) ,即 \(\displaystyle p\; |\; \left(u+v\right)\left(u-v\right)\) 。
由于 \(\displaystyle p\; \not |\; \left(u-v\right)\),所以 \(\displaystyle p\; |\; \left(u+v\right)\) ,即 \(\displaystyle u+v \equiv 0 \pmod p\) 。
所以 \((0,p-1]\) 中 \(\displaystyle x^2 \equiv \left(p-x\right)^2 \pmod p\),有且只有 \(\displaystyle \frac{p-1}{2}\) 对这样的数,所以有且只有 \(\displaystyle \frac{p-1}{2}\) 种不同的 \(n\),得证。
\(Lemma.2\):定义 勒让德符号 (Legendre symbol) \(\displaystyle \left(\frac{a}{p}\right)=\begin{cases}1,&a\text{在模$p$意义下是二次剩余}\\-1,&a\text{在模$p$意义下是非二次剩余}\\0,&a\equiv0\pmod p\end{cases}\)
然后有 \(\displaystyle \left(\frac{a}{p}\right)\equiv a^{\frac{p-1}{2}}\pmod p\) (即 欧拉判别法 (Euler’s criterion))。
\(Proof\):当 \(a\) 在模 \(p\) 意义下是二次剩余,即 \(\displaystyle \left(\frac{a}{p}\right)\equiv a^{\frac{p-1}{2}} \equiv 1\pmod p\) 时:设 \(\displaystyle x^2 \equiv a\pmod p\) ,所以有 \(\displaystyle x^{p-1} \equiv a^{\frac{p-1}{2}} \equiv 1\pmod p\) 。
根据费马小定理可知这样的 \(x\) 存在,故原命题得证。
当 \(a\) 在模 \(p\) 意义下不是二次剩余,即 \(\displaystyle \left(\frac{a}{p}\right)\equiv a^{\frac{p-1}{2}} \equiv -1\pmod p\) 时:设 \(\displaystyle x^2 \equiv a\pmod p\) ,所以有 \(\displaystyle x^{p-1} \equiv a^{\frac{p-1}{2}} \equiv -1\pmod p\) 。
根据费马小定理可知这样的 \(x\) 不存在,故原命题得证。
当 \(\displaystyle \left(\frac{a}{p}\right)\equiv a^{\frac{p-1}{2}} \equiv 0\pmod p\) 时:显然。
由于这个算法构造巧妙,所以先介绍其实现方法。
1,给定 \(n\) ,若 \(\displaystyle \left(\frac{n}{p}\right) \neq 1\) ,则可以直接给出答案。
2,随机 \(a\) 直到 \(\displaystyle \left(\frac{a^2-n}{p}\right) = -1\) 。由 引理1 可知,期望随机次数为 \(2\)。
3,扩展至 \(\displaystyle \mathbb{F}_{p^2}\),并把 \(\displaystyle \sqrt{a^2-n}\) 作为“第二条数轴”的单位长度(可类比复数)。这样,此域中的数可以表示为 \(\displaystyle a+b\omega\) ,其中 \(\displaystyle \omega^2 = a^2-n\) ,并且此数域依然满足交换律、结合律以及分配律。
Q:为什么贴图? A:因为我看不懂
4,将 \(x\) 也扩展至 \(\displaystyle \mathbb{F}_{p^2}\),\(\displaystyle x \equiv \left(a+\omega\right)^{\frac{p+1}{2}} \pmod p \) ,取 \(x\)“第一条数轴”中的部分即为答案。
要证明或者搞懂这个是在说什么,还需要再证明一些引理。
\(Lemma.3\):\(\displaystyle \omega^p \equiv – \omega \pmod p\)
\(Proof\):\(\displaystyle {\omega ^p} \equiv \omega \cdot {\omega ^{p – 1}} \equiv \omega {\left( {{\omega ^2}} \right)^{\frac{{p – 1}}{2}}} \equiv \omega {\left( {{a^2} – n} \right)^{\frac{{p – 1}}{2}}} \equiv – \omega \pmod p\)
\(Lemma.4\):\(\displaystyle \left(a+b\right)^p \equiv a^p+b^p \pmod p\)
\(Proof\):根据二项式定理展开,有 \(\displaystyle \left(a+b\right)^p \equiv \sum\limits_{i=0}^pC_p^ia^ib^{p-i} \pmod p\) ,由于当 \(0<i<p\) 时 \(\displaystyle C_p^i \equiv 0 \pmod p\),所以得证。
现在已经可以来证明「步骤4」了。
根据叙述,已知 \(\displaystyle x \equiv \left(a+\omega\right)^{\frac{p+1}{2}} \pmod p\),即 \(\displaystyle x^2 \equiv \left(a+\omega\right)^{p+1} \pmod p\) 。
根据 引理4 ,有 \(\displaystyle \left(a+\omega\right)^{p+1} \equiv \left(a+\omega\right)\left(a+\omega\right)^p \equiv \left(a+\omega\right)\left(a^p+\omega^p\right) \pmod p\) 。
根据费马小定理 \(\displaystyle a^p \equiv a \pmod p\) 以及 引理3,可得 \(\displaystyle \left(a+\omega\right)\left(a^p+\omega^p\right) \equiv \left(a+\omega\right)\left(a-\omega\right) \pmod p\) 。
所以有 \(\displaystyle x^2 \equiv \left(a+\omega\right)\left(a-\omega\right) \equiv a^2-\omega^2 \equiv a^2-\left(a^2-n\right) \equiv n \pmod p\) ,可证 \(x\) 为此方程的一个合法解。
等等,刚刚「步骤4」里不是说「取 \(x\)“第一条数轴”中的部分即为答案」吗?上面的证明都是在 \(\displaystyle \mathbb{F}_{p^2}\) 中进行的,怎么能直接取其中的一部分(保证解可以被压缩到 \(\displaystyle \mathbb{F}_{p}\))?
这里就要提到神奇的 拉格朗日定理 (Lagrange’s theorem) ,都知道 \(n\) 次方程最多有 \(n\) 个解,同时这个结论在模意义下仍是成立的。
所以在 \(\displaystyle \mathbb{F}_{p}\) 中,且 \(n\) 为模 \(p\) 意义下的二次剩余时,方程 \(\displaystyle x^2 \equiv n \pmod p\) 即 \(\displaystyle x^2-n \equiv 0 \pmod p\) 最多存在两根 \(\displaystyle x_0, x_1 \in \mathbb{F}_{p}\) 。
由拉格朗日定理可知即使扩展到 \(\displaystyle \mathbb{F}_{p^2}\) ,此方程仍然最多只会有两个根,并且这两根就“相当于是”原来的 \(x_0, x_1\),所以这两根一定可以被压缩到 \(\displaystyle \mathbb{F}_{p}\) 中,即 \(\displaystyle a+b\omega\) 中 \(b=0\) 。
(可类比复数运算/开根)
最后,分析过程可知此算法时间复杂度可近似于 \(O(\log p)\) 。
-
Section ?
本题总时间复杂度为:\(\displaystyle O\left( \sqrt p + \log p \right)\)。
涉及算法:BSGS + Cipolla + Time constant
// std::map 或者多个 log 或者哪里常数比较大都会被卡掉
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 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 |
#include <cstdio> #include <cstdlib> #include <cstring> #include <climits> #include <ctime> #define _min(_a_,_b_) ((_a_)<(_b_)?(_a_):(_b_)) #define MOD 1000000009 #define SQRT5MOD 383008016 // sqrt(5) #define INV2MOD 500000005 // 1/2 (mod) #define TMOD 691504013 // (sqrt(5)+1)/2 struct _pair { int first,second; _pair(int f,int s) :first(f),second(s) {}; }; inline _pair _make_pair(register int f,register int s) { return (_pair){f,s}; } inline int _pow(register int base,register int b) { register int ans=1; while(b) { if(b&1) ans=(long long)ans * base %MOD; base=(long long)base * base %MOD; b>>=1; } return ans; } int _cm_w; // unit of Fp^2 struct _cm { int a,b; // a+bw, w <-> _cm_w _cm(int a=0,int b=0) :a(a),b(b) {} inline _cm operator * (const _cm&o)const { _cm ans; ans.a=((long long)a * o.a %MOD + (long long)b * o.b %MOD * _cm_w %MOD) %MOD; ans.b=((long long)a * o.b %MOD + (long long)b * o.a %MOD) %MOD; return ans; } }; inline _cm _cm_pow(_cm base,register int b) { _cm ans=1; while(b) { if(b&1) ans=ans * base; base=base * base; b>>=1; } return ans; } inline _pair GetRoot(register int n) { if(!n) return _make_pair(0,0); if(_pow(n, (MOD-1)>>1) == MOD-1) return _make_pair(-1,-1); register int a; while(true) { a=rand()%MOD; _cm_w=((long long)a * a %MOD - n + MOD) %MOD; if(_pow(_cm_w, (MOD-1)>>1) == MOD-1) break; } _cm ans=_cm_pow((_cm){a,1},(MOD+1)>>1); return _make_pair(ans.a, (MOD - ans.a) %MOD); } struct _hash { #define _MAX_LEN 31624 #define _HASH_MOD 176503 int hcnt,hhead[_HASH_MOD+1],hnext[_MAX_LEN+1],harr[_MAX_LEN+1],hval[_MAX_LEN+1]; inline void clear() { hcnt=0,memset(hhead,0,sizeof(hhead)); } inline int& operator [](const int&arr) { register int idx=arr%_HASH_MOD; ++hcnt,hnext[hcnt]=hhead[idx],hhead[idx]=hcnt,harr[hcnt]=arr; return hval[hcnt]; } inline int at(register int arr)const { register int idx=arr%_HASH_MOD; for(register int j=hhead[idx];j;j=hnext[j]) if(harr[j]==arr) return hval[j]; return 0; } }hash; inline int BSGS(register int a,register int c) { hash.clear(); register int sp=__builtin_sqrt(MOD)+1; for(register int j=0,pwaj=1;j<sp;j++) hash[(long long)c * pwaj %MOD]=j+1, pwaj=(long long)pwaj * a %MOD; register int ans=0; for(register int i=1,aisp=(long long)_pow(a,sp),pwaisp=aisp;i<=sp; i++,pwaisp=(long long)pwaisp * aisp %MOD) if((ans=hash.at(pwaisp))) return i*sp-(ans-1); return -1; } inline int GetAns(register int c,register int sd,register bool left2) { register int t=TMOD,sc=(long long)c * SQRT5MOD %MOD; register int x=(long long)((sc + sd) %MOD + MOD) %MOD * INV2MOD %MOD; register int n=BSGS(t,x); //printf("%d^(%d)=%d\n",t,n,x); if(!n || n==-1 || n%2!=left2) return INT_MAX; return n; } int main() { //freopen("test.in","r",stdin); srand(time(NULL)); int _T;scanf("%d",&_T); while(_T--) { int c;scanf("%d",&c); register int ans=INT_MAX,tans; register int n=((long long)c * c %MOD * 5 %MOD - 4 + MOD) %MOD; _pair root = GetRoot(n); if(root.first != -1) tans=GetAns(c,root.first,1),ans=_min(ans,tans); if(root.second!=-1 && root.second!=root.first) tans=GetAns(c,root.second,1),ans=_min(ans,tans); n=((long long)c * c %MOD * 5 %MOD + 4) %MOD; root = GetRoot(n); if(root.first != -1) tans=GetAns(c,root.first,0),ans=_min(ans,tans); if(root.second!=-1 && root.second!=root.first) tans=GetAns(c,root.second,0),ans=_min(ans,tans); printf("%d\n",(ans==INT_MAX ? -1 : ans)); } return 0; } |