SPOJ3871 GCDEX

MATWETU posted @ 2011年8月15日 22:27 in 未分类 with tags 欧拉函数 约数 筛法 公约数 数论 SPOJ 积性函数 , 1553 阅读

这题的描述对我这种语言能力残缺的人来说,复述起来太难:

 

3871. GCD Extreme
Problem code: GCDEX
Given the value of N, you will have to find the value of G. The meaning of G is given in the following code
 
G=0;
for(k=i;k< N;k++)
for(j=i+1;j<=N;j++)
{
G+=gcd(k,j);
}
 
/*Here gcd() is a function that finds the greatest common divisor of the two input numbers*/
 
Input
 
The input file contains at most 20000 lines of inputs. Each line contains an integer N (1Output for each line of input produce one line of output. This line contains the value of G for the corresponding N. The value of G will fit in a 64-bit signed integer.
 
实质上是求这么一个东西:   G(n)=\sum_{k=1}^{n}\sum_{j=1}^{n}gcd(k,j)-\sum_{i=1}^{n}i
求G(n)的话,可以先求T(n)=\sum_{k=1}^{n}gcd(k,n),然后累加一下...额累加时第i项还要减个i~
 
如果学过数论,应该有印象,小于n的数里gcd(i,n)=d的i一共有\phi(\frac{n}{d})个~这意味着:
 
T(n)=\sum_{d|n}\frac{n}{d}\phi(d);
 

现在证明,T(n)是一个积性函数

证明:

设 m , n 满足 (m,n)=1:

T(n)=\sum_{d_n|n}\frac{n}{d_n}\phi(d_n)

T(m)=\sum_{d_m|m}\frac{m}{d_m}\phi(d_m)

则:

T(m)T(n)=(\sum_{d_m|m}\frac{m}{d_m}\phi(d_m))(\sum_{d_n|n}\frac{n}{d_n}\phi(d_n))

=\sum_{d_m|m}\sum_{d_n|n}\frac{mn}{d_md_n}\phi(d_m)\phi(d_n)

由于(m,n)=1,故dn与dm互素,所有dn,dm的组合正好取遍mn的约数

T(m)T(n)=\sum_{d|mn}\frac{mn}{d}\phi(d)=T(mn)

到了这里,只需考虑n等于素数的幂时T(n)的值:

根据T(n)的定义,显然有:

T(p^k)=\begin{Bmatrix}2p-1 & k=1\\p\cdot T(p^{k-1})+\phi(p^k) & k>1\end{matrix}

有了这个性质,就有可以在O(n)时间内用这个  牛X的筛法  把表打出来,然后求和还是O(n),查询什么的显然O(1)咯lol...

代码:

 

#include	<stdlib.h>
#include	<stdio.h>
#include	<string.h>
#include	<math.h>
#define FF(i,b) for(i=0;i<(b);i++)
#define FOR(i,a,b) for(i=(a);i<=(b);i++)
#define ZERO(x) do{memset(x,0,sizeof(x));}while(0)
#define SWAP(a,b) do{int _tmp=a;a=b;b=_tmp;}while(0)
#define INFI 2000000000
#define INFF 1e100
#define maxn 1000000
typedef long long  LL;
int p[maxn+1],ph[maxn+1],pk[maxn+1],T[maxn+1];
LL G[maxn+1];
void sevie()
{
	memset(p,0,sizeof(p));
	int i,j;
	p[0]=ph[0]=0;ph[1]=1;T[1]=1;
	for(i=2;i<=maxn;i++)
	{
		if(!p[i]){p[++p[0]]=i;ph[i]=i-1;pk[i]=i;T[i]=(i<<1)-1;}
		for(j=1;p[j]<=maxn/i;j++)
		{
			p[p[j]*i]=1;
			if(i%p[j]){
				ph[i*p[j]]=ph[i]*p[j];
				pk[i*p[j]]=p[j];
				T[i*p[j]]=T[i]*T[p[j]];
			}
			else{
				ph[i*p[j]]=ph[i]*p[j];
				pk[i*p[j]]=pk[i]*p[j];
				if(i==pk[i])T[i*p[j]]=p[j]*T[i]+ph[i*p[j]];
				else T[i*p[j]]=T[i/pk[i]]*T[pk[i]*p[j]];
				break;
			}
		}
	}
	G[0]=0;
	for(i=1;i<=maxn;i++)G[i]=G[i-1]+T[i]-i;
}
int main ( int argc, char *argv[] )
{
	sevie();
	int n;
	while(scanf("%d",&n)!=EOF)
	{
		if(!n)break;
		printf("%lld\n",G[n]);
	}
	return EXIT_SUCCESS;
}

登录 *


loading captcha image...
(输入验证码)
or Ctrl+Enter