SP7001 【VLATTICE - Visible Lattice Points】Solution

题目描述

一共T组数据,每组数据一个整数n,表示有一个大小为 n 3 n^3 n3的正方体。求:在(0,0,0)位置能够看到多少个不被遮挡的点(i,j,k)

特别地,i,j,k均为整数。

思路

不妨从二维的情况开始考虑。事实上,这道题目的确有二维形式。

POJ3090 Visible Lattice Points

不妨这样想,点(i,j)为什么会被看到?是因为它和(0,0)的连线上面没有点了。我们知道这条直线的斜率为 j − 0 i − 0 = j i \frac{j-0}{i-0}=\frac{j}{i} i0j0=ij
这说明j和i的最大公因数为1,i,j互质。

那么这道题就成了对于每一个 s &lt; = n s&lt;=n s<=n,求得所有小于等于s的与s互质的数字的个数,而这就是求欧拉函数
φ ( s ) \varphi(s) φ(s)的值,乘以2,再去除被计算两次的点(1,1).
所以 a n s = 2 ∗ ∑ s = 1 n φ ( s ) − 1 ans=2*\sum_{s=1}^{n}\varphi(s)-1 ans=2s=1nφ(s)1,利用线性筛求得n以内所有欧拉函数的值即可。

升维打击!

多了一维,思路保持不变。对于所有的i,j,k,我们要找出一共有多少种情况,使得 g c d ( i , j , k ) = 1 gcd(i,j,k)=1 gcd(i,j,k)=1.我们可以这样表达这个公式:

  • g ( n ) = ∑ i = 1 n ∑ j = 1 n ∑ k = 1 n p d ( i , j , k ) g(n)=\sum_{i=1}^{n}\sum_{j=1}^{n}\sum_{k=1}^{n}pd(i,j,k) g(n)=i=1nj=1nk=1npd(i,j,k)
  • 其中当且仅当 g c d ( i , j , k ) = 1 gcd(i,j,k)=1 gcd(i,j,k)=1 p d ( i , j , k ) = 1 pd(i,j,k)=1 pd(i,j,k)=1,否则为0.

显然,我们可以利用三重循环暴力枚举ijk,再用 l o g 2 n log_{2}n log2n的时间完成gcd,但是在 n = 1000000 n=1000000 n=1000000多组数据的情况下分分钟就T飞了!于是我们来考虑一下新的算法:

变换大法+数学乱搞

如果现在有一个函数 f ( x ) f(x) f(x),满足如下条件:

  • f ( n ) = ∑ i = 1 N ∑ j = 1 N ∑ k = 1 N p d ( i , j , k ) f(n)=\sum_{i=1}^{N}\sum_{j=1}^{N}\sum_{k=1}^{N}pd(i,j,k) f(n)=i=1Nj=1Nk=1Npd(i,j,k)
  • 其中当且仅当 g c d ( i , j , k ) = n gcd(i,j,k)=n gcd(i,j,k)=n p d ( i , j , k ) = 1 pd(i,j,k)=1 pd(i,j,k)=1,否则为0.

很容易看出这个函数表示在0-N范围内有多少对公因数为n的三数对。
很显然 f ( 1 ) f(1) f(1)就是我们要求的数值。
但是这个函数看起来不好求啊??

请看:


莫比乌斯反演变换

定义:如果有 F ( n ) = ∑ d ∣ n f ( d ) ( d ∣ n F(n)=\sum_{d|n}f(d)(d|n F(n)=dnf(d)(dn表示d整除n ) ) )

那么必然能够得到 f ( n ) = ∑ d ∣ n μ ( n d ) F ( d ) f(n)=\sum_{d|n}μ(\frac{n}{d})F(d) f(n)=dnμ(dn)F(d)

这被称作莫比乌斯反演公式,其中 μ ( ) μ() μ()函数代表莫比乌斯函数。关于这个函数,我们需要知道的是:
如果给出一个合数p的标准分解式: p = p 1 α 1 p 2 α 2 . . . p m α m p=p_1^{α_1}p_2^{α_2}...p_m^{α_m} p=p1α1p2α2...pmαm
那么就有

μ ( p ) = { 1 p = 1 ( − 1 ) m e a c h α = 1 0 o t h e r w i s e μ(p)=\begin{cases} 1 &amp;&amp; {p=1}\\ (-1)^m &amp;&amp; each α=1\\ 0 &amp;&amp;otherwise\end{cases} μ(p)=1(1)m0p=1eachα=1otherwise

这是莫比乌斯函数的定义公式


定义函数,进行反演

经验告诉我们, f ( n ) f(n) f(n)的选取常常与等于有关,而 F ( n ) F(n) F(n)的选取通常与整除有关。
这里,我们定义:

  • f ( n ) = ∑ i = 1 N ∑ j = 1 N ∑ k = 1 N p d ( i , j , k ) f(n)=\sum_{i=1}^{N}\sum_{j=1}^{N}\sum_{k=1}^{N}pd(i,j,k) f(n)=i=1Nj=1Nk=1Npd(i,j,k)
  • 其中当且仅当 g c d ( i , j , k ) = n gcd(i,j,k)=n gcd(i,j,k)=n p d ( i , j , k ) = 1 pd(i,j,k)=1 pd(i,j,k)=1,否则为0.

那么相应地,我们这样定义大F(n)

  • F ( n ) = ∑ i = 1 N ∑ j = 1 N ∑ k = 1 N p d ( i , j , k ) F(n)=\sum_{i=1}^{N}\sum_{j=1}^{N}\sum_{k=1}^{N}pd(i,j,k) F(n)=i=1Nj=1Nk=1Npd(i,j,k)
  • 其中当且仅当${n}\equiv 0\mod{gcd(i,j,k)} 时 时 pd(i,j,k)=1 , 否 则 为 0. 那 么 我 们 知 道 ,否则为0.那么我们知道 ,0.gcd(i,j,k)$必须是n的因数。

好了,很容易推导出: F ( n ) = ∑ d ∣ n f ( d ) F(n)=\sum_{d|n}f(d) F(n)=dnf(d)

莫比乌斯反演变换后得到 f ( n ) = ∑ d ∣ n μ ( n d ) F ( d ) . f(n)=\sum_{d|n}μ(\frac n d)F(d). f(n)=dnμ(dn)F(d).

接下来我们对这个式子做一点处理。因为要求得 f ( 1 ) f(1) f(1)的值,我们令 n = 1 n=1 n=1,得到 f ( 1 ) = ∑ d = 1 N μ ( d ) F ( d ) = ∑ d μ ( d ) ∑ i = 1 N ∑ j = 1 N ∑ k = 1 N p d ( i , j , k ) . f(1)=\sum_{d=1}^Nμ(d)F(d)=\sum_{d}μ(d)\sum_{i=1}^{N}\sum_{j=1}^{N}\sum_{k=1}^{N}pd(i,j,k). f(1)=d=1Nμ(d)F(d)=dμ(d)i=1Nj=1Nk=1Npd(i,j,k).

对右边部分处理。不妨设 i = a d , j = b d , k = c d i=ad,j=bd,k=cd i=ad,j=bd,k=cd,那么d一定是从1到N的。对于每一个d,abc都只能是1-[N/d]中的某一个,因此可得 F ( d ) = [ N d ] 3 F(d)={[\frac N d]}^3 F(d)=[dN]3.如果是二维,同理可得 F ( d ) = [ N d ] 2 F(d)={[\frac N d]}^2 F(d)=[dN]2

考虑到N的因数只有 N \sqrt N N 个,我们可以在线性筛莫比乌斯函数的时候预处理莫比乌斯函数的前缀和。利用整除分块,对 [ N d ] [\frac N d] [dN]相同的d进行组合处理.至此,本题得解。最终的答案应该是立体处理+三个平面的处理(分别是 x = 0 , y = 0 , z = 0 x=0, y=0, z=0 x=0,y=0,z=0)+ ( 0 , 1 , 1 ) , ( 1 , 1 , 0 ) , ( 1 , 0 , 1 ) (0,1,1),(1,1,0),(1,0,1) (0,1,1),(1,1,0),(1,0,1)三个点.

a n s = ∑ d = 1 N μ ( d ) ( [ N d ] 2 ( [ N d ] + 1 ) ) + 3 ans=\sum_{d=1}^Nμ(d)({[\frac N d]}^2([\frac N d]+1))+3 ans=d=1Nμ(d)([dN]2([dN]+1))+3

总的时间复杂度为 O ( T N ) O(T\sqrt N) O(TN ),不用整除分块是 O ( T N ) O(TN) O(TN),本来是过不掉的,但是数据似乎太水了

#include<bits/stdc++.h>
using namespace std;
inline int re()
{
    int n=0,k=1;
    char c=getchar();
    while(c>'9'||c<'0'){if(c=='-') k=-1; c=getchar();	}
    while(c<='9'&&c>='0'){n=(n<<1)+(n<<3)+(c^48); c=getchar();	}
    return n*k;
}
bool pd[10010010];
int prime[10010010];
int miu[10010010];
int he[10010010];
inline int shai(int n)
{
    int tot=0;
    memset(pd,true,sizeof(pd));
    miu[1]=1;
    he[1]=1;
    for(int i=2;i<=n;i++)
    {
        if(pd[i])
        {
            prime[++tot]=i;
            miu[i]=-1;
        }
        for(int j=1;j<=tot&&i*prime[j]<=n;j++)
        {
            pd[i*prime[j]]=false;
            if(i%prime[j]==0){miu[i*prime[j]]=0;break;}
            else miu[i*prime[j]]=-miu[i];
        }
        he[i]=he[i-1]+miu[i];
    }
    return tot;
}
int n,T;
int main()
{
    T=re();
    shai(1001000);
    long long ans=0;
    while(T--)
    {
        n=re();
        ans=3ll;
        long long l=1,r;
        while(l<=n)
        {
            long long k=n/l;
            r=n/k;
            ans+=(he[r]-he[l-1])*k*k*(k+3);
            l=r+1;
        }		
        printf("%lld\n",ans);
    }
    return 0;
}
  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值