注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

z55250825

一只蒟蒻

 
 
 

日志

 
 

【Spoj7001】【莫比乌斯反演】【Visible Lattice Points】  

2014-04-30 16:38:57|  分类: AC题目 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |
   题目大意:【仗仪队的升级版本】三维空间的一个n*n*n的方阵,左下角为(0,0,0),右上角为(n-1,n-1,n-1),求从(0,0,0)能看到多少个点?【能看到一个点当且仅当这个点到(0,0,0)的连线上没有其他的点】

   实际上考虑向量(x,y,z),从(0,0,0)沿(x,y,z)方向看(x,y,z)能看到的条件显然就是 (x,y,z)互质,否则的话会有一个更小的点挡住它们。所以问题实际上就是求1<=x,y,z<=n的互质数的个数。
   然后就是什么是莫比乌斯反演了...咱还是喜欢做 直接拉链接君了 恩= =【传送门
   上面那篇教程讲的非常不错0.0表示咱总算是理解了莫比乌斯反演的一点点额...
   其实咱们可以看看一个简单点的问题,解决了那个,这个其实也可以解决了0.0

   考虑求 互质数对(x,y)的个数,这里1<=x<=a,1<=y<=b。
   首先咱们定义 F[n]表示 在限制范围内 gcd(x,y)=n 的数对的个数。然后咱们再定义 G[n]表示 满足gcd(x,y)%n=0的数对的个数。
   而咱们知道gcd(x,y)的上限是U=max(a,b)
   则咱们可以通过G[]的定义得到 G[n]=∑F[i](i%n=0,i<=U)
   咱们又可以直接求出G[n],实际上gcd(x,y)%n=0的话,实际上要求x,y都包含n即可,则个数为 G[n]=(a/n)*(b/n)
   而咱们通过莫比乌斯反演,G[n]=∑F[i](i%n=0,i<=U)
   则F[n]=∑M[i/n]*G[i]=∑U[i/n]*(a/i)*(b/i) (i%n=0,i<=U)
   这个弄出来后,只需要求莫比乌斯函数M[]即可,这个由于是积性函数,所以可以类似欧拉函数一样在线性筛法的时候求出来。那么这道题就结束了。

   然后考虑这道题gcd(x,y,z)实际上是类似的...咱们可以同样得到 答案为 ∑M[i]*(n/i)*(n/i)*(n/i),但是!注意咱们没有考虑其中一位为0的情况,实际上咱们还得套用其中一位为0的情况,这个只需要把其中一个(n/i+3)即可【+3也就可以计算二维的答案了】。然后考虑还有两个0的情况,这个是3无疑,所以答案初始化的时候应该为3而不是0开始。

   然后WA一次,是没有考虑数据范围可以达到int64【话说每一次都是这样啊啊啊啊啊】
   然后TLE一次...加了个常数优化,只有miu[i]<>0的时候才计算,然后每一次除法只计算一次,用一个int64保存下来,然后自乘,然后每一个答案最后一口气输出。终于AC了额...
==========================

program spoj7001;
const maxn=1000010;
var miu,p:array[0..maxn]of longint;
    v:array[0..maxn]of boolean;
    a:array[0..100]of longint;
    ans:array[0..100]of int64;
    n:longint;
    q:int64;

procedure prepare;
var i,j:longint;
begin
  fillchar(v,sizeof(v),true);
  p[0]:=0;miu[1]:=1;
  for i:=2 to n do
   begin
     if v[i] then
      begin
        inc(p[0]);
        p[p[0]]:=i;
        miu[i]:=-1;
      end;
     for j:=1 to p[0] do
      begin
        if i*p[j]<=n then v[i*p[j]]:=false;
        if i*p[j]>n then break;
        if i mod p[j]=0 then
         begin
           miu[i*p[j]]:=0;
           break;
         end
        else miu[i*p[j]]:=-miu[i];
      end;
   end;
end;

procedure init;
var i,j,t:longint;
begin
  read(t);
  n:=-1;
  for i:=1 to t do
   begin
     read(a[i]);
     if n<a[i] then n:=a[i];
   end;
  prepare;
  for i:=1 to t do
   begin
     ans[i]:=3;
     for j:=1 to a[i] do
      if miu[j]<>0 then
       begin
         q:=a[i] div j;
         if miu[j]=1 then ans[i]:=ans[i]+q*q*(q+3)
                     else ans[i]:=ans[i]-q*q*(q+3);
       end;
   end;
  for i:=1 to t do writeln(ans[i]);
end;

begin
  init;
end.

   

   
   
   
  评论这张
 
阅读(51)| 评论(0)
推荐 转载

历史上的今天

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2017