The Boss on Mars
Problem's Link
Mean:
给定一个整数n，求1~n中所有与n互质的数的四次方的和.(1<=n<=1e8)
analyse:
看似简单，倘若自己手动推公式的话，还是需要一定的数学基础.
总的思路：先求出sum1=(1^4)+(2^4)+...(n^4)，再求出sum2=（1~n中与n不互质的数的四次方的和），answer=sum1-sum2.
如何求sum1呢？
有两种方法：
1.数列差分.由于A={Sn}={a1^4+a2^4+...an^4}对应一个五阶线性差分方程，只需要求出这个五阶线性差分方程的系数即可.
有关数列差分求幂数和通项的知识，click here.
2.利用低次幂数和来递推高次幂数和公式.
最终求得的公式为：Sn=(n*(n+1)*(2n+1)*(3*n*n+3*n-1))/30.
注意，上式中最后有除法，而我们的最终答案要对1e9+7取余，所以需要求30对1e9+7的模逆元.
由于1e9+7是质数，所以可以直接使用结论：
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;
-------------------------------------------------------------------------
求sum2时需要用容斥，当然直接容斥暴力统计的话也会超时.
注意到：
2^4+4^4+6^4+8^4 = 2^4*(1^4+2^4+3^4+4^4) .
所以再求sum2时仍然可以使用幂数求和公式，这样一来时间复杂度就非常低了.
Time complexity: O(logn)
view code
* this code is made by crazyacking
* Verdict: Accepted
* Submission Date: 2015-10-10-22.47
* Time: 0MS
* Memory: 137KB
*/
#include <queue>
#include <cstdio>
#include <set>
#include <string>
#include <stack>
#include <cmath>
#include <climits>
#include <map>
#include <cstdlib>
#include <iostream>
#include <vector>
#include <algorithm>
#include <cstring>
#define max(a,b) (a>b?a:b)
using namespace std;
typedef long long( LL);
typedef unsigned long long( ULL);
const double eps( 1e-8);
const int MOD = 1000000007;
typedef long long LL;
int cnt = 0;
int a [ 50 ];
LL n , inv;
// prime factor decomposition
void primeFactorization( int n)
{
cnt = 0;
for( int i = 2; i * i <=n; i ++)
{
if(n % i == 0)
{
a [ cnt ++ ] = i;
while(n % i == 0)
n /= i;
}
}
if(n > 1) a [ cnt ++ ] =n;
}
LL mulMod( LL a , LL b , LL m)
{
LL ans = 0;
while(b)
{
if(b & 1)
{
ans = ( ans + a) % m;
b --;
}
b >>= 1;
a = a * 2;
if( a >n) a %= m;
}
return ans;
}
LL quickPow( LL a , LL b , LL m)
{
LL ans = 1;
while(b)
{
if(b & 1)
{
ans = mulMod( ans , a , m);
b --;
}
b >>= 1;
a = mulMod( a , a , m);
}
return ans;
}
// Exponential sum
LL f( LL n)
{
LL ans =n;
ans =( ans *(n + 1)) % MOD;
ans =( ans *( 2 *n + 1)) % MOD;
ans =( ans *(( 3 *n *n + 3 *n - 1) % MOD)) % MOD;
ans =( ans * inv) % MOD;
return ans;
}
LL solve( LL n)
{
primeFactorization(n);
LL ans = 0;
for( int i = 1; i <( 1 << cnt); i ++)
{
LL tmp = 1;
int tt = 0;
for( int j = 0; j < cnt; j ++)
{
if(( 1 << j) & i)
{
tmp = tmp * a [ j ];
tt ++;
}
}
LL q =n / tmp;
LL t = mulMod( mulMod( tmp , tmp , MOD ), mulMod( tmp , tmp , MOD ), MOD);
if( tt & 1)
ans =( ans + mulMod( f( q ), t , MOD)) % MOD;
else
ans =( ans - mulMod( f( q ), t , MOD) + MOD) % MOD;
}
return ans;
}
int main()
{
int t;
scanf( "%d" , & t);
while( t --)
{
scanf( "%I64d" , &n);
if(n == 1)
{
puts( "0");
continue;
}
inv = quickPow( 30 , MOD - 2 , MOD);
LL tmp = f(n);
LL ans = solve(n);
ans =( tmp - ans + MOD) % MOD;
printf( "%I64d \n " , ans);
}
return 0;
}
