洛谷 P4714 「数学」约数个数和 解题报告

P4714 「数学」约数个数和

题意(假):每个数向自己的约数连边,给出(n,k(le 10^{18})),询问(n)的约数形成的图中以(n)为起点长为(k)的链有多少条(注意每个点都有自环)

这样想是做不出来题的。

正常的题意是:询问(n)的约数的约数的....(共(k)次复读后)约数个数和。


考虑(f_k(n))表示答案。

显然有(f_{k}(n)=sum_{d|n}f_{k-1}(d))

注意到用数论卷积的形式可以表示为

[mathtt f_k=mathtt f_{k-1}*mathtt {Id_0} ]

因为(mathtt f_0=mathtt d),即约数个数,又因为积性函数乘积性函数为积性函数

所以(mathtt f_k)是积性函数。

考虑求(mathtt f_k(p^c))

我们有(mathtt f_k(p^c)=sum_{i=0}^cf_{k-1}(p^i))

我们注意到每次(k-1),实际上是对(p)有一个划分,我们可以形象的理解为,(c)个物品,中间插(k+1)个板子,板子之间可以没有物品的方案数,其中物品无序,板子有序(因为每次只能从后往前放置)

根据插板法强制钦定板子的物品,可以得到

[mathtt f_k(p^c)=inom{c+k+1}{k+1} ]

这个(c)很小,直接算算就可以了。

分解大数可以直接使用Pollard_Rho算法


Code:

#include <cstdio>
#include <cctype>
#include <cstdlib>
#include <algorithm>
#define ll long long
using std::min;
const int SIZE=1<<21;
char ibuf[SIZE],*iS=ibuf,*iT=ibuf;
//#define gc() (iS==iT?(iT=(iS=ibuf)+fread(ibuf,1,SIZE,stdin),iS==iT?EOF:*iS++):*iT++)
#define gc() getchar()
template <class T>
void read(T &x)
{
    x=0;char c=gc();
    while(!isdigit(c)) c=gc();
    while(isdigit(c)) x=x*10+c-'0',c=gc();
}
const int pri[]={2,3,5,7,11,13,17,19,23,29,31,37};
void add(ll &a,ll b,ll p){a=a+b>=p?a+b-p:a+b;}
ll mul(ll d,ll k,ll p)
{
    ll f=0;
    while(k)
    {
        if(k&1) add(f,d,p);
        add(d,d,p);
        k>>=1;
    }
    return f;
}
ll qp(ll d,ll k,ll p)
{
    ll f=1;
    while(k)
    {
        if(k&1) f=mul(f,d,p);
        d=mul(d,d,p);
        k>>=1;
    }
    return f;
}
bool Miller_Rabin(ll n)
{
    if(n==1) return false;
    for(int i=0;i<12;i++) if(n%pri[i]==0) return n==pri[i];
    ll res=n-1;int k=0;
    while(!(res&1)) res>>=1,++k;
    for(int i=0;i<12;i++)
    {
        ll x=qp(pri[i],res,n);
        for(int j=0;j<k&&x>1;j++)
        {
            ll y=mul(x,x,n);
            if(y==1&&x!=n-1) return false;
            x=y;
        }
        if(x!=1) return false;
    }
    return true;
}
ll F(ll x,ll c,ll p) {return (mul(x,x,p)+c)%p;}
ll gcd(ll a,ll b){return b?gcd(b,a%b):a;}
ll Find(ll n)
{
    ll x,y=rand()%n,c=rand()%n;
    int w=1<<9;
    for(int l=1;;l<<=1)
    {
        x=y;
        for(int i=0;i<l;i++) y=F(y,c,n);
        for(int i=0;i<l;i++)
        {
            int le=min(l-i,w);
            ll g=1,las=y;
            for(int j=0;j<le;j++) y=F(y,c,n),g=mul(g,(y+n-x)%n,n);
            g=gcd(g,n);
            if(g==1) continue;
            if(g==n)
            {
                g=1,y=las;
                while(g==1) y=F(y,c,n),g=gcd((y+n-x)%n,n);
            }
            return g;
        }
    }
}
ll s[1<<19],n,k,ans=1;
int tot;
void Pollard_Rho(ll n)
{
    if(n==1) return;
    if(Miller_Rabin(n)) {s[++tot]=n;return;}
    ll d=Find(n);
    while(d==n) d=Find(n);
    Pollard_Rho(d),Pollard_Rho(n/d);
}
const ll mod=998244353;
ll C(ll m,ll n)
{
    ll f=1;
    for(ll i=1;i<=n;i++) f=mul(f,i,mod);
    f=qp(f,mod-2,mod);
    for(ll i=m;i>m-n;i--) f=mul(f,i,mod);
    return f;
}
int main()
{
    read(n),read(k);
    Pollard_Rho(n);
    std::sort(s+1,s+1+tot);
    tot=std::unique(s+1,s+1+tot)-s-1;
    for(int i=1;i<=tot;i++)
    {
        ll c=0;
        while(n%s[i]==0) n/=s[i],++c;
        ans=mul(ans,C(c+k+1,c),mod);
    }
    printf("%lld
",ans);
    return 0;
}

2019.4.29

原文地址:https://www.cnblogs.com/butterflydew/p/10789725.html