powerful number求积性函数前缀和

orzzzq

类似杜教筛的构造思想。

构造容易求的东西求出f的前缀和

不同之处是，杜教筛构造出两个可以快速求前缀和
但是powerful number（简称pn）构造方法有所不同。

 

对于这个题

数据范围，min_25筛显然是过不去的。最多跑到1e11

$f(p^e)=p^k$启发

构造$g(x)=x^k$和相似一些

再构造$h(x)$，使得$f=h*g$这里表示狄利克雷卷积

那么$\sum_i f(i)=\sum_i\sum_{d|i}h(i)\times g(i/d)=\sum_d h(d)\sum_{i=1}^{n/d}g(i)$

g的前缀和好算 。拉格朗日差值。

h？

先考虑h是什么

 

首先，f是积性函数，g也是，不妨考虑能不能构造h也是

要使得对于$f=h*g$

只需要对于$f(p^e)$都是成立的。

首先$h(1)=1$

$f(p)=h(1)\times g(p)+h(p)\times g(1)=p$故$h(p)=0$

每个质因子次数都>=2的才有值，

这些数就叫做powerful number

每个数都可以写成$a^2b^3$的形式。

 

那么发现，powerful number只有$\sqrt(n)$个

可以暴力枚举！

h有值情况是什么？

$f(p^e)=\sum_{i=0}^{e} h(p^i)\times g(p^{e-i})$

$h(p^e)=f(p^e)-\sum_{i=0}^{e-1} h(p^i)\times g(p^{e-i})$

已经可以前缀和优化求出$h(p^e)$了

根据本题的g的定义，归纳可以得到$h(p^e)=p^k-p^{2k}$

 

然后可以通过本题。

具体的，预处理$<=sqrt(n)$的质数。

预处理拉格朗日差值需要的东西。

枚举h的话，

法一：枚举$a^2b^3$再进行质因数分解。

法二：枚举每个质因子的次数。从2开始。dfs枚举，或者类似min_25筛的方法枚举也一样

（注意枚举方法：

枚举下一个质因子>=2的位置！不要什么dfs(d,x,h)->dfs(d+1,x,h)表示不选这样。会多余很多不选择的情况。

而枚举下一个，可以O(总个数)枚举所有的。见代码

）

 

$O(sqrt(n)*k)$

加强一下，

预处理$<=sqrt(n)$的g的前缀和（可以用线性筛）

可以做到：$O(sqrt(n)+sqrt(n)+sqrt(sqrt(n))\times k)$

k=5000,n=1e14也可以3s内跑过！

#include<bits/stdc++.h>
#define reg register int
#define il inline
#define fi first
#define se second
#define mk(a,b) make_pair(a,b)
#define numb (ch^'0')
#define pb push_back
#define solid const auto &
#define enter cout<<endl
#define pii pair<int,int>
//#define int long long  
using namespace std;
typedef long long ll;
template<class T>il void rd(T &x){
    char ch;x=0;bool fl=false;while(!isdigit(ch=getchar()))(ch=='-')&&(fl=true);
    for(x=numb;isdigit(ch=getchar());x=x*10+numb);(fl==true)&&(x=-x);}
template<class T>il void output(T x){if(x/10)output(x/10);putchar(x%10+'0');}
template<class T>il void ot(T x){if(x<0) putchar('-'),x=-x;output(x);putchar(' ');}
template<class T>il void prt(T a[],int st,int nd){for(reg i=st;i<=nd;++i) ot(a[i]);putchar('\n');}
namespace Modulo{
const int mod=1e9+7;
il int ad(int x,int y){return x+y>=mod?x+y-mod:x+y;}
il int sub(int x,int y){return ad(x,mod-y);}
il int mul(int x,int y){return (ll)x*y%mod;}
il void inc(int &x,int y){x=ad(x,y);}
il void inc2(int &x,int y){x=mul(x,y);}
il int qm(int x,int y=mod-2){int ret=1;while(y){if(y&1) ret=mul(x,ret);x=mul(x,x);y>>=1;}return ret;}
template<class ...Args>il int ad(const int a,const int b,const Args &...args) {return ad(ad(a,b),args...);}
template<class ...Args>il int mul(const int a,const int b,const Args &...args) {return mul(mul(a,b),args...);}
}
using namespace Modulo;
namespace Miracle{
const int N=1e7+4;
const int K=5005;
ll n,k;
int mi[N];
int pri[N],cnt;
int vis[N];
int sum[N];
int sqr;
void sieve(ll n){
    sum[1]=1;
    for(reg i=2;i<=n;++i){
        if(!vis[i]){
            pri[++cnt]=i;
            sum[i]=mi[cnt]=qm(i,k);
        }
        for(reg j=1;j<=cnt;++j){
            if(i*pri[j]>n) break;
            vis[i*pri[j]]=1;
            sum[i*pri[j]]=mul(sum[i],sum[pri[j]]);
            if(i%pri[j]==0){
                break;
            }
        }
    }
    for(reg i=1;i<=n;++i){
        inc(sum[i],sum[i-1]);
    }
}
int Y[K];
int zh[K],fu[K];
void pre(){
    zh[0]=1;fu[0]=1;
    for(reg i=1;i<=k+2;++i) {
        inc(Y[i],ad(Y[i-1],qm(i,k)));
        zh[i]=mul(zh[i-1],qm(i));
        fu[i]=mul(fu[i-1],qm(mod-i));
    }
//    prt(Y,1,k+2);
//    prt(zh,1,k+2);
//    prt(fu,1,k+2);
}
int pr[K],bc[K];
int num;
int calc(ll n){
//    cout<<" calc "<<n<<endl;
    if(n<=sqr) {
//        cout<<" sum "<<sum[n]<<endl;
        return sum[n];
    }
    if(n<=k+2) {
//        cout<<" Y "<<Y[n]<<endl;
        return Y[n];
    }
    n%=mod;
//    ++num;
    pr[0]=1;bc[k+3]=1;
    for(reg i=1;i<=k+2;++i){
        pr[i]=mul(pr[i-1],n-i);
    }
    for(reg i=k+2;i>=1;--i){
        bc[i]=mul(bc[i+1],n-i);
    }
    int ret=0;
    for(reg i=1;i<=k+2;++i){
        inc(ret,mul(Y[i],pr[i-1],bc[i+1],zh[i-1],fu[k+2-i]));
    }
//    cout<<" ret "<<ret<<endl;
    return ret;
}
int ans;
int tot=0;
void dfs(int d,ll x,int h){
    if(d>cnt) return;
    for(reg y=d;y<=cnt&&(ll)x<=(n/pri[y]/pri[y]);++y){
        int tmp=mul(h,sub(mi[y],mul(mi[y],mi[y])));
        ll now=x*pri[y]*pri[y];
        for(reg e=2;;++e){
//            cout<<" now "<<now<<endl;
            inc(ans,mul(tmp,calc(n/now)));
            dfs(y+1,now,tmp);
            if(now>n/pri[y]) break;
            now*=pri[y];
        }
    }
}
int main(){
    rd(n);rd(k);
    sqr=sqrt(n);
    sieve(sqr);
    pre();
//    cout<<" nd "<<endl;
    ans=calc(n);
    dfs(1,1,1);
    cout<<ans;//<<" time "<<tot<<" cha "<<num<<endl;
    return 0;
}

}
signed main(){
    Miracle::main();
    return 0;
}

/*
   Author: *Miracle*
*/

 

 

对于LOJ简单的函数——Min_25筛

也可以用powerful number做的

令$G(x)=\phi$，p=2需要特殊处理。

这里h就没有O(1)公式了，还不能前缀和优化，只能暴力卷积

h[p][e]表示h(p^e)那么实际上计算h的总复杂度也不是很大。

好吧其实感觉挺大的，但是看不懂zzq在写什么神仙操作。。。。

 

 

提供了一种处理前缀和的思路。

也是构造，但是复杂度降低的原因是，g前缀和好求，h有值的位置不多。且$h(p)=0$

一般地，如果可以找到合适的g去拟合（大概至少使得g(p)=f(p)?且g为积性函数），使得能构造出积性函数h使得$h(p)=0$，那么都是可以这样做的！

如果没有好的办法，$h(p^e)$理论上可以暴力计算。