HDU 6706 huntian oy

数学算法解析:求解特定条件下的数论问题

HDU  6706  huntian oyHDU\ \ 6706\ \ huntian\ oyHDU  6706  huntian oy

http://acm.hdu.edu.cn/showproblem.php?pid=6706

哎,开始就觉得,对于,i>j,gcd(i,j)=1,i>j,gcd(i,j)=1,i>j,gcd(i,j)=1

gcd(ia−ja,ib−jb)=igcd(a,b)−jgcd(a,b)gcd(i^a-j^a,i^b-j^b)=i^{gcd(a,b)}-j^{gcd(a,b)}gcd(iaja,ibjb)=igcd(a,b)jgcd(a,b)

手残,好久没刷题了 ,忘了很多东西 哈哈
借着题目又复习了一下。
其实对于:

(ia−ja)=(i−j)∑k=0a−1ia−1(ji)k(i^a-j^a)=(i-j)\sum_{k=0}^{a-1}i^{a-1}\Big(\frac{j}{i}\Big)^k(iaja)=(ij)k=0a1ia1(ij)k

令b>a,b=ta+r令b>a, b=ta+rb>a,b=ta+r

(ib−jb)=(ita+r−jta+r)=(i−j)∑k=0ta+r−1ita+r−1(ji)k=(i−j)∑d=0t−1i(t−1−d)a+rjda∑k=0a−1ia−1(ji)k+(i−j)jta∑k=0r−1ir−1(ji)k={∑d=0t−1i(t−1−d)a+rjda}(ia−ja)+jta(ir−jr)(i^b-j^b)=(i^{ta+r}-j^{ta+r})=(i-j)\sum_{k=0}^{ta+r-1}i^{ta+r-1}\Big(\frac{j}{i}\Big)^k\\ =(i-j)\sum_{d=0}^{t-1}i^{^{(t-1-d)a+r}}j^{da}\sum_{k=0}^{a-1}i^{a-1}\Big(\frac{j}{i}\Big)^k+(i-j)j^{ta}\sum_{k=0}^{r-1}i^{r-1}\Big(\frac{j}{i}\Big)^k\\=\Big\{\sum_{d=0}^{t-1}i^{^{(t-1-d)a+r}}j^{da}\Big\}(i^a-j^a)+j^{ta}(i^r-j^r)(ibjb)=(ita+rjta+r)=(ij)k=0ta+r1ita+r1(ij)k=(ij)d=0t1i(t1d)a+rjdak=0a1ia1(ij)k+(ij)jtak=0r1ir1(ij)k={d=0t1i(t1d)a+rjda}(iaja)+jta(irjr)
又因为gcd(i,j)=1gcd(i,j)=1gcd(i,j)=1.
(ib−jb)% (ia−ja)=(ib%a−jb%a)(i^b-j^b)\%\ (i^a-j^a)=(i^{b\%a}-j^{b\%a})(ibjb)% (iaja)=(ib%ajb%a)

所以:

ans=∑i=1n∑j=1i(i−j)[gcd(i,j)=1]ans = \sum_{i=1}^{n}\sum_{j=1}^i(i-j)[gcd(i,j)=1]ans=i=1nj=1i(ij)[gcd(i,j)=1]
ans=∑i=1niφ(i)−∑i=1n1+iφ(i)2=∑i=2niφ(i)2ans = \sum_{i=1}^ni\varphi(i)-\sum_{i=1}^n\frac{1+i\varphi(i)}{2}\\=\sum_{i=2}^n\frac{i\varphi(i)}{2}ans=i=1niφ(i)i=1n21+iφ(i)=i=2n2iφ(i)

这种计算方便统计小于10610^6106的答案。

自己用了一个很蠢的方法。。。。。。。。

如果前面使用莫比乌斯反演的话有:

令:F(n)=∑i=1n∑j=1i(i−j)F(n)=\sum_{i=1}^{n}\sum_{j=1}^i(i-j)F(n)=i=1nj=1i(ij)
ans=∑k=1nμ(k)k∑i=1⌊nk⌋∑j=1i(i−j)=∑k=1nμ(k)kF(⌊nk⌋)ans =\sum_{k=1}^n\mu(k)k\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{j=1}^i(i-j)=\sum_{k=1}^n\mu(k)kF(\Big\lfloor\frac{n}{k}\Big\rfloor)ans=k=1nμ(k)ki=1knj=1i(ij)=k=1nμ(k)kF(kn)

这个姿态有点慢。需要用前面的方法打标。后面这种方法需要计算 mu(k)kmu(k)kmu(k)k前缀和MMM,然后查表。

ans=∑L=1mF(L)(M(⌊nL⌋−M(⌊nL+1⌋))+∑k=1⌊nm+1⌋μ(k)kF(⌊nk⌋)ans = \sum_{L=1}^{m}F(L)(M(\Big\lfloor\frac{n}{L}\Big\rfloor-M(\Big\lfloor\frac{n}{L+1}\Big\rfloor))+\sum_{k=1}^{\lfloor\frac{n}{m+1}\rfloor}\mu(k)kF(\Big\lfloor\frac{n}{k}\Big\rfloor)ans=L=1mF(L)(M(LnM(L+1n))+k=1m+1nμ(k)kF(kn)
#pragma GCC diagnostic error "-std=c++11"
#pragma GCC target("avx")
#pragma GCC optimize(3)
#pragma GCC optimize("Ofast")
#pragma GCC optimize("inline")
#pragma GCC optimize("-fgcse")
#pragma GCC optimize("-fgcse-lm")
#pragma GCC optimize("-fipa-sra")
#pragma GCC optimize("-ftree-pre")
#pragma GCC optimize("-ftree-vrp")
#pragma GCC optimize("-fpeephole2")
#pragma GCC optimize("-ffast-math")
#pragma GCC optimize("-fsched-spec")
#pragma GCC optimize("unroll-loops")
#pragma GCC optimize("-falign-jumps")
#pragma GCC optimize("-falign-loops")
#pragma GCC optimize("-falign-labels")
#pragma GCC optimize("-fdevirtualize")
#pragma GCC optimize("-fcaller-saves")
#pragma GCC optimize("-fcrossjumping")
#pragma GCC optimize("-fthread-jumps")
#pragma GCC optimize("-funroll-loops")
#pragma GCC optimize("-fwhole-program")
#pragma GCC optimize("-freorder-blocks")
#pragma GCC optimize("-fschedule-insns")
#pragma GCC optimize("inline-functions")
#pragma GCC optimize("-ftree-tail-merge")
#pragma GCC optimize("-fschedule-insns2")
#pragma GCC optimize("-fstrict-aliasing")
#pragma GCC optimize("-fstrict-overflow")
#pragma GCC optimize("-falign-functions")
#pragma GCC optimize("-fcse-skip-blocks")
#pragma GCC optimize("-fcse-follow-jumps")
#pragma GCC optimize("-fsched-interblock")
#pragma GCC optimize("-fpartial-inlining")
#pragma GCC optimize("no-stack-protector")
#pragma GCC optimize("-freorder-functions")
#pragma GCC optimize("-findirect-inlining")
#pragma GCC optimize("-fhoist-adjacent-loads")
#pragma GCC optimize("-frerun-cse-after-loop")
#pragma GCC optimize("inline-small-functions")
#pragma GCC optimize("-finline-small-functions")
#pragma GCC optimize("-ftree-switch-conversion")
#pragma GCC optimize("-foptimize-sibling-calls")
#pragma GCC optimize("-fexpensive-optimizations")
#pragma GCC optimize("-funsafe-loop-optimizations")
#pragma GCC optimize("inline-functions-called-once")
#pragma GCC optimize("-fdelete-null-pointer-checks")
#pragma GCC optimize(3)
#include <stdio.h>
#include <cmath>
#define MAXN 1000006
using namespace std;
typedef long long LL;

int mod = 1e9 + 7;
int mu[MAXN + 5] = { 0,-1 };
int M[MAXN + 5];
int tmp[MAXN + 5];
int In[100];

inline int Sn(int n)
{
	if (n & 1)
		return (LL)n*((n + 1) >> 1) % mod;
	else
		return (LL)(n >> 1)*(n + 1) % mod;
}

void clat(int n, int d)
{
	int ans = 0;
	int m = (int)sqrt(n) + 1;
	for (int L = 1; L < m; L++)
	{
		ans += (LL)M[L] * (Sn(n / L) - Sn(n / (L + 1)) + mod) % mod;
		if (ans >= mod)ans -= mod;
	}
	for (int i = (int)(n / m); i > 1; i--)
	{
		int u = n / i;
		if (u < MAXN)
			ans += (LL)i*M[u] % mod;
		else
			ans += (LL)i*tmp[i*d] % mod;
		if (ans >= mod)ans -= mod;
	}
	tmp[d] = (1 - ans + mod) % mod;
}

int slove(int n)
{
	if (n < MAXN)return M[n];
	for (int i = (n / (MAXN - 1)); i; i--) clat(n / i, i);
	return tmp[1];
}

int phi[MAXN + 10];

void init()
{
	for (int i = 1; i <= MAXN; i++)  phi[i] = i;
	for (int i = 2; i <= MAXN; i += 2) phi[i] /= 2;
	for (int i = 3; i <= MAXN; i += 2)
		if (phi[i] == i) {
			for (int j = i; j <= MAXN; j += i) {
				phi[j] = phi[j] / i*(i - 1);
			}
		}
}

inline int F(int n)
{
	LL tm = LL(n)*(n + 1) % mod;
	int ret = (tm*(2 * n + 1) % mod)*In[12] % mod;
	ret -= tm*In[4] % mod;
	if (ret < 0)ret += mod;
	return ret;
}

int S[MAXN + 5];

inline int getM(int u, int d)// u = n/d
{
	if (u < MAXN)return M[u];
	return tmp[d];
}

int main()
{
	
	init();
	for (int i = 1; i < MAXN; i++)
	{
		mu[i] = -mu[i];
		for (int j = (i << 1); j < MAXN; j += i)mu[j] += mu[i];
		if (mu[i] < 0)mu[i] += mod;
		M[i] = M[i - 1] + LL(mu[i]) * i%mod;
		if (M[i] >= mod)M[i] -= mod;
	}
	In[1] = 1;
	for (int i = 2; i < 100; i++)
	{
		In[i] = LL(mod - (mod / i))*In[mod%i] % mod;
	}
	for (int i = 2; i < MAXN; i++)
	{
		S[i] = S[i - 1] + ((LL)phi[i] * i%mod)*In[2] % mod;
		if (S[i] >= mod)S[i] -= mod;
	}

	int T = 10;

	scanf("%d", &T);

	while (T--)
	{
		int n, a, b;
		scanf("%d %d %d", &n, &a, &b);
		if (n < MAXN)
		{
			printf("%d\n", S[n]);
			continue;
		}
		slove(n);
		int m = sqrt(n) + 1;
		int ans = 0;
		for (int L = 1; L < m; L++)
		{
			int mm = getM(n / L, L) - getM(n / (L + 1), L + 1);
			if (mm < 0)mm += mod;
			ans += (LL)F(L)*mm%mod;
			if (ans >= mod)ans %= mod;
		}
		int size = n / m;
		for (int k = 1; k <= size; k++)
		{
			ans += ((LL)k*mu[k] % mod) * F(n / k) % mod;
			if (ans >= mod)ans -= mod;
		}
		printf("%d\n", ans);
	}
	return 0;
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值