荔园在线
荔园之美，在春之萌芽，在夏之绽放，在秋之收获，在冬之沉淀
[回到开始]
[上一篇][下一篇]
发信人: kaman (Devil Aspire), 信区: ACMICPC
标 题: Number Theory Std Algorithms
发信站: 荔园晨风BBS站 (Tue May 9 10:54:18 2006), 站内
/*==========================================*\
*
* 数论算法函数库
*
* copyright starfish
* 2000/10/25
*
\*==========================================*/
/*********************************************
欧几里德算法求最大公约数
copyright starfish
2000/10/24
*********************************************/
//classic euclid algorithm to calculate the gcd(a,b),
//ie. greatest common divisor of a and b
int euclid(int a,int b)
{
if (b==0) return a;
else return(euclid(b,a%b));
}
/*********************************************
扩展欧几里德算法求gcd(a,b)=ax+by
copyright starfish
2000/10/24
*********************************************/
//extended euclid algorithm to calculate the gcd(a,b),
//as well as the integer x and y where gcd(a,b)=a*x+b*y
int ext_euclid(int a,int b,int &x,int &y)
{
int t,d;
if (b==0) {x=1;y=0;return a;}
d=ext_euclid(b,a %b,x,y);
t=x;
x=y;
y=t-a/b*y;
return d;
}
/*********************************************
Miller_Rabin随机素数测试算法
说明：这种算法可以快速地测试一个数是否
满足素数的必要条件，但不是充分条件。不
过也可以用它来测试素数，出错概率很小，
对于任意奇数n>2和正整数s,该算法出错概率
至多为2^(-s)，因此，增大s可以减小出错概
率，一般取s=50就足够了。
*********************************************/
int Witness(int a,int n)
{
int i,d=1,x;
for (i=ceil(log((float)n-1)/log(2.0))-1;i>=0;i--)
{ x=d;
d=(d*d)%n;
if ((d==1)&&(x!=1)&&(x!=n-1)) return 1;
if (((n-1)&(1<<i))>0) d=(d*a)%n;
}
return (d==1?0:1);*
}
int Miller_Rabin(int n,int s)
{
int j,a;
for (j=0;j<s;j++)
{ a=rand()*(n-2)/RAND_MAX+1;
if (Witness(a,n)) return 0;
*}
return 1;*
}
/*********************************************
*
**返回x的二进制长度
*********************************************/
int BitLength(int x)
{
*int d = 0;
*while (x > 0) {
**x >>= 1;
**d++;
*}
*return d;
}
/*********************************************
*
**返回x的二进制表示中从低到高的第i位
**，最低位为第一位
*********************************************/
int BitAt(int x, int i)
{
*return ( x & (1 << (i-1)) );
}
/*********************************************
模取幂运算 计算a^b mod n
*********************************************/
int Modular_Expoent(int a,int b,int n)
{
int i, y=1;
for (i = BitLength(b); i > 0; i--)
{
y = (y*y)%n;
if (BitAt(b,i) > 0)
*** y = (y*a)%n;
}
return y;
}
/********************************************
求解模线性方程 ax=b (mod n) ,n>0
*********************************************/
void modular_linear_equation_solver(int a,int b,int n)
{
int e,i,d;
int x,y;
d=ext_euclid(a,n,x,y);
if (b%d>0) printf("No answer!\n");
else
{ e=(x*(b/d))%n;
for (i=0;i<d;i++) //notice! here x maybe <0
printf("The %dth answer is : %ld\n",i+1,(e+i*(n/d))%n);
}
}
/*********************************************
求解模线性方程组 (中国余数定理)
a=B[1] (mod W[1])
a=B[2] (mod W[2])
........
a=B[n] (mod W[n])
其中W,B已知，W[i]>0且W[i]与W[j]互质, 求a
*********************************************/
int China(int B[],int W[],int k)
{ int i;
int d,x,y,a=0,m,n=1;
for (i=0;i<k;i++)
n*=W[i];
for (i=0;i<k;i++)
{ m=n/W[i];
d=ext_euclid(W[i],m,x,y);
a=(a+y*m*B[i])%n;
}
if (a>0) return a;
else return(a+n);
}
--
Long long ago,there is a hero stand at the edge of the land,a sword in
hand.......
︵～` ╭) ▁▁
︸ヾ) (╯ ▕天外▏
ゞ7( ︵ )) ▕飞仙▏
╰ ︶へ︶╯ ▔▔
ヅ╯ ジ〉
※ 来源:·荔园晨风BBS站 bbs.szu.edu.cn·[FROM: 192.168.20.225]
[回到开始]
[上一篇][下一篇]
荔园在线首页 友情链接：深圳大学 深大招生 荔园晨风BBS S-Term软件 网络书店