/* L07.c 
中国剰余アルゴリズム
*/

#include<stdio.h>
#include<stdlib.h>

/* 最大公約数 */
int gcd(int a, int b)
{
	return (b == 0) ? abs(a) : gcd(b,a%b);
}

/* ユークリッドのアルゴリズム拡張版 */
/* d = gcd(a,b) と d = ax + by を満たす x,y を返す */
int euclid(int a, int b, int *x, int *y)
{
	int r0=a,r1=b,r2,x0=1,x1=0,x2,y0=0,y1=1,y2,q;
	while(r1!=0){
		q=r0/r1;
		r2=r0-q*r1;
		x2=x0-q*x1;
		y2=y0-q*y1;
		// 番号を振り替える
		r0=r1;
		r1=r2;
		x0=x1;
		x1=x2;
		y0=y1;
		y1=y2;
	}
	if(r0<0){
		r0=-r0;
		x0=-x0;
		y0=-y0;
	}
	*x=x0;
	*y=y0;
	return(r0);
}

/* a を法 n の数に取り直す */
int mod(int a, int n)
{
	int b=a%n;
	return b<0 ? b+n : b;
	// C言語では a<0, n>0 のときの符号は処理系に依存
}

/* 中国剰余アルゴリズム */
/* 連立合同式  x ≡ a mod m, x ≡ b mod n の解 x を求める */
int chinese(int a, int m, int b, int n)
{
	int d,u,v,mn=m*n;
	d=euclid(m,n,&u,&v);
	if(d != 1){
		printf("法が互いに素ではありません。\n");
		exit(0);
	}
	return mod(n*mod(a*v,m)+m*mod(b*u,n),mn);
}

/* ３式以上の中国剰余アルゴリズム */
/* 連立合同式  x ≡ a[i] mod n[i] ( i = 0,...,k-1 ) の解 x を求める */
int multichinese(int k, int *a, int *n)
{
	int x,m,i;
	x=mod(a[0],n[0]);
	m=n[0];
	for(i=1;i<k;i++){
		x=chinese(x,m,a[i],n[i]);
		m*=n[i];
	}
	return x;
}

main()
{
	int x,a[8],n[8],i,k,m;

	printf("中国剰余アルゴリズムを用いて連立合同式\n");
	printf("   x ≡ a[i] mod n[i] ( i = 1,...,k )\n");
	printf("の解を求めます。\n\n");

	printf("式は幾つありますか。\n");
	scanf("%d",&k);
	if(k <= 0) exit(1);

	for(i=0;i<k;i++){
			printf("a[%d] = ",i+1);
			scanf("%d",a+i);
			printf("n[%d] = ",i+1);
			scanf("%d",n+i);
	}

	x=multichinese(k,a,n);
	printf("\n解 x = %d\n\n",x);
	printf("検算 :\n");
	for(i=0;i<k;i++) printf("x ≡ %d mod %d\n",mod(x,n[i]),n[i]);
}

/* 実行例

中国剰余アルゴリズムを用いて連立合同式
   x ≡ a[i] mod n[i] ( i = 1,...,k )
の解を求めます。

式は幾つありますか。
3
a[1] = 1
n[1] = 3
a[2] = 2
n[2] = 5
a[3] = 3
n[3] = 7

解 x = 52

検算 :
x ≡ 1 mod 3
x ≡ 2 mod 5
x ≡ 3 mod 7

*/