题目
设g(x)g(x)g(x)为xxx的可重质因子数目,f(x)=2g(x)f(x)=2^{g(x)}f(x)=2g(x)。
求∑i=1nf(i)\sum_{i=1}^nf(i)∑i=1nf(i)
正解
奇怪的数论知识增加了。
介绍一个叫power number的东西。power number是所有的质因子的指数都大于等于222的数。
每个power number都可以表示成a2b3a^2b^3a2b3的形式,其中a,b∈Z+a,b\in Z^+a,b∈Z+
power number的数量是比较少的,所以有没有什么办法将值挂在power number上优化求值呢?
构造数论函数uuu满足u∗I∗I=fu*I*I=fu∗I∗I=f,即u=f∗μ∗μu=f*\mu*\muu=f∗μ∗μ。容易得知uuu也是个积性函数,并且手玩可以玩出u(pk)u(p^k)u(pk)(ppp为质数)的表达式。然后就可以发现uuu只有在power number处有值。
接下来推一下式子:∑i=1nf(i)=∑i=1n∑ju(j)∑a,b[jab=i]=∑ju(j)∑a,b[jab≤n]=∑ju(j)∑a,b[ab≤nj]\sum_{i=1}^nf(i)=\sum_{i=1}^n\sum_j u(j)\sum_{a,b} [jab=i]=\sum_j u(j)\sum_{a,b} [jab\le n]=\sum_j u(j)\sum_{a,b} [ab\le \frac{n}{j}]∑i=1nf(i)=∑i=1n∑ju(j)∑a,b[jab=i]=∑ju(j)∑a,b[jab≤n]=∑ju(j)∑a,b[ab≤jn]
记d(n)=∑a,b[ab≤n]d(n)=\sum_{a,b}[ab\le n]d(n)=∑a,b[ab≤n]。
这是一个经典问题。朴素求法是变成∑i⌊ni⌋\sum_i \lfloor \frac{n}{i}\rfloor∑i⌊in⌋,整除分块。
然后发现时间卡不过去……
有个小优化:画出函数图像xy=nxy=nxy=n,我们要求函数图像和xxx正半轴和yyy正半轴之间的整点个数。这个东西可以看成一个正方形和两个斜边是弯的三角形。正方形的整点可以O(1)O(1)O(1)算,然后求弯的三角形内的整点就可以直接枚举每一行,这里行数是O(n)O(\sqrt n)O(n)的。
形式化地说,就是2∑i=1⌊n⌋⌊ni⌋+⌊n⌋22\sum_{i=1}^{\lfloor \sqrt n \rfloor} \lfloor \frac{n}{i} \rfloor+\lfloor \sqrt n \rfloor ^22∑i=1⌊n⌋⌊in⌋+⌊n⌋2
左边的这个东西直接暴力算(不要整除分块!!!)。
(当然还有个更厉害的做法是在SB树上二分来拟合凸包,时间复杂度理论上可以达到O(n13lgn)O(n^{\frac{1}{3}}\lg n)O(n31lgn))
代码
using namespace std;
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <cassert>
#define N 10000000
#define ll long long
ll n,sq,mo;
int p[N+1],np;
bool inp[N+1];
void init(int n){
for (int i=2;i<=n;++i){
if (!inp[i])
p[++np]=i;
for (int j=1;j<=np && i*p[j]<=n;++j){
inp[i*p[j]]=1;
if (i%p[j]==0)
break;
}
}
}
ll calc(int k){return (1ll<<k-2)%mo;}
int d1[N+1],d2[N+1];
#define d(i) (i<=sq?d1[i]:d2[n/i])
int Sd(ll k){
if (d(k)!=-1) return d(k);
ll sk=sqrt(k),res=0;
for (register int i=1;i<=sk;++i)
(res+=k/i)%=mo;
return d(k)=((res*2-sk*sk)%mo+mo)%mo;;
}
ll ans;
void dfs(int x,ll k,ll uk){
(ans+=uk*Sd(k))%=mo;
for (;x<=np && k>=(ll)p[x]*p[x];++x){
ll _k=k;
k/=(ll)p[x]*p[x];
for (int i=2;k;k/=p[x],++i)
dfs(x+1,k,uk*calc(i)%mo);
k=_k;
}
}
int main(){
freopen("remapping.in","r",stdin);
freopen("remapping.out","w",stdout);
scanf("%lld%d",&n,&mo);
sq=sqrt(n);
for (;(sq+1)*(sq+1)<=n;++sq);
init(sq);
memset(d1,255,sizeof(int)*(sq+1));
memset(d2,255,sizeof(int)*(sq+1));
dfs(1,n,1);
printf("%lld\n",ans);
return 0;
}