Solution -「洛谷 P5325」Min_25 筛

(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;
}