bzoj-3529 数表

141 篇文章 0 订阅
88 篇文章 0 订阅

题意:

令x的约数之和为g[x];

多组数据,输入n,m,a;

求∑g[gcd(i,j)],g[gcd(i,j)]<=a,1<=i<=n,1<=j<=m;


题解:

首先对于取值在非负整数集合的数论函数,有以下结论成立;

此定理即莫比乌斯反演定理;

有了这个结论,本题中我们定义两个函数;

f [n]为gcd(i,j)==n的(i,j)对数;

F[n]为n|gcd(i,j)的(i,j)对数;

显然F[n]=∑f[d],n|d;

那么莫比乌斯反演的②式成立;

而倘若不考虑题中小于等于a的影响,那么答案为:

∑ f[d]*g[d],1<=d<=min(n,m);//g[d]为d的约数和;

带入∑ g[d]*∑ μ(x/d)*(n/x)*(m/x);

1<=d<=min(n,m),d|x;

化简∑ (n/d)*(m/d)∑ μ(d/x)*g[x];

1<=d<=min(n,m),x|d;

(如果是单个询问)将所有满足g[x]<=a的∑ μ(d/x)*g[x]维护一个对d的前缀和;

预处理线性筛出一些函数之后;

可以在O(nlogn)的时间内解决(logn是在插入时一个d要由多个x更新而来造成的,这样均摊为log);

但是多组询问在线的O(Q*nlogn)无法解决问题;

所以为了防止多次维护这个前缀和,将询问离线并按a排序;

排序之后将当前g[x]<=a的东西全部插入,自然也要将g[x]排序;

插入用树状数组

然后这样再处理所有询问,复杂度在O(nlog^2n+Q√nlogn)的级别;

故可以解决此题;


代码:


#include<stdio.h>
#include<string.h>
#include<algorithm>
#define N 100100
#define M 22000
#define pr pair<int,int>
using namespace std;
typedef long long ll;
struct node
{
	int n,m,k,no;
}a[M];
int pri[N>>1],ans[M],tot;
int miu[N],popo[N],qqq[N],sum[N];
pr g[N];
bool vis[N];
int cmp(node a,node b)
{
	return a.k<b.k;
}
int lowbit(int x)
{
	return x&(-x);
}
void update(int x,int val)
{
	while(x<N)
	{
		sum[x]+=val;
		x+=lowbit(x);
	}
}
int query(int x)
{
	int ret=0;
	while(x)
	{
		ret+=sum[x];
		x-=lowbit(x);
	}
	return ret;
}
void init()
{
	miu[1]=g[1].first=g[1].second=popo[1]=qqq[1]=1;
	for(int i=2;i<N;i++)
	{
		if(!vis[i])
		{
			pri[++tot]=i;		//素数表
			miu[i]=-1;			//莫比乌斯函数
			popo[i]=i;			//最小质因子的最大p^t满足p^t|i
			qqq[i]=i+1;			//∑p^t
			g[i].first=i+1;		//i的约数和
			g[i].second=i;		//id(i),就是i本身
		}
		for(int j=1;j<=tot&&i*pri[j]<N;j++)
		{
			vis[i*pri[j]]=1;
			g[i*pri[j]].second=i*pri[j];
			if(i%pri[j]==0)
			{
				miu[i*pri[j]]=0;
				popo[i*pri[j]]=popo[i]*pri[j];
				qqq[i*pri[j]]=qqq[i]+popo[i*pri[j]];
				g[i*pri[j]].first=g[i].first/qqq[i]*qqq[i*pri[j]];
				break;
			}
			miu[i*pri[j]]=-miu[i];
			popo[i*pri[j]]=pri[j];
			qqq[i*pri[j]]=pri[j]+1;
			g[i*pri[j]].first=g[i].first*(pri[j]+1);
		}
	}
}
int Inv(int n,int m)
{
	int ret=0,i,last;
	for(i=1;i<=n;i=last+1)
	{
		last=min(n/(n/i),m/(m/i));
		ret+=(n/i)*(m/i)*(query(last)-query(i-1));
	}
	return ret&2147483647;
}
int main()
{
	init();
	int T,n,m,i,j,k,p;
	scanf("%d",&T);
	for(i=1;i<=T;i++)
	{
		scanf("%d%d%d",&a[i].n,&a[i].m,&a[i].k);
		a[i].no=i;
		if(a[i].n>a[i].m)	swap(a[i].n,a[i].m);
	}
	sort(a+1,a+1+T,cmp);
	sort(g+1,g+N);
	for(i=p=1;i<=T;i++)
	{
		while(g[p].first<=a[i].k)
		{
			for(j=1;g[p].second*j<N;j++)
				update(g[p].second*j,g[p].first*miu[j]);
			p++;
		}
		ans[a[i].no]=Inv(a[i].n,a[i].m);
	}
	for(i=1;i<=T;i++)
		printf("%d\n",ans[i]);
	return 0;
}





  • 0
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值