XTU 1440 Lonely prime number(Fibonacci)

时间：2010-03-17 来源：chhaya

http://acm.xtu.edu.cn/OnlineJudge/index.php/problem/read/id/1440

求（Fib[a]p+ Fib[a+1]p+ Fib[a+2]p……….. + Fib[b]p ） mod p

[a,b]是个区间，p是个素数。通过此题，算是把矩阵二分乘法快速求幂学会了，以后不会再想到矩阵就怕了。

还有两个公式： (x+y)^p mod p = (x^p+y^p) mod p
Fib[n+2]-1 = Fib[1]+Fib[2]+...Fib[n]   (前n项Fib和)

（Fib[a]p+ Fib[a+1]p+ Fib[a+2]p……….. + Fib[b]p ） mod p
= (Fib[a]+ Fib[a+1]+ Fib[a+2]……….. + Fib[b]）^p mod p
= ((Fib[b+2]-1) - (Fib[a+1]-1) ) ^ p mod p

而Fibonacci 数列跟矩阵的关系：

| Fib(n+1) Fib(n)   |        |1   1 |
|                           | =    |        |   ^n
| Fib(n)      Fib(n-1)|      |1   0 |

为了快速求幂，用到了二分乘法：来自http://blog.sina.com.cn/s/blog_51cea4040100g8e8.html

先摆代码：

long long qpow(int a, int b)

{

long long c, d;

c = 1; //存a^b

d = a; //存a的倍幂

while (b > 0)

{

if (b & 1) //或 if (b % 2 == 1)

c *= d;

b = b >> 1; //或 b = b / 2

d = d * d;

}

return c;

}

这里，我们举个列子，比如，a^22

（22）10进制 == (10110)2进制

(a^22)

(a^16)*(a^6)

(a^4)*(a^2)

(a^2)*1

(10110)%2=0 c=1,d=a^2;

(1011)%2=1 c=(a^2)*1,d=a^4;

(101)%2=1 c=(a^4)*(a^2)*1,d=a^8;

(10)%2=1 c=(a^4)*(a^2)*1,d=a^16;

(1)%2=1 c=(a^16)*(a^4)*(a^2)*1;

这样子矩阵求幂跟后面和求幂都可以用到二分，终于算是明白了。。

#include <stdio.h>
long long a, b, p;

/*二分求矩阵幂，次函数求出的是Fib(a)*/
long long fib(long long a)
{
    int i, j;
    long long r[2][2] = {1LL, 1LL, 1LL, 0LL}; //存放矩阵的倍幂
    long long c[2][2] = {1LL, 0LL, 0LL, 1LL}; //存放结果
    long long tmp[2][2]; //放临时结果

    while(a > 0)
    {
        if(a & 1)
        {

             /* c[2][2] = r[2][2] * c[2][2] % p */
            tmp[0][0] = r[0][0] * c[0][0] + r[0][1] * c[1][0];
            tmp[0][1] = r[0][0] * c[0][1] + r[0][1] * c[1][1];
            tmp[1][0] = r[1][0] * c[0][0] + r[1][1] * c[1][0];
            tmp[1][1] = r[1][0] * c[0][1] + r[1][1] * c[1][1];
            for(i = 0; i < 2; i++)
                for(j = 0; j < 2; j++)
                    c[i][j] = tmp[i][j] % p;
        }

             /*r[2][2] = r[2][2] * r[2][2] % p*/
            tmp[0][0] = r[0][0] * r[0][0] + r[0][1] * r[1][0];
            tmp[0][1] = r[0][0] * r[0][1] + r[0][1] * r[1][1];
            tmp[1][0] = r[1][0] * r[0][0] + r[1][1] * r[1][0];
            tmp[1][1] = r[1][0] * r[0][1] + r[1][1] * r[1][1];
            for(i = 0; i < 2; i++)
                for(j = 0; j < 2; j++)
                    r[i][j] = tmp[i][j] % p;
            a = a >> 1;
    }
    return c[0][1];
}

/*二分乘法求幂*/
long long solve(long long n)
{
    long m = p;
    long long c = 1LL;
    while(m > 0)
    {
        if(m & 1)
        {
            c = c * n;
            c = c % p;
        }
        n = n * n;
        n = n % p;
        m = m >> 1;
    }
    return c;
}

int main()
{
    long long n1, n2;
    freopen("1440.in", "r", stdin);
    while(scanf("%lld %lld %lld", &a, &b, &p))
    {
        if(a == 0 && b == 0 && p == 0)
            break;
        n1 = fib(a + 1LL) - 1LL;   //前a-1项和
        n2 = fib(b + 2LL) - 1LL;   //前b项和
        printf("%lld\n", solve((n2 - n1 + p) % p));
    }
    return 0;

}