Sumdiv （ 求因子和+逆元）

题目

Consider two natural numbers A and B. Let S be the sum of all natural divisors of A^B. Determine S modulo 9901 (the rest of the division of S by 9901).

Input

The only line contains the two natural numbers A and B, (0 <= A,B <= 50000000)separated by blanks.

Output

The only line of the output will contain S modulo 9901.

题意

求a^b的因子数和%9901

具体分析

首先我们第一反应出

m为它的素因子个数

则就有

t为它的素因子个数（就是上一句的m）

重点来了，我们需要知道：

这个公式可以求n的因子和

那么，等比数列怎么求呢？

百度一下就可知道：

好的，现在似乎我们就可以求出结果了，我们只需要把n改成n*b即可，再算快速幂

但是，你认为vjudge题的心机有这么简单吗

首先，我们如果直接求先不管时间

除了还要模，可以得到正确答案吗？

这里，很多大佬都用了

二分求等比数列和

可是，我才不会用这么拗口的名字求呢

竟然有除法和模运算，用逆元不就行了？？

懒懒的我就是这样做的~~

信心满满地提交后：

然后，我就是这样的：

  

连极限数据都对了的，为什么呢？

我才不会告诉你们我是偷偷看了题解的呢

我们知道，只有

a*b%c=1 && gcd（a,c）= 1

时，才有逆元，但是题目并不能保证这个

万一有a是c的倍数呢

所以我们需要用到另一个公式：

    mod c = a mod b*c   / b

证明：设a/b=tm+x

                 a=tbm+bx

                 a mod(bm) = bx

                 a mod(bm) /b =x

               所以 (a/b)mod m = amod(bm)/b

本以为逆元掌握得差不多的我竟然不知道！！！

你以为这样就完了吗，当然不可能！

我说过，vjudge是把我们编程社坑惨了滴

所以，我们在做前面爆longlong了怎么办，如果n很大，确实有可能

所以我们还需要引进一个新知识：

快速乘

怎么做呢，其实二进制是满足乘法分配律的

举个栗子：

10 =（1010）2  ，4*13等于4*（1010）2 ，展开得到4*13 == 4*（1000+10）2

所以就有了：

ll qc( ll x , ll y , ll m ){
    ll sum = 0;
    while( y ){
        if( y % 2 ) sum = ( sum + x ) % m;
        x = ( x + x ) % m;
        y /= 2;
    }
    return sum;
}

这是logn的做法，当然还有o(1)的

ll qc(ll a, ll b, ll m){
    ll c = a * b - ( ll ) ( ( long double ) a * b / m + 0.5 ) * m;
    if( c < 0 )
        c += m;
    return c;
}

卡一次啊

代码

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
using namespace std;
const int MAXN = 7900;//可以算一下，只需要根号a的最大值以内的质数就行了
#define ll long long
const ll mod = 9901;
int pn;
ll a , b , prime[MAXN+3] , p[MAXN + 2];
bool vis[MAXN + 3];
ll sum[MAXN+2];
void pr(){//还是习惯用欧拉筛
    for( int i = 2 ; i <= MAXN ; i ++ ){
        if( !vis[i] ){
            prime[++pn] = i;
        }
        for( int j = 1 ;  j <= pn && 1ll * i * prime[j] <= MAXN ; j ++ ){
            vis[i*prime[j]] = 1;
            if( i % prime[j] == 0 )
                break;
        }
    }
}
ll qc(ll a, ll b, ll m){//快速乘
    ll c = a * b - ( ll ) ( ( long double ) a * b / m + 0.5 ) * m;
    if( c < 0 )
        c += m;
    return c;
}
ll qpow( ll x , ll y , ll m )//快速幂
{
    x %= m;
    ll sum = 1;
    while( y )
    {
        if( y % 2 )
            sum = qc( sum , x , m ) % m;
        x = qc( x , x , m ) % m;
        y /= 2;
    }
    return sum;
}
int main()
{
    pr();
    while( ~scanf( "%lld%lld" , &a , &b ) )
    {
        memset( sum , 0 , sizeof( sum ) );//习惯要养好
        memset( p , 0 , sizeof( p ) );
        ll cnt = 0;
        if( a == 0 )
        {
            printf( "0" );
            return 0;
        }//其实没啥用，防止万一
        for( int i = 1 ; prime[i] <= sqrt( 1.0 * a ) && i <= pn ; i++ )
        {//这里就是分解因数了
            if( a % prime[i] == 0 )
            {
                cnt ++;
                p[cnt] = prime[i];
                while( a % prime[i] == 0 )
                {
                    sum[cnt] ++;
                    a /= prime[i];
                }
            }
        }
        if( a > 1 ){//注意，这里不能打掉
            p[++cnt] = a;
            sum[cnt] ++;
        }
        ll ans = 1;
        for( int i = 1 ; i <= cnt ; i ++ ){
            ll m = ( p[i] - 1 ) * mod;//根据公式可得，我们把mod看做模数，p[i]-1就是b，而qpow( p[i] , sum[i] * b + 1 , m )-1就是a
            ans = ans * ( ( qpow( p[i] , sum[i] * b + 1 , m ) + m - 1 ) / ( p[i] - 1 ) ) % mod ;
        }
        printf( "%I64d\n" , ( ans + mod ) % mod );//防止模后的答案改变
    }
    return 0;
}

其实倒数第六行那个加上m-1的m可以不要