【数论】【矩阵乘法】【NOI2011】兔农

Description农夫栋栋近年收入不景气，正在他发愁如何能多赚点钱时，他听到隔壁的小朋友在讨论兔子繁殖的问题。问题是这样的：第一个月初有一对刚出生的小兔子，经过两个月长大后，这对兔子从第三个月开始，每个月初生一对小兔子。新出生的小兔子生长两个月后又能每个月生出一对小兔子。问第n个月有多少只兔子？聪明的你可能已经发现，第n个月的兔子数正好是第n个Fibonacci(斐波那契)数。栋栋不懂什么是Fibonacci数，但他也发现了规律：第i+2个月的兔子数等于第i个月的兔子数加上第i+1个月的兔子数。前几个月的兔子数依次为：1 1 2 3 5 8 13 21 34 …栋栋发现越到后面兔子数增长的越快，期待养兔子一定能赚大钱，于是栋栋在第一个月初买了一对小兔子开始饲养。每天，栋栋都要给兔子们喂食，兔子们吃食时非常特别，总是每k对兔子围成一圈，最后剩下的不足k对的围成一圈，由于兔子特别害怕孤独，从第三个月开始，如果吃食时围成某一个圈的只有一对兔子，这对兔子就会很快死掉。我们假设死去的总是刚出生的兔子，那么每个月的兔子数仍然是可以计算的。例如，当k=7时，前几个月的兔子数依次为：1 1 2 3 5 7 12 19 31 49 80 …给定n，你能帮助栋栋计算第n个月他有多少对兔子么？由于答案可能非常大，你只需要告诉栋栋第n个月的兔子对数除p的余数即可。Input输入一行，包含三个正整数n, k, p。Output输出一行，包含一个整数，表示栋栋第n个月的兔子对数除p的余数。Sample Input6 7 100Sample Output7HINT1<=N<=10^182<=K<=10^62<=P<=10^9

为了方便看，也给一段Libreoffice Math脚本，请用Libreoffice Math打开：

"设"lbrace a_i rbrace"为斐波拉契数列"1,` 1,` 2,` dotsaxis "且"a_0 `=` 0"，"lbrace b_i rbrace"为目标数列．"newline"将"lbrace b_i rbrace"从0处分段，那么每一段的最开始两个数一定是一样的（设为"x"），"newline"所以某一段（记作"lbrace c_i rbrace"）可以看成"c_i `=` x a_i"，而该段的最后一个数（记作"c_l"）"newline"一定满足" c_l `equiv` 1 (mod K)"即"x a_l `equiv` 1 (mod K)"．也就是说，知道了每一段的"newline"第一项也就能确定该段的最后一个数并进而确定该段的长度．"newline "令" boldA `=` left(matrix{0 # 1 # 0 ##1 # 1 # 0 ##0 # 0 # 1}right)"，" boldB `=` left(matrix{0 # 0 # 0 ##0 # 0 # 0 ##1 # 1 # 0}right)newline"那么"{bold A}^i `=` left(alignl matrix{a_{i `-` 1}# a_i# 0 ##a_i# a_{i `+` 1}# 0 ##0 # 0# 1}right)newline"而在实际的递推数列中，若第"i"项被减掉了"1"，那么对后面第"i `+` k newline"项的贡献为"a_{k `+` 1}"（即后面会依次减去"1,` 2,` 3,` 5,` dotsaxis "）"newline"于是可以通过以下方式构造出矩阵：" newline"对每一段都构造出一个矩阵"{bold M}_t `=` {bold A}^{l_t} `+` bold B "（"l_t"为该段长度）"newline "那么"{bold M}_t `times` {bold M}_{t `+` 1} `=` left(alignl matrix{a_{l_t `-` 1}# a_{l_t}# 0 ##a_{l_t}# a_{l_t `+` 1}# 0 ##1# 1# 1}right) `times`left(alignl matrix{a_{l_{t `+` 1} `-` 1}# a_{l_{t `+` 1}}# 0 ##a_{l_{t `+` 1}} # a_{l_{t `+` 1} `+` 1}# 0 ##1# 1# 1}right)newline alignl ~~~~~~~~~~~~~~~~~~~`=` left(alignl matrix{a_{l_t `+` l_{t `+` 1} `-` 1}# a_{l_t `+` l_{t `+` 1}}# 0 ##a_{l_t `+` l_{t `+` 1}} # a_{l_t `+` l_{t `+` 1} `+` 1}# 0 ##a_{l_{t `+` 1} `+` 1}# a_{l_{t `+` 1} `+` 2}# 1}right) newline "仔细观察，发现这个矩阵的第"2"行第"1"列的数恰好就是"newline"正常斐波拉契数列的某一项，而它下面那个数就正好是"newline"由于该段开头被减一而需要被减去的数﹒对每一段都构造出"newline"这样的一个矩阵，最后将第一个矩阵一直乘到最后一个矩阵"newline"所得到的第"2"行第"1"列的数减去第"3"行第"1"列的数就得到了最终答案！"newline"但这样做的代价仍然是很大的……但是我们不怕——"newline"因为将"b"数列分段后，仍然存在循环节！"newline"于是又可以把有循环节那一块的矩阵乘到一起，用快速幂来加速！"newline"另外，如果在给"b"数列分段的时候，某一处找不到乘法逆元，"newline"则可以直接用"bold A"的幂来直接表示出剩下的所有矩阵"newline"（也即剩下的所有项都满足斐波拉契递推关系）．"

代码：

/************************************\ * @prob: NOI2011 兔农               * * @auth: Wang Junji                * * @stat: Accepted.                 * * @date: July. 1st, 2012           * * @memo: 数论、欧拉定理、矩阵乘法      *\************************************/#include <cstdio>#include <cstdlib>#include <algorithm>#include <cstring>#include <string>typedef long long int64;const int maxK = 1000010;struct Matrix{    int ele[3][3]; static int Mod;    Matrix() {memset(ele, 0, sizeof ele);} /* Matrix */    Matrix(const Matrix& b) {*this = b;} /* Matrix */    Matrix& operator=(const Matrix& b)    {memcpy(ele, b.ele, sizeof ele); return *this;}    /* operator= */    int* const operator[](const int& Ind) {return ele[Ind];}  /* operator[] */    const int* const operator[](const int& Ind) const {return ele[Ind];} /* operator[] */    Matrix& operator*=(const Matrix& b)    {        static Matrix res;        for (int i = 0; i < 3; ++i)        for (int j = 0; j < 3; ++j)        {            int64 tmp = 0;            for (int k = 0; k < 3; ++k)                tmp += (int64)ele[i][k] * b[k][j];            res[i][j] = tmp % Mod;        } /* for */        return *this = res;    } /* operator*= */} matr[maxK], A, B, C;int Fib[maxK << 3], Ind_Fib[maxK], Ind_Cyc[maxK], len[maxK], p, K, Matrix::Mod = 0;int64 n;Matrix pow(const Matrix& a, int64 n){    Matrix res, tmp(a);    res[0][0] = res[1][1] = res[2][2] = 1;    for (; n; n >>= 1, tmp *= tmp) if (n & 1) res *= tmp;    return res;} /* pow */int pow(int a, int64 n, int Mod){    int64 res = 1, tmp = a;    for (; n; n >>= 1, (tmp *= tmp) %= Mod) if (n & 1) (res *= tmp) %= Mod;    return res;} /* pow */int gcd(int n, int m){    for (int r; m; n = m, m = r) r = n % m;    return n;} /* gcd */int main(){    freopen("rabbit.in" , "r", stdin );    freopen("rabbit.out", "w", stdout);    scanf("%lld%d%d", &n, &K, &p);    Matrix::Mod = p;    Fib[1] = Fib[2] = 1;    for (int i = 3; ; ++i)    {        Fib[i] = (Fib[i - 1] + Fib[i - 2]) % K;        if (!Ind_Fib[Fib[i]]) Ind_Fib[Fib[i]] = i;        if (Fib[i] == 1 && Fib[i - 1] == 1) break;    } /* for */    int tmp = K, mul_Inv = 1;    for (int i = 2; i * i <= tmp; ++i)    if (tmp % i == 0)    {        mul_Inv *= i - 1;        while ((tmp /= i) % i == 0) mul_Inv *= i;    } /* if */    if (tmp > 1) mul_Inv *= tmp - 1;    A[0][1] = A[1][0] = A[1][1] = A[2][2] = 1;    B[0][0] = B[1][1] = B[2][2] = 1;    C[0][0] = C[1][1] = C[2][2] = 1;    memset(Ind_Cyc, 0xff, sizeof Ind_Cyc);    int tot = 0, ths = 1, ans = 0; int64 totlen = 0;    for (; Ind_Cyc[ths] == -1; ++tot)    {        if (gcd(ths, K) > 1 || !Ind_Fib[pow(ths, mul_Inv - 1, K)])        {            for (int i = 0; i < tot; ++i) B *= matr[i];            B *= pow(A, n - totlen);            ans = ((B[1][0] - B[2][0]) % p + p) % p;            printf("%d\n", ans);            return 0;        } /* if */        Ind_Cyc[ths] = tot;        int nxt = pow(ths, mul_Inv - 1, K);        len[tot] = Ind_Fib[nxt];        if (n < totlen + len[tot])        {            for (int i = 0; i < tot; ++i) B *= matr[i];            B *= pow(A, n - totlen);            ans = ((B[1][0] - B[2][0]) % p + p) % p;            printf("%d\n", ans);            return 0;        } /* if */        totlen += len[tot];        matr[tot] = pow(A, len[tot]);        ++matr[tot][2][0] %= p;        ++matr[tot][2][1] %= p;        ths = ((int64)ths * Fib[len[tot] - 1]) % K;    } /* for */    int start = Ind_Cyc[ths];    for (int i = 0; i < start; ++i) B *= matr[i], n -= len[i];    totlen = 0;    for (int i = start; i < tot; ++i) C *= matr[i], totlen += len[i];    B *= pow(C, n / totlen);    n %= totlen, totlen = 0;    for (int i = start; i < tot; ++i)    {        if (n < totlen + len[i])        {            B *= pow(A, n - totlen);            ans = ((B[1][0] - B[2][0]) % p + p) % p;            printf("%d\n", ans);            return 0;        } /* if */        B *= matr[i];        totlen += len[i];    } /* for */    return 0;} /* main */

本文章来自CSDN博客，转载请注明出处：http://blog.csdn.net/whjpji/article/details/7709301