-
Notifications
You must be signed in to change notification settings - Fork 0
/
multiplicative.hpp
47 lines (42 loc) · 1.37 KB
/
multiplicative.hpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
#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
*/