数论及其应用——同余问题

  写在前面：《在数论及其应用——素数问题》一文中，笔者只是简单的介绍了一下拓展欧几里得算法可以解决同余问题，但没有详细的展开论述。这篇文章，就是针对数论中的同余问题，展开较为详细的探讨。

  我们通过一个简单的题目，来进入关于同余问题的学习。(Problem source : pku 2769)

Reduced ID Numbers

Time Limit: 2000MS		Memory Limit: 65536K
Total Submissions: 9527		Accepted: 3824

Description

T. Chur teaches various groups of students at university U. Every U-student has a unique Student Identification Number (SIN). A SIN s is an integer in the range 0 ≤ s ≤ MaxSIN with MaxSIN = 10⁶-1. T. Chur finds this range of SINs too large for identification within her groups. For each group, she wants to find the smallest positive integer m, such that within the group all SINs reduced modulo m are unique.

Input

On the first line of the input is a single positive integer N, telling the number of test cases (groups) to follow. Each case starts with one line containing the integer G (1 ≤ G ≤ 300): the number of students in the group. The following G lines each contain one SIN. The SINs within a group are distinct, though not necessarily sorted.

Output

For each test case, output one line containing the smallest modulus m, such that all SINs reduced modulo m are distinct.

    题目大意：给出n个数，你需要找出一个最小的m，使得这n个数中对m取模有n个不同的结果(任意两个数对m取模的结果不等)。
  这道题不论从数理分析还是编程实现，都很基础，我们只是通过这个问题来初步了解一下什么是同余问题。
  数理分析：这其实就是最基本的同余问题——a = b (mod m)。但是题目要求任意的两个数字不可同余，那么我们只需要排除掉所有同余情况即可。
   编程实现：我们只需设置两层循环，对m由小到大进行穷举即可。穷举过程中，我们需要开设一个记录余数的布尔数组p，当穷举时算出余数为i，如果p[i] = 1，表明之前出现过此余数，那么就跳出循环寻找新的m。
   参考代码如下。

#include<cstdlib>
#include<iostream>
#include<string.h>
using namespace std;
bool p[100001];
int d[1000001];
int main()
{
    int ncases;
    int num , i , j , k;
    cin >> ncases;
    while(ncases--)
    {
         cin >> num;

          for(i = 0;i < num;i++)
            cin >> d[i];
        bool find;

        for(i = 1;;i++)
        {
            memset(p , 0 , sizeof(p));
            find = 1;
              for(j = 0;j < num;j++)
              {
                   if(p[d[j]%i])
                   {
                         find = 0;
                         break;
                   }
                   p[d[j]%i] = 1;
              }

                 if(find)
                     break;
        }
        cout << i << endl;

    }
    return 0;
}

  我们再来看一道简单的同余问题。(Problem source : hdu 1021)

Problem Description

There are another kind of Fibonacci numbers: F(0) = 7, F(1) = 11, F(n) = F(n-1) + F(n-2) (n>=2).

 

Input

Input consists of a sequence of lines, each containing an integer n. (n < 1,000,000).

 

Output

Print the word "yes" if 3 divide evenly into F(n).
Print the word "no" if not.

  这道题目在数理分析上没有什么难度，困难点就是n的数据的大小。随着n的增大，我们得到的这个数列的值会变大，由于其递推式是前两项的加和，我们可以考虑只记录这个数列每一项余3的结果，由于是加法，结果不会发生改变。
  这样我们容易得到下面的代码。

#include<stdio.h>
using namespace std;
const int maxn = 1000;
int main()
{
     int F[maxn];
     int i , n;
     F[0] = 7%3, F[1] = 11%3;

     for(i = 2;i < maxn;i++)
     {
          F[i] = (F[i-1] + F[i-2])%3;
     }
     while(scanf("%d",&n) != EOF)
     {
          if(F[n] == 0)   printf("yes\n");
          else            printf("no\n");
     }

}

  但是我们可以看到，这里n的最大值可以取到1000000，如果规定这个数为数组长度，会溢出的。所以我们还需要进行优化。此时我们想到在解决素数问题的时候，面对大数据，我们往往避开直接运算而采用寻找分布规律的方法。这里同样也是适用的，我们通过简单的打表会发现，n = 2,6,10,14,18……的时候会输出yes，也就是说，n%4 = 2时，输出yes。
  这样就不会出现空间溢出的情况了，参考代码如下。
 

#include<stdio.h>
using namespace std;

int main()
{
   int n;
     while(scanf("%d",&n) != EOF)
     {
          if(n%4 == 2)   printf("yes\n");
          else            printf("no\n");
     }

}

 

  承接在《数论及其应用应用——整除性问题》一文中对拓展欧几里得算法的学习，我们知道，它能够解同余方程，准确的说，是一元线性同余方程。在这里我们将更加进一步的讨论解同余方程的问题。    

  解同余方程组。  

  我们先就最简单的方程组分析。  

  x = r1(mod a1)   ①  

  x = r2(mod a2)   ②  

  这个方程组有如下等价形式。  

  x = r1 + y1a1  ③  

  x = r2 + y2a2  ④  

  我们可得r1+y1a1 = r2+y2a2。  

  即y2a2 - y1a1 = r2 - r1  ⑤   这里y1和y2是未知的，但却可以通过我们学习过的拓展欧几里得算法可以对方程⑤进行求解，求出y2后带入④，即可求出x。

  但是如果对含有n个式子的方程组，我们该如何解决呢？我们既然会了两两方程求解，那么对于含n个式子的方程组也不成问题。但是对解空间却要求更加严格。因为，我们整体来看，一个同余方程可以确定一个解空间，而解两个同余方程的过程其实就是在找两个解空间的交集，得到一个新的解空间之后再与第三个同余方程进行匹配，所以在刚开始我们解两个同余方程的时候，显而易见，我们必须得到最小解，这将使得新的解空间更大，否则我们会面临方程组有解而我们却没有找到的窘境。   那么现在问题来了，如何找到两个同余方程的最小解x呢？这个问题我们在《数论及其应用——整除性问题》一文当中已经有所探讨，这里就不再累述。

  那么还剩下最后一个问题，我们成功的得到了两个同余方程的最小解x1，接下来该怎么办呢？即如何用这个解构造一个新的同余方程来继续和第三个同余方程进行联立求解。

  我们构造的这个新的方程要有这样的特点  

  1.最小解是x1。  

  2.这个解必须同时满足前两个同余方程。  

  我们给出一个具体的例子。  

  x = r1+y1a1   ①'  

  x = r2+y2a2   ②'  

  x = r3+y3a3   ③'  

  我们对①'②'进行求解后，得到了最小解x1，基于新方程需要满足的条件，我们可以构造出方程——x = x1 + y4(lcm(a1,a2))。   即x = x1 (mod lcm(a1,a2))  ④'   然后将③'④'继续联立，即可求解。

  以上便是求解同余方程组的方法。   那么我么不妨找一个具体的题目实现以下算法。(Problem source :pku 2891 )

Description

 

Elina is reading a book written by Rujia Liu, which introduces a strange way to express non-negative integers. The way is described as following:

Choose k different positive integers a₁, a₂, …, a_k. For some non-negative m, divide it by every a_i (1 ≤ i ≤ k) to find the remainder r_i. If a₁, a₂, …, a_k are properly chosen, m can be determined, then the pairs (a_i, r_i) can be used to express m.

“It is easy to calculate the pairs from m, ” said Elina. “But how can I find m from the pairs?”

Since Elina is new to programming, this problem is too difficult for her. Can you help her?

Input

The input contains multiple test cases. Each test cases consists of some lines.

• Line 1: Contains the integer k.
• Lines 2 ~ k + 1: Each contains a pair of integers a_i, r_i (1 ≤ i ≤ k).

Output

Output the non-negative integer m on a separate line for each test case. If there are multiple possible values, output the smallest one. If there are no possible values, output -1.

  题目大意：典型的求解同余方程的题目，这里不再累述。

  编程实现：基于对上文模型的分析，编程上实现这一算法并不困难。  

  参考代码如下。

#include<cstdlib>
#include<stdio.h>
#include<iostream>
using namespace std;
void Ex_gcd(long long a,long long b,long long &d,long long&x,long long &y)
{
     if(!b)
     {
         x = 1 , y = 0 , d = a;
     }
     else
        {
            Ex_gcd(b,a%b,d,y,x) , y -= x*(a/b);
        }
}
int main()
{
    long long i,n,a1,r1,a2,r2,a,b,c,d,x0,y0;
    while(scanf("%I64d",&n) != EOF)
    {
         bool ifhave = 1;
         scanf("%I64d%I64d",&a1,&r1);
         for(i = 1;i < n;i++)
         {
              scanf("%I64d%I64d",&a2,&r2);
              a = a1  ,b = a2 , c = r2 - r1;
              Ex_gcd(a,b,d,x0,y0);
              if(c%d != 0)
                  ifhave = 0;
              int t = b/d;
              x0 = (x0*(c/d)%t + t)%t;
              r1 = a1*x0 + r1;
              a1 = a1*(a2/d);
         }
         if(!ifhave)
             printf("-1\n");
         else
             printf("%I64d\n",r1);


    }
    return 0;
}

  求逆元问题：
  在线性代数这门课中，读者可能知道逆元这个概念，在这里我们讨论模的逆元，结合其定义，我们将a在模n下的逆元记做x，则有ax ≡ 1(mod n)，容易看到我们下面需要做的就是求解一个同余方程，那么这正是拓展欧几里得算法能够做的事情。
  当且仅当gcd(a,n) = 1，可以实现如下转化，ax ≡ 1(mod n)   <=>   ax + ny = gcd(a,n) ， 而x即为a在模n下的逆元。
  我们结合一个具体问题来讲其代码实现。(Problem source : hdu 1576)
  

  Problem Description

  要求(A/B)%9973，但由于A很大，我们只给出n(n=A%9973)(我们给定的A必能被B整除，且gcd(B,9973) = 1)。

  由同余运算的基本性质，(A/B)%MOD = A%MOD * B'%MOD = n * B' % MOD，其中B'是A在模9973下的逆元。
  简单的参考代码如下。

#include<cstdio>
const int MOD = 9973;

//ax + by = gcd(a,b)，返回gcd(a,b)
long long extend_gcd(long long a , long long b, long long &x , long long &y)
{
     if(a == 0 && b == 0)  return -1;
     if(b == 0){x = 1 , y = 0;return a;}
     long long d = extend_gcd(b,a%b,y,x);
     y -= a/b*x;
     return d;
}

//ax = 1(mod n)，返回x，即为a在n下的逆元
long long mod_reverse(long long a , long long n)
{
    long long x , y;
    long long d = extend_gcd(a,n,x,y);
    if(d == 1)  return (x%n+n)%n;
    else        return -1;
}

int main()
{
    int T;
    int n , B;
    scanf("%d",&T);
    while(T--)
    {
         scanf("%d%d",&n,&B);
         int x = mod_reverse(B,MOD);
         printf("%d\n",n*x % MOD);
    }
}

  

   接下来我们讨论两种快速幂算法：整数快速幂和矩阵快速幂。   

  我们先来看整数快速幂，它基于的模型问题很简单，求解A^B(mod m)的值。  

  如果B过大，显然我们就无法用计算机直接进行计算，在之前的探讨中，我们知道数论中我们应对大数的常用方法就是在计算过程中进行求余运算。  

  举一个最简单的例子。  

  (a * b) % m = ((a % m) * b) % m   ①  

  这个结论的证明也是很明显的，我们知道求余是不影响+ -运算的结果的，那么我们就可以把乘法编程加法。左式 = (a + a + a +……)%m，右式 = (a%m + a%m +……)%m，显然左式等于右式。  

  那么我们如何利用上面的方法利用到幂的运算中呢？很简单，因为幂可以转化成乘法，然后再将乘法转化成上文中的加法就可以了。   我们将B进行二进制展开，即B = b(n) *2^n + b(n-1)*2^(n-1)+……+b1*2^1 + b0*2^0。  

  则A^B = A^(b(n) *2^n + b(n-1)*2^(n-1)+……+b1*2^1 + b0*2^0)              

            = A^(b(n)*2^n) * A^(b(n-1)*2^(n-1)) * …… * A^(b1*2^1) * A^(b0*2^0)                                    ②  

  这里的序列{bn}是数字B二进制表示形式各个位数上的数字。这里我们就完成了将幂运算转化成乘法运算的过程。   那么这里我们可以看到，bi = 0时，A^(bi*2^i) = 1，不需要进行计算。bi = 1时，这需要计算。   那么我们按照二进制数B的位数，从右往左开始计算，如果该位数上为1，我们进入计算。基于A^(2^i) = (A^(2^(i-1)))^2递推关系式，所以我们在运算到A^(bi * 2^i) % m 的时候，可以利用(A ^ (b(i-1) * 2^(i-1)))^2 % m，这样就从右往左形成了递推，从右端开始，计算结果又可以作为下一位计算的条件。  

  然后再基于一开始我们给出的等式①，恰好服务于等式②乘积的形式，这样我们就可以在计算的过程中进行取余操作了。

  我们来通过一个具体的问题来完成整数快速幂的编程实现。(Problem source : pku 1995)  

  

Description

People are different. Some secretly read magazines full of interesting girls' pictures, others create an A-bomb in their cellar, others like using Windows, and some like difficult mathematical games. Latest marketing research shows, that this market segment was so far underestimated and that there is lack of such games. This kind of game was thus included into the KOKODáKH. The rules follow: 
Each player chooses two numbers Ai and Bi and writes them on a slip of paper. Others cannot see the numbers. In a given moment all players show their numbers to the others. The goal is to determine the sum of all expressions Ai^Bi from all players including oneself and determine the remainder after division by a given number M. The winner is the one who first determines the correct result. According to the players' experience it is possible to increase the difficulty by choosing higher numbers. 
You should write a program that calculates the result and is able to find out who won the game. 

Input

The input consists of Z assignments. The number of them is given by the single positive integer Z appearing on the first line of input. Then the assignements follow. Each assignement begins with line containing an integer M (1 <= M <= 45000). The sum will be divided by this number. Next line contains number of players H (1 <= H <= 45000). Next exactly H lines follow. On each line, there are exactly two numbers Ai and Bi separated by space. Both numbers cannot be equal zero at the same time.

Output

For each assingnement there is the only one line of output. On this line, there is a number, the result of expression 

(A₁^B₁+A₂^B₂+ ... +A_H^B_H)mod M.

  

  题目大意：给你n组(Ai，Bi)，让你计算(∑Ai^Bi )%m的结果。  

  很典型的整数快速幂，只需将多组Ai，Bi结果加和然后再对m取余即可。  

  参考代码如下。

#include<iostream>
#include<stdio.h>
using namespace std;
long long quick_mod(long long a,long long b,long long m)
{
    long long ans = 1;
    while(b)
    {
         if(b&1)
         {
              ans = (ans * a) % m;
              b--;
         }
         b /= 2;
         a = a * a % m;
    }
    return ans;
}

int main()
{
     int ncases , i;
     long long n , m , ans , a , b;
     cin>>ncases;
     while(ncases--)
     {
          cin >> m >> n;
          ans = 0;
          for(i = 1;i <= n;i++)
          {
               cin >> a >> b;
               ans = (ans + quick_mod(a,b,m))%m;
          }
          cout << ans << endl;
     }
}

   

  我们再来学习一道关于乘法快速幂的问题。(Problem source : hdu 1452)

  

Problem Description

Consider a positive integer X,and let S be the sum of all positive integer divisors of 2004^X. Your job is to determine S modulo 29 (the rest of the division of S by 29).
Take X = 1 for an example. The positive integer divisors of 2004^1 are 1, 2, 3, 4, 6, 12, 167, 334, 501, 668, 1002 and 2004. Therefore S = 4704 and S modulo 29 is equal to 6.

 

Input

The input consists of several test cases. Each test case contains a line with the integer X (1 <= X <= 10000000).
A test case of X = 0 indicates the end of input, and should not be processed.

 

Output

For each test case, in a separate line, please output the result of S modulo 29.

  题目大意：给出整数X，计算2004^X的因子和模29的结果。

  数理分析：这道题目既然提到因子和了，那么我们这里引出因子和函数的一些性质。如果百度一下“积性函数”，我们会发现，在常见的积性函数中式有因子和函数σ(n)，并且当且仅当gcd(a,b) = 1的时候，有积性函数的性质σ(a*b) = σ(a)*σ(b)。

  关于这个函数该性质的证明，笔者将会在以后的拓展阅读中给出。

  而当p是素数的时候，有性质2：σ(p^n) = 1+p+p^2+……+p^n = (p^(n+1)-1)/(p-1)。这条性质很显然，无需证明。

  那么我们再次回到这个问题到中来。本题需要我们计算σ(2004^X)，我们开始利用上面的性质进行化简运算。

  σ(2004^X) = σ((4^X)  * (3^X) * (167^X)) 且gcd(3,4,167) = 1,167与22关于29同余，再基于上文给出的性质2，进行化简。

  所以原式  = σ(2^2X) * σ(3^X) * σ(22^X)          

               = (2^(2X+1)-1) *   (3^(X+1)-1)/2  *  (22^(X+1)-1)/21   ①   

  最终我们即求①式mod 29的结果。

  我们再结合取余方程的性质——(a*b)%m = (a%m) * (b%m)，将①式化作三部分来看。

   a =  (2^(2X+1)  - 1) mod 29

   b =  (3^(X+1)-1)/2  mod 29

    c =   (22^(X+1)-1)/21 mod 29   观察到计算b、c的时候有除法，这种形式是无法利用快速幂进行连续求余的，因此我们只需通过拓展欧几里得算法求出2和21的逆元Inv(2)和Inv(21)即可，然后再通过乘法快速幂得到结果。

   Ps：逆元的概念，就是针对一个元素，它的逆元素取消了该元素的某个运算操作。在这里，a*b = 1(mod m)，a、b即互为逆元素。这里Inc(2) = 15 。

   即原式 x / 2 = m  (mod 29)  ②

   变为了 x * 15 = m(mod 29) ③

   可以看到③ = ② * 2 * 15 ，而2 * 15 mod 29 = 1，②③两式本质没有发生改变，但是我们却将含除法的式子转化成了乘法式子，这样便有利于我们快速幂取模运算了。

   参考代码如下。

#include<iostream>
#include<stdio.h>
using namespace std;
long long quick_mod(long long a,long long b,long long m)
{
    long long ans = 1;
    while(b)
    {
         if(b&1)
         {
              ans = (ans * a) % m;
              b--;
         }
         b /= 2;
         a = a * a % m;
    }
    return ans;
}

int main()
{

     long long a , b , c;
     long long x;
     while(cin>>x)
     {
          if(x == 0)  break;
        a = (quick_mod(2,2*x+1,29)-1)%29;
        b = (quick_mod(3,x+1,29)-1)*15%29;
        c = (quick_mod(22,x+1,29)-1)*18%29;
        int ans = ((a * b % 29) *c) % 29;
        printf("%d\n",ans);

     }
}

  基于上面对所谓因子和函数和快速幂的结合，我们再来看一道类似的问题。(Problem source : hdu 1852)
 

Problem Description

As we all know, the next Olympic Games will be held in Beijing in 2008. So the year 2008 seems a little special somehow. You are looking forward to it, too, aren't you? Unfortunately there still are months to go. Take it easy. Luckily you meet me. I have a problem for you to solve. Enjoy your time.
Now given a positive integer N, get the sum S of all positive integer divisors of 2008^N. Oh no, the result may be much larger than you can think. But it is OK to determine the rest of the division of S by K. The result is kept as M.
Pay attention! M is not the answer we want. If you can get 2008^M, that will be wonderful. If it is larger than K, leave it modulo K to the output. See the example for N = 1，K = 10000: The positive integer divisors of 20081 are 1、2、4、8、251、502、1004、2008，S = 3780, M = 3780, 2008^M % K = 5776.

 

Input

The input consists of several test cases. Each test case contains a line with two integers N and K (1 ≤ N ≤ 10000000, 500 ≤ K ≤ 10000). N = K = 0 ends the input file and should not be processed.

 

Output

For each test case, in a separate line, please output the result.

  题目大意：给出数字n和k，m = σ(2008^n)  mod  k， 计算2008^m % k的值。
  数理分析：这道题目在上面的问题上增加了一些难度，其主要体现在取模的数变成了变量，并且多了一步计算2008^m % k的操作。
  基于上面的讨论，我们还是要先计算m。
  由于2008 = 8 * 251，因此m = σ(2008^n)
                                                    = σ(2^3n) * σ(251^n)
                                                    = (2^(3n+1) - 1) * (251^(n+1) - 1) /250   (mod k)
  在上题题当中，我们在遇到除法而影响快速幂取模的问题的时候，采取了计算对应元素的逆元来将等式做等价转换，但是在这里由于k是变量，因此我们无法保证gcd(251,k) = 1，就是无法保证逆元一定存在。
  因此这里介绍另外的处理方法，对于求余运算，有如下性质：
  x / q = p (mod a)    等价于  (x mod aq) / q = p
  基于这条性质，就可以很好的处理含有除式的快速取幂问题了。而关于这条性质的证明，需要参考相关的数论书籍，笔者将会在以后整理资料的时候给出，在这里先不展开论述。
  有了以上分析，编程实现就很容易了，参考代码如下。
   

#include<iostream>
#include<stdio.h>
using namespace std;
long long quick_mod(long long a,long long b,long long m)
{
    long long ans = 1;
    while(b)
    {
         if(b&1)
         {
              ans = (ans * a) % m;
              b--;
         }
         b /= 2;
         a = a * a % m;
    }
    return ans;
}

int main()
{

     long long a , b , c;
     long long n , k;
     while(cin>>n>>k)
     {
          if(n == 0 && k == 0)  break;
        a = quick_mod(2,3*n+1,250*k)-1;
        b = quick_mod(251,n+1,250*k)- 1;

        int ans = (a * b)%(250*k);
        ans /= 250;
        printf("%d\n",quick_mod(2008,ans,k));

     }
}

  下面我们开始介绍矩阵快速幂。
  我们可以通过一个具体的题目来了解一些相关的概念。(Problem source : pku 3070)
 

Description

In the Fibonacci integer sequence, F₀ = 0, F₁ = 1, and F_n = F_{n − 1} + F_{n − 2} for n ≥ 2. For example, the first ten terms of the Fibonacci sequence are:

0, 1, 1, 2, 3, 5, 8, 13, 21, 34, …

An alternative formula for the Fibonacci sequence is

.

Given an integer n, your goal is to compute the last 4 digits of F_n.

Input

The input test file will contain multiple test cases. Each test case consists of a single line containing n (where 0 ≤ n ≤ 1,000,000,000). The end-of-file is denoted by a single line containing the number −1.

Output

For each test case, print the last four digits of F_n. If the last four digits of F_n are all zeros, print ‘0’; otherwise, omit any leading zeros (i.e., print F_n mod 10000).

Sample Input

0
9
999999999
1000000000
-1

Sample Output

0
34
626
6875

Hint

As a reminder, matrix multiplication is associative, and the product of two 2 × 2 matrices is given by

.

Also, note that raising any 2 × 2 matrix to the 0th power gives the identity matrix:

.

  题目大意：给出第n项斐波那契数列的后四位数字，注意不需要前缀0.
  数理分析：题目直接给出了斐波那契数列用矩阵乘方计算的形式，若想证明该等式，可查询线性代数方面的书籍资料，在这里不再展开论述。同时后面也给出里的矩阵乘方的相关概念。那么在这里，我们可以原封不动的迁移我们曾经学习过的乘法快速幂，只不过这里的只需把两个整数的相乘写成两个矩阵的相乘即可。
  相关的代码会结合后面的题目一并给出。

 

  让我们再来看一道极其类似的题目。(Problem source : UVA 10689)

  Description
  Let's define another number sequence ,given by the follow function:
      f(0) = a
      f(1) = b
      f(n) = f(n - 1) + f(n - 2)  ( n > 2)
  When a = 0 and b = 1,this sequence gives the Fibonacci Sequence.Changing the values of a and b,you can get many different sequences.Given the values of a,b,you have to find last m digits of f(n).
  题目大意：根据斐波那契数列的递推式，键盘输入f(0)(用变量a记录)和f(1)(用变量b记录)的值，那么你就会得到一个类斐波那契数列，再给出n、m的值，让你给出f(n)的后m位数字。
  数理分析：f(0) = a , f(1) = b,那么我们不妨继续写几项。
  f(3) = a+b
  f(4) = a + 2b
  f(5) = 2a + 3b
  f(6) = 3a + 5b
  f(7) = 5a + 8b
  ……
  写到这里，我们不难看出这个类斐波那契数列中其实还是隐藏着斐波那契数列，很容易归纳出f(n) = Fib(n-3) * a + Fib(n-2)*b ,而题目要求的输出这个数字的后m位，其实就是一个取余运算，而根据求余运算的规律，x * y (mod m) = (x mod m)*(y mod m)，我们就成功的将这个问题与我们上文所探讨的利用矩阵快速幂对很大的斐波那契数取余的模型当中去。
  参考代码如下。
  

#include <stdio.h>
#include <string.h>
#include <iostream>
#include <algorithm>
#include <math.h>
using namespace std;
const int mod[5]= {0,10,100,1000,10000};
const int MOD =1e9+7;
#define LL long long
LL X,Y,N,M,i,j;
struct Matrlc
{
    int  mapp[2][2];
} ans,base;
Matrlc unit= {1,0,0,1};
Matrlc mult(Matrlc a,Matrlc b)  //两个矩阵的相乘
{
    Matrlc c;
    for(int i=0; i<2; i++)
        for(int j=0; j<2; j++)
        {
            c.mapp[i][j]=0;
            for(int k=0; k<2; k++)
                c.mapp[i][j]+=(a.mapp[i][k]*b.mapp[k][j])%mod[M];
            c.mapp[i][j]%=mod[M];
        }
    return c;
}
void  pow1(int  n)   //矩阵快速幂
{
    base.mapp[0][0] =base.mapp[0][1]=base.mapp[1][0]=1;
    base.mapp[1][1]=0;
    ans.mapp[0][0] = ans.mapp[1][1] = 1;// ans 初始化矩阵
    ans.mapp[0][1] = ans.mapp[1][0] = 0;
    while(n)
    {
        if(n&1)   ans=mult(ans,base);
        base=mult(base,base);
        n>>=1;
    }

}
int main()
{
    int t;
    scanf("%d",&t);
    while(t--)
    {
        scanf("%lld%lld%lld%lld",&X,&Y,&N,&M);
        if(N==0) return X;
        else if(N==1) return Y;
        else
        {
            pow1(N-1);
            LL  result=(ans.mapp[0][0]*Y+ans.mapp[0][1]*X)%mod[M];
                  printf("%lld\n",result);
        }
    }
    return 0;
}

  基于上文中对2*2矩阵快速幂实现对斐波那契数列的求余，我们会问，如何实现n*n的矩阵的快速幂取模呢？(Problem source : pku 3233)
  

Description

Given a n × n matrix A and a positive integer k, find the sum S = A + A² + A³ + … + A^k.

Input

The input contains exactly one test case. The first line of input contains three positive integers n (n ≤ 30), k (k ≤ 10⁹) and m (m < 10⁴). Then follow n lines each containing n nonnegative integers below 32,768, giving A’s elements in row-major order.

Output

Output the elements of S modulo m in the same way as A is given.

  题目大意：给出n*n的矩阵A，整数k，计算S = ∑A^i   (i∈[1,k])。
  编程实现：我们这里需要做改动是将上文中2*2的矩阵相乘在这里编程n*n的矩阵相乘，然后遍历指数i再求和即可。由于这里只涉及了加法运算和乘法运算，因此在计算过程中连续取模是合法的操作。
  参考代码如下。

#include<cstdlib>
#include<stdio.h>
#include<iostream>
#define maxn 101
using namespace std;
typedef struct
{
    int m[maxn][maxn];
}matrax;
matrax a,per;
int n , m;
void init()
{
    int i , j;
    for(i = 0;i < n;i++)
          for(j = 0;j < n;j++)
    {
         scanf("%d",&a.m[i][j]);
         a.m[i][j] %= m;
         per.m[i][j] = (i == j);
    }
}
matrax add(matrax a, matrax b)
{
    matrax c;
    int i , j;
    for(i = 0; i < n;i++)
    for(j = 0;j < n;j++)
    {
        c.m[i][j] = (a.m[i][j] + b.m[i][j])%m;
    }
    return c;
}
matrax multi(matrax a , matrax b)
{
    matrax c;
    int i , j , k;
    for(i = 0;i < n;i++)
        for(j = 0;j < n;j++)
    {
        c.m[i][j] = 0;
        for(k = 0;k < n;k++)
        {
              c.m[i][j] += a.m[i][k] * b.m[k][j];
        }
    c.m[i][j] %= m;
    }
    return c;
}

matrax power(int k)
{
  matrax c , p , ans = per;
  p = a;
  while(k)
  {
       if(k&1)
       {
            ans = multi(ans,p);
            k--;
       }
       else
       {
           k/=2;
           p = multi(p,p);
       }
  }
  return ans;
}

matrax MatraxSum(int k)
{
     if(k == 1)
         return a;
     matrax temp , b;
     temp = MatraxSum(k/2);
     if(k&1)
     {
          b = power(k/2 + 1);
          temp = add(temp , multi(temp,b));
          temp = add(temp , b);
     }
     else
        {
            b = power(k/2);
           temp = add(temp,multi(temp,b));
        }
        return temp;
}

int main()
{
      int k , i , j;
      while(scanf("%d%d%d",&n,&k,&m) != EOF)
      {


            init();
            matrax ans = MatraxSum(k);
            for(i = 0;i < n;i++)
            {
                for(j = 0;j < n-1;j++)
                      printf("%d ",ans.m[i][j]);
                printf("%d\n",ans.m[i][j]);
            }

      }
      return 0;
}

 

  基于之前对解决一元线性同余方程、线性同余方程、高次同余方程的探讨，下面我们来讨论它们有一定关联的一个定理——中国剩余定理。

  给出中国剩余定理：给出如下的线性同余方程组。

 

Description

Some people believe that there are three cycles in a person's life that start the day he or she is born. These three cycles are the physical, emotional, and intellectual cycles, and they have periods of lengths 23, 28, and 33 days, respectively. There is one peak in each period of a cycle. At the peak of a cycle, a person performs at his or her best in the corresponding field (physical, emotional or mental). For example, if it is the mental curve, thought processes will be sharper and concentration will be easier.  Since the three cycles have different periods, the peaks of the three cycles generally occur at different times. We would like to determine when a triple peak occurs (the peaks of all three cycles occur in the same day) for any person. For each cycle, you will be given the number of days from the beginning of the current year at which one of its peaks (not necessarily the first) occurs. You will also be given a date expressed as the number of days from the beginning of the current year. You task is to determine the number of days from the given date to the next triple peak. The given date is not counted. For example, if the given date is 10 and the next triple peak occurs on day 12, the answer is 2, not 3. If a triple peak occurs on the given date, you should give the number of days to the next occurrence of a triple peak.

Input

You will be given a number of cases. The input for each case consists of one line of four integers p, e, i, and d. The values p, e, and i are the number of days from the beginning of the current year at which the physical, emotional, and intellectual cycles peak, respectively. The value d is the given date and may be smaller than any of p, e, or i. All values are non-negative and at most 365, and you may assume that a triple peak will occur within 21252 days of the given date. The end of input is indicated by a line in which p = e = i = d = -1.

Output

For each test case, print the case number followed by a message indicating the number of days to the next triple peak, in the form:
Case 1: the next triple peak occurs in 1234 days.
Use the plural form ``days'' even if the answer is 1.

 

 

 #include<iostream>
#include<stdio.h>
using namespace std;
int m[4],a[4],M;
void exgcd(int a,int b,int &d,int &x,int &y)
{
    if(b==0)
    {
        x = 1 , y = 0;
        d = a;
        return;
    }
    else
        {
             exgcd(b,a%b,d,x,y);
             int temp = x;
             x = y;
             y = temp - (a/b)*y;
        }
}
int China(int r)
{
     M = 1;
     int i,Mi ,x0,y0,d,ans = 0;
     for(i = 1;i <= r;i++)
          M *= m[i];
     for(i = 1;i <= r;i++)
     {
          Mi = M/m[i];
          exgcd(Mi,m[i],d,x0,y0);
          ans = (ans + Mi*x0*a[i])%M;
     }
     if(ans < 0)
         ans += M;
     return ans;
}
int main()
{
    int ncases = 0,p,e,i,d,temp;
    while(scanf("%d%d%d%d",&p,&e,&i,&d) != EOF)
    {
         if(p == -1 && e == -1 && i == -1 && d == -1)
              break;
         ncases++;
         a[1] = p,a[2] = e , a[3] = i;
         m[1] = 23,m[2] = 28 , m[3] = 33;
         int ans = China(3);
         while(ans <= d)
              ans += M;
         printf("Case %d: the next triple peak occurs in %d days.\n",ncases,ans - d);

    }
    return 0;
}