1. 程式人生 > >POJ 1845(Sumdiv) 數論好題

POJ 1845(Sumdiv) 數論好題

Time Limit: 1000MS Memory Limit: 30000K
Total Submissions: 19733 Accepted: 4987

Description

Consider two natural numbers A and B. Let S be the sum of all natural divisors of A^B. Determine S modulo 9901 (the rest of the division of S by 9901).

Input

The only line contains the two natural numbers A and B, (0 <= A,B <= 50000000)separated by blanks.

Output

The only line of the output will contain S modulo 9901.

Sample Input

2 3

Sample Output

15

Hint

2^3 = 8. 
The natural divisors of 8 are: 1,2,4,8. Their sum is 15. 
15 modulo 9901 is 15 (that should be output). 

分析(轉自某部落格):

這是一道數論的好題,需要較好的數學基礎

題意: 給定A,B,求A^B的所有因數的和,再MOD 9901

這裡用到了數論當中相當一部分知識

a. 唯一分解定理

任何一個整數都可以分解為若干個素數的冪的乘積的形式

A = ( p1 ^ q1 * p2 ^ q2 * ..... * pn ^ qn ) p為素數

A^B = ( p1 ^ (q1*B) * p2 ^ (q2*B) * ..... * pn ^ (qn*B) )

b. 約數個數公式 約數和公式

ys( A ) = ( 1+q1)(1+q2)(1+q3).....(1+qn)

Sum( A ) = ( 1 + p1 + p1^2 + .... + p1^q1 ) * ( 1 + p2 + p2^2 + .... + p2^q2 ) * ...... * ( 1 + pn + pn^2 +.....+ pn^qn )

Sum( A^B ) = (1 + p1 + p1^2 + .... + p1^(q1*B) ) * ( 1 + p2 + p2^2 + .... + p2^(q2*B) ) * ...... * ( 1 + pn + pn^2 +.....+ pn^(qn*B) )

一個數首先我們對他進行質因數分解.由上面的唯一分解定理得到

A = ( p1^q1 ) * (p2^q2) * (p3^q3) * ......*(pn^qn)

對於每一個pi ,可以取的次數為 0 - qi, 由乘法公式

所以一個數的約數個數為ys( A )  = (1+q1)(1+q2)(1+q3)....(1+qn)  

約數和公式  

A = ( p1^q1 ) * (p2^q2) * (p3^q3) * ......*(pn^qn)

每一個約數,可以取0個到q1個,那麼就是這個約數對總和的貢獻....

然後,乘法原理。

c. 快速冪取模

二分的思想遞迴求解A^B%C

d. 簡單的模計算公式 

(a*b)%M = (a%M)*(b%M)%M

e. 遞迴求解等比數列前n項和 (或者是用乘法逆元來求,還不會,馬上補上)

 ( 1 + p1 + p1^2 + .... + p1^n )    分奇偶進行討論:

1.)    n 為 奇數 ( n & 1 ), 總共有偶數項 ,假設n == 5

 ( 1 + p1 + p1^2 + p1^3 + p1^4 + p1^5  ) = ( 1 + p1 + p1^2 + p1^3 (1 + p1 + p1^2 ) )

=( 1 + p1 + p1^2 + p^(5/2+1) (1 + p1 + p1^2) )

=( 1+ p^(5/2+1) ) * (1 + p1 + p1^2 ) 

轉化為一般式應該為

 ( 1 + p1 + p1^2 + .... + p1^n )

= ( 1 + p1 + p1^2 + .... + p1^(n/2) + p^(n/2+1) * (1 + p1 + p1^2 + ..... p1^(n/2)) )

=(1 + p^(n/2+1)) * (1 + p1 + p1^2 + ..... + p1^(n/2))

=(1 + fast_mod(p1,n/2+1) ) * sum(p1,n/2) 

2.)    n 為 偶數 總共有奇數項 ,假設n == 6

 ( 1 + p1 + p1^2 + p1^3 + p1^4  + p1^5 + p1^6 ) = ( 1 + p1 + p1^2 + p1^3 + p1^4 (1 + p1 + p1^2) )

=( 1 + p1 + p2 + p1^4( 1 + p1 + p2 ) + p1^3 )

=( (1 + p^4)( 1 + p1 + p2 ) + p1^3  )

轉化為一般式子為

 ( 1 + p1 + p1^2 + .... + p1^n )

 =( 1 + p1^(n/2+1) ) * ( 1 + p1 + p2 + ... + p^(n/2-1) ) + p1^(n/2)

=( 1 + fast_mod(p1,n/2+1) )* sum( p1, n/2-1 )  + fast_mod(p1,n/2)

f.乘法逆元

首先我們需要知道等比數列的前N項和 

然後除以(1-q)轉化為乘法逆元來做

這裡應該是求的前n+1項和。。看前面的公式就知道了

Sum( A^B ) = (1 + p1 + p1^2 + .... + p1^(q1*B) ) * ( 1 + p2 + p2^2 + .... + p2^(q2*B) ) * ...... * ( 1 + pn + pn^2 +.....+ pn^(qn*B) )

單拆開這一項來說 (1 + p1 + p1^2 + .... + p1^(q1*B) ) 

有B*q1 +1項, Sn上下兩邊去負號

Sn = (q^(B*q1 + 1) - 1) / ( q - 1)

對於乘法逆元:在mod m的操作下(即Zm中),

a存在乘法逆元當且僅當a與m互質。不定方程ab+mx=1的任意一組整數解(b,x),b就是a的乘法逆元。

所以這題還有一個特判 如果 num % MOD == 1, 那就是說明 (num-1)%MOD == 0

這個情況下,不能num-1沒有乘法逆元。要另外計算。。

具體怎麼計算,也是搞了老半天不懂,,後來還是要自己推。

如果 a % MOD == 1

那麼( 1 + a^1 + a^2 + a^3 + ... +a^n )%MOD

=( 1 + a-1 + 1 + (a-1 + 1)^2 + (a-1 + 1)^3 + ..... + (a-1 +1)^n ) % MOD

=( 1 + 1 + 1 +1 ..... + 1 + (a-1)*K ) % MOD    (這裡暫且不管K是多少,但是因為(a-1) % MOD == 0,所以可以直接全部約去 )

=( 1 + 1 + 1 + 1.... + 1 ) % MOD

= n + 1

程式碼:
#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <algorithm>
#include <cmath>
using namespace std;

typedef long long int ll;

const int mod = 9901;



struct Node{
	int Num;
	int Count;
};
Node factor[50000];
int Count_factor = -1; //合數分解N(1/4)的時間複雜度 
void Init(int N){
	int sum=0;
	for(int i=2; i<=N; i++){
		if(N%i==0){
			sum = 0;
			while(N%i==0){
				sum++;
				N /= i;
			}
			factor[++Count_factor].Num = i;
			factor[Count_factor].Count = sum;
		}	
	}
}


ll quick_pow(ll a,ll m){
	ll ans = 1;
	while(m){
		if(m%2)
			ans = (ans*a)%mod;
		a = (a*a)%mod;
		m /= 2;
	}
	return ans;
}


ll F(int a,int b){
	
	ll tmp1 = quick_pow(a,b+1);
	tmp1 = (tmp1-1+mod)%mod;
	ll tmp2 = quick_pow(a-1,mod-2);
	return (tmp1*tmp2)%mod;

}

int main(){
		
	int A,B;
	scanf("%d%d",&A,&B);
	if(A==0){
		printf("0\n");
		return 0;
	}
	Init(A);
	ll ans=1;
	for(int i=0; i<=Count_factor; i++){
		factor[i].Count *= B;
		if( (factor[i].Num-1)%mod==0 )
			ans = (ans*(factor[i].Count+1))%mod;
		else
			ans = (ans*F(factor[i].Num,factor[i].Count))%mod;
	}
	printf("%I64d\n",ans);
	
	return 0;
}