题目大意
给你 (n,m)，求
[sum_{i=1}^nsum_{x_1,x_2,ldots,x_m=1}^ioperatorname{lcm}(gcd(i,x_1),gcd(i,x_2),ldots,gcd(i,x_m))
]
对 ({10}^9+7) 取模。
(nmleq {10}^9)
题解
先推一下式子：
[ans=sum_{i=1}^nsum_{x_1,x_2,ldots,x_m=1}^ioperatorname{lcm}(gcd(i,x_1),gcd(i,x_2),ldots,gcd(i,x_m))\
=sum_{i=1}^nsum_{x_1,x_2,ldots,x_m=1}^ifrac{i}{gcd(frac{i}{gcd(i,x_1)},frac{i}{gcd(i,x_2)},ldots,frac{i}{gcd(i,x_m)})}\
]
然后枚举 (y_j=frac{i}{gcd(i,x_j)})。对于一个 (y_j)，有多少个满足要求的 (x_j) 呢？
[frac{i}{gcd(i,x_j)}=y_j\
gcd(i,x_j)=frac{i}{y_j}\
gcd(y_j,frac{x_j}{y_j})=1\
]
所以有 (varphi(y_j)) 个。
[ans=sum_{i=1}^nsum_{y_1,y_2,ldots,y_mmid i}varphi(y_1)varphi(y_2)cdotsvarphi(y_m)frac{i}{gcd(y_1,y_2,ldots,y_m)}
]
记
[f(n)=sum_{y_1,y_2,ldots,y_mmid n}varphi(y_1)varphi(y_2)cdotsvarphi(y_m)frac{n}{gcd(y_1,y_2,ldots,y_m)}
]
容易发现 (f(n)) 是积性函数。
[f(p^k)=sum_{i=0}^kp^{k-i}sum_{j=i}^kmu(frac{j}{i})g(j)\
=sum_{j=0}^kg(j)sum_{i=0}^jp^{k-i}mu(frac{j}{i})\
g(j)=egin{cases}
p^{km}&,j=1\
(p^k-p^{j-1})^m&,2leq jleq k
end{cases}\
]
[f(p^k)=sum_{x_1,x_2,ldots,x_m=1}^{p^k}lcm(gcd(p^k,x_1),gcd(p^k,x_2),ldots,gcd(p^k,x_m))\
=sum_{x_1,x_2,ldots,x_m=0}^{k}p^{max(x_1,x_2,ldots,x_m)}prod_{j=1}^mvarphi(p^{k-x_j})\
=sum_{i=0}^kp^i(g(i)-g(i-1))\
g(i)=sum_{x_1,x_2,dots,x_m=0}^iprod_{j=1}^mvarphi(p^{k-x_j})\
=prod_{j=1}^m(sum_{l=0}^ivarphi(p^{k-l}))\=
egin{cases}
p^{km}&,i=k\
(p^k-p^{k-1-i})^m&,0leq i<k
end{cases}\
f(p^{k+1})=p^mf(p^k)-p^k(p^{(k+1)m}-(p^{k+1}-1)^m)+p^{k+1}(p^{(k+1)m}-(p^{k+1}-1)^m)\
=p^mf(p^k)+p^k(p-1)(p^{(k+1)m}-(p^{k+1}-1)^m)
]
当 (nleq {10}^7) 的时候可以线性筛算出 (f(1)sim f(n))。
时间复杂度：(O(n(1+frac{log m}{log n})))
当 (n) 大的时候：
[f(p)=g(0)+pg(1)-pg(0)\
=(1-p)(p-1)^m+p^{m+1}\
=p^{m+1}-(p-1)^{m+1}\
=-sum_{i=0}^m{(-1)}^{m+1-i}p^iinom{m+1}{i}\
=sum_{i=0}^m{(-1)}^{m-i}p^iinom{m+1}{i}
]
就可以上 min_25 筛了。
时间复杂度：(O(frac{n^{frac{3}{4}}}{log n}m+sqrt nm^2))
代码
这个代码有些东西没有预处理，时间复杂度可能比上面写的高一些。
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cstdlib>
#include<ctime>
#include<functional>
#include<cmath>
#include<vector>
#include<assert.h>
//using namespace std;
using std::min;
using std::max;
using std::swap;
using std::sort;
using std::reverse;
using std::random_shuffle;
using std::lower_bound;
using std::upper_bound;
using std::unique;
using std::vector;
typedef long long ll;
typedef unsigned long long ull;
typedef double db;
typedef std::pair<int,int> pii;
typedef std::pair<ll,ll> pll;
void open(const char *s){
#ifndef ONLINE_JUDGE
char str[100];sprintf(str,"%s.in",s);freopen(str,"r",stdin);sprintf(str,"%s.out",s);freopen(str,"w",stdout);
#endif
}
void open2(const char *s){
#ifdef DEBUG
char str[100];sprintf(str,"%s.in",s);freopen(str,"r",stdin);sprintf(str,"%s.out",s);freopen(str,"w",stdout);
#endif
}
int rd(){int s=0,c,b=0;while(((c=getchar())<'0'||c>'9')&&c!='-');if(c=='-'){c=getchar();b=1;}do{s=s*10+c-'0';}while((c=getchar())>='0'&&c<='9');return b?-s:s;}
void put(int x){if(!x){putchar('0');return;}static int c[20];int t=0;while(x){c[++t]=x%10;x/=10;}while(t)putchar(c[t--]+'0');}
int upmin(int &a,int b){if(b<a){a=b;return 1;}return 0;}
int upmax(int &a,int b){if(b>a){a=b;return 1;}return 0;}
const ll p=1000000007;
ll fp(ll a,ll b)
{
ll s=1;
for(;b;b>>=1,a=a*a%p)
if(b&1)
s=s*a%p;
return s;
}
int n,m;
namespace gao1
{
const int N=10000010;
int f[N];
int pw[N];
bool b[N];
int pri[N];
int c[N];
int cnt;
void gao()
{
f[1]=1;
pw[1]=1;
for(int i=2;i<=n;i++)
{
if(!b[i])
{
pri[++cnt]=i;
pw[i]=fp(i,m);
}
for(int j=1;j<=cnt&&i*pri[j]<=n;j++)
{
b[i*pri[j]]=1;
pw[i*pri[j]]=(ll)pw[i]*pw[pri[j]]%p;
if(i%pri[j]==0)
break;
}
}
for(int i=2;i<=n;i++)
{
if(!b[i])
{
c[i]=i;
f[i]=fp(i,m+1)-fp(i-1,m+1);
}
for(int j=1;j<=cnt&&pri[j]<=i&&i*pri[j]<=n;j++)
{
int v=i*pri[j];
if(i%pri[j]==0)
{
c[v]=c[i]*pri[j];
if(c[v]==v)
f[v]=((ll)pw[pri[j]]*f[i]+(ll)i*(pri[j]-1)*(pw[v]-pw[v-1]))%p;
else
f[v]=(ll)f[c[v]]*f[v/c[v]]%p;
break;
}
c[v]=pri[j];
f[v]=(ll)f[i]*f[pri[j]]%p;
}
}
ll ans=0;
for(int i=1;i<=n;i++)
ans+=f[i];
ans=(ans%p+p)%p;
printf("%lld
",ans);
}
}
namespace gao2
{
const int N=100010;
ll inv[N],fac[N],ifac[N];
ll binom(int x,int y)
{
return x>=y&&y>=0?fac[x]*ifac[y]%p*ifac[x-y]%p:0;
}
int _;
int b[N],pri[N],cnt;
void init()
{
inv[1]=fac[0]=fac[1]=ifac[0]=ifac[1]=1;
for(int i=2;i<=100000;i++)
{
inv[i]=-p/i*inv[p%i]%p;
fac[i]=fac[i-1]*i%p;
ifac[i]=ifac[i-1]*inv[i]%p;
}
for(int i=2;i<=_;i++)
{
if(!b[i])
pri[++cnt]=i;
for(int j=1;j<=cnt&&i*pri[j]<=_;j++)
{
b[i*pri[j]]=1;
if(i%pri[j]==0)
break;
}
}
}
ll pw[N],s[N];
void init(int x)
{
pw[1]=1;
for(int i=2;i<=x+2;i++)
{
if(!b[i])
pw[i]=fp(i,x);
for(int j=1;j<=cnt&&i*pri[j]<=n&&pri[j]<=i;j++)
{
pw[i*pri[j]]=pw[i]*pw[pri[j]]%p;
if(i%pri[j]==0)
break;
}
}
for(int i=1;i<=x+2;i++)
s[i]=(s[i-1]+pw[i])%p;
}
ll pre[N],suf[N];
ll calc(int x,int y)
{
if(x<=y+2)
return s[x];
pre[0]=1;
for(int i=1;i<=y+2;i++)
pre[i]=pre[i-1]*(x-i)%p;
suf[y+3]=1;
for(int i=y+2;i>=1;i--)
suf[i]=suf[i+1]*(x-i)%p;
ll res=0;
for(int i=1;i<=y+2;i++)
res=(res+s[i]*pre[i-1]%p*suf[i+1]%p*ifac[i-1]%p*ifac[y+2-i]%p*((y+2-i)&1?-1:1))%p;
return res;
}
ll s1[N],s2[N],ans1[N],ans2[N];
void gao(int x)
{
init(x);
calc(3,x);
for(int i=1;i<=_;i++)
s1[i]=(fp(i,x)+s1[i-1])%p;
for(int i=1;n/i>_;i++)
s2[i]=calc(n/i,x);
for(int i=1;i<=cnt;i++)
{
ll v=fp(pri[i],x);
int j;
for(j=1;n/j/pri[i]>_&&n/j>=(ll)pri[i]*pri[i];j++)
s2[j]=(s2[j]-v*(s2[j*pri[i]]-s1[pri[i]-1]))%p;
for(;n/j>_&&n/j>=(ll)pri[i]*pri[i];j++)
s2[j]=(s2[j]-v*(s1[n/j/pri[i]]-s1[pri[i]-1]))%p;
for(j=_;j>=(ll)pri[i]*pri[i];j--)
s1[j]=(s1[j]-v*(s1[j/pri[i]]-s1[pri[i]-1]))%p;
}
}
ll g(int z,int x,int y)
{
if(z<0)
return 0;
return z==y?fp(x,(ll)y*m):fp(fp(x,y)-fp(x,y-z-1),m);
}
ll get(int x,int y)
{
ll res=0;
for(int i=0;i<=y;i++)
res=(res+fp(x,i)*(g(i,x,y)-g(i-1,x,y)))%p;
return res;
}
ll solve(ll x,int y)
{
if(x<=1||x<pri[y])
return 0;
if(y>cnt)
return (x<=_?ans1[x]:ans2[n/x])-ans1[_];
ll res=(x<=_?ans1[x]:ans2[n/x])-ans1[pri[y]-1];
for(int i=y;i<=cnt&&(ll)pri[i]*pri[i]<=x;i++)
{
ll x1=x/pri[i];
for(int j=1;x1>=pri[i];j++,x1/=pri[i])
res=(res+solve(x1,i+1)*get(pri[i],j)+get(pri[i],j+1))%p;
}
return res;
}
void gao()
{
_=sqrt(n)+0.5;
init();
for(int i=0;i<=m;i++)
{
gao(i);
ll v=binom(m+1,i)*((m-i)&1?-1:1);
for(int i=1;i<=_;i++)
ans1[i]=(ans1[i]+v*s1[i])%p;
for(int i=1;n/i>_;i++)
ans2[i]=(ans2[i]+v*s2[i])%p;
}
ll ans=solve(n,1);
ans++;
ans=(ans%p+p)%p;
printf("%lld
",ans);
}
}
int main()
{
open("c");
scanf("%d%d",&n,&m);
if(n<=1e7)
gao1::gao();
else
gao2::gao();
return 0;
}