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
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.
2 4 5
82 354HintCase1: sum=1+3*3*3*3=82 Case2: sum=1+2*2*2*2+3*3*3*3+4*4*4*4=354
这道题需要求的是在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;
}此外，还有非递归写法，用位运算来实现，这里我们需要考虑数字有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^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();
}
}
183
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
