一本通1633【例 3】Sumdiv
1633:【例 3】Sumdiv
时间限制: 1000 ms 内存限制: 524288 KB【题目描述】
原题来自:Romania OI 2002
求 ABAB 的所有约数之和 mod9901。
【输入】
输入两个整数 A,B。
【输出】
输出答案 mod9901。
【输入样例】
2 3
【输出样例】
15
【提示】
样例说明
23=8,8 的所有约数为 1,2,4,8,1+2+4+8=15,15mod9901=15,因此输出 15。
数据范围与提示:
对于全部数据,0≤A,B≤5×107。
sol:写了半天几乎吐血终于水过这道板子题
首先有个很显然的东西
对于一个数
如果他是p1a1*p1a2*p3a3*~~*pnan
那么他的因数和就是 (p10+p11+p12+...+p1a1)*(p20+p21+...+p2a2)*...*(pn1+pn2+...+pnan)
于是这道题就是要把所有质因数的幂次的乘积,先来推一下等比公式
令 S = p0+p1+p2+...+pa, (1)
则p*S = p1+p2+...+pa+pa+1, (2)
(2)-(1)可以得到 S*(p-1) = pa+1-p0,所以S= (pa+1-p0)/ (p-1) ,前面一项显然Ksm就可以解决,考虑后一项
需要用到p-1的逆元,如何求呢??
定义一个数a(就求a的逆元)
对于任意一个数 X/a ≡ X*Inva (%Mod)
易知 a*Inva ≡ 1 (%Mod)
随便推导下
原式: a*Inva = 1 (%Mod)
--> a*Inva-k*Mod = 1 (%Mod)
--> a*Inva+P*Mod = 1 (%Mod) (类ax+by=c的形式)
配上巨丑无比的代码
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
/* 原式: a*Inva = 1 (%Mod) --> a*Inva-k*Mod = 1 (%Mod) --> a*Inva+P*Mod = 1 (%Mod) (类ax+by=c的形式) */ #include <bits/stdc++.h> using namespace std; typedef long long ll; inline ll read() { ll s=0; bool f=0; char ch=' '; while(!isdigit(ch)) { f|=(ch=='-'); ch=getchar(); } while(isdigit(ch)) { s=(s<<3)+(s<<1)+(ch^48); ch=getchar(); } return (f)?(-s):(s); } #define R(x) x=read() inline void write(ll x) { if(x<0) { putchar('-'); x=-x; } if(x<10) { putchar(x+'0'); return; } write(x/10); putchar((x%10)+'0'); return; } #define W(x) write(x),putchar(' ') #define Wl(x) write(x),putchar('\n') const ll Mod=9901; ll A,B; int cnt=0; struct Prime { ll Shuz,Xis; }Prim[50]; inline ll Ksm(ll x,ll y) { ll ans=1; while(y) { if(y&1) { ans=ans*x%Mod; } x=x*x%Mod; y>>=1; } return ans; } ll gcd(ll x,ll y) { return (!y)?(x):(gcd(y,x%y)); } inline void Exgcd(ll a,ll b,ll &X,ll &Y) { if(b==0) { X=1; Y=0; return; } Exgcd(b,a%b,X,Y); ll XX=X,YY=Y; X=YY; Y=XX-a/b*YY; return; } //a*Inva+Mod*k = 1 (%Mod) inline ll Inv(ll Num) { ll a,b,c,r,X,Y; a=Num; b=Mod; c=1; r=gcd(a,b); Exgcd(a,b,X=0,Y=0); X=X*r/c; ll tmp=b/r; X=(X>=0)?(X%tmp):(X%tmp+tmp); // printf("X=%lld Test=%lld\n",X,Num*X%Mod); return X; } inline void Solve(ll Num) { int i; for(i=2;i<=sqrt(Num);i++) { if(Num%i==0) { Prim[++cnt].Shuz=i; Prim[cnt].Xis=0; while(Num%i==0) { Prim[cnt].Xis++; Num/=i; } Prim[cnt].Xis*=B; } } if(Num>1) { Prim[++cnt]=(Prime){Num,B}; } return; } int main() { int i; ll ans=1ll; R(A); R(B); Solve(A); for(i=1;i<=cnt;i++) { // printf("%lld %lld\n",Prim[i].Shuz,Prim[i].Xis); ll tmp=(Ksm(Prim[i].Shuz,Prim[i].Xis+1)-1+Mod)%Mod; tmp=tmp*Inv(Prim[i].Shuz-1)%Mod; ans=ans*tmp%Mod; } Wl(ans); return 0; } /* input 2 3 output 15 input 10 6 output 5187 input 45279441 39876543 output 6060 input 9 8 output 5660 */
河田は河田、赤木は赤木……。
私は誰ですか。教えてください、私は誰ですか。
そうだ、俺はあきらめない男、三井寿だ!