非常恶心的一道数学题,推式子推到吐血。
光是(gcd)求和我还是会的,但是多了个(ij)是什么鬼东西。
很套路的把后面的(d)提出来:
利用(sum_{d|n}mu(d)=[n=1])替换掉([gcd(i,j)=1])就有:
然后第一个难点来了,我们把(x)弄到前面去枚举,然后利用和式的变换,枚举(i,j)是(x)的多少倍,然后用分配率凑在一起就有:
PS:如果上面理解了那么下面就不难了,我们设(T=dx)带进去就有:
后面那个东西可以直接提出来就得到:
把(T^2)提到前面去:
后面那个(sum_{d|T}dmu(frac{T}{d}))显然就是狄利克雷卷积的形式啊,这就等于((idastmu)(T))
然后回忆一下,((idastmu)(T))不就是(phi)么,因此带进去就有:
把(sum_{i=1}^{lfloorfrac{n}{T} floor}i)用等差数列求和的公式展开就变成:
如果我们对(lfloorfrac{n}{T} floor)做除法分块,那么现在要求的就是(T^2phi(T))的前缀和了
考虑用杜教筛的方法化式子,我们设(f(i)=i^2phi(i)),这个时候我们要找一个(g)和(f)卷起来可以方便求。
由于(i^2)很烦,我们考虑先把它消去。令(g(i)=i^2),则有:
(h(n)=(fast g)(n)=sum_{d|n} d^2phi(d)(frac{n}{d})^2=n^2ast(phiast epsilon)(n)=n^3)
所以用杜教筛的套路式就是(g(1)f(n)=h(n)-sum_{d=2}^ng(d)f(lfloorfrac{n}{d} floor))
即(f(n)=sum_{i=1}^ni^3-sum_{d=2}^ng(d)f(lfloorfrac{n}{d} floor)),我们知道(sum_{i=1}^ni^3=(frac{n(n+1)}{2})^2,sum_{i=1}^ni^2=frac{n(n+1)(2n+1)}{6}),所以都可以(O(1))求啦。
CODE
#include<cstdio>
#include<map>
#define RI register int
#define RL register LL
using namespace std;
typedef long long LL;
const int P=10000000;
map <LL,int> sum_f; LL n; int prime[P+5],phi[P+5],cnt,f[P+5],mod,inv2,inv6,ans; bool vis[P+5];
inline void inc(int &x,int y)
{
if ((x+=y)>=mod) x-=mod;
}
inline void dec(int &x,int y)
{
if ((x-=y)<0) x+=mod;
}
inline int quick_pow(int x,int p,int mul=1)
{
for (;p;p>>=1,x=1LL*x*x%mod) if (p&1) mul=1LL*mul*x%mod; return mul;
}
inline int sum(int x)
{
return 1LL*x*(x+1)%mod*inv2%mod;
}
inline int sqr(int l,int r)
{
int t=1LL*r*(r+1)%mod*((r<<1LL)+1)%mod*inv6%mod;
dec(t,1LL*l*(l-1)%mod*((l<<1LL)-1)%mod*inv6%mod); return t;
}
#define Pi prime[j]
inline void Euler(void)
{
vis[1]=phi[1]=1; RI i,j; for (i=2;i<=P;++i)
{
if (!vis[i]) prime[++cnt]=i,phi[i]=i-1;
for (j=1;j<=cnt&&i*Pi<=P;++j)
{
vis[i*Pi]=1; if (i%Pi) phi[i*Pi]=phi[i]*(Pi-1);
else { phi[i*Pi]=phi[i]*Pi; break; }
}
}
for (i=1;i<=P;++i) f[i]=f[i-1],inc(f[i],1LL*phi[i]*i%mod*i%mod);
}
#undef Pi
inline int F(LL x)
{
if (x<=P) return f[x]; if (sum_f.count(x)) return sum_f[x];
int ans=sum(x%mod); ans=1LL*ans*ans%mod; for (RL l=2,r;l<=x;l=r+1)
{
r=x/(x/l);
int t=sqr(l%mod,r%mod);
dec(ans,1LL*t*F(x/l)%mod);
}
return sum_f[x]=ans;
}
int main()
{
//freopen("CODE.in","r",stdin); freopen("CODE.out","w",stdout);
scanf("%d%lld",&mod,&n); Euler(); inv2=quick_pow(2,mod-2);
inv6=quick_pow(6,mod-2); for (RL l=1,r;l<=n;l=r+1)
{
r=n/(n/l); int t=sum(n/l%mod); t=1LL*t*t%mod;
int Fs=F(r); dec(Fs,F(l-1)); inc(ans,1LL*t*Fs%mod);
}
return printf("%d",ans),0;
}