[POJ1845]Sumdiv(数论+矩乘)

该博客主要介绍了一道数论题目POJ1845,求解两个数AB的约数和并模9901。博主受到线性筛约数和的启发,通过质因数分解和矩乘优化,实现了O(30^3)的时间复杂度解决方案。内容包含题目描述、题解思路及代码实现。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

题目描述

传送门
题意:求 AB 的约数和,对9901取模。

题解

做这道题的时候受到线性筛约数和的启发
线性筛的方法是,令f(i)表示i的约数和,p为质数,那么f(i*p)=f(i)*p+f(?),其中?表示i除去所有质因子p剩下的数

那么对A分解质因数并且记录质因子次数,那么 AB 的质因子次数应该在原先的基础上扩大B倍
令f(i)表示将i~n所有质因子(包括次数)表示的数的约数和,那么实际上使x=x*prime(i)+f(i+1)做cnt遍就求出了f(i)
做cnt遍的过程用矩乘优化,这样复杂度就是 O(30309)

代码

#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
using namespace std;
#define N 30
#define Mod 9901

int a,b;
int prime[N],cnt[N],f[N];
struct data{int a[5][5];}unit,st,mi;

void calc(int x)
{
    for (int i=2;x>1&&i*i<=x;++i)
        if (x%i==0)
        {
            prime[++prime[0]]=i%Mod;
            while (x%i==0)
            {
                ++cnt[prime[0]];
                x/=i;
            }
            cnt[prime[0]]*=b;
        }
    if (x>1) prime[++prime[0]]=x%Mod,cnt[prime[0]]=b;
}
data cheng(data a,data b)
{
    data ans;memset(ans.a,0,sizeof(ans.a));
    for (int i=1;i<=3;++i)
        for (int j=1;j<=3;++j)
        {
            for (int k=1;k<=3;++k)
                ans.a[i][j]=ans.a[i][j]+a.a[i][k]*b.a[k][j];
            ans.a[i][j]%=Mod;
        }

    return ans;
}
data fast_pow(data a,int p)
{
    data ans=unit;
    for (;p;p>>=1,a=cheng(a,a))
        if (p&1)
            ans=cheng(ans,a);
    return ans;
}
int main()
{
    scanf("%d%d",&a,&b);
    if (b==0)
    {
        puts("1");
        return 0;
    }
    calc(a);
    f[prime[0]+1]=1;
    for (int i=1;i<=3;++i) unit.a[i][i]=1;
    for (int i=prime[0];i>=1;--i)
    {
        int f1=(f[i+1]*prime[i]+f[i+1])%Mod,f2=(f1*prime[i]+f[i+1])%Mod;
        if (cnt[i]==1)
        {
            f[i]=f1;
            continue;
        }
        memset(st.a,0,sizeof(st.a));
        st.a[1][1]=f2,st.a[1][2]=f1,st.a[1][3]=1;
        memset(mi.a,0,sizeof(mi.a));
        mi.a[1][1]=prime[i],mi.a[3][1]=f[i+1],mi.a[1][2]=1,mi.a[3][3]=1;

        mi=fast_pow(mi,cnt[i]-2);
        st=cheng(st,mi);
        f[i]=st.a[1][1];
    }
    printf("%d\n",f[1]);
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值