【noi.ac #2030】数学题

题目

http://noi.ac/problem/2030

思路

如果你没有发现标算1的神奇结论,也可以直接上套路来做
显然可以枚举g,对于每个g计算P(g) = ∑{i=1…n, j=1…m, gcd(i,j)=g} floor(i/j)(g=1…m)
欲计算P,仅需要对每个g计算Q然后反演得到,其中Q(d) = ∑{i=1…n, d|i, j=1…m, d|j} floor(i/j)(d=1…m)
注意m最大有10^7,O(m log m)可能会超时,建议用O(m log log m)的反演
而Q(d) = S(floor(n/d), floor(m/d)),回到了标算1中求解的问题
总复杂度为O(n^0.25 (n^0.5 + m^0.5) + m log log m)

代码

#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
const int P = 1e7;
const int PR = 1e6;
int n, m;
int np[P + 10], pr[PR], cnt = 0;
ll phi[P + 10];
/*void update(int u)
{
	int ls=u<<1,rs=u<<1|1;
	if(tr[ls].id) tr[u].id=1,tr[u].v=max(tr[ls].v,tr[u].v);
	if(tr[rs].id) tr[u].id=1,tr[u].v=max(tr[rs].v,tr[u].v);
	if(tr[u].id==0) tr[u].bz=0,tr[u].v=0;
}
void pushdown(int u,int st,int ed)
{
	int ls=u<<1,rs=u<<1|1;
	if(tr[ls].id) tr[ls].bz+=tr[u].bz,tr[ls].v+=tr[u].bz;
	if(tr[rs].id) tr[rs].bz+=tr[u].bz,tr[rs].v+=tr[u].bz;
	tr[u].bz=0;
}
void add(int t)
{
	for(int i=t; i<=n; i+=i&-i) T[i]+=1;
}
int total(int t)
{
	int ret=0;
	for(int i=t; i; i-=i&-i) ret+=T[i];
	return ret;
}
void ins(int u,int st,int ed,int x)
{
	if(st==ed)
	{
		tr[u].bz=0; tr[u].v=0; tr[u].id=0; return;
	}
	pushdown(u,st,ed);
	int mid=st+ed>>1;
	if(mid>=x) ins(u<<1,st,mid,x);else ins(u<<1|1,mid+1,ed,x);
	update(x);
}
void modify(int u,int st,int ed,int l,int r,int v)
{
	if(l<=st&&ed<=r)
	{
		tr[u].bz+=v; tr[u].v+=v; return;
	}
	pushdown(u,st,ed);
	int mid=st+ed>>1;
	if(mid>=l) modify(u<<1,st,mid,l,r,v);
	if(mid<r) modify(u<<1|1,mid+1,ed,l,r,v);
	update(u);
}
int query(int u,int st,int ed,int l,int r)
{
	if(l<=st&&ed<=r)
	{
		return tr[u].v;
	}
	pushdown(u,st,ed);
	int mid=st+ed>>1,aii=0;
	if(mid>=l) aii=max(aii,query(u<<1,st,mid,l,r));
	if(mid<r) aii=max(aii,query(u<<1|1,mid+1,ed,l,r));
	return aii;
}*/
void sieve() {
	phi[1] = 1;
	for (int i = 2; i <= P; i++) {
		if (!np[i]) pr[++cnt] = i, phi[i] = i - 1;
		for (int j = 1; j <= cnt && (ll)i * pr[j] <= P; j++) {
			np[i * pr[j]] = 1;
			if (i % pr[j] == 0) {
				phi[i * pr[j]] = phi[i] * pr[j];
				break;
			} else phi[i * pr[j]] = phi[i] * (pr[j] - 1);
		}
	}
	for (int i = 1; i <= P; i++) phi[i] += phi[i - 1];
}

ll calc(int n, int m) {
	ll ret = 0;
	for (int j = 1, r; j <= m; j = r + 1) {
		if (j > n) break;
		r = min(m, n / (n / j));
		ll s = (ll)(r + j) * (r - j + 1) / 2;
		ll x = n / j, p = x * (x + 1) / 2;
		ret += (ll)(n + 1) * (r - j + 1) * x - s * p;
	}
	return ret;
}

void put128(__int128 x) {
	if (x >= 10) put128(x / 10);
	putchar('0' + (x % 10));
}

int main() {
	sieve();
	scanf("%d%d", &n, &m);
	ll ans = 0;
	for (int l = 1, r; l <= m; l = r + 1) {
		r = min(n / (n / l), m / (m / l));
		ans += (phi[r] - phi[l - 1]) * calc(n / l, m / l);
		ans %= 1LL << 60;
	}
	printf("%lld\n", ans);
	return 0;
}
  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值