The Boss on Mars
Time Limit: 2000/1000 MS (Java/Others) Memory Limit: 32768/32768 K (Java/Others)
Total Submission(s): 2109 Accepted Submission(s): 637
Problem Description On Mars, there is a huge company called ACM (A huge Company on Mars), and it’s owned by a younger boss.
Due to no moons around Mars, the employees can only get the salaries per-year. There are n employees in ACM, and it’s time for them to get salaries from their boss. All employees are numbered from 1 to n. With the unknown reasons, if the employee’s work number is k, he can get k^4 Mars dollars this year. So the employees working for the ACM are very rich.
Because the number of employees is so large that the boss of ACM must distribute too much money, he wants to fire the people whose work number is co-prime with n next year. Now the boss wants to know how much he will save after the dismissal.
Input The first line contains an integer T indicating the number of test cases. (1 ≤ T ≤ 1000) Each test case, there is only one integer n, indicating the number of employees in ACM. (1 ≤ n ≤ 10^8)
Output For each test case, output an integer indicating the money the boss can save. Because the answer is so large, please module the answer with 1,000,000,007.
Sample Input
2
4
5
Sample Output
82
354
Hint
Case1: sum=1+3*3*3*3=82
Case2: sum=1+2*2*2*2+3*3*3*3+4*4*4*4=354
Author ZHANG, Chao
Source 2011 Asia Dalian Regional Contest
Recommend lcy | We have carefully selected several similar problems for you: 4056 4053 4057 4052 4054
這道題需要求的是在1~n的範圍内與n互質的數的4次方之和。我們單純地去周遊1~n判斷數是否與n的公約數為1是最容易想到的做法。但是這樣的複雜度太高，沒有可行性。因而我們需要從反面考慮，我們可以求出1~n之間不與n互質的所有數之和，然後用1^4+2^4+……+n^4之和去減。求1~n之間不與n互質的所有數之和需要用到容斥原理。
容斥原理是離散數學集合論中的一個重要結論。相當于是求集合的并集中有多少個元素。是以我們需要先對n進行分解質因數，然後再進行容斥原理的操作。先加上1~n中能被一個質因數整除的數的4次方之和，再減去能被兩個質因數之積整除的數的4次方之和，再加上能被三個質因數之積整除的數的4次方之和……
關于容斥原理，有兩種寫法，一種是遞歸寫法，另一種是非遞歸寫法。
遞歸寫法如下：
ll dfs(int indx,ll n) //表示從indx這個位置開始,不與n互質的所有數的四次方之和
{
ll res=0;
for(int i=indx;i<s;i++) //s表示質因子的個數
{
ll cnt=factor[i]; //表示目前枚舉的質因子
res=(res+sum(n/cnt)*pow_mod(cnt,4))%mod; n/cnt表示在1~n的範圍内能被cnt整除的數的個數
res=(res-dfs(i+1,n/cnt)*pow_mod(cnt,4))%mod;
}
return res;
}
當cnt取第一個位置的時候，與後邊的(s-1)個相乘形成(s-1)個兩個質因子，cnt取第二個位置時， 與後邊的(s-1)個相乘形成(s-1)個兩個質因子……這樣兩個質因子之積的個數就是1+2+……+(s-1)=s*(s-1)/2個，這就是C(s,2)個。是以是減号。當進入第二層遞歸時，這時暫不考慮第一層遞歸與之相乘的數是哪個，隻考慮目前的遞歸。把目前遞歸層的結果取為正号，再下一層的取為負号，這樣因為遞歸的層層調用，可以保證正負号的正确性。同時代碼比較簡潔。但是了解起來有些困難。不過仔細思考該遞歸的過程，這種遞歸版的容斥原理寫法還是很優雅的。
此外，還有非遞歸寫法，用位運算來實作，這裡我們需要考慮數字有s個質因子，相當于s個盒子，把s個數放進盒子中，不一定要保證s個數都放進盒子裡，每個盒子可放可不放。某個盒子放了數字記為1，沒放記為0。則相當于0~(1<<s)-1這些數字的二進制變化。二進制表示中含有0個1的數有C(s,0)個，含有1個1的數有C(s,1)個，含有2個1的數有C(s,2)個,以此類推。這樣了解以後代碼也比較好寫。
非遞歸寫法如下：
ll sum1(ll k)
{
int t=(1<<s);
ll res=0;
for(int i=1;i<t;i++)
{
ll num=0,p=1;
for(int j=0;j<s;j++)
{
if(i&(1<<j))
{
num++;
p=(p*factor[j])%mod;
}
}
if(num%2) res=(res+sum(k/p)*pow_mod(p,4))%mod;
else res=((res-sum(k/p)*pow_mod(p,4))%mod+mod)%mod;
}
return res;
}
下面需要說的是1^4+2^4+……+n^4的公式。這個可以根據二項式定理進行推導。
1^3=(1+0)^3=1+0+0+0
2^3=(1+1)^3=1+3*1+3*1^2+1*1^3
3^3=(1+2)^3=1+3*2+3*2^2+1*2^3
4^3=(1+3)^3=1+3*3+3*3^2+1*3^3
……
n^3=(1+n-1)^3=1+3*(n-1)+3*(n-1)^2+1*(n-1)^3
(n+1)^3=(1+n)^3=1+3*n+3*n^2+1*n^3
設S(1)=1+2+3+……+n=n*(n+1)/2,
S(2)=1^2+2^2+……+n^2
将(n+1)個等式相加可得(n+1)^3=(n+1)+3*(1+2+3+……+n)+3*(1^2+2^2+……+n^2)=(n+1)+3*n*(n+1)/2+3*S(2);
等式兩邊乘以2，得到2*(n+1)^3=2*(n+1)+3*n*(n+1)+6*S(2)
S(2)=n*(n+1)*(2*n+1)/6
同理可得，
1^4=(1+0)^4=1+0+0+0+0
2^4=(1+1)^4=1+4*1+6*1^2+4*1^3+1*1^4
3^4=(1+2)^4=1+4*2+6*2^2+4*2^3+1*2^4
4^4=(1+3)^4=1+4*3+6*3^2+4*3^3+1*3^4
……
n^4=(1+n-1)^4=1+4*(n-1)+6*(n-1)^2+4*(n-1)^3+1*(n-1)^4
(n+1)^4=(1+n)^4=1+4*n+6*n^2+4*n^3+1*n^4
設S(1)=1+2+3+……+n=n*(n+1)/2,
S(2)=1^2+2^2+……+n^2
S(3)=1^3+2^3+……+n^3
将這(n+1)個式子相加并将S(1),S(2),S(3)代入得，S(3)=n*n*(n+1)*(n+1)/4
同理可得，
1^5=(1+0)^5=1+0+0+0+0+0
2^5=(1+1)^5=1+5*1+10*1^2+10*1^3+5*1^4+1*1^5
3^5=(1+2)^5=1+5*2+10*2^2+10*2^3+5*2^4+1*2^5
4^5=(1+3)^5=1+5*3+10*3^2+10*3^3+5*3^4+1*3^5
……
n^4=(1+n-1)^5=1+5*(n-1)+10*(n-1)^2+10*(n-1)^3+5*(n-1)^4+1*(n-1)^5
(n+1)^4=(1+n)^5=1+5*n+10*n^2+10*n^3+5*n^4+1*n^5
設S(1)=1+2+3+……+n=n*(n+1)/2,
S(2)=1^2+2^2+……+n^2
S(3)=1^3+2^3+……+n^3
S(4)=1^4+2^4+……+n^4
将這(n+1)個式子相加并将S(1),S(2),S(3),S(4)代入得，S(4)=n*(n+1)*(2*n+1)*(3*n*n+3*n-1)/30
這裡要對mod取模，而該式除以了30，除法沒有取模運算，是以要将30變為在MOD下的逆元。
根據費馬小定理可得，b^(MOD-1)=1mod MOD(MOD為質數)。
則等式兩邊同時除以b可得， b^(MOD-2)=(1/b)mod MOD
是以用快速幂取模即可算出30在1000000007下的模逆元。
然後正常計算出1^4+2^4+……+n^4即可。本題比較坑的細節在取模運算上。當進行減法運算時，要防止出現負數。
是以利用該式((a-b)%mod+mod)%mod将負數變為正數。
還有在計算乘法時也要防止溢出，利用(a*(b%mod))%mod可以防止溢出。當a和b相乘之積不會溢出時，直接用(a*b)%mod即可。
參考代碼：
#include<cstdio>
#include<iostream>
#include<cmath>
#include<cstring>
#include<algorithm>
#include<string>
#include<vector>
#include<map>
#include<set>
#include<stack>
#include<queue>
#include<ctime>
#include<cstdlib>
#include<iomanip>
#include<utility>
#define pb push_back
#define mp make_pair
#define CLR(x) memset(x,0,sizeof(x))
#define _CLR(x) memset(x,-1,sizeof(x))
#define REP(i,n) for(int i=0;i<n;i++)
#define Debug(x) cout<<#x<<"="<<x<<" "<<endl
#define REP(i,l,r) for(int i=l;i<=r;i++)
#define rep(i,l,r) for(int i=l;i<r;i++)
#define RREP(i,l,r) for(int i=l;i>=r;i--)
#define rrep(i,l,r) for(int i=1;i>r;i--)
#define read(x) scanf("%d",&x)
#define put(x) printf("%d\n",x)
#define ll long long
#define lson l,m,rt<<1
#define rson m+1,r,rt<<11
using namespace std;
const int mod=1000000000+7;
int t;
ll n,inx; //inx是30在1000000007為mod下的模逆元
int s,factor[110]; //factor數組為n的質因數數組，s表示n的質因數的個數
ll pow_mod(ll x,ll n) //快速幂取模
{
ll ans=1;
while(n)
{
if(n&1) ans=(ans*x)%mod;
x=(x*x)%mod;
n>>=1;
}
return ans;
}
ll sum(ll n) //利用公式求1^4+2^4+……+n^4的值并取模
{
ll ans=n;
ans=(ans*(n+1))%mod;
ans=(ans*(2*n+1))%mod;
ans=(ans*((3*n*n+3*n-1)%mod))%mod; //這裡一定要先對裡邊的3*n*n-3*n+1取模再乘，否則會導緻溢出。這裡雖然減去了1，但是不會出現負數，無須+mod,再模mod
ans=(ans*inx)%mod; //乘以30的乘法逆元
return ans;
}
ll dfs(int indx,ll n) //遞歸版容斥原理
{
ll res=0;
for(int i=indx;i<s;i++)
{
ll cnt=factor[i];
res=(res+sum(n/cnt)*pow_mod(cnt,4))%mod;
res=(res-dfs(i+1,n/cnt)*pow_mod(cnt,4))%mod; //這裡雖然做減法運算，但是可以肯定1~n中被一個質因數整除的數之和必定大于被兩個質因數之積整除的數之和大，比三個質因數之積整除的數之和大。是以不會出現負數，不用先加mod再模mod
}
return res;
}
ll sum1(ll k) //非遞歸版容斥原理
{
int t=(1<<s);
ll res=0;
for(int i=1;i<t;i++)
{
ll num=0,p=1;
for(int j=0;j<s;j++)
{
if(i&(1<<j))
{
num++;
p=(p*factor[j])%mod;
}
}
if(num%2) res=(res+sum(k/p)*pow_mod(p,4))%mod;
else res=((res-sum(k/p)*pow_mod(p,4))%mod+mod)%mod;
}
return res;
}
void get_factor(ll n) //求出n的質因數
{
s=0;
for(ll i=2;i*i<=n;i++)
{
if(n%i==0)
{
factor[s++]=i;
while(n%i==0)
n/=i;
}
}
if(n!=1) factor[s++]=n;
}
void solve()
{
inx=pow_mod(30,mod-2);
get_factor(n);
ll cnt=sum(n),res=dfs(0,n);
cnt=((cnt-res)%mod+mod)%mod; //這裡一定要注意，cnt-res可能為負數，是以要先加上mod,然後模mod
printf("%I64d\n",cnt);
}
int main()
{
read(t);
while(t--)
{
scanf("%I64d",&n);
solve();
}
}
![各種二分查找[圖]](data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACwAAAAAAQABAAACAkQBADs=)