[NOI2016]bzoj 4652 循环之美 - 数论 - 杜教筛

27 篇文章 0 订阅
5 篇文章 0 订阅

首先题目是在求:

ans=i=1nj=1m[ij][jk] a n s = ∑ i = 1 n ∑ j = 1 m [ i ⊥ j ] [ j ⊥ k ]

注意到j和k互质来得简单,因为k只有2000并且是常数.
ans=j=1m[jk]i=1n[ij] a n s = ∑ j = 1 m [ j ⊥ k ] ∑ i = 1 n [ i ⊥ j ]

ans=j=1m[jk]i=1nd|i,d|jμ(d) a n s = ∑ j = 1 m [ j ⊥ k ] ∑ i = 1 n ∑ d | i , d | j μ ( d )

ans=d=1min(n,m)μ(d)ndd|j,jm[jk] a n s = ∑ d = 1 m i n ( n , m ) μ ( d ) ⌊ n d ⌋ ∑ d | j , j ≤ m [ j ⊥ k ]

注意到若 j=td j = t d 那么 [jk]=[tdk]=[dk][tk] [ j ⊥ k ] = [ t d ⊥ k ] = [ d ⊥ k ] [ t ⊥ k ]
ans=d=1nμ(d)[dk]ndj=1md[jk] a n s = ∑ d = 1 n μ ( d ) [ d ⊥ k ] ⌊ n d ⌋ ∑ j = 1 ⌊ m d ⌋ [ j ⊥ k ]

=d=1nμ(d)[dk]ndf(md) = ∑ d = 1 n μ ( d ) [ d ⊥ k ] ⌊ n d ⌋ f ( ⌊ m d ⌋ )

其中 f(n)=ni=1[ik]=nkϕ(k)+f(n  mod  k) f ( n ) = ∑ i = 1 n [ i ⊥ k ] = ⌊ n k ⌋ ϕ ( k ) + f ( n     m o d     k )
因此 f(n) f ( n ) 可以在 O(klgk)O(1) O ( k l g k ) − O ( 1 ) 内解决
因此对最后的式子做数论分块,那么我们需要求:
g(n,k)=nd=1μ(d)[dk] g ( n , k ) = ∑ d = 1 n μ ( d ) [ d ⊥ k ]
这个怎么做呢?设 n=pcq n = p c q ,其中p是n的某个质因子,并且 pq p ⊥ q
那么考虑 g(n,q) g ( n , q ) 中的所有数字,显然除了答案的部分还有和q互质但是不和p互质的部分要减去;
记这些中的数字为 y=tp,tnp,tq y = t p , t ≤ ⌊ n p ⌋ , t ⊥ q (显然因为p是质数所以不和p互质意味着是p的倍数)
考虑我们求的是这些 y y μ和,那么如果t不和p互质的话 μ(y)=0 μ ( y ) = 0 ,因此可以认为t也和p互质,因此t就和k互质。
g(n,k)=g(n,q)t=1np[tk]μ(tp) g ( n , k ) = g ( n , q ) − ∑ t = 1 ⌊ n p ⌋ [ t ⊥ k ] μ ( t p )

注意到莫比乌斯函数是积性函数,而 tp t ⊥ p ,因此 μ(tp)=μ(t)μ(p) μ ( t p ) = μ ( t ) μ ( p )
g(n,k)=g(n,q)μ(p)t=1np[tk]μ(t) g ( n , k ) = g ( n , q ) − μ ( p ) ∑ t = 1 ⌊ n p ⌋ [ t ⊥ k ] μ ( t )

注意到p是质数,所以 μ(p)=1 μ ( p ) = − 1
g(n,k)=g(n,q)+t=1np[tk]μ(t)=g(n,q)+g(np,k) g ( n , k ) = g ( n , q ) + ∑ t = 1 ⌊ n p ⌋ [ t ⊥ k ] μ ( t ) = g ( n , q ) + g ( ⌊ n p ⌋ , k )

按照固定的顺序除以p,那么有效的状态只有 d(k)2 d ( k ) 2 ,其中 d(k) d ( k ) 表示k的因子个数,实现可以记忆化,因此只传k此时的质因子个数,会好些一些。
边界条件,当 n=0 n = 0 的时候 g(n,k)=0 g ( n , k ) = 0 ,当k=1的时候需要杜教筛。综上总复杂度 O(d(k)2+n23) O ( d ( k ) 2 + n 2 3 ) .

#include<iostream>
#include<cstring>
#include<cstdio>
#include<algorithm>
#include<map>
#define lint long long
#define N 6000000
#define hv(n,c) (n*10000ll+c)
#define debug(x) cerr<<#x<<"="<<x
#define sp <<" "
#define ln <<endl
using namespace std;
int ms[N],m[N],pri[N];bool np[N];
map<int,int> sav;
map<lint,int> gs;
int v[50],f[2010];
inline int prelude(int n)
{
    ms[1]=m[1]=1;
    for(int i=2,c=0;i<=n;i++)
    {
        if(!np[i]) pri[++c]=i,m[i]=-1;ms[i]=ms[i-1]+m[i];
        for(int j=1;j<=c&&(lint)i*pri[j]<=n;j++)
        {
            int x=i*pri[j];np[x]=true;
            if(i%pri[j]) m[x]=-m[i];
            else { m[x]=0;break; }
        }
    }
    return 0;
}
inline int S(int n)
{
    if(n<N) return ms[n];int ans=0ll;
    if(sav.count(n)) return sav[n];
    for(int s=2,t;s<=n;s=++t)
        t=n/(n/s),ans+=(t-s+1)*S(n/t);
    return sav[n]=1-ans;
}
inline int g(int n,int c)
{
    lint nc=hv(n,c);if(gs.count(nc)) return gs[nc];
    if(!n) return 0;if(!c) return gs[nc]=S(n);
    return gs[nc]=g(n,c-1)+g(n/v[c],c);
}
inline int g(int a,int b,int c) { return g(b,c)-g(a-1,c); }
inline int F(int n,int k) { return n<=k?f[n]:(n/k)*f[k]+f[n%k]; }
inline int gcd(int a,int b) { return a?gcd(b%a,a):b; }
int main()
{
    int n,m,k,c=0;scanf("%d%d%d",&n,&m,&k);
    prelude(N-1);lint ans=0ll;int ks=k;
    for(int i=2;i<=ks;i++) if(ks%i==0)
    {   v[++c]=i;while(ks%i==0) ks/=i;  }
    for(int i=1;i<=k;i++) f[i]=f[i-1]+(gcd(i,k)==1);
    for(int s=1,t;s<=min(n,m);s=++t)
        t=min(n/(n/s),m/(m/s)),ans+=(lint)F(m/s,k)*(n/t)*g(s,t,c);
    return !printf("%lld\n",ans);
}

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值