2019南京网络赛 K Sum

Define function

fn(k)=n∑l1=1n∑l2=1...n∑lk=1(gcd(l1,l2,...,lk))2fn(k)=∑l1=1n∑l2=1n...∑lk=1n(gcd(l1,l2,...,lk))2.

Given n(1≤n≤109),k(2≤k≤10105)n(1≤n≤109),k(2≤k≤10105), please calculate

k∑i=2fn(i)∑i=2kfn(i)

modulo 109+7109+7.

Input

There are multiple test cases, the first line contains an integer T(1≤T≤10)T(1≤T≤10), denoting the number of test cases.

For the next TT lines, each line contains two integers n(1≤n≤109)n(1≤n≤109) and k(2≤k≤10105)k(2≤k≤10105).

Output

For each test case, print one line contains an integer, which is the value of ∑ki=2fn(i)∑i=2kfn(i).

Sample Input

2
2 2
100 3

Sample Output

7
4331084

题意:给定函数

[f_n(k)=sum_{l_1=1}^nsum_{l_2=1}^ncdotssum_{l_k=1}^ngcd(l_1,l_2,cdots,l_k)^2 ]

求:

[sum_{i=2}^kf_n(k) ]

的值。

(其中,(nleqslant 10^9),(kleqslant 10^{10^5}))。

题解:

究极变态数论缝合怪题目,我们一点一点的去剖析。

先搞搞(f_n(k)),发现

[f_n(k)=sum_{l_1=1}^nsum_{l_2=1}^ncdotssum_{l_k=1}^ngcd(l_1,l_2,cdots,l_k)^2 ]

[=sum_{l_1=1}^nsum_{l_2=1}^ncdotssum_{l_k=1}^nsum_{d|l_1,l_2,cdots,l_k}d^2[gcd(l_1,l_2,cdots,l_k)=d] ]

[=sum_{d=1}^nd^2sum_{l_1=1}^{n/d}sum_{l_2=1}^{n/d}cdotssum_{l_k=1}^{n/d}[gcd(l_1,l_2,cdots,l_k)=1] ]

[=sum_{d=1}^nd^2sum_{l_1=1}^{n/d}sum_{l_2=1}^{n/d}cdotssum_{l_k=1}^{n/d}sum_{p|l_1,l_2,cdots,l_k}mu(p) ]

[=sum_{d=1}^nd^2sum_{p=1}^{n/d}mu(p)lfloorfrac{n}{dp} floor^k ]

老套路了,我们令(T=dp),则:

[=sum_{T=1}^nlfloorfrac nT floor^ksum_{d|T}d^2mu(frac Td) ]

再令

[F(T)=sum_{d|T}d^2mu(frac Td) ]

显然(F)是一个积性函数,我们可以写成卷积的样子:

[F=id_2*mu (id_2(T)=T^2) ]

这个式子很重要,待会用杜教筛的时候会用到。

我们先去推(F)的线性筛递推公式,考虑

[F(p_i^{a_i})=sum_{d|p_i^{a_i}}mu(d)left(frac {p_i^{a_i}}{d} ight)^2 ]

发现只有当(d=1,p_i)时,函数值不为(0),其他的时候函数的值都为(0)

[F(p_i^{a_i})=p_i^{2a_i}-p_i^{2(a_i-1)} ]

[F(p_i^{a_i+1})=p_i^2F(p_i^{a_i}) ]

所以我们现在推出了函数(F)的线性筛递推式:

[F(icdot p_j)=F(i)cdot p_j^2qquad (iperp p_j) ]

初始条件

[F(p)=p^2-1 ]

但是我们这还是无法处理高于一亿的数据(线性筛的复杂度所致),所以对于更高的数据范围我们要用杜教筛,还是推一下杜教筛的式子,我们利用上面说的那个卷积式:

[F=id_2*mu ]

两边同时卷积上(I)得:

[F*I=id_2*mu*I=id_2*varepsilon=id_2 ]

所以我们选定

[h=id_2,g=I ]

这两个函数都可以快速计算,于是杜教筛部分也就成了。

但是,这个问题还没结束,我们仅仅是找到了(f_n(k))的计算方法,但是我们实际是要算

[sum_{i=2}^kf_n(i) ]

但其实也还好,我们再把式子重新写一遍

[sum_{i=2}^ksum_{T=1}^nlfloorfrac nT floor^kF(T) ]

交换求和顺序

[sum_{T=1}^nsum_{i=2}^klfloorfrac nT floor^kF(T) ]

发现

[sum_{i=2}^klfloorfrac nT floor^k ]

可以直接套用等比数列求和公式计算,但是要特别注意公比为1的情况,此时要化为等差数列求和公式计算,不再赘述。

另外由于我们在用等比数列求和公式的时候难免会碰到计算形如(q^k)的形式,其中(k)又特别大,大到需要用高精度处理,这时候我们就必须用欧拉定理去处理高精度幂取模,这里的欧拉定理指的是

[a^{phi(m)}equiv1quad( ext{mod} m) ]

我们一边输入一边取模就可以了。

代码:

#include<bits/stdc++.h>
#include<tr1/unordered_map>
using namespace std;
typedef long long ll;//simplify long long
typedef unsigned long long ull;
typedef double db;
typedef long double ldb;
#define inf 2147483647
#define pi 3.14159265358979
#define rep(i, l, r) for(int i = l; i <= r; i ++)
#define lop(i, r, l) for(int i = r; i >= l; i --)
#define step(i, l, r, __step) for(int i = l; i <= r; i += __step)
#define revp(i, r, l, __step) for(int i = r; i >= l; i -= __step)
#define reg regsiter
#define RI regsiter int
#define RL regsiter long long
#define rerep(i, l, r) for(regsiter int i = l; i <= r; i ++)
#define relop(i, r, l) for(regsiter int i = r; i >= l; i --)
#define i8 __int128
#define __int128 ll//don't forget delete it in contest!
inline int read()//fast read
{
	int ret = 0, sgn = 1;
	char chr = getchar();
	while(chr < '0' || chr > '9')
	{if(chr == '-') sgn = -1; chr = getchar();}
	while(chr >= '0' && chr <= '9')
	{ret = ret * 10 + chr - '0'; chr = getchar();}
	return ret * sgn;
}
i8 write(i8 x)//int128's output
{
	if(x < 0) putchar('-'), x = -x;
	if(x > 9) write(x / 10);
	putchar(x % 10 + '0');
}
const int N = 6e6 + 50;
const int M = 9e5 + 5;
int prime[N], cnt;
bool vis[N];
const ll mod = 1e9 + 7;
ll f[N], sumf[N];
tr1::unordered_map<ll, ll> wf;
inline ll ksm(ll x, ll y)
{
	ll ret = 1;
	while(y)
	{
		if(y & 1) ret *= x, ret %= mod;
		x *= x, x %= mod;
		y >>= 1;
	}
	return ret;
}
ll inv6;
void shai(int x)
{
	f[1] = 1;
	rep(i, 2, x)
	{
		if(!vis[i])
		{
			prime[++ cnt] = i;
			f[i] = (1ll * i * i - 1) % mod;
		}
		rep(j, 1, cnt)
		{
			if(i * prime[j] > x) break;
			vis[i * prime[j]] = 1;
			if(!(i % prime[j]))
			{
				f[i * prime[j]] = f[i] * prime[j] % mod * prime[j] % mod;
				break;
			}
			f[i * prime[j]] = f[i] * f[prime[j]] % mod;
		}
	}
	rep(i, 1, x)
		sumf[i] = sumf[i - 1] + f[i], sumf[i] %= mod;
}
ll djsf(int x)
{
	if(x <= N - 10) return sumf[x];
	if(wf[x]) return wf[x];
	ll ret = 1ll * x * (x + 1) % mod * (2ll * x + 1) % mod * inv6 % mod;
	for(int l = 2, r; l <= x; l = r + 1)
	{
		r = x / (x / l);
		ret += mod - djsf(x / l) * (r - l + 1) % mod;
		ret = (ret % mod + mod) % mod;
	}
	ret %= mod;
	return wf[x] = ret;
}
inline ll inv(ll x)
{
	return ksm(x, mod - 2);
}
ll K, P;
inline ll calc(ll x)
{
	if(x == 1) return (P - 1 + mod) % mod;
	else return x * (ksm(x, K) - x + mod) % mod * inv(x - 1) % mod;
}
int main()
{
	int T;
	T = read();
	inv6 = ksm(6, mod - 2);
	shai(N - 10);
	while(T --)
	{
		int n;
		scanf("%d", &n);
		char t = getchar();
		K = 0, P = 0;
		while(1)
		{
			t = getchar();
			if(t == '
') break;
			K = K * 10 + t - '0';
			P = P * 10 + t - '0';
			K %= mod - 1;
			P %= mod;
		}
		ll ans = 0, d, S;
		for(int l = 1, r; l <= n; l = r + 1)
		{
			r = n / (n / l);
			ans += (djsf(r) - djsf(l - 1) + mod) % mod * calc(n / l) % mod;
			ans %= mod;
		}
		printf("%lld
", ans);
	}
	return 0;
}
原文地址:https://www.cnblogs.com/loi-frank/p/13955465.html