/*********************************************

欧几里德算法求最大公约数

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    :    %ldn",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);
}