HDU 4059 The Boss on Mars

题目求1-n中与n互质的数的4次方之和，即S=a1^4+a2^4+……; a1,a2……均小于等于n且与n互质。

先求出1^4+2^4+……n^4然后减去与n不互质的数的4次方。

必然要先要用到4次方的求和公式。接下来简单的证明一下，这里前提是你知道3次方的公式，如果不会照下面的模式可以利用2次公式推出3次公式

（x+1)^5=x^5+5*x^4+10*x^3+10*x^2+5*x+1;

则 1=1;

2^5=(1+1)^5=1^5+5*1^4+10*1^3+10*1^2+5*1^1+1;

3^5=(2+1)^5=2^5+5*2^4+10*2^3+10*2^2+5*2^1+1;

……

……

(n+1)^5=(n+1)^5=n^5+5*n^4+10*n^3+10*n^2+5*n^1+1;

全部叠加起来，则(n+1)^5=5*(1^4+2^4+……n^4)+10*(1^3+2^3+……+n^3)+10*(1^2+2^2+……+n^2)+5*(1+2+……+n)+n+1;

然后将(1^3+2^3+……n^4)=(n+1)^2*n^2/4; (1^2+2^2+……n^2)=(n*(n+1)*(2*n+1))/6; 代入。

化简后得到(1^4+2^4+……+n^4)=(n*(n+1)*(2n+1)*(3*n*n+3*n-1))/30;

这里还用到一个公式 a%m = (b/c)%m ==> a % m = b * c ^(m-2)%m ( m为素数 )；

证明：

 b = a * c % m  ===> b = a % m * c %m;  因为根据费马小定理 a^(p-1)= 1 %p;(p是素数)

所以 a % m = a*1%m = a * c^(m-1)%m = a*c*c^(m-2)%m = b*c^(m-2)%m;

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<algorithm>
#include<cmath>
#include<queue>
#include<set>
#include<map>
#include<cstring>
#include<vector>
#include<string>
#define LL long long
using namespace std;
const LL MOD = 1000000007;
LL factor[64],res;

int prime[2000],cnt=0;
void PowMode(  )
{
    res = 1;
    LL a = 30LL;
    LL b = MOD  - 2;
    while( b )
    {
         if( b&1 ) res = (res*a)%MOD;
         a = ( a * a )%MOD;
         b >>= 1;
    }    
}
void Prime( )
{
    bool hash[5024]={0};
    for( int i = 3 ; i <= 101 ; i +=2 ){
         if( !hash[i>>1] ){
             int x = i <<1;
             for( int j = i * i ; j <= 10000; j += x )
                  hash[j>>1] = true;
         }    
    }    
    prime[cnt++] = 2;
    for( int i = 1 ; i <= 5000 ; i++ ){
         if( !hash[i] ){
                prime[cnt++] = ( i << 1 ) + 1;
         }    
    }
}
int Get_factor( LL num )
{
    int count = 0;
    for( int i = 0 ; i < cnt ; i++ ){
            if( num < prime[i] ) break;
            if( num % prime[i] == 0 ){
                factor[count++] =(LL) prime[i];
                while( num % prime[i] == 0 )
                       num /= prime[i];
                }
    }    
    if( num != 1 ) factor[count++] = num;
    return count;
}
LL Get_sum( LL N )
{
   LL ans = N;
   ans = ( ans * ( N + 1 ) )%MOD;
   ans = ( ans * ( 2 *N + 1 ) )%MOD;
   ans = ( ans * ( ( 3 * N * N )%MOD + (3*N)%MOD - 1 + MOD)%MOD )%MOD;
   return (ans * res)%MOD;    
}
LL POW( LL num )
{
   LL ans = (((((num*num)%MOD)*num)%MOD)*num)%MOD;
   return ans; 
}
LL DFS( LL d , int start , int count )
{
   LL ans = 0;
   for( int i = start ; i < count ; i ++ )
   {
       LL temp =  POW( factor[i] ); 
       ans = (ans + ( temp  * Get_sum( d / factor[i] ))%MOD - (temp * DFS( d / factor[i] , i + 1 , count ))%MOD + MOD)%MOD;        
   }
   return ans;
}
int main(  ){
     Prime(  );
     PowMode(  );
     LL N;
     int T;
     while( scanf( "%d",&T )==1 )
     {
         for( int i = 1 ; i <= T ; i ++ ){
                scanf( "%I64d",&N );
                int count = Get_factor( N );
                LL ans = 0;
                ans = ((Get_sum( N   ) - DFS( N  , 0 , count ))%MOD+MOD)%MOD;
                printf( "%I64d\n",ans );
                }        
     }
    //system( "pause" );
    return 0;
}