XTU 1440 Lonely prime number(Fibonacci)

求（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;

}