codeforce 17D Notepad

 

这题实际求的是\(a^b\ mod\ n\),但是\(gcd(a,n)\)[下面记作(a,n)]不一定等于1.......

(a,n)等于1时直接费马小定理能推倒...

(a,n)不为1时,有如下事实(当然互素时显然成立):

\[a^{2\phi (n)}\ mod\ n = a^{\phi (n)}\ mod\ n\]

即

\[a^{\phi (n)}(a^{\phi(n)}-1)\ mod \ n =0\]

简单证明如下:

1.如果\(a^{\phi (n)}\ mod \ n =0\),等式显然成立.

2.否则,命题等价于

\[\frac{a^{\phi (n)} }{(a^{\phi (n)},n)} (a^{\phi (n)}-1) \ mod \ \frac{n}{(a^{\phi (n)} ,n)} = 0\]

即:

\[a^{\phi (n)}\ mod \ \frac{n}{(a^{\phi (n)} ,n)} = 1\]

到这里就能使用费马小定理了~

回到题目上:\(a^b\ mod\ n\)$实际上分两种情况:

1.\(b<\phi (n)\),此时直接算

2.else,\(a^b\ mod\ n=a^{b\ mod\ \phi (n)+\phi (n)}\ mod\ n\)

对任意(?)a,n成立

 

 

未分类 Comments(28) 2011年11月25日 14:00

hud4048

核武的思路:http://hi.baidu.com/aekdycoin/blog/item/69e075f5e200bcf47709d769.html

参考完核武报告,没完全懂~翻了翻组合数学书,这实质上就是求元素可重复的圆排列--许多组合数学书里讲莫比乌斯反演时唯一的例题~

先求 r 种不同元素,放在圆的n个点上 , 有几种放法, 旋转后相同算一种.

做法是:设周期为为d,d|n,的方案数是Td(r) 则 每个这样的圆排列,对应着d个线排列

线排列总数是r^n,所以:

r^n = \sum_{d|n} d\cdot T_d(r)

莫比乌斯反演后:

n\cdot T_n(r) = \sum_{d|n}\mu(n/d)\cdot r^{d}

总的方案数M_n就是各种周期的和:

M_n(r) = \sum_{k|n}\frac{1}{k}\sum_{d|k}\mu(k/d)\cdot r^{d}

因  d|k|n  , 设k'=k/d . n'=n/d

M_n(r) = \sum_{\frac{k}{d}|\frac{n}{d}}\sum_{d|n}\frac{1}{k}\mu(k/d)\cdot r^{d}

M_n(r) = \sum_{d|n}\frac{r^d}{n}\sum_{{k^,}|{n^,}}\mu(k^,)\cdot \frac{n^,}{k^,}

注意到后面一部分和等于n'的欧拉函数,故Mn可以写成:

M_n(r) = \frac{1}{n}\sum_{d|n}\phi(n/d)\cdot r^d

问题基本解决,剩下就是核武报告里说差不多,把最小公约数是i的倍数的方法数求出,容斥原理搞搞.

最后的答案应该是:

ANS=\sum_{n=1}^{maxWeight} \mu(n)\cdot \frac{1}{n}\sum_{d|n}\phi(n/d)\cdot r_n^d

其中r_n表示重量是n的倍数的石子种数

mu(n)  和 phi(n) 都可以在预处理时筛出来,

n的分解只需要搞一遍就够了~

遗憾的是,还是很慢,board的第一名时空复杂度都比我写的好~ym

未分类 Comments(0) 2011年9月21日 13:01

SPOJ3871 GCDEX

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

 

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;
}

未分类 Comments(0) 2011年8月15日 22:27

牛逼的线性时间筛法

这个牛B的筛法是在这里看到的

复杂度是线性的，扩展性超好～

自己YY了几个可以扩展的内容，权当模板～

 

#define N 1024
//niubee linear sevie method
int phi[N],//欧拉函数
    p[N],//素数表
    ld[N],//最小素因子
    ldt[N],//最小素因子的次数
    min_pk[N],//if p==ld[n],min_pk[n]=1+p+p^2+...+p^ldt[n];
    sig[N],//约数和
    nd[N];//因子数
//依赖性:min_pk<-sig;ldt<-nd;
void sevie()
{
	memset(p,0,sizeof(p));
	int i,j;

	phi[0]=0;phi[1]=1;//phi
	ld[0]=0;ld[1]=1;//ld
	ldt[0]=0;ldt[1]=1;//ldt
	min_pk[0]=0;min_pk[1]=1;//min_pk
	sig[0]=0;sig[1]=1;//sig
	nd[0]=0;nd[1]=1;//nd
	for(i=2;i<=N;i++)
	{
		if(!p[i]){
			p[++p[0]]=i;
			phi[i]=i-1;
			ld[i]=i;
			ldt[i]=1;
			min_pk[i]=i+1;
			sig[i]=i+1;
			nd[i]=2;
		}
		for(j=1;p[j]<=N/i;j++)
		{
			p[p[j]*i]=1;
			ld[p[j]*i]=p[j];
			if(i%p[j]){
				phi[p[j]*i]=phi[p[j]]*phi[i];
				ldt[p[j]*i]=1;
				min_pk[p[j]*i]=p[j]+1;
				sig[p[j]*i]=sig[i]*sig[p[j]];
				nd[p[j]*i]=nd[p[j]]*nd[i];
			}
			else{
				phi[p[j]*i]=phi[i]*p[j];
				ldt[p[j]*i]=ldt[i]+1;
				min_pk[p[j]*i]=p[j]*min_pk[i]+1;
				sig[p[j]*i]=sig[i]/min_pk[i]*min_pk[p[j]*i];
				nd[p[j]*i]=nd[i]/(ldt[i]+1)*(ldt[p[j]*i]+1);
				break;
			}
		}
	}
}

ALGO Comments(0) 2011年8月14日 15:35

数论相关算法(未完)

辗转相除法最大公约数:

//要求a,b不同时为零
int gcd(int a,int b)
{
    if(!b)return a;
    else return gcd(b,a%b);
}

利用最大公约数求最小公倍数:

int lcm(int a,int b)
{
    if(!(a*b))return 0;
    else return a*b/gcd(a,b);
}

求$a^b\mod n$:

int modExp(int a,int b,int n)
{
    int t=1,y=(a%n+n)%n;
    while(b)
    {
            if(b%2==1)t=t*y%n;
            y=y*y%n;
            b/=2;
    }
    return t;
}

扩展的Euclid算法:

//返回a，b的最大公约数d，并使ax+by=d
int ExtEuclid(int a,int b,int *x,int *y)
{
    int d,tmp;
    if(b==0)
    {
            *x=1;
            *y=0;
            return a;
    }
    d=ExtEuclid(b,a%b,x,y);
    tmp=*x;
    *x=*y;
    *y=tmp-a/b*(*x);
    return d;
}

解线性同余方程$ax\equiv b(\mod n):

//返回最小的x...通解是X=x(mod n/d) 
int modularLinerEquation(int a,int b,int n)
{
    int x,y,d;
    d=ExtEuclid(a,n,&x,&y);
    if(b%d!=0)return -1; // no solution
    x=b/d*x;
    x=(x%(n/d)+(n/d))%(n/d);
    return x; 
}

用中国剩余定理解同余方程组$a\equiv b_i(\mod m_i)$:

//用中国剩余定理解线性同与方程组$a\equiv b_i(\mod m_i)$
int solModularEquations(int b[],int m[],int z)
{
    int M=1,i,y,result=0;
    for(i=0;i<z;i++)M*=m[i];
    for(i=0;i<z;i++)
    {
        y=modularLinerEquation(M/m[i],1,m[i]);
        result=(result+y*b[i]*M/m[i])%M;
    }
    return result;
}

Miller-Rabin测试:

//Miller-Rabin素性测试
//令大奇数$n=2^lm+1$,其中l非负整数,m为正奇数.
//若$b^m\equiv 1(\mod n)$或者$b^{2^jm}\equiv -1(\mod n)$ ,0<=j<=l-1
//则称n通过以b为基的miller-rabin测试 
int passMillerRabinTest(int n) 
{
    int l,m;
    m=n-1;
    l=0;
    while(m%2==0)
    {
        l++;
        m/=2;
    }
    
    b=rand()%(n-1)+1;
    
    if(modExp(b,m,n)==1)return 1;
    
    int i,k;
    k=m;
    for(i=0;i<l;i++)
    {
        if(modExp(b,k,n)==n-1)return 1;
        k*=2;
    }
    
    return 0;
}

 

ALGO Comments(0) 2010年2月03日 00:25