zoukankan      html  css  js  c++  java
  • * SPOJ PGCD Primes in GCD Table (需要自己推线性筛函数,好题)

    题目大意:

    给定n,m,求有多少组(a,b) 0<a<=n , 0<b<=m , 使得gcd(a,b)= p , p是一个素数

    这里本来利用枚举一个个素数,然后利用莫比乌斯反演可以很方便得到答案,但是数据量过大,完全水不过去

    题目分析过程(从别人地方抄来的)

    ans = sigma(p, sigma(d, μ(d) * (n/pd) * (m/pd)))
    
    Let s = pd, then
    
    ans = sigma(s, sigma(p, μ(s/p) * (n/s) * (m/s)))
        = sigma(s, (n/s) * (m/s) * sigma(p, μ(s/p)))
    
    Let g(x) = sigma(p, μ(x/p)), then
    
    ans = sigma(s, (n/s) * (m/s) * g(s))

    如果我们能预处理g(x)的话就能和前面一样分块搞了。这个时候我们多么希望g(x)μ(x)一样是积性函数。看完题解后,发现有一个不是积性函数,胜似积性函数的性质。由于题解没有给证明,所以就意淫了一个证明。

    考虑质数p'g(p'x) = sigma(p | p'x, μ(p'x/p))

    • x % p' == 0,那么考虑sigma中的变量p的所有取值,它和g(x)p是相同的。而μ(x)这个函数,如果x有平方因子的话就等于0,因此当p != p'μ(p'x/p) = 0,当p == p'是,μ(p'x/p) = μ(x)。所以g(p'x) = μ(x)
    • x % p' != 0,同样考虑p,会发现它的取值只比g(x)中的p多出一个p'。同理按照p是否等于p'讨论,可以得到g(p'x) = -g(x) + μ(x)

    因此g(x)可以在线性筛素数的时候算出。剩下的就是前缀和、分块了。

    然后自己看了很久才表示看懂- - ,这里加上自己的理解更加便于阅读

    因为我们根据线性筛法得到的莫比乌斯函数,欧拉函数这样的积性函数

    所以我们也总是希望上方的g(x)也可以是积性函数方便在线性时间内计算出答案

    对于一个积性函数来说,总是不断往前判断它乘上一个素数和当前的联系

    列一个素数 p' 那么g(p'x)的结果分析

    线性筛的中间过程总是判断当前 i % prime[j] 是否等于 0 , 所以这里也是判断当前x和p'的整除关系

    1.x%p' = 0  , 那么说明x中必然有一个p'的因子 , 那么sigma(p | p'x, μ(p'x/p)) 中,如果列举的p!=p' , 说明此时 p'*p' | (p'*x/p) ,也就是有至少2个p'的素数因子,所以

    mu[p'*x/p] = 0 , 唯一就只要考虑 p 取到 p'的时候  mu[p'*x/p] = mu[x] , 因为其他mu[]都为0->g(p'x)=mu[x]

    2.x%p'!=0 , 那么对于先前所有的 x/p 来说,此时乘了p'  , 若p!=p' , 那么因为多了一个因子 mu[p'*x/p] = -mu[x/p] , 所以在p!=p'时,所有的情况相加为-g(x),

    在考虑枚举到的p'=p , mu[p'*x/p]  = mu[x] , 所以所有答案叠加所得就是 -g[x]+mu[x]

     1 #include <cstdio>
     2 #include <cstring>
     3 #include <iostream>
     4 using namespace std;
     5 const int MAXN = 10000000;
     6 #define ll long long
     7 int mu[MAXN+5] , prime[MAXN/10] , g[MAXN+5] , sum[MAXN+5];
     8 bool check[MAXN+5];
     9 int tot;
    10 
    11 void getMu(int n)
    12 {
    13     memset(check , 0 ,  sizeof(check));
    14     tot=0;
    15     check[1] = true , mu[1]=1 , g[1]=0;
    16     for(int i=2 ; i<=n ; i++){
    17         if(!check[i]) prime[tot++]=i , mu[i]=-1 , g[i]=1;
    18         for(int j=0 ; j<tot ; j++){
    19             if(i*prime[j]>n) break;
    20             check[i*prime[j]]=true;
    21             if(i%prime[j] == 0){
    22                 mu[i*prime[j]] = 0;
    23                 g[i*prime[j]]=mu[i];
    24                 break;
    25             }
    26             else mu[i*prime[j]] = -mu[i] , g[i*prime[j]]=-g[i]+mu[i];
    27         }
    28     }
    29     sum[1]=g[1];
    30     for(int i=2 ; i<=n ; i++) sum[i]=g[i]+sum[i-1];
    31 }
    32 
    33 void swap(int &x , int &y)
    34 {
    35     int t=x;
    36     x=y , y=t;
    37 }
    38 
    39 int main()
    40 {
    41    // freopen("a.in" , "r" , stdin);
    42     getMu(MAXN);
    43     int T , n , m;
    44     scanf("%d" , &T);
    45     while(T--)
    46     {
    47         scanf("%d%d" , &n , &m);
    48         if(n>m) swap(n , m);
    49         ll ans = 0;
    50         int la=0;
    51         for(int s=1 ; s<=n ; s=la+1){
    52             //这里分块的意思很明显,因为在一个区间内,n/s ,m/s 的取值是不会变的,所以我们不用连续计算
    53             la = min(n/(n/s) , m/(m/s));
    54             ans = ans + (ll)(n/s)*(m/s)*(sum[la]-sum[s-1]);
    55         }
    56         printf("%lld
    " , ans);
    57     }
    58     return 0;
    59 }
  • 相关阅读:
    Linux常用的命令
    【练习】分区
    【测试】RAC搭建(裸设备)
    【练习】使用事务和锁定语句
    【练习】使用事务控制语句
    【练习】设置数据类型
    【练习】显示MySQLadmin 库户籍选项
    【练习】显示MYSQL客户机选项
    【练习】如何显示本地主机上的MySQL客户机版本
    【测试】切换保护模式,最大性能到最大可用
  • 原文地址:https://www.cnblogs.com/CSU3901130321/p/4426334.html
Copyright © 2011-2022 走看看