POJ 1845 Sumdiv(因子和 大数模 费马小定理)

POJ 1845

题目大意

给你A和B(0≤A,B≤50000000)​,定义S为AB​的所有因子的和，让你输出S mod​ 9901

分析

AB可以写成(pa11pa22..pakk)B

因子的和S=∑(pa1b11pa2b22...pakbkk),1≤bi≤B

S=(p01+p11+...pa1∗B1)(p02+p12+...pa2∗B2)...(p0k+p1k+...pak∗Bk)

=pa1∗B+11−1p1−1⋅pa2∗B+12−1p2−1...pak∗B+1k−1pk−1

带有除法的取模运算可以用乘法逆元来解决，分子对幂运算取模可以用二分法取模法做。

一个优化是通过费马小定理把幂次方缩小到9901以内：

xy中的y=y%9901;

费马小定理:ap−1≡1(modp)(a和p互素)

​

这道题思路不难但涉及的知识点多：扩展欧几里得、欧拉函数、费马小定理、大数取模、同余

思路想好后就在看书上代码做了（细节问题太多）。需要注意的一个细节是：在扩展欧几里得之后需要x=(x%MOD+MOD)%MOD; 防止x出现负数的情况

代码

#include<cstdio>#include<iostream>#include<cmath>#include<cstring>#include<cstdlib>#include<queue>using namespace std;const int MOD=9901;int p[50000];//保存A的因子int pnum[50000];//A的因子对应的个数int pcnt;void ExtGcd(int a,int b,int &x,int &y){    if(!b){x=1;y=0;return ;}    ExtGcd(b,a%b,x,y);    int t=x;x=y;    y=t-(a/b)*y;}void Find_Factor(int A)//将A分解为素数幂{      memset(pnum,0,sizeof(pnum));      pcnt=0;      for(int i=2;i*i<=A;i++)      {            if(A%i==0)            {                  p[++pcnt]=i;                  while(A%i==0)                  {                       pnum[pcnt]++;                       A=A/i;                  }            }      }      if(A>1)      {           p[++pcnt]=A;           pnum[pcnt]=1;      }}int Mul(int x,int y)//返回x^y%MOD{      x=x%MOD;      if(x==0)return MOD;      if(x==1)return 1;      int m=1;      y=y%(MOD-1);//费马小定理优化      while(y)//二进制取模      {            if(y&1)m=m*x%MOD;            x=x*x%MOD;            y=y>>1;      }      return m;}void Work(int A,int B){      long long int ans=1;      int x,y;      for(int i=1;i<=pcnt;i++)      {            if(p[i]%MOD==1){ans=ans*(pnum[i]*B+1)%MOD;continue;}            ans=ans*(Mul(p[i],pnum[i]*B+1)-1)%MOD;            ExtGcd(p[i]-1,MOD,x,y);//x即为p[i]-1的逆(%MOD)            x=(x%MOD+MOD)%MOD;            ans=ans*x%MOD;      }      cout<<ans<<endl;}int main(){     int A,B;     while(scanf("%d%d",&A,&B)!=EOF)     {          if(A==0){cout<<0<<endl;continue;}          if(A==1 || B==0){cout<<1<<endl;continue;}          Find_Factor(A);          Work(A,B);     }}