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}^n\sum_{l_2=1}^n\cdots\sum_{l_k=1}^ngcd(l_1,l_2,\cdots,l_k)^2 \]

求:

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

的值。

(其中,\(n\leqslant 10^9\),\(k\leqslant 10^{10^5}\))。

题解:

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

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

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

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

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

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

\[=\sum_{d=1}^nd^2\sum_{p=1}^{n/d}\mu(p)\lfloor\frac{n}{dp}\rfloor^k \]

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

\[=\sum_{T=1}^n\lfloor\frac nT\rfloor^k\sum_{d|T}d^2\mu(\frac Td) \]

再令

\[F(T)=\sum_{d|T}d^2\mu(\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}\right)^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(i\cdot p_j)=F(i)\cdot p_j^2\qquad (i\perp 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}^k\sum_{T=1}^n\lfloor\frac nT\rfloor^kF(T) \]

交换求和顺序

\[\sum_{T=1}^n\sum_{i=2}^k\lfloor\frac nT\rfloor^kF(T) \]

发现

\[\sum_{i=2}^k\lfloor\frac nT\rfloor^k \]

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

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

\[a^{\phi(m)}\equiv1\quad(\text{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 == '\n') 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\n", ans);
	}
	return 0;
}
posted @ 2020-11-10 19:07  Frank喵^_^  阅读(98)  评论(0编辑  收藏  举报