POJ 1845 求a^b的约数和

题目大意就是给定a和b，求a^b的约数和

f(n) = sigma(d) [d|n] 

这个学过莫比乌斯反演之后很容易看出这是一个积性函数

那么f(a*b) = f(a)*f(b)  (gcd(a,b)=1)

那么这道题就可以将a分解为每一个素数的k次方，求出相对应的f(p^k),将每一个乘在一起就行了

因为每一个素数得到的都是只有唯一的素数因子，那么f(n) 又变成了求 a^0+a^1+a^2....+a^k的值了 （n=a^k）

这里因为要取模，所以我用的是矩阵快速幂求的，网上别人用的都是二分分治求等比数列，反正这两种方法都不需要逆元，都能过

但我最开始写的要逆元的方法过不了也不知道为什么，希望路过大神知道的提个醒

矩阵快速幂的矩阵构造两维，第一维是等比关系，第二维保存前面所有值的和 

{

a , a

0 , 1

}

 

分治求等比数列的和：

3: 用递归二分求等比数列1+pi+pi^2+pi^3+...+pi^n：

（1）若n为奇数,一共有偶数项，则：
      1 + p + p^2 + p^3 +...+ p^n

      = (1+p^(n/2+1)) + p * (1+p^(n/2+1)) +...+ p^(n/2) * (1+p^(n/2+1))
      = (1 + p + p^2 +...+ p^(n/2)) * (1 + p^(n/2+1))

上式红色加粗的前半部分恰好就是原式的一半，那么只需要不断递归二分求和就可以了，后半部分为幂次式，将在下面第4点讲述计算方法。

 

（2）若n为偶数,一共有奇数项,则:
      1 + p + p^2 + p^3 +...+ p^n

      = (1+p^(n/2+1)) + p * (1+p^(n/2+1)) +...+ p^(n/2-1) * (1+p^(n/2+1)) + p^(n/2)
      = (1 + p + p^2 +...+ p^(n/2-1)) * (1+p^(n/2+1)) + p^(n/2);

   上式红色加粗的前半部分恰好就是原式的一半，依然递归求解

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
using namespace std;
#define ll long long
const int MOD = 9901;

struct Matrix{
    int m[2][2];
    void init(){memset(m , 0 , sizeof(m)); m[0][0] = m[1][1] = 1;}
    Matrix operator*(const Matrix &p)const {
        Matrix ans ;
        for(int i=0 ; i<2 ; i++)
            for(int j=0 ; j<2 ; j++){
                ans.m[i][j] = 0;
                for(int k=0 ; k<2 ; k++)
                    ans.m[i][j] = (ans.m[i][j]+m[i][k]*p.m[k][j]%MOD)%MOD;
            }
        return ans;
    }
    void out(){
        cout<<m[0][0]<<" "<<m[0][1]<<endl<<m[1][0]<<" "<<m[1][1]<<endl;
    }
};

Matrix q_pow(Matrix a , int b)
{
    Matrix ans;
    ans.init();
    while(b){
        if(b&1) ans = ans*a;
        a = a*a ; b>>=1;
    }
    return ans;
}

int get(int a , int b)
{
    //a^0+a^1+a^2...+a^b
    Matrix p;
    p.m[0][0] = a%MOD , p.m[0][1] = a%MOD , p.m[1][0] = 0 , p.m[1][1] = 1;
    p = q_pow(p , b);
    return (p.m[0][1]+p.m[1][1])%MOD;
}

void fenjie(int a , int b)
{
    int mx = (int)sqrt(a+0.5) , ret = 1;
    for(int i=2 ; i<=mx ; i++){
        if(i>a) break;
        int cnt = 0;
        while(a%i==0){
            cnt+=b , a/=i;
        }
        if(cnt) ret = (ret*get(i , cnt))%MOD;
    }
    if(a>1) ret = (ret*get(a , b))%MOD;
    printf("%d\n" , ret);
}

int main() {
   // freopen("a.in" , "r" , stdin);
   // freopen("compare.txt" , "w" , stdout);
    int a , b;
    while(~scanf("%d%d" , &a , &b)){
        if(a == 0) puts("0");
        else fenjie(a , b);
    }
    return 0;
}

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
using namespace std;
#define ll long long
const int MOD = 9901;

int q_pow(int a , int b)
{
    a = a%MOD;
    int ret = 1;
    while(b){
        if(b&1) ret = ret*a%MOD;
        a = a*a%MOD ; b>>=1;
    }
    return ret;
}

int inv(int a){return q_pow(a , MOD-2);}

int sum(int a , int b)
{
    //a^0+a^1+a^2...+a^b
    if(b==0) return 1;
    if(a==0) return 0;
    if(b&1) return ((1+q_pow(a , b/2+1))%MOD*(sum(a , b/2)%MOD))%MOD;
    else return ((1+q_pow(a , b/2+1))%MOD*(sum(a , b/2-1)%MOD)%MOD+(q_pow(a , b/2)%MOD))%MOD;
}

void fenjie(int a , int b)
{
    int mx = (int)sqrt(a+0.5) , ret = 1;
    for(int i=2 ; i<=mx ; i++){
        if(i>a) break;
        int cnt = 0;
        while(a%i==0){
            cnt+=b , a/=i;
        }
        if(cnt) ret = (ret*sum(i , cnt))%MOD;
    }
    if(a>1) ret = (ret*sum(a , b))%MOD;
    printf("%d\n" , ret);
}

int main() {
   // freopen("a.in" , "r" , stdin);
   // freopen("compare.txt" , "w" , stdout);
    int a , b;
    while(~scanf("%d%d" , &a , &b)){
        if(a == 0) puts("0");
        else fenjie(a , b);
    }
    return 0;
}