POJ 1845（Sumdiv）数论好题

最新推荐文章于 2021-10-20 20:51:18 发布

spring371327

最新推荐文章于 2021-10-20 20:51:18 发布

阅读量738

点赞数

分类专栏：数论文章标签：费马定理求逆元合数分解因子和快速幂取模 POJ 1845

本文链接：https://blog.csdn.net/spring371327/article/details/52725077

版权

数论专栏收录该内容

17 篇文章 0 订阅

订阅专栏

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

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 +１项， 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;
}