注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

vfleaking的博客

My name is VFlea King

 
 
 

日志

 
 

2432: [Noi2011]兔农  

2013-05-17 14:38:13|  分类: NOI |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |
【题目大意】
F[n] = F[n - 1] + F[n - 2]
if (n >= 3 && F[n] % G == 1)
    F[n]--;

给出N、G、P,求F[N] % P的值。

N <= 10^18, G <= 10^6, P <= 10^9

当年在考场上看到这道题时我还只会玩脚丫子……2432: [Noi2011]兔农 - vfleaking - vfleaking的博客

首先我们考察模G意义下的Fibonacci数列的循环节吧~
这个有个高端的名字叫Pisano period。
超好懂。

OEIS列出了前几项:http://oeis.org/A001175

现在我们来证,循环节P <= 6G
设G = 2^t * 5^s * p1^k1 * p2^k2 * ... * pa^ka
…………t = 0或s = 0的情况就不考虑了。
P <= lcm(3 * 2^(t - 1), 4 * 5^s, 2(p1 + 1) * p1^(k1 - 1), 2(p2 + 1) * p2^(k2 - 1), ..., 2(pa + 1) * pa^(ka - 1))
= lcm(3, 2^(t - 1), 4, 5^s, ((p1 + 1) / 2) * p1^(k1 - 1), ((p2 + 1) / 2) * p2^(k2 - 1), ..., ((pa + 1) / 2) * pa^(ka - 1))
<= 3 * 4 * (2^(t - 1) * 5^s * p1^k1 * p2^k2 * ... * pa^ka)
<= 12 * (G / 2)
= 6G

这个是不可改进的,因为形如2 * 5^s的数能达到这个上界。

2432: [Noi2011]兔农 - vfleaking - vfleaking的博客……好帅气的结论!!!

Fibonacci数列第n项对G取模的值为Fib[n]。

考虑题目中给的F与Fib的区别,就在于那个减一了~~

神马时候会减一呢?
考虑第一次减一,假设出现在k:
则Fib[k] = 1
那么此时:
...... Fib[k - 2] Fib[k - 1] 0 Fib[k - 1] Fib[k - 1] 2Fib[k - 1] 3Fib[k - 1] ....
考虑后面那一串,现在暂时不对G取模鸟~:
Fib[k - 1] Fib[k - 1] 2Fib[k - 1] 3Fib[k - 1] 5Fib[k - 1] 8Fib[k - 1] ...
它们都是Fib[k - 1]的倍数耶……
同时除以Fib[k - 1]就有:
1 1 2 3 5 8 ...
新的Fibonacci数列。

于是我们注意到每次减一之后一定会出现一个0咯……而0前面一个数如果是a,那么这一段的数列就是:
... a 0 Fib[1]a Fib[2]a Fib[3]a Fib[4]a ...
我们把“Fib[1]a Fib[2]a Fib[3]a Fib[4]a ....... Fib[L - 1]a 0”这样的一串称作数列b[a]。

那么数列b的长度是多少呢?就是问L等于多少,就是问使得Fib[L]a = 1 (mod G)的最小的L。
由于Fib循环节最长是6G <= 6 * 10^6,所以我们可以暴力记录下来Fib的一个循环节,用逆元乱搞下就能求出L。

我们还能求出next[a]数组,表示“Fib[L - 1]a”的值。
这样如果b[a]数列出现了,后面一定跟着b[next[a]]数列。
当然如果不存在a的逆元或者没有Fib[L] = a的逆元,就永远不会有下一个b……
此时反正不会减一,直接当Fibonacci数列处理。

用矩阵乘法算模P的答案。

暴力搞搞都行了。

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <climits>
using namespace std;

typedef long long s64;

const s64 S64_MAX = 9223372036854775807ll;

const int MaxNG = 1000000;

inline int rev(const int &a, const int &m)
{
int x1 = 1, x2 = 0, x3 = a;
int y1 = 0, y2 = 1, y3 = m;

while (y3 != 0)
{
int k = x3 / y3;
int t1 = x1 - y1 * k, t2 = x2 - y2 * k, t3 = x3 - y3 * k;
x1 = y1, x2 = y2, x3 = y3;
y1 = t1, y2 = t2, y3 = t3;
}
if (x3 != 1)
return 0;
return x1 >= 0 ? x1 : x1 + m;
}

s64 n;
int nG, mod;

int fib[MaxNG * 6 + 1];
int firstPos[MaxNG];

inline void calcFirstPos()
{
fill(firstPos, firstPos + nG, -1);

fib[0] = 1, fib[1] = 1;

int i = 2;
do
{
int c = fib[i - 1] + fib[i - 2];
if (c >= nG)
c -= nG;

fib[i] = c;
if (firstPos[c] == -1)
firstPos[c] = i;

i++;
}
while (!(fib[i - 1] == 1 && fib[i - 2] == 1));
}

struct matrix
{
int a[3][3];

inline int& operator()(const int &x, const int &y)
{
return a[x][y];
}
inline int operator()(const int &x, const int &y) const
{
return a[x][y];
}

inline matrix &operator*=(const matrix &rhs)
{
static s64 c[3][3];
for (int i = 0; i < 3; i++)
for (int j = 0; j < 3; j++)
{
c[i][j] = 0;
for (int k = 0; k < 3; k++)
c[i][j] += (s64)a[i][k] * rhs.a[k][j] % mod;
c[i][j] %= mod;
}
for (int i = 0; i < 3; i++)
for (int j = 0; j < 3; j++)
a[i][j] = c[i][j];
return *this;
}
};

inline matrix _genMatI()
{
matrix mat;
mat(0, 0) = 1, mat(0, 1) = 0, mat(0, 2) = 0;
mat(1, 0) = 0, mat(1, 1) = 1, mat(1, 2) = 0;
mat(2, 0) = 0, mat(2, 1) = 0, mat(2, 2) = 1;
return mat;
}

matrix matI = _genMatI();

inline matrix matpow(const matrix &b, const s64 &n)
{
matrix res = matI;
matrix t = b;

for (s64 i = n; i > 0; i >>= 1)
{
if (i & 1)
res *= t;
t *= t;
}
return res;
}

int main()
{
cin >> n >> nG >> mod;

calcFirstPos();

static int preRev[MaxNG];
for (int i = 0; i < nG; i++)
preRev[i] = rev(i, nG);

static int next[MaxNG];
static s64 aL[MaxNG];
for (int i = 0; i < nG; i++)
{
if (!preRev[i] || !firstPos[preRev[i]])
next[i] = -1, aL[i] = S64_MAX;
else
{
int k = firstPos[preRev[i]];
next[i] = (s64)fib[k - 1] * i % nG, aL[i] = k + 1;
}
}

matrix matA, matB;
matA(0, 0) = 1, matA(0, 1) = 1, matA(0, 2) = 0;
matA(1, 0) = 1, matA(1, 1) = 0, matA(1, 2) = 0;
matA(2, 0) = 0, matA(2, 1) = 0, matA(2, 2) = 1;

matB(0, 0) = 1, matB(0, 1) = 1, matB(0, 2) = 0;
matB(1, 0) = 1, matB(1, 1) = 0, matB(1, 2) = 0;
matB(2, 0) = mod - 1, matB(2, 1) = 0, matB(2, 2) = 1;

static bool book[MaxNG + 1];
int cirS;
for (cirS = 1; cirS != -1 && !book[cirS]; cirS = next[cirS])
book[cirS] = true;

matrix matR = matI;
for (int i = 1; i != cirS; i = next[i])
{
if (n < aL[i])
{
matR *= matpow(matA, n);
n = 0;
break;
}
else
{
matR *= matpow(matA, aL[i] - 1);
matR *= matB;
n -= aL[i];
if (n == 0)
break;
}
}
if (n > 0)
{
int cur = cirS;
s64 totL = 0;
matrix matC = matI;
do
{
matC *= matpow(matA, aL[cur] - 1);
matC *= matB;
totL += aL[cur];
cur = next[cur];
}
while (cur != cirS);

matR *= matpow(matC, n / totL);
n %= totL;

for (int i = cirS; ; i = next[i])
{
if (n < aL[i])
{
matR *= matpow(matA, n);
break;
}
else
{
matR *= matpow(matA, aL[i] - 1);
matR *= matB;
n -= aL[i];
if (n == 0)
break;
}
}
}

static int vec[3] = {0, 1, 1};
s64 res = 0;
for (int k = 0; k < 3; k++)
res += (s64)vec[k] * matR(k, 0);
res %= mod;
cout << res << endl;

return 0;
}

  评论这张
 
阅读(3778)| 评论(3)
推荐 转载

历史上的今天

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2017