Multiplicative Sum
(Math/multiplicative.hpp)

• View this file on GitHub
• Last update: 2023-01-16 20:41:46+09:00
• Include: #include "Math/multiplicative.hpp"

使い方

T MultiplicativeSum(ll N): 乗法的関数 $f$ について $\sum_{k=1}^N f(k)$ を計算する。テンプレートには

• T: 返り値の型
• T (*pe)(int p,int e): $f(p^e)$ を返す関数
• T (*psum)(ll x): $\sum_{p \leq x:\mbox{prime}} f(p)$ を返す関数 ( $x$ は $\lfloor N/d \rfloor$ の形に限られる)

T AdditiveSum(ll n): 上の $f$ が 加法的関数 である場合。

Depends on

• Prime Sieve (Math/sieve.hpp)

Verified with

• Verify/YUKI_1781.test.cpp

Code

#pragma once
#include "Math/sieve.hpp"


template<typename T,T (*pe)(int,int),T (*psum)(ll)>T MultiplicativeSum(ll N){
    ll SQ=sqrtl(N);
    auto ps=sieve(SQ);
    
    T ret=psum(N)+1;
    auto dfs=[&](auto& dfs,ll x,int i,int e,T cur,T pre)->void{
        T nxt=pre*pe(ps[i],e+1);
        ret+=cur*(psum(double(N)/x)-psum(ps[i]));
        ret+=nxt;
        ll L=sqrtl(double(N)/x);
        if(ps[i]<=L)dfs(dfs,x*ps[i],i,e+1,nxt,pre);
        rep(j,i+1,ps.size()){
            if(ps[j]>L)break;
            dfs(dfs,x*ps[j],j,1,cur*pe(ps[j],1),cur);
        }
    };
    rep(i,0,ps.size())dfs(dfs,ps[i],i,1,pe(ps[i],1),1);
    return ret;
}

template<typename T,T (*pe)(int,int),ll (*pcnt)(ll),T (*psum)(ll)>T AdditiveSum(ll N){
    ll SQ=sqrtl(N);
    auto ps=sieve(SQ);
    
    T ret=psum(N);
    auto dfs=[&](auto& dfs,ll x,int i,int e,T cur,T pre)->void{
        T nxt=pre+pe(ps[i],e+1);
        ret+=cur*(pcnt(double(N)/x)-pcnt(ps[i]))+(psum(double(N)/x)-psum(ps[i]));
        ret+=nxt;
        ll L=sqrtl(double(N)/x);
        if(ps[i]<=L)dfs(dfs,x*ps[i],i,e+1,nxt,pre);
        rep(j,i+1,ps.size()){
            if(ps[j]>L)break;
            dfs(dfs,x*ps[j],j,1,cur+pe(ps[j],1),cur);
        }
    };
    rep(i,0,ps.size())dfs(dfs,ps[i],i,1,pe(ps[i],1),0);
    return ret;
}

/**
 * @brief Multiplicative Sum
 * @docs docs/multiplicative.md
 */

#line 2 "Math/sieve.hpp"

template<int L=50101010>vector<int> sieve(int N){
    bitset<L> isp;
    int n,sq=ceil(sqrt(N));
    for(int z=1;z<=5;z+=4){
        for(int y=z;y<=sq;y+=6){
            for(int x=1;x<=sq and (n=4*x*x+y*y)<=N;++x){
                isp[n].flip();
            }
            for(int x=y+1;x<=sq and (n=3*x*x-y*y)<=N;x+=2){
                isp[n].flip();
            }
        }
    }
    for(int z=2;z<=4;z+=2){
        for(int y=z;y<=sq;y+=6){
            for (int x=1;x<=sq and (n=3*x*x+y*y)<=N;x+=2){
                isp[n].flip();
            }
            for(int x=y+1;x<=sq and (n=3*x*x-y*y)<=N;x+=2){
                isp[n].flip();
            }
        }
    }
    for(int y=3;y<=sq;y+=6){
        for(int z=1;z<=2;++z){
            for(int x=z;x<=sq and (n=4*x*x+y*y)<=N;x+=3){
                isp[n].flip();
            }
        }
    }
    for(int n=5;n<=sq;++n)if(isp[n]){
        for(int k=n*n;k<=N;k+=n*n){
            isp[k]=false;
        }
    }
    isp[2]=isp[3]=true;

    vector<int> ret;
    for(int i=2;i<=N;i++)if(isp[i]){
        ret.push_back(i);
    }
    return ret;
}

/**
 * @brief Prime Sieve
 */
#line 3 "Math/multiplicative.hpp"

template<typename T,T (*pe)(int,int),T (*psum)(ll)>T MultiplicativeSum(ll N){
    ll SQ=sqrtl(N);
    auto ps=sieve(SQ);
    
    T ret=psum(N)+1;
    auto dfs=[&](auto& dfs,ll x,int i,int e,T cur,T pre)->void{
        T nxt=pre*pe(ps[i],e+1);
        ret+=cur*(psum(double(N)/x)-psum(ps[i]));
        ret+=nxt;
        ll L=sqrtl(double(N)/x);
        if(ps[i]<=L)dfs(dfs,x*ps[i],i,e+1,nxt,pre);
        rep(j,i+1,ps.size()){
            if(ps[j]>L)break;
            dfs(dfs,x*ps[j],j,1,cur*pe(ps[j],1),cur);
        }
    };
    rep(i,0,ps.size())dfs(dfs,ps[i],i,1,pe(ps[i],1),1);
    return ret;
}

template<typename T,T (*pe)(int,int),ll (*pcnt)(ll),T (*psum)(ll)>T AdditiveSum(ll N){
    ll SQ=sqrtl(N);
    auto ps=sieve(SQ);
    
    T ret=psum(N);
    auto dfs=[&](auto& dfs,ll x,int i,int e,T cur,T pre)->void{
        T nxt=pre+pe(ps[i],e+1);
        ret+=cur*(pcnt(double(N)/x)-pcnt(ps[i]))+(psum(double(N)/x)-psum(ps[i]));
        ret+=nxt;
        ll L=sqrtl(double(N)/x);
        if(ps[i]<=L)dfs(dfs,x*ps[i],i,e+1,nxt,pre);
        rep(j,i+1,ps.size()){
            if(ps[j]>L)break;
            dfs(dfs,x*ps[j],j,1,cur+pe(ps[j],1),cur);
        }
    };
    rep(i,0,ps.size())dfs(dfs,ps[i],i,1,pe(ps[i],1),0);
    return ret;
}

/**
 * @brief Multiplicative Sum
 * @docs docs/multiplicative.md
 */

Back to top page