zoukankan      html  css  js  c++  java
  • 积性函数前缀和

    最近突然做到一些求积性函数前缀和的题,用到了各种筛,有一题用到 min_25 筛法,于是好好学习了一波,运用极不熟练。后来又遇到一道杜教筛的题,结果发现自己连 phi(x) 前缀和都不会推了?吓得我赶紧复习 + 写博客。

    前置技能 

    常见(完全)积性函数;整除分块; dirichlet 卷积;埃氏筛


    这里还是简单介绍一下 dirichlet 卷积

    我们用 * 表示两个函数的 dirichlet 卷积。

    定义函数 f 和 g 的 dirichlet 卷积形式为 $(f*g)(n)=sum_{d|n} f(d)gleft (frac{n}{d} ight )$ 。

    接下来,设 $f$ 为积性函数,我们要求 $sum_{i=1}^{n} f(i)$ , n 是 $10^9sim 10^{12}$ 左右的级别。本文忽略复杂度分析(需用到积分相关知识,本人太弱不会)。


    杜教筛

    (先扔链接:https://www.cnblogs.com/zzqsblog/p/5461392.html

    适用范围:存在另两个积性函数 $g,h$ ,且满足 $f*g=h$ , $g,h$ 函数的前缀和比较容易求。

    先来举几个例子:

    • 比如说要求 $mu$ 的前缀和,我们知道 $mu * 1=e$ ( $e(i)=[i=1]$ ), $e$ 的前缀和显然恒为 $[nge 1]$ , $1$ 的前缀和也很好求,所以可以使用杜教筛;
    • 又比如说要求 $varphi$ 的前缀和,我们知道 $varphi * 1=id$ , $1,id$ 的前缀和都很好求 ,所以也可以使用杜教筛。

    那么具体怎么求呢?

    设 $S(n)=sum_{i=1}^{n} f(i)$ ,那么有

    $$egin{aligned} sum_{i=1}^{n} (f*g)(i) &= sum_{i=1}^{n} sum_{d|i} f(d)gleft (frac{i}{d} ight ) \ S(n)=sum_{i=1}^n f(i) &= sum_{i=1}^{n} (f*g)(i) - sum_{i=1}^{n} sum_{d|i,d<i} f(d)gleft (frac{i}{d} ight ) \ &=sum_{i=1}^n (f*g)(i) - sum_{i=2}^n g(i) sum_{d=1}^{left lfloor frac{n}{i} ight floor} f(d)\ &=sum_{i=1}^n (f*g)(i) - sum_{i=2}^n g(i)Sleft (left lfloor frac{n}{i} ight floor ight ) end{aligned}$$

    可以发现我们得到了一个 $S(n)$ 的子问题,可以递归下去求解。

    至于边界情况,考虑如果 $nleq B$ 的时候我们直接用线性筛线性求解 $S(n)$ ,当 B 取 $n^{frac 23}$ 的时候复杂度最优,大概为 $mathcal{O}(n^{frac23})$ 。足以解决这个问题。

    贴代码 (一百年前写的)

     1 #include<bits/stdc++.h>
     2 #define rep(i,x,y) for (int i=(x);i<=(y);i++)
     3 #define ll long long
     4 using namespace std;
     5 const int N=1e6+5;
     6 int p[N/10],cnt,vis[N],mu[N]; map<ll,ll> mp;
     7 void init(int n){
     8     mu[1]=1;
     9     rep (i,2,n){
    10         if (!vis[i]) p[++cnt]=i,mu[i]=-1;
    11         for (int j=1;j<=cnt&&i*p[j]<=n;j++){
    12             vis[i*p[j]]=1;
    13             if (i%p[j]==0) break;
    14             mu[i*p[j]]=-mu[i];
    15         }
    16     }
    17     rep (i,2,n) mu[i]+=mu[i-1];
    18 }
    19 ll S(ll n){
    20     if (n<N) return mu[n];
    21     if (mp.count(n)) return mp[n];
    22     ll ret=1;
    23     for (ll i=2,last;i<=n;i=last+1){
    24         last=n/(n/i);
    25         ret-=(last-i+1)*S(n/i);
    26     }
    27     return mp[n]=ret;
    28 }
    29 int main(){
    30     init(N-1);
    31     ll a,b; cin>>a>>b; cout<<S(b)-S(a-1);
    32     return 0;
    33 }
    51nod1244 莫比乌斯函数之和

    min_25 筛

    (先扔链接:https://www.cnblogs.com/cjyyb/p/9185093.html

    一些吐槽:这个杜教筛真的是筛??好吧其实我也不知道它为什么叫杜教筛,这确实只是个递归求和的方法。( yhx : 应该叫杜氏求和法)

    不管这些,下面介绍一种真正的筛—— min_25 筛。它本质和洲阁筛相同,但小一倍常数。众所周知数论题非常容易被卡常数,所以学好 min_25 筛还是很重要的 ~

    适用范围:函数 f 需要满足 $f(p)$ 是一个低阶多项式, $f(p^c)$ 有较为容易的求法,最好是 $mathcal{O}(1)$ 。这里 p 代表任意素数。

    首先 min_25 筛和洲阁筛都是基于埃氏筛的。

    要求 $sum_{i=1}^{n} f(i)$ ,我们把 $1sim n$ 中的素数和合数分开考虑。

    假设素数集合为 $P$ , $P_i$ 为第 i 个素数, $minp(i)$ 表示 i 的最小素因子。 

    素数

    以下假设 $f(p)$ 为完全积性函数(但 $f(i)$ 不一定。原因下面会说)。

    要求 $sum_{i=1}^n [i in P]f(i)$ 。

    为了方便计算,这里所有数的函数值都定义为 $f(p)$ ,即都和 $f(p)$ 使用一样的表达式。

    设 $g(n,j)=sum_{i=1}^n [iin P or minp(i)>P_j]f(i)$ ,即筛到第 j 个素数为止剩下的数和所有素数的函数值之和。

    考虑 $g(n,j)$ 的转移:

    •  $P_{j}^2>n$ :显然 $g(n,j)=g(n,j-1)$ 。
    •  $P_{j}^2leq n$ :此时需要减去一些贡献。这个贡献怎么求呢?考虑筛 $P_j$ 的过程:找到剩下的数中所有 $P_j$ 的倍数删去,减去它们的函数值。由于 $f$ 为完全积性函数(上面说的原因在这里,如果不是完全积性就无法提出 $f(P_j)$ ,因而无法进行下面的变换),可以提出 $f(P_j)$ ,即 $f(P_j) gleft (left lfloor frac{n}{P_j} ight floor,j-1 ight )$ 。注意到若直接减去这个值,我们会把质数部分减去,而事实上我们需要把质数保留,因此这个贡献实际上是 $$f(P_j) [gleft (left lfloor frac{n}{P_j} ight floor,j-1 ight ) - sum_{i=1}^{j-1} f(P_i)]\ =f(P_j) [gleft (left lfloor frac{n}{P_j} ight floor,j-1 ight ) - g(P_j-1,j-1)]$$ 
    • 故 $$g(n,j)=egin{cases} g(n,j-1) & P_j^2>n \ g(n,j-1)-f(P_j) [gleft (lfloor frac{n}{P_j} floor,j-1 ight) - g(P_j-1,j-1)] & P_j^2leq n end{cases}$$

    那么 $sum_{i=1}^n [i in P]f(i)=g(n,|P|)$ ,这里 $P$ 是前 $sqrt{n}$ 个素数的集合。

    【稍微解释一下这个过程】

    由于我们最终用到的函数值只有 $g(n,|P|)$ ,即我们这步求的是素数的函数值,按照我们的筛法,实际上是先将所有数都按照 $f(p)$ 的式子计算函数值,然后把合数筛掉,最后只剩下素数,那么就能得到正确的函数值之和。但是中间筛的过程中,如果 $f(x)$ 的表达式和 $f(p)$ 不一样的话, $g$ 的值并不是真实的函数值。所以注意这里只求出了素数函数值之和

    合数

    这里没有素数那栏里的若干限制。

    设 $S(n,j)=sum_{i=1}^n [minp(i)geq P_j]f(i)$ 。和 $g$ 类似的定义,但注意这里取到了等号,且不包含素数。

    $S$ 的值分 2 部分计算:质数,合数。 1 的函数值放到最后计算。

    质数的函数值即 $g(n,|P|)-sum_{i=1}^{j-1} f(P_i)$ 。至于合数,由于 $f(x)$ 是积性函数,那么可以把最小质因子全部提出来(它和剩下部分互质),所以枚举它的最小质因子 $P_k$ 及幂次 $e$ ,得 $$S(n,j)=g(n,|P|)-sum_{i=1}^{j-1}f(P_i)+sum_{kgeq j}sum_{e>0}^{P_k^{e+1}leq n}(f(P_k^e)Sleft (left lfloor frac{n}{P_k^e} ight floor,k+1 ight )+f(P_k^{e+1}))$$

    后面的 $f(P_k^{e+1})$ 是单独处理的特殊情况,由于前面无法统计到。

    这样总复杂度为 $mathcal{O}(frac{n^{frac34}}{log n})$ 。

    【为什么要求 $f(p)$ 为低阶多项式呢】

    我们在处理素数部分的时候说过假设 $f(p)$ 为完全积性函数。那么低阶多项式一定能拆成 $sum_{i} a_ix^i$ ,这里每一项都是完全积性的,所以只需拆项使用上述方法求解即可。复杂度 $ imes 阶数$ 。

    栗子

    loj6053 简单的函数

    这题的函数比较奇怪,但好在它是积性的; $f(p)$ 除 $p=2$ 时均为 $p-1$ ,是一个低阶多项式; $f(p^c)$ 可以 $mathcal{O}(1)$ 求解。所以可以用 min_25 筛求解。

    把 $f(p)$ 拆掉, $p=2$ 的时候特判一下。别的直接套用上面的公式即可。

    code:

     1 #include<bits/stdc++.h>
     2 #define rep(i,x,y) for (int i=(x);i<=(y);i++)
     3 #define ll long long
     4 
     5 using namespace std;
     6 
     7 const int N=2e5+10,mod=1e9+7;
     8 ll n,w[N]; int m,sq,id1[N],id2[N],h[N],g[N],p[N],cnt,pre[N],vis[N];
     9 
    10 inline int get_id(ll x){return x<=sq?id1[x]:id2[n/x];}
    11 
    12 void sieve(int n){
    13     rep (i,2,n){
    14         if (!vis[i]) p[++cnt]=i,pre[cnt]=(pre[cnt-1]+i)%mod;
    15         for (int j=1;j<=cnt&&i*p[j]<=n;j++){
    16             vis[i*p[j]]=1;
    17             if (i%p[j]==0) break;
    18         }
    19     }
    20 }
    21 
    22 int S(ll n,int i){
    23     if (n<=1||p[i]>n) return 0;
    24     int k=get_id(n);
    25     int res=((ll)g[k]-h[k]-pre[i-1]+i-1)%mod; res=(res+mod)%mod;
    26     if (i==1) res+=2;
    27     for (int j=i;j<=cnt&&(ll)p[j]*p[j]<=n;j++){
    28         ll pw=p[j];
    29         for (int e=1;pw<=n/p[j];e++)
    30             res=(res+((ll)(p[j]^e)*S(n/pw,j+1)%mod+(p[j]^(e+1)))%mod)%mod,pw*=p[j];
    31     }
    32     return res;
    33 }
    34 
    35 int main(){
    36     scanf("%lld",&n); sq=sqrtl(n); sieve(sq);
    37     for (ll i=1,j;i<=n;i=j+1){
    38         ll x=n/i; j=n/x; w[++m]=x;
    39         if (x<=sq) id1[x]=m; else id2[j]=m;
    40         h[m]=(x-1)%mod,g[m]=(x+2)%mod*((x-1)%mod)%mod;
    41         (g[m]&1)?g[m]+=mod:0; g[m]>>=1;
    42     }
    43     rep (j,1,cnt)
    44         for (int i=1;i<=m&&(ll)p[j]*p[j]<=w[i];i++){
    45             int k=get_id(w[i]/p[j]);
    46             h[i]=((ll)h[i]-h[k]+j-1)%mod,h[i]=(h[i]+mod)%mod;
    47             g[i]=((ll)g[i]+mod-(ll)p[j]*(g[k]+mod-pre[j-1])%mod)%mod;
    48         }
    49     printf("%d
    ",S(n,1)+1);
    50     return 0;
    51 }
    loj6053 - 简单的函数
  • 相关阅读:
    代理模式
    spring aop
    mybatis从入门到精通(其他章节)
    mybatis从入门到精通(第6章)
    Java中Compareable和Comparator两种比较器的区别
    Java的equals方法的使用技巧
    Dubbo的配置过程,实现原理及架构详解
    什么是IPFS?IPFS与区块链有什么关系
    leetCode242 有效的字母异位词
    需要多个参数输入时-----------------考虑使用变种的Builder模式
  • 原文地址:https://www.cnblogs.com/bestFy/p/10701494.html
Copyright © 2011-2022 走看看