zoukankan      html  css  js  c++  java
  • 51nod1222 最小公倍数计数 莫比乌斯反演 数学

    求$sum_{i = 1}^{n} sum_{j = 1}^{i} [lcm(i, j) le n]$
    因为这样不好求,我们改成求$sum_{i = 1}^{n} sum_{j = 1}^{n} [lcm(i, j) le n]$.
    这样求出来的值把除了(i, i)这样的点对以外所有点对都重复统计了一次。因此$ans = frac{rnt + n}{2}$(先加上没有重复统计的点对个数,使得所有点对都重复统计了一次,然后再除2就是不重复统计的点对个数)
    接下来就是化式子了...
    $$sum_{i = 1}^{n} sum_{j = 1}^{n} [lcm(i, j) le n]$$
    $$sum_{i = 1}^{n} sum_{j = 1}^{n} [frac{ij}{(i, j)} le n]$$
    枚举$(i, j) = d$
    $$sum_{d = 1}^{n} sum_{i = 1}^{n} sum_{j = 1}^{n} [frac{ij}{(i, j)} le n] [(i, j) == d]$$
    枚举$i = frac{i}{d}$
    $$sum_{d = 1}^{n} sum_{i = 1}^{lfloor{ frac{n}{d} } floor} sum_{j = 1}^{lfloor{ frac{n}{d} } floor} [ijd le n][(i, j) == 1]$$
    上反演
    $$sum_{d = 1}^{n} sum_{i = 1}^{lfloor{ frac{n}{d} } floor} sum_{j = 1}^{lfloor{ frac{n}{d} } floor}[ijd le n] sum_{k | (i, j)} mu(k)$$
    $$sum_{d = 1}^{n} sum_{i = 1}^{lfloor{ frac{n}{d} } floor} sum_{j = 1}^{lfloor{ frac{n}{d} } floor} sum_{k | (i, j)} mu(k)[ijd le n]$$
    把$k$提到前面来,然后枚举$k$的倍数
    $$sum_{d = 1}^{n} sum_{k = 1}^{lfloor{ frac{n}{d} } floor} mu(k) sum_{i = 1 }^{lfloor{ frac{n}{dk} } floor} sum_{j = 1}^{lfloor{ frac{n}{dk} } floor}[i j k^2 d le n]$$
    可以发现,原式中对$d,k$的限制实际上是$dk le n$,因此我们再次将$k$向前提,那么此时$k$已经是第一个枚举的了,因此我们对中括号内的式子进行移项得到:
    $$sum_{k = 1}^{n} mu(k) sum_{d = 1}^{lfloor{ frac{n}{k} } floor} sum_{i = 1}^{lfloor{ frac{n}{dk} } floor} sum_{j = 1}^{lfloor{ frac{n}{dk} } floor}[ijd le frac{n}{k^2}]$$
    那么为了保证$[ijd le frac{n}{k^2}]$,对现在枚举的变量有如下限制:
    $k le sqrt{n}, d le lfloor{ frac{n}{k^2} } floor, x le lfloor{ frac{n}{k^2d} } floor, y le lfloor{ frac{n}{k^2dx} } floor$
    因此原式变为:
    $$sum_{k = 1}^{sqrt{n}} mu(k) sum_{d = 1}^{lfloor{ frac{n}{k^2} } floor} sum_{i = 1}^{lfloor{ frac{n}{k^2d} } floor} sum_{j = 1}^{lfloor{ frac{n}{k^2dx} } floor} 1$$
    设$H(n) = sum_{a = 1}^{n} sum_{b = 1}^{lfloor{ frac{n}{a} } floor} sum_{c = 1}^{lfloor{ frac{n}{ab} } floor} 1$,则原式为:
    $$sum_{k = 1}^{sqrt{n}} mu(k) H(lfloor{ frac{n}{k^2} } floor)$$
    因此我们的目的就是快速求出$H(n)$。
    $H(n)$可以看做是求满足如下关系的三元组的个数$abc le n$.
    因此我们不妨假设$a < b < c$,那么有$a le n ^ {frac{1}{3}}, b le sqrt{lfloor{ frac{n}{a} } floor}$(因为除去$a$后剩余数为$lfloor{ frac{n}{a} } floor$,而$b$最大也要满足$b le c$,因此$b$最大也就是开根)
    那么此时$c$的范围就为$[b + 1, lfloor frac{n}{ab} floor]$,因此$c$的个数可以$O(1)$算出
    因为我们假定了$a < b < c$,而原式中没有这样的限制,所以这样会少算,比如$123$会被计算一次,但实际上应该被计算$123, 132, 213, 231, 321, 312$共$6$次。其他情况也是类似的。
    因此我们枚举$a, b, c$,分别求3个均不相同,前2个相同,后2个相同,3个都相同的方案数,
    然后对于这些方案数,分别乘上对应系数$6, 3, 3, 1$.最后加起来就等于$H(n)$.

     1 #include<bits/stdc++.h>
     2 using namespace std;
     3 #define AC 401000
     4 #define LL long long
     5 #define R register int
     6 #define RL register LL 
     7 
     8 LL l, r, tot, maxn;
     9 LL pri[AC], mu[AC];
    10 bool z[AC];
    11 
    12 void pre()
    13 {
    14     scanf("%lld%lld", &l, &r), maxn = sqrt(r);
    15     mu[1] = 1;
    16     for(R i = 2; i <= maxn; i ++)
    17     {
    18         if(!z[i]) pri[++ tot] = i, mu[i] = -1;
    19         for(R j = 1; j <= tot; j ++)
    20         {
    21             int now = pri[j];
    22             if(i * now > maxn) break;
    23             z[i * now] = true;
    24             if(!(i % now)) break;
    25             mu[i * now] = - mu[i];
    26         }
    27     }
    28 }
    29 
    30 LL get(LL n)
    31 {
    32     LL A = 0, B = 0, C = 0, D = 0;
    33     for(RL i = 1; i * i * i <= n; i ++)//枚举3个不相同的
    34     {
    35         LL m2 = sqrt(n / i), to = n / i;//因为k要从j + 1 枚举到 n / i / j,所以相减就是个数
    36         for(RL j = i + 1; j <= m2; j ++) A += (LL)(to / j - j);//除法如此之慢。。。。
    37     }
    38     for(RL i = 1; i * i * i <= n; i ++) B += (LL)(n / i / i - i);//加上只有前2个相同的
    39     for(RL i = 1; i * i * i <= n; i ++) C += ((LL)sqrt(n / i) - i);//枚举后2个相同的,因为后2个相同,所以只需要知道j的个数即可
    40     for(RL i = 1; i * i * i <= n; i ++) ++ D;//因为诡异精度误差,,,,就算是加减,也要先转成LL....
    41 //    D += m1;//因为前面是double类型,所以一遇到乘法,就可能造成进位。。。
    42     return A * 6 + B * 3 + C * 3 + D;
    43 } 
    44 
    45 LL cal(LL n)
    46 {
    47     int block = sqrt(n); LL rnt = 0;
    48     for(R i = 1; i <= block; i ++) rnt += mu[i] * get(n / i / i);
    49     return (rnt + n) / 2;
    50 }
    51 
    52 int main()
    53 {
    54     freopen("in.in", "r", stdin);
    55     pre();
    56     printf("%lld
    ", cal(r) - cal(l - 1));
    57     fclose(stdin);
    58     return 0;
    59 }
    View Code
  • 相关阅读:
    《ActionScript 3 CookBook 简体中文完整版》下载
    打开组件服务超慢,打不开属性窗口。
    无法引用Microsoft.Office.Interop.Excel的解决
    HttpWebResponse类
    反射性能优化 标记个
    配置文件入门 WebConfig.config常用配置节点介绍
    配置文件的读写
    HTTP权威指南阅读记录 第一章
    IIS
    锁机制与原子操作 <第四篇>
  • 原文地址:https://www.cnblogs.com/ww3113306/p/10176059.html
Copyright © 2011-2022 走看看