程序師世界是廣大編程愛好者互助、分享、學習的平台,程序師世界有你更精彩!
首頁
編程語言
C語言|JAVA編程
Python編程
網頁編程
ASP編程|PHP編程
JSP編程
數據庫知識
MYSQL數據庫|SqlServer數據庫
Oracle數據庫|DB2數據庫
 程式師世界 >> 編程語言 >> C語言 >> C++ >> C++入門知識 >> POJ 1845 Sumdiv (因子和)

POJ 1845 Sumdiv (因子和)

編輯:C++入門知識

POJ 1845 Sumdiv (因子和)


Sumdiv Time Limit: 1000MS   Memory Limit: 30000K Total Submissions: 15404   Accepted: 3800

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
這個題和HDU的1452基本上一樣,只不過這裡是HDU的推廣。
先說這個題的解題步驟: ① 將A分解A=(p1^a1)*(p2^a2)*.....*(pk^ak),那麼A^B就是A^B=[p1^(a1*B)] * [p2^(a2*B)] *.....*[pk^(ak*B)] ② 根據約數和的公式求和S並mod9901 別急,下面還有一些細節問題:
這個題我一開始就是先分解,然後按約數和的公式算.但是一直WA,我改了許多地方都不行。 如果一個數n=(p1^a1)*(p2^a2)*.....*(pk^ak) 約數和s=(p1^(a1+1)-1)/(p1-1) * (p2^(a2+1)-1)/(p2-1) *.....* (pk^(ak+1)-1)/(pk-1) 我就是按這個求s的公式算的,但是一直不對,肯定是因為這個公式中涉及到了除法,雖然在模運算中出現除法可以借助逆元來解決,但是因為求逆元是有條件限制的(a在模m意義下存在逆元的充要條件是a和m互素),我想是肯定是在處理這個地方的時候一直沒處理好,所以一直WA。
其實是我忽略了另一個公式,我上面給的求s公式是變形之後的,變形之前是: s=(1+p1+p1^2+p1^3+...p1^a1) * (1+p2+p2^2+p2^3+….p2^a2) * (1+p3+ p3^3+…+ p3^a3) * .... * (1+pk+pk^2+pk^3+...pk^ak) 這個公式是不涉及除法的,用起來會很舒服。(而對這個公式每個括號裡都是一個等比數列,分別對其求和就是上面的那個帶除法的公式)。 至此問題就明了了。
那麼就拋開這個題看下面的知識: 類似於這樣的一個問題:S(k)=A^0+A^1+A^2+....+A^k的快速求和。 方法:二分+遞歸求解。 下面看個簡單例子: ① k=4為偶數 A^0 + A^1 + A^2 + A^3 + A^4 =(A^0 + A^1)+A^2 + A^3(A^0 + A^1)=(1+A^3)*[A^0 + A^1]+A^2=(1+A^(k/2+1))*S(k/2-1)+A^(n/2) ② k=5為奇數 A^0 + A^1 + A^2 + A^3 + A^4 + A^5
=[A^0 + A^1 + A^2] + A^3*[A^0 + A^1 + A^2] = (1+A^3)*[A^0 + A^1 + A^2] = (1+A^(k/2+1))*S(k/2) 上面式子中,藍色部分用快速冪,紅色部分繼續遞歸,這樣就可以快速求出S(k)了。 遞歸出口n==0 return 1;
代碼中divi[i][0]代表第i個素因子,divi[i][0]代表這個素數的次數。
#include 
#include 
#include 
#include 
using namespace std;
typedef __int64 ll;

const int MOD=9901;
ll divi[10000][2],tot;

ll quick_mod(ll a,ll b){		//a^b%MOD
	ll res=1;
	a%=MOD;
	while(b){
		if(b&1) res=(res*a)%MOD;
		b/=2;
		a=(a*a)%MOD;
	}
	return res;
}

ll cal(int p,int n){
	if(n==0) return 1;
	if(n&1)
		return (1+quick_mod(p,n/2+1))*cal(p,n/2)%MOD;
	else
		return ((1+quick_mod(p,n/2+1))*cal(p,n/2-1)+quick_mod(p,n/2))%MOD;
}

void pre_solve(ll n){
	ll i;
	tot=0;
	for(i=2;i*i<=n;){
		if(n%i==0){
			divi[tot][0]=i;
			divi[tot][1]=0;
			do{
				n/=i;divi[tot][1]++;
			}while(n%i==0);
			tot++;
		}
		if(i==2) i++;
		else i+=2;
	}
	if(n>1){
		divi[tot][0]=n;divi[tot][1]=1;tot++;
	}
}

int main()
{
	ll A,B,i,res;
	while(scanf("%I64d%I64d",&A,&B)!=EOF){
		if(A==0){			//別忘了特判
			printf("0\n");continue;
		}
		if(A==1 || B==0){
			printf("1\n");continue;
		}
		pre_solve(A);
		res=1;
		for(i=0;i


  1. 上一頁:
  2. 下一頁:
Copyright © 程式師世界 All Rights Reserved