HDU_3875 Euclidean Algorithm 数论

http://acm.hdu.edu.cn/showproblem.php?pid=3875

题意：

给一个n = p * q ，一个c ，问 ( sum{ lcm(i , n) } - sum{ gcd(i , n) } ) % c == 0 是否成立。

思路：
          这是一个数论的综合题，整整花了3个小时啊。我们还是把问题分开处理， 就是先处理ans1 = sum{ gcd(i, n ) } ; 我们令f(x) = sum{ gcd(i , x ) }  ; 那么f(x) 就是一个积性函数。那么我们就可以将n分解质因数之后求解，也就是这样n = p1^k1 * p2 ^ k2 * p3 ^k3 *... * pn^kn ，那么上面的式子就可以变成：f(n) = f( p1^k1 ) * f( p2^k2 ) * f(p3 ^k3 ) * ..... * f( pn^kn) ，下面我们只需要求出每个f就可以了，对于每个f( pi^ki ) = sum{ gcd( i , pi^ki ) } , gcd(i , pi^ki) 只可能是1、pi、pi^2 ....  pi^ki 中的一个，并且对于某个特定的pi^a，一共有phi( pi^(ki-a) )个， 这样最后的和就会变成这样：f( pi^ki )  = 1 * phi( pi^ki ) + pi * phi( pi^( ki-1) ) + .... + pi^ki * phi(1)  = (ki + 1) *pi^ki - ki * pi^(ki-1) 。到这里第一个问题就可以在n分解质因数之后解决了。

         接下去就来解决第二个问题，求ans2 = sum{ lcm(i , n ) } ，ans2 = sum{ i * n / gcd(i , n) }  = n * sum{ i / gcd(i , n) }  = n * sum{ n / d *phi(n / d)  / 2 }   = n / 2 * sum{ n/d * phi(n/d) } = n / 2 * sum{ d * phi(d) } + n / 2 = n / 2 * ( sum( d *phi(d) ) + 1 ) ，问题就变成了求 f(n) = sum{ d * phi(d) } 1<=d<=n && d|n  ，同样这样函数也是积性的，那么问题也同样又变成了将n分解质因数，再求解，到这里我们就不能得出f(pi^ki) = ( pi^*(2*ki+1) + 1 ) / (pi + 1) 这个公式了。 最后值得注意的是，因为n最大可以有2^50，所以我们要用pollard_rho在O(n^1/4)的时间内去分解质因数。 这样本题就可以完美解决了。 

#include <stdio.h>
#include <string.h>
#include <vector>
#include <algorithm>
#include <math.h>
#include <stdlib.h>
#include <time.h>
using namespace std;

typedef __int64 LL ;
const int S = 10 ;
const int CC = 20 ;
LL N, C ;
LL _min ;
vector< LL > prime;
LL ans1 , ans2 ;
vector< pair<LL , LL> > fac ;

LL mul(LL a, LL b ,LL n){
    LL res = 0 , add = a ;
    while( b ){
        if(b&1) res = (res + add) % n ;
        add = (add + add) % n ;
        b >>= 1;
    }
    return res ;
}

LL pow( LL a , LL b , LL n ){
    LL res = 1 , add = a;
     while( b ){
        if( b&1 )   res = mul( res ,add , n) ;
        add = mul( add , add , n ) ;
        b >>= 1 ;
     }
     return res ;
}

int witness(LL a, LL n){
    LL u = n-1 , t = 0 ;
    while( !(u&1) ) u >>=1 , t ++ ;
    LL x = pow( a, u , n ) , y;
    for(int i=1;i<=t;i++){
        y = mul( x , x , n ) ;
        if(y==1 && x!=1 && x!=n-1 ) return 1 ;
        x = y ;
    }
    if( x!=1 )  return 1 ;
    return 0 ;
}


int Miller_Rabin( LL n ){
    if(n==2)  return 1 ;
    if( n<1 || !(n&1) ) return 0 ;
    for(int i=2;i<100;i++)
        if( n%i == 0 )   if( n==i ) return 1 ;
                        else        return 0 ;
    srand( (LL)time( NULL ) )  ;
    for(int i=0;i<S;i++){
        LL a = (LL)rand()%(n-2) + 2 ;
        if( witness(a , n) )    return 0 ;
    }
    return 1 ;
}

LL gcd(LL a, LL b){
    while( b ){
        LL c = a;
        a = b ;
        b = c % b ;
    }
    return a;
}

LL pollard_rho(LL n ,LL c){
    LL x, y, z , i = 1 ,k = 2 ,d ;
    srand( (LL)time(NULL) )  ;
    x = ( (LL)rand() ) % (n-1) + 1 ;
    y = x;
    while( true ){
        i ++ ;
        x = ( mul(x , x , n) + c ) % n ;
        d = gcd( y-x+n , n ) ;
        if( d!=1 && d!=n )  return d ;
        if(y == x)  return n ;
        if(i == k)  y = x , k <<= 1 ;
    }
}

void find( LL n , LL c) {
    if( n<=1 )  return ;
    if( Miller_Rabin(n) ) {
        prime.push_back( n );
        return ;
    }
    LL r = pollard_rho( n , c-- ) ;
    find( n/r , c ) ;
    find( r , c )  ;
}

void get_prime( LL n){
    LL tmp = n  ;
    if( n==1 )  return ;
    if( Miller_Rabin( tmp ) ){
        prime.push_back( tmp )  ;
        return ;
    }
    find( n , CC );
}
int main(){
    int T ; scanf("%d",&T);
    while( T-- ){
        scanf("%I64d %I64d",&N,&C);
        if( N==1 ){
            printf("yes\n");    continue ;
        }
        fac.clear(); prime.clear() ;
        get_prime( N );
        sort( prime.begin() , prime.end() ) ;
        LL front = prime[0] , c = 1 ;
        for(int i=1;i<prime.size();i++){
            if( front != prime[i] ){
                fac.push_back( make_pair(front , c) ) ;
                front = prime[i] ; c = 1 ;
            }
            else    c++ ;
        }
        fac.push_back( make_pair(front, c) ) ;
        ans1 = ans2 = 1 ;
        for(int i=0;i<fac.size();i++){
            LL p = fac[i].first , cc = fac[i].second;
            LL res = ( pow(p , cc , C) - pow(p , cc-1 ,C ) + C ) % C ;
            res = mul( res , cc , C ) ;
            res = (res + pow(p, cc , C)) % C ;
            ans1 = mul( ans1 , res , C ) ;
            res = 1 ;
            LL kk = ( p + 1 ) * 2 * C ; //这里之所以可以这么乘是因为有n=p*q这个条件
            res = ( pow( p , 2*cc+1  ,kk ) + 1 ) % kk ;
            res = ( res + 1 ) / (p + 1) ;
            ans2 = mul( ans2 , res , 2*C ) ;
        }
        ans2 = (ans2 + 1) % (2*C) ;
        ans2 = mul( ans2 , N%(2*C) , 2*C) ;
        ans2 = ( ans2 / 2 ) % C ;
        LL rr = (ans2 - ans1 + C) % C ;
        if(rr == 0 )    printf("yes\n") ;
        else            printf("no\n") ;
    }
    return 0 ;
}