P1593 因子和

P1593 因子和
新算法：
#define ni 逆元
先质因数分解，
（1+p1^1+p1^2...p1^x）*(1+p2^1+p2^2...p2^x)
然后套等比数列公式就可以了。

  1 #include<iostream>
  2 #include<cstdio>
  3 #include<queue>
  4 #include<algorithm>
  5 #include<cmath>
  6 #include<ctime>
  7 #include<set>
  8 #define mod 9901
  9 #include<cstring>
 10 #define inf long long_MAX
 11 #define For(i,a,b) for(register long long i=a;i<=b;i++)
 12 #define p(a) putchar(a)
 13 #define g() getchar()
 14 //by war
 15 //2017.10.27
 16 using namespace std;
 17 long long num,b;
 18 struct node
 19 {
 20     long long cnt;
 21     long long p;
 22 }a[10000];//40M
 23 long long prime[10000];//40M
 24 long long cnt;
 25 long long tot;
 26 bool vis[10000];
 27 long long ans;
 28 
 29 void in(long long &x)
 30 {
 31     long long y=1;
 32     char c=g();x=0;
 33     while(c<'0'||c>'9')
 34     {
 35     if(c=='-')
 36     y=-1;
 37     c=g();
 38     }
 39     while(c<='9'&&c>='0')x=x*10+c-'0',c=g();
 40     x*=y;
 41 }
 42 void o(long long x)
 43 {
 44     if(x<0)
 45     {
 46         p('-');
 47         x=-x;
 48     }
 49     if(x>9)o(x/10);
 50     p(x%10+'0');
 51 }
 52 
 53 void Euler(long long x)
 54 {
 55     For(i,2,x)
 56     {
 57         if(!vis[i])prime[++cnt]=i;
 58         for(register long long j=1;j<=cnt&&prime[j]*i<=x;j++)
 59         {
 60             vis[prime[j]*i]=true;
 61             if(i%prime[j]==0)
 62             break;
 63         }
 64     }
 65 }
 66 
 67 void resolve(long long x)
 68 {
 69     For(i,1,cnt)
 70     {
 71         if(x%prime[i]==0)
 72         {
 73             a[++tot].p=prime[i];
 74             while(x%prime[i]==0)
 75             {
 76             a[tot].cnt++;
 77             x/=prime[i];
 78             }
 79         }
 80     }
 81     if(x>1)
 82     {
 83     a[++tot].p=x;
 84     a[tot].cnt++;    
 85     }
 86 }
 87 
 88 long long ksm(long long a,long long b)
 89 {
 90     if(b==0)
 91     return 1;
 92     while(b%2==0)
 93     {
 94         a=(a*a)%mod;
 95         b>>=1;
 96     }
 97     long long r=1;
 98     while(b>0)
 99     {
100         if(b%2==1)
101         r=(r*a)%mod;
102         a=(a*a)%mod;
103         b>>=1;
104     }
105     return r%mod;
106 }
107 
108 void exgcd(long long a,long long b,long long &x,long long &y)
109 {
110     if(!b)
111     {
112         x=1;
113         y=0;
114         return;
115     }
116     exgcd(b,a%b,x,y);
117     long long temp=x;
118     x=y;
119     y=temp-(a/b)*y;
120 }
121 
122 long long series(long long q,long long n)
123 {
124     long long fz=ksm(q,n)-1;
125     long long x,y,b;
126     exgcd(q-1,mod,x,y);
127     long long ni=x;
128     ni=(ni%mod+mod)%mod;
129     return (fz*ni%mod+mod)%mod;
130 }
131 
132 int main()
133 {
134     in(num),in(b);
135     Euler(sqrt(num));
136     resolve(num);
137     ans=1;
138     For(i,1,tot)
139     ans=ans*series(a[i].p,a[i].cnt*b+1)%mod;
140     o(ans%mod);
141     return 0;
142 }

 

posted @ 2017-10-27 19:39  WeiAR  阅读(335)  评论(0)  编辑  收藏  举报