hdu 4910 Problem about GCD （miller rabbin 大质数判定）

Problem about GCD
Time Limit: 2000/1000 MS (Java/Others)    Memory Limit: 32768/32768 K (Java/Others)
Total Submission(s): 815    Accepted Submission(s): 169

 

Problem Description

Given integer m. Find multiplication of all 1<=a<=m such gcd(a, m)=1 (coprime with m) modulo m.

 

Input

Input contains multiple tests,  one test per line.
Last line contains -1, it should be skipped.

[Technical Specification]
m <= 10^18

 

Output

For each test please output result. One case per line. Less than 160 test cases.

 

Sample Input

1
2
3
4
5
-1

 

Sample Output

0
1
2
3
4

 

Source

 BestCoder Round #3 

 

Recommend

We have carefully selected several similar problems for you:  5654 5653 5652 5651 5650

题目大意：给定一个数字n，求出1 - n之间与他互质的数的乘积mod n

题解：这道题需要打标找规律，但是傻逼的我并没有找出来，这看出了4的倍数（不包括4）答案一定是1，于是只能从网上查了一下规律。

规律转载自：http://www.cnblogs.com/hundundm/p/3889927.html

首先 m = 1 时 ans = 0
对于 m > 1 的 情况
　　由于 1 到 m-1 中所有和m互质的数字，在 对m的乘法取模 运算上形成了群
　　ai =  ( 1<=a<m && gcd(a,m) == 1 )
　　所以 对于 a 必然存在b = a^(-1) = inv(a) 使得 a * b = 1 (mod m)
　　这里存在两种情况
　　a != b 那么最后的连乘式中a b均出现一次，相乘得1
　　a == b 那么最后的连乘式中只出现一个a
　　实际上所有 a = inv(a) 的 ai 连乘就是答案
　　　　继续考虑假如 gcd(a,m) == 1 则 gcd(m - a, m) == 1
　　　　记m - a = -a (mod m)
　　　　那么 a * (-a) = - (a*a) = -1 (mod m)
　　　　　　m != 2时, m - a != a (否则 a = m/2 , gcd(m, m/2) = m/2 != 1)
　　　　　　　　所以a 和 -a 总是成对出现
　　　　　　　　所以a^2 = 1 (mod m)的解的个数/2 为奇数时，答案为-1，为偶数时 答案为1
　　　　　　m == 2时，求得答案为1(由于此时1和-1等价，出现了特殊性)
　　　　　　
　　　　　　所以对于m > 2的情况，只需求a^2 = 1 (mod m)的解的个数是不是4的倍数

 

a^2 = 1 (mod m) 等价变换
(a + 1)(a - 1) = 0 (mod m)
假设 m = p0^k0 * p1^k1 * ... * pi^ki (pi为素数)
那么根据中国剩余定理 原方程等价于
方程组 (a + 1)(a - 1) = 0 (mod pi^ki)
　　先考虑单个方程:
　　　　pi > 2 时，(a + 1) 和 (a - 1) 必定有一个和pi互质(否则 pi % 2 == 0)
　　　　所以该条方程的解为 ±1 (mod pi^ki)
　　
　　　　pi == 2时，
　　　　　　k == 1时 方程解为 1 (mod 2)
　　　　　　k == 2时 方程解为 ±1 (mod 4)
　　　　　　k > 2 方程解为 ±1, (2^(k-1)+1), (2^(k-1)-1) (mod 2^k)
　　当方程组只有一条方程时，情况如上所示
　　然后考虑多条方程，合并的情况
　　　　根据中国剩余定理，各个式子的各个取值，所有情况在范围内均有且只有一个解
　　　　所以方程组解的个数，就是各个方程解的个数的乘积
m > 2时，解的个数不是4的倍数情况（也就是2）只有以下几种
m = 2^2 = 4
m = p^k (p != 2, 且为素数)
m = 2 * p^k (p != 2, 且为素数)

 

 

简而言之：就是1，2，4，以及只有一个质因子(奇数)或者它的1/2只有一个质因子(偶数)，答案是N-1.其余情况答案均是 1。

 

知道这些之后发现n的范围是10^18，如果分解质因子的话肯定会TLE。

于是接触到了Miller_rabin算法，神奇的随机性素数测试算法

解析转载自：http://blog.csdn.net/z690933166/article/details/9860937

普通的素数测试我们有O(√ n)的试除算法。事实上，我们有O(slog³n)的算法。

定理一：假如p是质数，且(a,p)=1，那么a^(p-1)≡1(mod p)。即假如p是质数，且a,p互质，那么a的(p-1)次方除以p的余数恒等于1。（费马小定理）

该定理的逆命题是不一定成立的，但是令人可喜的是大多数情况是成立的。

于是我们就得到了一个定理的直接应用，对于待验证的数p，我们不断取a∈［1，p-1]且a∈Z，验证a^(p-1) mod p是否等于1，不是则p果断不是素数，共取s次。其中a^(p-1) mod p可以通过把p-1写成二进制，由(a*b)mod c=(a mod c)*b mod c，可以在t=log(p-1)的时间内计算出解，如考虑整数相乘的复杂度，则一次计算的总复杂度为log³(p-1)。这个方法叫快速幂取模。

为了提高算法的准确性，我们又有一个可以利用的定理。
定理二：对于0<x<p，x^2 mod p =1 => x=1或p-1。

我们令p-1=(2^t)*u，即p-1为u二进制表示后面跟t个0。我们先计算出x[0]=a^u mod p ，再平方t次并在每一次模p，每一次的结果记为x[i]，最后也可以计算出a^(p-1) mod p。若发现x[i]=1而x[i-1]不等于1也不等于p-1，则发现p果断不是素数。

可以证明，使用以上两个定理以后，检验s次出错的概率至多为2^(-s)，所以这个算法是很可靠的。

需要注意的是，为了防止溢出（特别大的数据），a*b mod c 也应用类似快速幂取模的方法计算。当然，数据不是很大就可以免了。

 

 

然后就可以做这道题了。先去试10^6以内的素数，判断能否由一个质因子一直除到x==1，如果10^6范围内找不到可以整除的数，再利用miller rabin判下是否是素数，如果不是的话，把这个数字开根在平方，判断是不是完全平方数(一个数如果是另一个整数的完全平方，那么我们就称这个数为完全平方数)，这样做的原因是数字最大10^18,如果没有10^6以内的质因子，又不是质数的话，那么他最多只能包含2个质因子了，因为如果包含3个10^6以外的质因子，那么这个数一定不在10^18范围之内，那么如果他不是一个完全平方数的话，那么就肯定不是了。

 

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<cstdlib>
#include<ctime>
#define ll long long
using namespace std;
ll x;
ll prime[1000003],p[1000003];
void calc()
{
 for (ll i=2;i<=1000000;i++)//筛法求素数 
  {
      if (!p[i]) prime[++prime[0]]=i;
      for (ll j=1;j<=prime[0];j++)
      {
        if (i*prime[j]>1000000) break;
        p[i*prime[j]]=1;
        if (i%prime[j]==0) break;
      }
  }
}
ll mul(ll a, ll b, ll mod)  //快速求a * b % n 的答案 ，例如: b = 1011101那么a * b mod n = (a * 1000000 mod n + a * 10000 mod n + a * 1000 mod n + a * 100 mod n + a * 1 mod n) mod n  
{  
    ll res=0;  
    while (b)  
    {  
        if (b&1) res = (res+a)%mod;  
        a = (a+a)%mod;  
        b >>= 1;  
    }  
    return res;  
}  
ll make_pow(ll x,ll num,ll p)//快速幂 
{
  ll sum=1;  ll k=x%p;
  while (num!=0)
  {
   if (num&1)  sum=mul(sum,k,p);
   num=num>>1;
   k=mul(k,k,p);
  }
 return sum;
}
bool miller_rabbin(ll n)//miller rabbin 算法，随机性素数测试 
{
 if (n==2) return true;
 if (n<2||!(n&1)) return false;
 ll t=0;
 ll a,x,y,u=n-1;
 while (!(u&1))  t++,u>>=1;//如果u为偶数则u右移，用k记录移位数 
 for (int i=0;i<=1000;i++)
 {
  a=rand()*rand()%(n-1)+1;
  x=make_pow(a,u,n);
  for (int j=0;j<t;j++)//把移位减掉的量补上，并在这里加上二次探测
  {
      y=mul(x,x,n);
      if (y==1&&x!=1&&x!=n-1) return false;//二次探测定理，这里如果y = 1则x 必须等于 1，或则 n-1否则可以判断不是素数 
      x=y;
  }
  if (x!=1) return false;//利用费马小定理 
 }
 return true;
}
bool check(long long  x)
{
 long long num=0;
 for (ll i=1;prime[i]<=x&&i<=prime[0];i++)
 {
  if (x%prime[i]==0)
  {
    num++;
    while (x%prime[i]==0) x/=prime[i];
    if (num>1) return false;
  }
  if (x==1)  break;
 }
 if (x>1) num++;
 if (num>1)  return false;
 if (miller_rabbin(x))  return true;
 ll m=(ll)sqrt(x*1.0);
 if (m*m==x) return true;
 return false;
}
int main()
{
  calc();
  scanf("%I64d",&x);
  while (x!=-1)
  {
   if (x!=4&&x%4==0)
    {
      printf("1\n");
      scanf("%I64d",&x);
      continue;
    }
   if (x==1||x==2||x==4)
    {
     printf("%I64d\n",x-1);
     scanf("%I64d",&x);
     continue;
    }
   if (x%2)
    {
     if (check(x))
      printf("%I64d\n",x-1);
    else printf("1\n");
    }
   else
    {
     if (check(x/2))
      printf("%I64d\n",x-1);
     else printf("1\n");
    }
   scanf("%I64d",&x);
  }
}