线性筛约数个数、约数和的新方法

最近本人脑洞大开,发现了一种线性筛约数个数和约数和的一种神奇方法。

目前网上的方法基本都是利用num[i]数组记录i最小的质因子的个数,然后进行转移。

其实可以省去num[i]数组,直接进行递推。

n的标准分解式为:

n=p1r1p2r2pkrk

那么n的约数个数及约数和分别为:

d(n)=i=1k(ri+1)

σ(n)=i=1k(1+pi+pi2++piri)

要想线性筛一个积性函数f(i),只需要知道4个东西。

1:f(1)=?

  很显然,d(1)=1,σ(1)=1;

2:f(p)=?,其中p为质数。

  同样很显然,d(p)=2,σ(p)=p+1;

34:f(ipj)=?ipj互质或不互质。

先说线性筛约数个数的方法,

ipj互质,则我们可以利用积性函数的性质直接推出:

d(ipj)=d(i)d(pj)=2d(i)

最后思考当ipj不互质的时候,考虑iipjipj的关系(由线性筛过程可知,ipj不互质即pj|i)

这里不妨设pj=p1,rj=r1

i=p1r1p2r2pkrk

ip1=p1r1+1p2r2pkrk

ip1=p1r11p2r2pkrk

则有:

d(i)=(r1+1)(r2+1)(rk+1)

d(ip1)=(r1+2)(r2+1)(rk+1)

d(ip1)=r1(r2+1)(rk+1)

d(i)r1+1后面的一大坨为T,即可表示为:

d(i)=(r1+1)T

d(ip1)=(r1+2)T=d(i)+T

d(ip1)=r1T=d(i)T

23式相加,整理得

d(ip1)=2d(i)d(ip1)

补上线性筛约数个数的代码:

复制代码
int pr[N],vis[N],d[N],cnt;
void make(int n)
{
    d[1] = 1;
    for(int i = 2;i<=n;i++)
    {
        if(!vis[i])pr[++cnt] = i,d[i] = 2;
        for(int j = 1;i*pr[j]<=n&&j<=cnt;j++)
        {
            vis[i*pr[j]] = 1;
            if(i%pr[j]==0)
            {
                d[i*pr[j]] = d[i]*2-d[i/pr[j]];
                break;
            }d[i*pr[j]] = d[i]*2;
        }
    }
}
复制代码

 

相比于使用num数组,这段代码就显得更加简洁,更重要的是它节省了内存.

同样,约数和也可以像这样筛出来.

ipj互质时,

σ(ipj)=σ(i)σ(pj)=σ(i)(pj+1)

不互质时,同样考虑σ(i)σ(ip1)σ(ip1)的关系.

设:

i=p1r1p2r2pkrk

ip1=p1r1+1p2r2pkrk

ip1=p1r11p2r2pkrk

则有:

σ(i)=(1+p1++p1r1)(1+p2++p2r2)(1+pk++pkrk)

σ(ip1)=(1+p1++p1r1+1)(1+p2++p2r2)(1+pk++pkrk)

σ(ip1)=(1+p1++p1r11)(1+p2++p2r2)(1+pk++pkrk)

同理,设后面的那一大串为T.

σ(i)=(1+p1++p1r1)T

σ(ip1)=(1+p1++p1r1+1)T=σ(i)+p1r1+1T

σ(ip1)=(1+p1++p1r11)T=σ(i)p1r1T

3式乘上p1再与2式相加,得到

σ(ip1)+p1σ(ip1)=(p1+1)σ(i)

整理,即:

σ(ip1)=(p1+1)σ(i)p1σ(ip1)

最后补上线性筛约数和的代码

复制代码
int pr[N],vis[N],sigma[N],cnt;
void make(int n)
{
    sigma[1] = 1;
    for(int i = 2;i<=n;i++)
    {
        if(!vis[i])pr[++cnt] = i,sigma[i] = i+1;
        for(int j = 1;i*pr[j]<=n&&j<=cnt;j++)
        {
            vis[i*pr[j]] = 1;
            if(i%pr[j]==0)
            {
                sigma[i*pr[j]] = sigma[i]*(pr[j]+1)-pr[j]*sigma[i/pr[j]];
                break;
            }sigma[i*pr[j]] = sigma[i]*(pr[j]+1);
        }
    }
}
复制代码

事实上它还可以再短一点(附上约数个数和约数和放在一起的版本):

复制代码
int n,pr[N],vis[N],d[N],sigma[N],cnt;
void make(int n)
{
    d[1] = sigma[1] = 1;
    for(int i = 2;i<=n;i++)
    {
        if(!vis[i])pr[++cnt] = i,d[i] = 2,sigma[i] = i+1;
        for(int j = 1;i*pr[j]<=n&j<=cnt;j++)
        {
            vis[i*pr[j]] = 1;
            d[i*pr[j]] = d[i]<<1;
            sigma[i*pr[j]] = sigma[i]*(pr[j]+1);
            if(i%pr[j]==0)
            {
                d[i*pr[j]]-=d[i/pr[j]];
                sigma[i*pr[j]]-=pr[j]*sigma[i/pr[j]];
                break;
            }
        }
    }
}
View Code
复制代码

 

posted @   ldysy2102  阅读(773)  评论(2编辑  收藏  举报
编辑推荐:
· go语言实现终端里的倒计时
· 如何编写易于单元测试的代码
· 10年+ .NET Coder 心语,封装的思维:从隐藏、稳定开始理解其本质意义
· .NET Core 中如何实现缓存的预热?
· 从 HTTP 原因短语缺失研究 HTTP/2 和 HTTP/3 的设计差异
阅读排行:
· 分享一个免费、快速、无限量使用的满血 DeepSeek R1 模型,支持深度思考和联网搜索!
· 基于 Docker 搭建 FRP 内网穿透开源项目(很简单哒)
· ollama系列01:轻松3步本地部署deepseek,普通电脑可用
· 按钮权限的设计及实现
· 25岁的心里话
点击右上角即可分享
微信分享提示