hdu/hoj 3875 Euclidean Algorithm

首先，不得不说这是一道综合性较强的数论题。

 

题意：给你n，c,其中n = pq,1<=p<=2^25,1<=q<=2^25,1<=c<=10^9,求

     

 是否是c的倍数，是的话输出yes，否则输出no。

 

思路：对于gcd的求和，可以推出

 

可以证明gcd的和是积性函数，所以对n分解因子，得到

 只要将每个的gcd之和 求出，相乘即可得到答案。

令gcd之和为函数F(x),则F() = .

于是gcd之和解决。

 

接下来，求lcm之和。对于lcm的求和，可以推出如下公式

 

其中可以证明，

是积性函数。令其值为h(n).

于是，可以将n分解质因子，得到 的形式后，

求出h( )的值相乘即可。而

 

 

于是  的值可以求出，进而可以求出

的值。于是lcm的和得到解决。

对n质因子分解，由于p,q较大，所以采用Pollard_rho 算法分解，时间复杂度为 

。注意除法取模的处理。

源代码：

#include <stdio.h>
#include <ctime>
#include <algorithm>
#include <stdlib.h>
#include <math.h>
const int TIME = 12;
long long factor[100],fac_num[100];
int count;
int cmp(const void *a,const void *b)
{
    long long temp;
    temp = *(long long*)a-*(long long *)b;
    if(temp==0)
     return 0;
    else if(temp<0)
     return -1;
    else return 1;
}
long long gcd(long long a,long long b)
{
    if(b==0) return a;
    return gcd(b,a%b);
}
long long multi_mod(long long a,long long b,long long n)
{
    long long res = 0;
    a%=n;
    for(;b;b>>=1)
    {
      if(b&1)
        res = (res+a)%n;
      a<<=1;
      a%=n;
    }
    return res;
}
long long multi_exp(long long a,long long b,long long n)
{
    long long res = 1;
    for(;b;b>>=1)
    {
       if(b&1)
         res = multi_mod(res,a,n);
       a = multi_mod(a,a,n);
    }
    return res;
}
bool witness(long long a,long long n)
{
   long long  m,x,y;
   int i,j = 0;
   m = n-1;
   while(m%2==0)
   {
       j++;
       m>>=1;
   }
   x = multi_exp(a,m,n);
   for(i = 1;i<=j;i++)
   {
       y = multi_exp(x,2,n);
       if(y==1&&x!=1&&x!=n-1)
         return false;
       x = y;
   }
   if(y!=1)
     return false;
   return true;
}
bool Miller_Rabin(long long n,int times)
{
    int i;
    long long a;
    if(n==1)
      return false;
    if(n==2)
      return true;
    if(n%2==0)
      return false;
    for(i = 1;i<=times;i++)
    {
        a = rand()%(n-2)+2;
        if(!witness(a,n))
          return false;
    }
    return true;
}
long long pollard_rho(long long n)
{
    long long i = 1, x = (long long)fabs(1.0*(rand() % (n - 1))) + 1, y, k, d, c;
    y = x, k = 2;
    do { c = (long long)fabs(1.0*(rand() % (n - 1))) + 1; } while (c == 0 || c == 2);
    do
    {
        i++;
        d = gcd(n + y - x, n);
        if (d != 1 && d != n) return d;
        if (i == k) y = x, k <<= 1;
        x = (multi_mod(x, x, n) + n - c) % n;
    } while (y != x);
    return n;
}
void get_factor(long long n)
{
    long long temp;
    if(n==1)
     return ;
    if(Miller_Rabin(n,TIME)==true)
    {
     factor[count++] = n;
     return ;
    }
    temp = n;
    while(temp>=n)
    temp = pollard_rho(temp);
    get_factor(temp);
    get_factor(n/temp);
}
int main()
{
    long long n,c,temp1,temp2,temp3,keep;
    int count1,i,count2,T;
    srand(time(NULL));
    scanf("%d",&T);
    while(T--)
    {
        scanf("%I64d%I64d",&n,&c);
        count = 0;
        get_factor(n);
        qsort(factor,count,sizeof(factor[0]),cmp);
        count1 = 0;
        count2 = 1;
        for(i = 1;i<count;i++)
        {
            if(factor[i]!=factor[i-1])
            {
                factor[count1] = factor[i-1];
                fac_num[count1++] = count2;
                count2 = 1;
            }
            else
                count2++;
        }
        factor[count1] = factor[count-1];
        fac_num[count1++] = count2;
        temp1 = 1;
        for(i=0;i<count1;i++)
        {
            temp3 = (fac_num[i]+1)*multi_exp(factor[i],fac_num[i],c)-fac_num[i]*multi_exp(factor[i],fac_num[i]-1,c);
            temp3 = (temp3+c)%c;
            temp1*=temp3;
            temp1%=c;
        }
        temp2 = 1;
        for(i = 0;i<count1;i++)
        {
            keep = (factor[i]+1)*2*c;
            temp3 = (multi_exp(factor[i],2*fac_num[i]+1,keep)+1)%keep/(factor[i]+1);
            temp2*=temp3;
            temp2%=2*c;
        }
        temp3 = (temp2-1+2*c)%(2*c);
        temp2 = ((((n%(2*c))*temp3)%(2*c)/2)+n%c)%c;
        temp3 = (temp2-temp1+c)%c;
          if(temp3==0)
         printf("yes\n");
        else
         printf("no\n");
    }
    return 0;
}