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

vfleaking的博客

My name is VFlea King

 
 
 

日志

 
 

SRM 582 Div1 1000pts 题解  

2013-06-20 16:12:48|  分类: Topcoder |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |
【题目大意】
一个数n被称为半完美数当且仅当存在a、b、c三个整数使得b > 1, 1 <= c < a, n = c * a ^ b。
给你L,R,求区间[L, R]中半完美数的个数。

1 <= L <= R <= 8 * 10 ^ 16

【题解】
显然把问题转化成[1, N]上半完美数的个数。

首先b > 1是坑你的。b只可能为2、3。
因为如果b >= 4,且b是偶数,那么设b = 2k。如果n = c * a ^ (2k),那么必有n = c * (a ^ 2) ^ k。显然k > 1, 1 <= c < a ^ 2。
如果b是奇数,那么设b = 2k + 1。如果n = c * a ^ (2k + 1),那么必有n = (c * a) * (a ^ 2) ^ k。显然k > 1, 1 <= c * a < a ^ 2。
于是只要n能表示成b >= 4,那么n一定可以表示成b = 2或b = 3。
所以考虑b = 2和b = 3的情况就好了。

那么能表示成a * x ^ 2的n有多少个呢?
如果对于一个能写成这种形式的n,我们取其中一个形式,且a有平方因子,一定能把平方因子移到x中去。这样a变小,x变大,仍然有1 <= a < x,仍是合法的表示方式。
无论取哪种形式最后都能转化为a是无平方因子的数,而且把n写成a * x ^ 2, a为无平方因子的数,只有唯一一种。
这样统计n与有(a, x)一一对应关系。a是无平方因子的数,1 <= a < x。
好……那么假设已经知道了a:
a * x ^ 2 <= N
1 <= a < x
解得:
a < x <= (N / a) ^ (1 / 2)
a < x <= floor((N / a) ^ (1 / 2))
所以对于给定的a,如果a < floor((N / a) ^ (1 / 2)),那么合法x的个数为floor((N / a) ^ (1 / 2)) - a个。
所以能表示成a * x ^ 2的n的个数为
--------------------
 \
  \
   |                  [a为无平方因子的数] floor((N / a) ^ (1 / 2)) - a
  /
 /
--------------------
1 <= a < N ^ (1 / 3)

同样,能表示成b * y ^ 3的n个个数为
--------------------
 \
  \
   |                  [b为无立方因子的数] floor((N / b) ^ (1 / 3)) - b
  /
 /
--------------------
1 <= b < N ^ (1 / 4)

然后发现布吉岛怎么做了,因为有些数既可以表示成a * x ^ 2也可以表示成b * y ^ 3。

考虑一个数n,n = b * y ^ 3, 1 <= b < y
那么n = (b * y) * y ^ 2
再设最大的k满足(k ^ 2) \ (b * y)。这里a \ b表示a是b的约数。
那么可以写成:n = (b * y / (k ^ 2)) * (y * k) ^ 2。
显然如果此时b * y / (k ^ 2) >= y * k,那么n就肯定不能表示成a * x ^ 2的形式了吧。
搞一搞这个不等式发现:k <= b ^ (1 / 3)
所以k要满足什么条件呢?
k >= 1, (k ^ 2) \ (b * y), (b * y) / (k ^ 2)是无平方因子的数, k <= b ^ (1 / 3)。

好像很可捉了。
我们首先求出能表示成a * x ^ 2的数的个数。
这一步是O(N ^ (1 / 3) * log N)的。因为开根号有个log。
然后再来求能表示成b * y ^ 3但不能表示成a * x ^ 2的数的个数。
简记P2(n)表示“n为无平方因子的数”,P3(n)表示“n为无立方因子的数”。
就是求:
-------------------- ----------- -----------
 \                    \           \
  \                    \           \
   |                    |           |        [P3(b)][b * y ^ 3 <= N][(k ^ 2) \ (b * y)]
   |                    |           |        [P2((b * y) / (k ^ 2))][k <= b ^ (1 / 3)]
  /                    /           /
 /                    /           /
-------------------- ----------- -----------
1 <= b < N ^ (1 / 4)    b < y      1 <= k
(后面写成两行了主要是因为一行排不下)
然后交换后面两个求和式
-------------------- --------------------- -----------
 \                    \                     \
  \                    \                     \
   |                    |                     |        [P3(b)][b * y ^ 3 <= N]
   |                    |                     |        [(k ^ 2) \ (b * y)][P2((b * y) / (k ^ 2))]
  /                    /                     /
 /                    /                     /
-------------------- --------------------- -----------
1 <= b < N ^ (1 / 4) 1 <= k <= b ^ (1 / 3)   b < y
然后再玩一下:
--------------------        --------------------- --------------------------
 \                           \                     \
  \                           \                     \
   |                 [P3(b)]   |                     |                       [(k ^ 2) \ (b * y)]
   |                           |                     |                       [P2((b * y) / (k ^ 2))]
  /                           /                     /
 /                           /                     /
--------------------        --------------------- --------------------------
1 <= b < N ^ (1 / 4)        1 <= k <= b ^ (1 / 3) b < y <= (N / b) ^ (1 / 3)

首先看[(k ^ 2) \ (b * y)]
设:
k' = k ^ 2 / gcd(k ^ 2, b)
b' = b / gcd(k ^ 2, b)
那么[(k ^ 2) \ (b * y)] = [k' \ (b' * y)]
然后k'与b'互质嘛,所以[(k ^ 2) \ (b * y)] = [k' \ y]

然后再看[P2((b * y) / (k ^ 2))]
发现
[P2((b * y) / (k ^ 2))] = [P2((b' * y) / k')]

于是继续玩玩吧……
--------------------        --------------------- --------------------------------
 \                           \                     \                          
  \                           \                     \                         
   |                 [P3(b)]   |                     |                             [P2(b' * y')]
  /                           /                     /
 /                           /                     /
--------------------        --------------------- --------------------------------
1 <= b < N ^ (1 / 4)        1 <= k <= b ^ (1 / 3) b < k' * y' <= (N / b) ^ (1 / 3)

--------------------        --------------------- -------------------------------------
 \                           \                     \                          
  \                           \                     \                         
   |                 [P3(b)]   |                     |                                  [P2(b' * y')]
  /                           /                     /
 /                           /                     /
--------------------        --------------------- -------------------------------------
1 <= b < N ^ (1 / 4)        1 <= k <= b ^ (1 / 3) b / k' < y' <= (N / b) ^ (1 / 3) / k'

啊啊啊!!!屏幕不够用了!!!
SRM 582 Div1 1000pts 题解 - vfleaking - vfleaking的博客

再来玩一下[P2(b' * y')]
[P2(b' * y')] = [P2(b')][P2(y')][gcd(b', y') = 1]
很好丽洁吧……
然后自行脑补前面的sigma吧……SRM 582 Div1 1000pts 题解 - vfleaking - vfleaking的博客
-------------------------------------
 \                          
  \                         
   |                                  [P2(b')][P2(y')][gcd(b', y') = 1]
  /
 /
-------------------------------------
b / k' < y' <= (N / b) ^ (1 / 3) / k'

-------------------------------------                 ---------------
 \                                                     \
  \                                                     \
   |                                  [P2(b')][P2(y')]   |            mu(d)
  /                                                     /
 /                                                     /
-------------------------------------                 ---------------
b / k' < y' <= (N / b) ^ (1 / 3) / k'                 d \ gcd(b', y')
其中mu是莫比乌斯函数。

        ---------------      -------------------------------------
         \                    \                                    
          \                    \                                   
[P2(b')]   |            mu(d)   |                                 [d \ y'][P2(y')]
          /                    /                                   
         /                    /                                    
        ---------------      -------------------------------------
             d \ b'          b / k' < y' <= (N / b) ^ (1 / 3) / k'

然后就好捉了……
y'不会超过N ^ (1 / 3)……
于是我们处理出来cntP2[d][m],表示[1, dm]中有多少个既是d的倍数又是无完全平方因子的数。
这个数组大小显然O(N ^ (1 / 3) log (N ^ (1 / 3))) = O(N ^ (1 / 3) log N)吧……

其中涉及枚举b'的约数。
我sb,先用筛法求出对于一个数n,哪些素数能整除它……用个vector<int>存下来。由于b'没有平方因子,所以可以直接枚举每个素数的次幂是0还是1了……然后再O(log N)把这些素数乘起来……
于是……时间复杂度应该是:
   --------------------        ---------------------         ---------------
    \                           \                             \             
     \                           \                             \            
      |                 [P3(b)]   |                  [P2(b')]   |             O(log N)
     /                           /                             /            
    /                           /                             /                 
   --------------------        ---------------------         ---------------
   1 <= b < N ^ (1 / 4)        1 <= k <= b ^ (1 / 3)              d \ b'

   -------------------- ---------------------- ---------------
    \                    \                      \             
     \                    \                      \            
<=    |                    |                      |             O(log N)
     /                    /                      /            
    /                    /                      /                 
   -------------------- ---------------------- ---------------
   1 <= b < N ^ (1 / 4) 1 <= k <= N ^ (1 / 12)       d \ b'

               -------------------- ---------------
                \                    \             
                 \                    \            
<= N ^ (1 / 12)   |                    |            O(log N)
                 /                    /            
                /                    /             
               -------------------- ---------------
               1 <= b < N ^ (1 / 4)      d \ b

注意最后一步把b'替换成了b是因为b'是b的约数,其约数个数肯定没有b多。
而1~n的约数的个数之和是O(n log n)的。

所以求“能表示成b * y ^ 3但不能表示成a * x ^ 2的数的个数”的时间复杂度大致是
O(N ^ (1 / 12) * N ^ (1 / 4) log (N ^ (1 / 4)) log N) = O(N ^ (1 / 3) log N log N)

至于预处理嘛……用埃拉托斯特尼筛法乱搞。应该是O(N ^ (1 / 3) log N)。
这样总时间复杂度就是O(N ^ (1 / 3) log N log N)了。足以通过此题。

(好像公式还是超过屏幕了,我假定你们复制下来到奇怪的地方看那些公式好了……)

代码:

#include <iostream>
#include <cstdio>
#include <vector>
#include <algorithm>
using namespace std;

typedef long long s64;

inline int gcd(int a, int b)
{
while (b != 0)
{
a %= b;
swap(a, b);
}
return a;
}

class SemiPerfectPower
{
private:
static const s64 MaxN = 80000000000000000ll;
static const int MaxN_2 = 282842712;
static const int MaxN_3 = 430886;
static const int MaxN_4 = 16817;

// n = c * (a ** b)

bool isP2[1 + MaxN_3], isP3[1 + MaxN_4];
int *cntP2[1 + MaxN_3];
vector<int> fac[1 + MaxN_3];

inline int floor_sqrt(const s64 &n)
{
int le = 1, ri = MaxN_2 + 1;
while (le != ri)
{
int mid = (le + ri) >> 1;
if ((s64)mid * mid <= n)
le = mid + 1;
else
ri = mid;
}
return le - 1;
}
inline int floor_cubert(const s64 &n)
{
int le = 1, ri = MaxN_3 + 1;
while (le != ri)
{
int mid = (le + ri) >> 1;
if ((s64)mid * mid * mid <= n)
le = mid + 1;
else
ri = mid;
}
return le - 1;
}

inline s64 calc(const s64 &n)
{
s64 res = 0;

// n = c * (a ** 2)
for (int c = 1; (s64)c * c * c < n; c++)
if (isP2[c])
res += floor_sqrt(n / c) - c;

// n = b * (y ** 3) && n != a * (x ** 2)
for (int b = 1; (s64)b * b * b * b < n; b++)
if (isP3[b])
for (int k = 1; k * k * k <= b; k++)
{
int g = gcd(b, k * k);
int sb = b / g;
int sk = k * k / g;

if (!isP2[sb])
continue;
int minY = b, maxY = floor_cubert(n / b);
// y belongs to (minY, maxY]

for (int ds = 0; ds < (1 << fac[sb].size()); ds++)
{
int d = 1;
int mu = 1;
for (int i = 0; i < (int)fac[sb].size(); i++)
if (ds >> i & 1)
d *= fac[sb][i], mu *= -1;
res += mu * (cntP2[d][maxY / sk / d] - cntP2[d][minY / sk / d]);
}
}
return res;
}
public:
SemiPerfectPower()
{
fill(isP2 + 1, isP2 + MaxN_3 + 1, true);
for (int i = 2; i * i <= MaxN_3; i++)
for (int j = i * i; j <= MaxN_3; j += i * i)
isP2[j] = false;

fill(isP3 + 1, isP3 + MaxN_4 + 1, true);
for (int i = 2; i * i * i <= MaxN_4; i++)
for (int j = i * i * i; j <= MaxN_4; j += i * i * i)
isP3[j] = false;

for (int i = 1; i <= MaxN_3; i++)
{
int l = MaxN_3 / i;
cntP2[i] = new int[1 + l];

cntP2[i][0] = 0;
for (int j = 1; j <= l; j++)
cntP2[i][j] = cntP2[i][j - 1] + isP2[i * j];
}

for (int i = 2; i <= MaxN_3; i++)
if (fac[i].empty())
for (int j = i; j <= MaxN_3; j += i)
fac[j].push_back(i);
}
~SemiPerfectPower()
{
for (int i = 1; i <= MaxN_3; i++)
delete []cntP2[i];
}

s64 count(s64 numL, s64 numR)
{
return calc(numR) - calc(numL - 1);
}
};


  评论这张
 
阅读(1729)| 评论(12)
推荐 转载

历史上的今天

评论

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

页脚

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