(mathcal{Description})
Link.
对于积性函数 (f(x))，有 (f(p^k)=p^k(p^k-1)~(pinmathbb P,kinmathbb N_+))。求 (sum_{i=1}^nf(i)mod(10^9+7))。
(nle10^{10})。
(mathcal{Solution})
Min_25 筛是不可能的。
Powerful Number 三步走咯！考虑素数点值：
[f(p)=p^2-p
]
那么令 (g=operatorname{id}cdotvarphi)（点乘号即数值相乘），就有 (g(p)=p^2-p)。积性函数的点乘亦为积性函数。
求 (g) 的前缀和，杜教筛基础操作，卷上一个 (operatorname{id})：
[egin{aligned}
lbrack(operatorname{id}cdotvarphi)staroperatorname{id}
brack(n)&=sum_{imid n}(operatorname{id}cdotvarphi)(i)cdotfrac{n}{i}\
&=sum_{imid n}nvarphi(i)\
&=n^2
end{aligned}
]
自然数平方和易求，丢到杜教筛的式子里，推导后得出
[S(n)=frac{n(n+1)(2n+1)}{6}-sum_{i=2}^niSleft(lfloorfrac{n}{i}
floor
ight)
]
其中 (S(n)) 即为 (sum_{i=1}^ng(i))。
求 (h(p^k))，可以用 Bell 级数推导。令 (F_p,G_p,H_p) 分别为 (f,g,h) 在某一素数 (p) 的 Bell 级数，则
[egin{cases}
F_p=operatorname{OGF}langle1,p(p-1),p^2(p^2-1),cdots
angle=frac{1}{1-p^2z}-frac{1}{1-pz}+1\
G_p=operatorname{OGF}langle1,p(p-1),p^3(p-1),cdots
angle=frac{1-pz}{1-p^2z}
end{cases}
]
应用“两函数 Bell 级数的乘法卷积”为“原函数 Dirichlet 卷积之 Bell 级数”的性质，得到
[egin{aligned}
H_p&=frac{F_p}{G_p}\
&=frac{frac{1}{1-p^2z}-frac{1}{1-pz}+1}{frac{1-pz}{1-p^2z}}\
&=frac{1-frac{1-p^2z}{1-pz}+1-p^2z}{1-pz}\
&=frac{1}{1-pz}-frac{1-p^2z}{(1-pz)^2}+frac{1-p^2z}{1-pz}\
end{aligned}
]
我们仅仅想求 (h(p^k))，即 ([z^k]H_p)，那么
[egin{aligned}
lbrack z^k
brack H_p&=[z^k]frac{1}{1-pz}-[z^k]frac{1-p^2z}{(1-pz)^2}-[z^k]frac{1-p^2z}{1-pz}\
&=p^k-[(k+1)p^k-kp^{k+1}]+(p^k-p^{k+1})\
&=(k-1)(p^{k+1}-p^k)
end{aligned}
]
最终，(mathcal O(n^{frac{2}{3}})) 就能求出答案啦。
(mathcal{Code})
/* Clearink */
#include <cstdio>
#include <unordered_map>
#define rep( i, l, r ) for ( int i = l, repEnd##i = r; i <= repEnd##i; ++i )
#define per( i, r, l ) for ( int i = r, repEnd##i = l; i >= repEnd##i; --i )
typedef long long LL;
const int MOD = 1e9 + 7, MAXSN = 1e7, INV2 = 500000004, INV6 = 166666668;
int pn, pr[MAXSN + 5], gs[MAXSN + 5], phi[MAXSN + 5];
bool npr[MAXSN + 5];
inline int mul( const long long a, const int b ) { return a * b % MOD; }
inline int sub( int a, const int b ) { return ( a -= b ) < 0 ? a + MOD : a; }
inline void subeq( int& a, const int b ) { ( a -= b ) < 0 && ( a += MOD ); }
inline int add( int a, const int b ) { return ( a += b ) < MOD ? a : a - MOD; }
inline void addeq( int& a, const int b ) { ( a += b ) >= MOD && ( a -= MOD ); }
inline void sieve() {
phi[1] = gs[1] = 1;
rep ( i, 2, MAXSN ) {
if ( !npr[i] ) phi[pr[++pn] = i] = i - 1;
for ( int j = 1, t; j <= pn && ( t = i * pr[j] ) <= MAXSN; ++j ) {
npr[t] = true;
if ( !( i % pr[j] ) ) { phi[t] = phi[i] * pr[j]; break; }
phi[t] = phi[i] * ( pr[j] - 1 );
}
gs[i] = add( gs[i - 1], mul( i, phi[i] ) );
}
}
inline int gSum( const LL n ) {
static std::unordered_map<LL, int> mem;
if ( n <= MAXSN ) return gs[n];
if ( mem.count( n ) ) return mem[n];
int ret = mul( n % MOD, mul( mul( ( n + 1 ) % MOD,
( n << 1 | 1 ) % MOD ), INV6 ) );
for ( LL l = 2, r; l <= n; l = r + 1 ) {
r = n / ( n / l );
subeq( ret, mul(
mul( mul( ( l + r ) % MOD, ( r - l + 1 ) % MOD ), INV2 ),
gSum( n / l ) ) );
}
return mem[n] = ret;
}
LL n;
inline int powerSum( const int pid, LL x, const LL v ) {
if ( !v ) return 0;
int ret = 0, p = pr[pid];
if ( pid == 1 || !( x % pr[pid - 1] ) ) ret = mul( v, gSum( n / x ) );
if ( pid > pn ) return ret;
if ( ( x *= p ) > n ) return ret;
if ( ( x *= p ) > n ) return ret;
LL pwr = 1ll * p * p;
if ( pid < pn ) addeq( ret, powerSum( pid + 1, x / pwr, v ) );
for ( int j = 2; x <= n; ++j, x *= p, pwr *= p ) {
addeq( ret, powerSum( pid + 1, x,
mul( v, mul( j - 1, pwr % MOD * ( p - 1 ) % MOD ) ) ) );
}
return ret;
}
int main() {
sieve();
scanf( "%lld", &n );
printf( "%d
", powerSum( 1, 1, 1 ) );
return 0;
}