bzoj2154和luogu3768
内容
推荐先看完这个
bzoj2154
题意
求
$$\sum_{i=1}^n\sum_{j=1}^mlcm(i,j)\ \ n,m \leq 10^7$$
题解
来治治公式恐惧症吧：
我们设$n < m$
显然：
$$\sum_{i=1}^n\sum_{j=1}^mlcm(i,j)=\sum_{i=1}^n\sum_{j=1}^m{ij\over gcd(i,j)}$$
我们设$gcd(i,j)=d$
现在想要枚举$d$，就确定$ij$如何取。那么就可以先把$i,j$中的$d$因子去掉就可以了。也就是$gcd(i/d,j/d)=1$。
设
$$F(x, y) = \sum_{i=1}^{x} \sum_{j=1}^{y} ij[gcd(i,j)=1]$$
那么原式中还需要乘上$d^2$，因为除掉了两个：
$$原式=\sum_{d=1}^{n} \frac{d^2 F( \frac{n}{d} , \frac{m}{d} )}{d} = \sum_{d=1}^{n} d F(\frac{n}{d} , \frac{m}{d} )$$
考虑化简一下
$$F(x, y)= \sum_{i=1}^{x} \sum_{j=1}^{y} ij[gcd(i,j)=1]$$
反演一下
$$= \sum_{i=1}^{x} \sum_{j=1}^{y} ij \sum_{d|(i,j)} \mu (d) $$
改成枚举$d$，计算一下贡献
$$= \sum_{d=1}^{x} \mu (d) \sum_{d|i}^{x} i \sum_{d|j}^{y} j$$
换个元
$$= \sum_{d=1}^{x} \mu (d) d^2 \sum_{i=1}^{ \frac{x}{d}} i \sum_{j=1}^{\frac{y}{d}} j$$
等差数列公式
$$= \sum_{d=1}^{x} \mu (d) d^2 \frac{\frac{x}{d} ( \frac{x}{d}+1)}{2} \frac{\frac{y}{d} (\frac{y}{d} +1)}{2} $$
带回原式
$$原式=\sum_{d=1}^{n} d F( \frac{n}{d} , \frac{m}{d} )$$
带入
$$\large = \sum_{d=1}^{n} d \sum_{i=1}^{ \frac{n}{d} } \mu (i) i^2 \frac{ \frac{ \frac{n}{d} }{i} ( \frac{ \frac{n}{d} }{i} +1)}{2} \frac{ \frac{ \frac{m}{d} }{i} ( \frac{ \frac{m}{d} }{i} +1)}{2} $$
把$d$乘下来
$$= \sum_{d=1}^{n} d \sum_{i=1}^{ \frac{n}{d} } \mu (i) i^2 \frac{ \frac{n}{di} ( \frac{n}{di} +1)}{2} \frac{ \frac{m}{di} ( \frac{m}{di} +1)}{2}$$
现在已经$O(\sqrt n \times \sqrt{n})=O(n)$的查询了。以及足够通过了。
但是我们还想优化：
令$T=di$，则$i|T,d=T/i$：
$$\sum_{d=1}^{n} d \sum_{i=1}^{ \frac{n}{d} } \mu (i) i^2 \frac{ \frac{n}{di} ( \frac{n}{di} +1)}{2} \frac{ \frac{m}{di} ( \frac{m}{di} +1)}{2} $$
换个元
$$= \sum_{T=1}^{n} \frac{ \frac{n}{T} ( \frac{n}{T} +1)}{2} \frac{ \frac{m}{T} ( \frac{m}{T} +1)}{2} \sum_{i|T} \frac{T}{i} \mu (i) i^2 $$
把$T$带出来
$$= \sum_{T=1}^{n} \frac{ \frac{n}{T} ( \frac{n}{T} +1)}{2} \frac{ \frac{m}{T} ( \frac{m}{T} +1)}{2} T \sum_{i|T} \mu (i) i $$
设$g(T)=T \sum_{i|T} \mu (i) i$，我们再设$f(T)=\sum_{i|T} \mu (i) i$那么$g(T)=Tf(T)$，考虑求$f(T)$
在线性筛中，外层为$k$，内层为$p$，所以求$f(kp)=\sum_{i|kp} \mu(i) i$
当$p|k$时
当$i$取的数的因子中不包含新加入的$p$时，答案就是$f(k)$
当$i$取包含新加入的因子$p$时，由于此时$p$指数已经$>=2$，所以$\mu (i)=0$，因此贡献为0
综上，当$p|k$时，答案为$f(k)$
当$p \nmid k$时当
$i$取的数的因子中不包含新加入的$p$时，同上，答案是$f(k)$
当$i$取的数的因子包含新加入的$p$时，由于指数为$1$，所以我们考虑$i=ap$，原式变为
$$\sum_{i|T} \mu(i) i$$
$$= \sum_{ap|kp} \mu (ap) ap $$
$$= p \sum_{a|k} \mu (a) \mu(p) a$$
$$= -p \sum_{a|k} \mu(a) a $$
$$= -p f(k) $$
综上，当$p_y \nmid k$时，答案为$(1-p)f(k)$
然后在线性筛的时候乱搞就可以了
#include <cstdio>
#define mod 20101009
using namespace std;
typedef long long ll;
const int N=1e7+10;
int prime[N],tot,n,m;
ll g[N];
bool not_prime[N];
int min(int x,int y){if (x>y) return y;return x;}
int max(int x,int y){if (x<y) return y;return x;}
void swap(int &x,int &y){int t=x;x=y;y=t;}
void init(int n)
{
g[1]=1;
for (int i=2;i<=n;i++)
{
if (!not_prime[i])
prime[++tot]=i,g[i]=1-i;
for (int j=1;j<=tot;j++)
{
int temp=prime[j]*i;
if (temp>n)
break;
not_prime[temp]=1;
if (i%prime[j]==0)
{
g[temp]=g[i];
break;
}
g[temp]=g[i]*(1-prime[j]);
}
}
for (int i=2;i<=n;i++)
g[i]*=i;
for (int i=1;i<=n;i++)
g[i]+=g[i-1],g[i]%=mod;
}
main()
{
//init(N-10);
ll ans=0;
scanf("%d%d",&n,&m);
if (n>m)
swap(n,m);
init(n);
for (int i=1,last;i<=n;i=last+1)
{
last=min(n/(n/i),m/(m/i));
ll t1=((ll)(n/i)*(n/i+1)/2)%mod;
ll t2=((ll)(m/i)*(m/i+1)/2)%mod;
ans=(ans+((g[last]-g[i-1])*((t1*t2)%mod))%mod)%mod;
}
printf("%lld\n",((ans%mod)+mod)%mod);
}
luogu3768
题意
求
$$\sum_{i=1}^n\sum_{j=1}^nijgcd(i,j) \bmod p$$
其中$n\leq 10^{10},p \in 质数$
题解
这个题和上个题惊人的类似啊
反演吧
我们设$gcd(i,j)=d$
现在想要枚举$d$，就确定$ij$如何取。那么就可以先把$i,j$中的$d$因子去掉就可以了。也就是$gcd(i/d,j/d)=1$。
设
$$F(x, y) = \sum_{i=1}^{x} \sum_{j=1}^{y} ij[gcd(i,j)=1]$$
那么原式中还需要乘上$d^2$，因为除掉了两个：
$$原式=\sum_{d=1}^{n} d^3 F( \frac{n}{d} , \frac{n}{d} )$$
考虑化简一下
$$F(x, y)= \sum_{i=1}^{x} \sum_{j=1}^{y} ij[gcd(i,j)=1]$$
反演一下
$$= \sum_{i=1}^{x} \sum_{j=1}^{y} ij \sum_{d|(i,j)} \mu (d) $$
改成枚举$d$，计算一下贡献
$$= \sum_{d=1}^{x} \mu (d) \sum_{d|i}^{x} i \sum_{d|j}^{y} j$$
换个元
$$= \sum_{d=1}^{x} \mu (d) d^2 \sum_{i=1}^{ \frac{x}{d}} i \sum_{j=1}^{\frac{y}{d}} j$$
等差数列公式
$$= \sum_{d=1}^{x} \mu (d) d^2 \frac{\frac{x}{d} ( \frac{x}{d}+1)}{2} \frac{\frac{y}{d} (\frac{y}{d} +1)}{2} $$
带回原式
$$原式=\sum_{d=1}^{n} d^3 F( \frac{n}{d} , \frac{n}{d} )$$
带入
$$\large = \sum_{d=1}^{n} d^3 \sum_{i=1}^{ \frac{n}{d} } \mu (i) i^2 (\frac{ \frac{ \frac{n}{d} }{i} ( \frac{ \frac{n}{d} }{i} +1)}{2} )^@ $$
把$d$乘下来
$$= \sum_{d=1}^{n} d^3 \sum_{i=1}^{ \frac{n}{d} } \mu (i) i^2( \frac{ \frac{n}{di} ( \frac{n}{di} +1)}{2})^2$$
令$T=di$，则$i|T,d=T/i$，换个元
$$= \sum_{T=1}^{n} (\frac{ \frac{n}{T} ( \frac{n}{T} +1)}{2})^2 \sum_{i|T} (\frac{T}{i})^3 \mu (i) i^2 $$
把$T$带出来
$$= \sum_{T=1}^{n} (\frac{ \frac{n}{T} ( \frac{n}{T} +1)}{2})^2 T^2 \sum_{i|T} \mu (i) \frac{T}{i} $$
所以$T \over i$也是一个因子，都是等价的，所以和下面是一样的
$$= \sum_{T=1}^{n} (\frac{ \frac{n}{T} ( \frac{n}{T} +1)}{2})^2 T^2 \sum_{i|T} \mu (\frac{T}{i}) i $$
我们现在就要研究这个函数
$$f(x)=x^2\sum_{i|x}\mu({x \over i})i$$
处理出它的前缀和就可以解决这个题了
我们再来看一下它的内函数
$$g(x)=\sum_{i|x}\mu({x \over i})i =\sum_{i|x}\mu(i){x \over i}$$
再用线性筛的思路试试
在线性筛中，外层为$k$，内层为$p$，所以求
$$g(kp)=\sum_{i|kp} \mu({kp \over i})i=\sum_{i|kp} \mu(i){kp \over i}$$
当$p|k$时
当$i$取的数的因子中不包含新加入的$p$时，答案就是$f(k)$
当$i$取包含新加入的因子$p$时，设$i=ap$
$$=\mu({kp \over ap})ap=\mu({k \over a})ap$$
也就是多了个$p$
综上，当$p|k$时，答案为$g(k)p$
当$p \nmid k$时当
$i$取的数的因子中不包含新加入的$p$时，同上，答案是$g(k)$
当$i$取的数的因子包含新加入的$p$时，由于指数为$1$，所以我们考虑$i=ap$，原式变为
$$\sum_{i|x} \mu({x \over i})i$$
$$= \sum_{ap|kp} \mu(ap){kp \over ap} $$
$$= \sum_{a|k} \mu (a) \mu(p) {k \over a}$$
$$= -\sum_{a|k} \mu(a) a $$
$$= – g(k) $$
综上，当$p_y \nmid k$时，答案为$(p-1)g(k)$
筛的时候发现这个函数竟然是$\phi$！
这个方法非常玄学。我们用更加正经的方法证明一下。
令$h(i)=i$，那么
$$g(x)=(h*\mu)(x)$$
我们给$f(x)$卷一个1试试
$$h*\mu*1=h*(\mu*1)=h*\epsilon=h$$
也就是$(h*\mu)*1=h$
又因为$\phi *1=h$
所以$g(x)=\phi(x),f(x)=x^2\phi(x)$
现在就是要用杜教筛来解决这个函数的前缀和。
令$c(x)=x^2$
那么：
$$(f*c)(x)$$
$$=\sum_{d|x}f(x)c({x\over d})$$
$$=\sum_{d|x}d^2\phi{x^2 \over d^2}$$
$$=x^2\sum_{d|x}\phi(x)=x^3$$
这样就可以杜教筛了。。。
令$s(n)=\sum_{i=1}^nf(x)$
那么
$$c(1)s(n)=\sum_{i=1}^n(f*c)(i)-\sum_{i=2}^nc(i)s({n \over i})$$
$$s(n)=\sum_{i=1}^ni^3-\sum_{i=2}^ni^2s({n\over i})$$
然后就可以做了
代码
#include <cstdio>
#include <algorithm>
#include <map>
#define maxn 4700000
#define N 4700050
using namespace std;
typedef long long ll;
ll p,n;
ll inv2,inv6;
ll phi[N],prime[N],tot;
bool not_prime[N];
map <ll,ll> ma;
ll pow(ll x,ll y,ll p)
{
ll ans=1;
for (;y;y>>=1,x=x*x%p)
if (y&1)
ans=ans*x%p;
return ans;
}
ll sqr(ll x) {x%=p;return x*x%p;}
ll g(ll x){x%=p;return sqr((x*(1+x)>>1)%p);}
ll get_h(ll x){x%=p;return x*(x+1)%p*(2*x+1)%p*inv6%p;}
ll min(ll a,ll b){if (a<b) return a;return b;}
void init(ll n)
{
phi[1]=1;
for (ll i=2;i<=n;i++)
{
if (!not_prime[i])
prime[++tot]=i,phi[i]=i-1;
for (ll j=1;j<=tot;j++)
{
ll temp=prime[j]*i;
if (temp>n)
break;
not_prime[temp]=1;
if (i%prime[j]==0)
{
phi[temp]=(phi[i]*prime[j])%p;
break;
}
phi[temp]=(phi[i]*(prime[j]-1))%p;
}
}
for(ll i=1;i<=n;i++)
phi[i]=(phi[i]*sqr(i)%p+phi[i-1])%p;
inv2=pow(2,p-2,p),inv6=pow(6,p-2,p);
}
ll get_f(ll n)
{
if (n<=maxn) return phi[n];
if (ma.count(n)) return ma[n];
ll temp=g(n);
for (ll i=2,last;i<=n;i=last+1)
{
last=n/(n/i);
temp=(temp-(get_h(last)-get_h(i-1))*get_f(n/i)%p)%p;
}
return ma[n]=temp;
}
main()
{
scanf("%lld%lld",&p,&n);
init(min(maxn,n));
ll ans=0;
for (ll i=1,last;i<=n;i=last+1)
{
last=n/(n/i);
ans=(ans+g(n/i)*((get_f(last)-get_f(i-1))%p)%p)%p;
}
printf("%lld\n",(ans+p)%p);
}