HDU 1573 X问题（中国剩余定理）

X问题

Time Limit: 1000/1000 MS (Java/Others)    Memory Limit: 32768/32768 K (Java/Others)
Total Submission(s): 7240    Accepted Submission(s): 2561

Problem Description

求在小于等于N的正整数中有多少个X满足：X mod a[0] = b[0], X mod a[1] = b[1], X mod a[2] = b[2], …, X mod a[i] = b[i], … (0 < a[i] <= 10)。

 

Input

输入数据的第一行为一个正整数T，表示有T组测试数据。每组测试数据的第一行为两个正整数N，M (0 < N <= 1000,000,000 , 0 < M <= 10)，表示X小于等于N，数组a和b中各有M个元素。接下来两行，每行各有M个正整数，分别为a和b中的元素。

 

Output

对应每一组输入，在独立一行中输出一个正整数，表示满足条件的X的个数。

 

Sample Input

310 31 2 30 1 2100 73 4 5 6 7 8 91 2 3 4 5 6 710000 101 2 3 4 5 6 7 8 9 100 1 2 3 4 5 6 7 8 9

 

Sample Output

103

 

解析:

同余式定义:

同余式是数论的基本概念之一，设m是给定的一个正整数，a、b是整数，若满足m|(a-b)，则称a与b对模m同余，记为a≡b(mod m)，或记为a≡b(m)。这个式子称为模m的同余式，若m∤ (a-b)，则称a、b对模m不同余，同余概念又常表达为：

1.a=b+km(k∈Z)；

2.a和b被m除时有相同的余数。

 不互质的中国剩余定理：

题意很简单，给出k组a r  每组代表x ≡ r (mod a) ;其中要注意的就是所有的a不一定互素，因为a不互素就不能直接用中国剩余定理来做，查了很多资料，感觉数学家的思维不是凡人可以理解的，还是自己写一下计算的过程

首先来计算两组 x ≡ r1 ( mod a1 ) ; x ≡  r2 ( mod a2 ) ;定义变量 k1 k2 得到 x = k1 * a1 + r1 ; x = k2 * a2 + r2 ;由上面的等式得到 k1 * a1 + r1 = k2 * a2 + r2 ;转化为 k1*a1 = (r2 - r1) + k2 *a2 ;对左右取模a2，因为 （k2*a2）%s2 = 0，所以等式转化为  k1 * a1 ≡ ( r2 - r1 ) (mod a2)；使用扩展欧几里得可以求解到 k1的值（判断是否存在k1的值），将k1带回到x1 = k1 * a1 + r1 ;得到同时满足于{ x = k1 * a1 + r1 ; x = k2 * a2 + r2 ; }的一个特解 ， 所以x ≡ x1 (mod lcm(a1,a2) ) ; 也就是x ≡ ( k1*a1+r1 ) ( mod ( a1*a2/d ) );这样也就将两个同余式转化为了一个，通过不断的转化，将k个等式合并为一个 ，用扩展欧几里得求出最小的正解x

对于x=(x%s+s)%s取到最小值的解释：

同余方程 ax≡b (mod n)对于未知数x 有解，当且仅当gcd(a,n) | b。且方程有解时，方程有gcd(a,n) 个解。

    求解方程ax≡b (mod n) 相当于求解方程ax+ ny= b, (x, y为整数)

    设d= gcd(a,n)，假如整数x 和y，满足d= ax+ ny(用扩展欧几里德得出)。如果d| b，则方程

   a* x0+ n* y0= d， 方程两边乘以b/ d，(因为d|b，所以能够整除)，得到a* x0* b/ d+ n* y0* b/ d= b。
    所以 x= x0* b/ d，y= y0* b/ d为 ax+ ny= b的一个解，所以 x= x0* b/ d为 ax= b (mod n )的解。

ax≡b (mod n)的一个解为x0= x* (b/ d ) mod n，且方程的d 个解分别为xi= (x0+ i* (n/ d ))mod n {i= 0... d-1}。

    设ans=x*(b/d),s=n/d;

方程ax≡b (mod n)的最小整数解为：(ans%s+s)%s;  //例如下面的例子,(11%3+3)%3=2

//if a = b * q + r (q > 0 and 0 <= r < q), then we have  r=(a%q+q)%q，(无论a是正还是负)

 

该式用于任意序列中的数都能求出最小的那个值

#include<cstdio>
#include<iostream>
#include<cstring>
#include<algorithm>
#define MAXN 15
using namespace std;
int flag;
long long int n;

long long int Extend(long long int a,long long int b,long long int &x,long long int &y)   //求p * a+q * b = Gcd(a, b)的一组解p0,q0
{
	long long int d;
	if(b==0)
	{
		x=1;
		y=0;
		return a;     //a是a,b的最大公约数
	}
	d=Extend(b,a%b,y,x);
	y-=(a/b)*x;            //x1=y2; y1=x2-(a/b)*y2;
	return d;
}

long long int China(long long int w[],long long int r[],long long int len,long long int w1,long long int r1)
{
	long long int i,x,y,d,c,w2,r2,s;
	for(i=2;i<=len;i++)
	{
		w2=w[i];
		r2=r[i];
		c=r2-r1;
		d=Extend(w1,w2,x,y);       // ?x * a1 ≡?( r2 - r1 ) (mod a2) 即x*a1 = (r2 - r1) + y *a2
		//求出了x*a1 = gcd(w1,w2) + y *a2的解
		if(c%d)    //d是w1,w2的最大公因数,所以c一定能被d整除
		{
			flag=1;  //k1 * a1 - k2 * a2 = r2-r1     (k1 * a1 - k2 * a2)%d = (r2-r1)%d =0 
			return 0;
		}
		else
		{
			x=x*c/d;       //扩展成x * a1 ≡?( r2 - r1 ) (mod a2)的解
			s=w2/d;        //根据Extend(w1,w2,x,y)得出，间隔
			x=(x%s+s)%s;   //将x转换成最小的整数解
			r1=x*w1+r1;     //将k1带回到 x1 = x * a1 + r1
			w1=w1*w2/d;    //w1,w2的最小公倍数

		}
	}
	d=Extend(1,w1,x,y);   //当所有式子化完后,1是a1的值(式子就变成x * 1 ≡?( r2 - r1 ) (mod w1))
	if(r1%d)
	{
		flag=1;
		return 0;
	}
	else
	{
		x=x*r1/d;
		s=w1/d;
		x=(x%s+s)%s;  //求出满足所有方程的最小值

		if(x>n)
		{
			flag=1;
			return 0;
		}
		else
		{
			if(x==0)return (n-x)/s;
			else  return (n-x)/s+1;
		}

	}
}

int main()
{
	long long int t,m,i,w1,r1,d,ans;
	long long int w[MAXN],r[MAXN];
	cin>>t;
	while(t--)
	{
		scanf("%lld%lld",&n,&m);
		for(i=1;i<=m;i++)
		{
			scanf("%lld",&w[i]);
		}
		for(i=1;i<=m;i++)
		{
			scanf("%lld",&r[i]);
		}
		flag=0;
		w1=w[1];
		r1=r[1];
		ans=China(w,r,m,w1,r1);
		if(flag)
		{
			printf("0\n");
		}
		else
		{
			printf("%lld\n",ans);
		}

	}
	return 0;

}

下面有一个只求最小的满足条件的X，不求个数的代码

#include <bits/stdc++.h>
typedef long long int lli;

const int MAXN = 1e5+10;

lli w[MAXN];
lli r[MAXN];
lli w1,r1;
int flag;

lli Extend(lli a,lli b,lli &x,lli &y)
{
    lli d;
    if(b==0)
    {
        x=1;
        y=0;
        return a;   //a是a,b的最大公约数
    }
    d=Extend(b,a%b,y,x);
    y-=(a/b)*x;     //x1=y2; y1=x2-(a/b)*y2;
    return d;
}

lli   China(int k)
{
    lli w1=w[1];
    lli r1=r[1];
    lli w2,r2,c,d,x,y,s;
    for(int i=2;i<=k;i++)
    {
        w2=w[i];
        r2=r[i];
        c=r2-r1;
        d=Extend(w1,w2,x,y);
        if(c%d)
        {
            flag=1;
            return -1;
        }
        else
        {
            x=x*c/d;
            s=w2/d;
            x=(x%s+s)%s;
            r1=x*w1+r1;
            w1=w1*w2/d;
        }
    }
    d=Extend(1,w1,x,y);
    if(r1%d)
    {
        flag=1;
        return -1;
    }
    else
    {
        x=x*r1/d;
        s=w1/d;
        x=(x%s+s)%s;
        if(x==0)
        {
            flag=1;
            return -1;
        }
        else
        {
            return x;
        }

    }
}

int main()
{
    int k;
    scanf("%d",&k);
    for(int i=1;i<=k;i++)
    {
        scanf("%lld%lld",&w[i],&r[i]);
    }
    flag=0;
    lli ans=China(k);
    if(flag) printf("-1\n");
    else printf("%lld\n",ans);

}