您的位置:首页 > 其它

POJ1845 Sumdiv 费马小定理+欧拉函数+素因子分解

2013-12-25 23:56 387 查看
题目有点难,过题率只有百分之二十,今天做了一个下午,做不出停了一下,在做别的题目以后,回来速度A了,原来费马小定理那里写错了,9900写成了9901,有点粗心,

题目意思是求A,B,,求出A^B的所有约数之和S,然后输出S%9901;

开始分析:

对A进行素因子分解,pi为素数序列,A=p1^m1*p2^m2……pn^mn,那么 A^B=p1^(m1*B)*p2^(m2*B)……pn^(mn*B),在因子分解里有一个基本的理论:如果f是乘性函数,那么f的和函数F(n)=∑f(d),(其中   d|n   ),也是乘性函数

那么接下来利用乘性函数,所有因子和就是 S=(1+p1+p1^2+……+p1^(n*B))*(1+p2+p2^2……+p2^(n*B))……(1+pn+pn^2+……+pn^(n*B)),然后利用等比数列求和  S=(1+p^(n*B+1))/(p1-1)*(1+p2^(n*B+1))/(p2-1)*……(1+p^(n*B+1))/(pn-1),

求和是 求每项的分子是用二分快速幂法,处理分母时,由于9901是一个素数,则会时候乘法逆元的思想就出现了,乘以这个数然后模9901的逆元就可以了

#include<iostream>
#include<cstdio>
#include<list>
#include<algorithm>
#include<cstring>
#include<string>
#include<queue>
#include<stack>
#include<map>
#include<vector>
#include<cmath>
#include<memory.h>
#include<set>

#define ll long long

#define eps 1e-8

#define inf 0xfffffff
const ll INF = 1ll<<61;

using namespace std;

//vector<pair<int,int> > G;
//typedef pair<int,int > P;
//vector<pair<int,int> > ::iterator iter;
//
//map<ll,int >mp;
//map<ll,int >::iterator p;
//

int exgcd(int a,int &x,int b,int &y)
{
if(b==0)
{
x=1;
y=0;
return a;
}
int r=exgcd(b,x,a%b,y);
int t=x;
x=y;
y=t-a/b*y;
}

int multi(int a,int b)
{
a%=9901;
if(a==0)return 9900;
if(a==1)return 0;
int ans=1;
b%=9900;//9901为素数,所以可以使用费马小定理
while(b)
{
if(b&1)
ans=ans*a%9901;
a=a*a%9901;
b>>=1;
}
return ans-1;
}

int main()
{
int a,b;
int prime[102][2];
while(scanf("%d %d",&a,&b)==2)
{
if(a == 0)
{
puts("0");
continue;
}
if(a == 1 || b == 0)
{
puts("1");
continue;
}
int cnt=0;
for(int i=2;i*i<=a;i++)
{
if(a%i == 0)
{
prime[cnt][0]=i;
prime[cnt][1]=1;
a/=i;
while(a%i == 0)
{
a/=i;
prime[cnt][1]++;
}
cnt++;
}
}
if(a > 1)
{
prime[cnt][0]=a;
prime[cnt][1]=1;
cnt++;
}
for(int i=0;i<cnt;i++)
prime[i][1]*=b;
ll ans=1;
int x,y;
for(int i=0;i<cnt;i++)
{
if(prime[i][0]%9901 == 0)continue;
if(prime[i][0]%9901 == 1)
{
ans=ans*(prime[i][1]+1)%9901;
continue;
}
ans=ans*multi(prime[i][0],prime[i][1]+1)%9901;
exgcd(prime[i][0]-1,x,9901,y);
x=(x%9901+9901)%9901;
ans=ans*x%9901;
}
printf("%I64d\n",ans);
}
return EXIT_SUCCESS;
}
内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
标签: