ckfan's Blog

这是在is-programmer的第一篇blog，拉宾米勒算法由一同学提供，谢谢他！

#define type int
type get_rand()
{
     type a;
     a  = rand();
     a* = rand();
     return abs(a);
}

/******a^b mod c *******/
type PowerMod(type a, type b, type c)
{
     type r =1 ;
     while(b)
     {
          if( b & 1 )
          {
               r = r * a % c;
          }
          b >>= 1;
          a = a * a % c;
     }
     return r;
}

/********* 素数随机判定算法 *******/
type mod_mul (type x , type y , type n)
{
     type T =(type)( sqrt((double)n) + 0.5);
     type t = T * T - n;
     type a = x / T;
     type b = x %T;
     type c = y / T;
     type d = y % T;
     type e = a * c / T;
     type f = a * c % T;
     type v = (( a * d + b * c) % n + e * t ) % n;
     type g = v / T;
     type h = v % T;
     type ans = ((( f + g) * t % n + b * d) % n + h * T) % n;
     while ( ans < 0)
     ans += n;
     return ans ;
}
/*
 米勒-拉宾法算法
*/
type Rabin_Miller(type n)
{
     type k = 0, i, j, m, a;
     if( n < 2)
          return 0;
     if(n == 2)
          return 1;
     if( !(n & 1) )
      return 0;

     m = n - 1;
 
     while( !( m & 1 ) )
          m >>= 1, k++;
 
     for(i = 0; i < 10; i++)
     {
          a = get_rand() % (n - 2) + 2;
          a = PowerMod(a, m, n);
          if(a == 1)
               continue;
          for(j = 0; j < k; j++)
          {
               if(a == n-1)
                    break;
               a = mod_mul(a, a, n);
          }
          if(j < k)
               continue;
      return 0;
     }
     return 1;
}

/*
获取特定长度的随机素数
*/
type Get_Rand_Prime(type bytes)
{
     type i;
     type prime;

     bytes = 0;
     prime = 0;
     i = 0;
     while (1)
     {
          prime = 0;
          i = bytes;
          while(i--)
          {
               prime = prime * 10 + rand() % 10;
               if (i ==( bytes - 1) && prime == 0)
               {
                    prime = 1;
               }
          }
          prime = prime | 0x01;
          if (Rabin_Miller(prime))
          {
               break;
          }
    }
    return prime;
}