1. 程式人生 > >尤拉函式程式設計實現

尤拉函式程式設計實現

 尤拉函式phi(x)是指不大於正整數x的與x互質的正整數的個數。例如phi(1)=1,phi(2)=1,phi(3)=2,phi(4)=2,phi(5)=4,phi(6)=2等等。很顯然,對每一個質數p,phi(p)=p-1。而對每一個質數的冪phi(p^n)=(p-1)×p^(n-1)。尤拉函式是積性函式,如果m、n互質,那麼phi(m×n)=phi(m)×phi(n)。
    尤拉函式的取值需要用到算術基本定理。將n寫成其質因子冪的連乘積的形式,n=p1^k1×p2^k2×…×pr^kr。那麼,phi(n)=n×(1-1/p1)×…×(1-1/pr)。利用這個公式可以寫出求尤拉函式的程式碼。當然,首先要求質數。如果問題範圍是MAXSIZE,那麼只需求出SIZE=sqrt(MAXSIZE)以內的質數即可。即保證任何一個數的最小質因子一定在這其中。

  1. int euler(int n){  
  2.     int ret = n;  
  3.     for(int i=0;P[i]*P[i]<=n;++i){ //P是儲存質數的陣列
  4.         if ( n % P[i] ) continue;  
  5.         ret -= ret / P[i];//計算尤拉函式
  6.         n /= P[i];  
  7.         while( 0 == n % P[i] ) n /= P[i];  
  8.         if ( 1 == n ) break;//這三行用於加速並有助於最後得到大的質因子
  9.     }  
  10.     if ( n != 1 ) ret -= ret / n;  //假如n含有超過SIZE的質因子
  11.     return ret;  
  12. }  


    利用尤拉函式的積性性質,可以將一定範圍內的尤拉函式值全部求出。任何一個合數n都可以寫作n=k×p^r。其中p是一個質數,通常就會取n的最小的質因子。此時可以保證k與p^r肯定是互質的。所以phi(n)=phi(k)×phi(p^r)=phi(k)×(p-1)×p^(r-1)。當r大於1的時候,該式可以寫作phi(n)=phi(k)×(p-1)×p^(r-2)×p=phi(k)×phi(p^(r-1))×p=phi(n/p)×p。而當r等於1的時候,明顯有phi(n)=phi(k)×(p-1)=phi(n/p)×(p-1)。因為是順序求值,所以在求n的尤拉函式值時,所有比n小的尤拉函式值均已求出。當然,這種方法同樣要求事先篩出質數。
  1. int Euler[SIZE] = {0};  
  2. void euler(){  
  3.     Euler[1] = 1;  
  4.     for(int i=2;i<SIZE;++i){  
  5.         if ( !isComp[i] ){//isComp[i]標識i是否為合數,可以由篩法得到
  6.             Euler[i] = i - 1;  
  7.             continue;  
  8.         }  
  9.         for(int j=0;j<PCnt;++j){  
  10.             if ( i % P[j] ) continue;  
  11.             if ( i / P[j] % P[j] ) Euler[i] = Euler[i/P[j]] * ( P[j] - 1 );  
  12.             else Euler[i] = Euler[i/P[j]] * P[j];  
  13.             break;//上述式子只需計算一次即可
  14.         }  
  15.     }      
  16. }  


    上述原則求尤拉函式只需要知道n的其中一個質因子即可,一般而言就用最小質因子。所以配合線性時間素數篩法,可以直接線上性時間內求出範圍內的所有質數以及每個數的最小質因子,計算尤拉函式值也是線性的。

自己的補充 ,覺得方法二的可讀性有點低,一些地方沒說清楚,我自己按照思路寫了如下方法二的程式碼

int Euler_Fomuler2(int n)
{
	int i,ret;
	ret = 1;
	for(i = 0;P[i] * P[i] < n;i++)
	{	
		if(n % P[i]) continue;
		while(n % P[i] == 0)
		{
			ret *= P[i];
			n /= P[i];
		}
		ret -= ret/P[i];		
	}
	if(n != 1)
	{
		ret *= (n-1);
	}
	return ret;
}

以下是生成P的程式碼

void make_p()
{
	
	int t[300];
	memset(t,0,sizeof(t));
	int i,j;
	t[0] = 1;
	t[1] = 1;
	for(i = 2;i < sqrt(double(300));i++)
	{
		j = i*2;
		while(j < 300)
		{
			t[j] = 1;
			j += i;
		}
	}
	for(i = 2;i < 300;i++)
	{
		if(t[i] == 0)
		{
			P[pos++] = i;
		}
	}
}