【POJ 1845】Sumdiv

首先可以把A质因数分解，然后用等比数列求和公式求出每个质因数可能的贡献。
现在问题来了，在模9901的情况下如何做除法。
一个显然的想法就是求出逆元，因为9901是一个质数，所以直接乘上逆元即可。
我一开始也是这么想的，但是疯狂WA啊。
后来网上找到了另一种方法：二分求等比数列和。大概的意思就是每次将这个等比数列分成两半，提出一个公因式，然后这个数列就缩短一半。如此递归，最后求解。这样做是对的。
两个做法拍了无数组数据，都找不到问题。
直到后来看到了这个博客：（果然自古正理出discuss）

http://blog.csdn.net/zhengnanlee/article/details/9886917

嗯对这个博客让我成功地睡了个好觉。。。
原因就在于这个模数9901太小了，的确比9901小的数都是和9901互质的，但是如果比9901大呢？这几组错的数据分子都是9901的倍数，这样求出来的逆元显然是错的。
经验教训！！！

错误的做法

#include<cmath>
#include<cstdio>
#include<vector>
#include <queue>
#include<cstring>
#include<iomanip>
#include<stdlib.h>
#include<iostream>
#include<algorithm>
#define ll long long
#define inf 1000000000
#define mo 9901
#define N 100000
#define fo(i,a,b) for(i=a;i<=b;i++)
#define fd(i,a,b) for(i=a;i>=b;i--)
using namespace std;
int i,tot;
ll x,y,p,q,res,num;
int prime[N],f[N];
void get_prime()
{
    int i,j;
    fo(i,2,N)
        {
            if (!f[i]) prime[++tot] = i; j = i;
            while (j <= N) {f[j] = 1; j += i;}
        }
}
ll quick_multi(ll x,ll p)
{
    ll res=1,t=x;
    while (p)
        {
            if (p&1) res = res * t % mo;
            t = (t * t); 
            t = t % mo;

            p /= 2;
        }
    return res;
}
void Exgcd(ll a,ll b,ll &x,ll &y) //ax+by=1
{
    if (b == 0) {x = 1; y = 0; return;}
    Exgcd(b,a%b,y,x); y = y - (a / b) * x;
}
ll inv(ll a,ll p)
{
    ll d,x,y; Exgcd(a,p,x,y);
    return (x + p) % p;
}
int main()
{
    get_prime();
    scanf("%lld%lld",&x,&y);
    res = 1;
    fo(i,1,tot)
        {
            if (prime[i] * prime[i] > x) break;
            if (x % prime[i]) continue;
            num = 0; while (!(x % prime[i])) {num++; x /= prime[i];}
            num *= y;
            p = quick_multi(prime[i],num+1); p = (p - 1 + mo) % mo;
            q = inv(prime[i]-1,mo);
            res = res * p * q % mo;
        }
    if (x > 1) 
        {
            num = y; 
            p = quick_multi(x,num+1); p = (p - 1 + mo - 1) % mo + 1;
            q = inv(x-1,mo);
            res = res * p * q % mo;
        }
    printf("%lld\n",res);
    return 0;   
}

正确的做法

#include<cmath>
#include<cstdio>
#include<vector>
#include <queue>
#include<cstring>
#include<iomanip>
#include<stdlib.h>
#include<iostream>
#include<algorithm>
#define ll long long
#define inf 1000000000
#define mo 9901
#define N 100000
#define fo(i,a,b) for(i=a;i<=b;i++)
#define fd(i,a,b) for(i=a;i>=b;i--)
using namespace std;
int x,y,i,num;
ll res;
ll pow(ll x,ll p)
{
    ll res=1,t=x;
    while (p)
        {
            if (p&1) res = res * t % mo;
            t = (t * t); 
            t = t % mo;

            p /= 2;
        }
    return res;
}
ll sum(ll p,ll n)
{
    if (n == 0) return 1;
    if (n & 1) return ((1+pow(p,n/2+1)) * sum(p,n/2)) % mo;
    return ((1+pow(p,n/2+1)) * sum(p,n/2-1) + pow(p,n/2)) % mo;  
}

int main()
{
    scanf("%d%d",&x,&y);
    res = 1;
    fo(i,2,x)
        {
            if (i * i > x) break;
            num = 0;
            while (!(x % i)) {x /= i; num++;}
            res = (res * sum(i,num*y)) % mo;
        }
    if (x-1) res = (res * sum(x,y)) % mo;
    printf("%lld",res);
    return 0;
}