poj 1845 Sumdiv (同余定理,快速幂取余)

21 篇文章 0 订阅

链接:poj 1845

题意:求A^B的所有因子的和对9901取余后的值

如:2^3=8,8的因子有 1,2,4,8,所有和为15,取余后也是15

应用定理主要有三个:

(1)整数的唯一分解定理:

     任意正整数都有且只有一种方式写出其素因子的乘积表达式。

      A=(p1^k1)*(p2^k2)*(p3^k3)*....*(pn^kn)   其中pi均为素数

(2)约数和公式

     对于已经分解的整数A=(p1^k1)*(p2^k2)*(p3^k3)*....*(pn^kn)

    有A的所有因子之和为

    S = (1+p1+p1^2+p1^3+...p1^k1) * (1+p2+p2^2+p2^3+….p2^k2)

      * (1+p3+ p3^3+…+ p3^k3) * .... * (1+pn+pn^2+pn^3+...pn^kn)

(3)同余定理

     (a+b)%m=(a%m+b%m)%m

       (a*b)%m=(a%m*b%m)%m

A^B的所有约数之和为:

   sum = [1+p1+p1^2+...+p1^(k1*B)] * [1+p2+p2^2+...+p2^(k2*B)] *...

         * [1+pn+pn^2+...+pn^(kn*B)]

必须得先求A的素因子及其个数,再求和

但是A,B都到5000万了,直接求和肯定会超时,我们再来看

递归二分求等比数列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))

 

(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);

由以上式子可知,前半部分恰好就是原式的一半,可以递归二分求解

其中p^k暴力求解也会超时,可以用快速幂取余

 

#include<stdio.h>
#include<string.h>
#include<math.h>
#define M 10000
#define N 9901
int p[M],cnt[M];
__int64 PowMod(__int64 a,__int64 b)  //a^b%N快速幂取余
{
    __int64 ans=1;
    a%=N;
    while(b){
        if(b%2)
            ans=ans*a%N;
        b/=2;
        a=a*a%N;
    }
    return ans;
}
__int64 sum(int p,int n)  //递归二分求和
{
    if(n==0)
        return 1;
    if(n%2)
        return sum(p,n/2)*(1+PowMod(p,n/2+1))%N;  //奇数时
    return (sum(p,n/2-1)*(1+PowMod(p,n/2+1))%N+PowMod(p,n/2))%N;  //偶数时
}
int main()
{
    int a,b,i,k;
    __int64 ans;
    while(scanf("%d%d",&a,&b)!=EOF){
        memset(cnt,0,sizeof(cnt));
        k=1;
        for(i=2;i*i<=a;i++){   //计算a的素因子和个数
            if(a%i==0){
                p[k]=i;
                while(a%i==0){
                    cnt[k]++;
                    a/=i;
                }
                k++;
            }
        }
        if(a!=1){   //记住最后判断一下
            p[k]=a;
            cnt[k++]++;
        }
        ans=1;
        for(i=1;i<k;i++)
            ans=ans*(sum(p[i],cnt[i]*b)%N)%N;
        printf("%I64d\n",ans);
    }
    return 0;
}

 

 

  • 0
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
以下是Java解决POJ3233—矩阵序列问题的代码和解释: ```java import java.util.Scanner; public class Main { static int n, k, m; static int[][] A, E; public static void main(String[] args) { Scanner sc = new Scanner(System.in); n = sc.nextInt(); k = sc.nextInt(); m = sc.nextInt(); A = new int[n][n]; E = new int[n][n]; for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) { A[i][j] = sc.nextInt() % m; E[i][j] = (i == j) ? 1 : 0; } } int[][] res = matrixPow(A, k); int[][] ans = matrixAdd(res, E); printMatrix(ans); } // 矩阵乘法 public static int[][] matrixMul(int[][] a, int[][] b) { int[][] c = new int[n][n]; for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) { for (int k = 0; k < n; k++) { c[i][j] = (c[i][j] + a[i][k] * b[k][j]) % m; } } } return c; } // 矩阵快速 public static int[][] matrixPow(int[][] a, int b) { int[][] res = E; while (b > 0) { if ((b & 1) == 1) { res = matrixMul(res, a); } a = matrixMul(a, a); b >>= 1; } return res; } // 矩阵加法 public static int[][] matrixAdd(int[][] a, int[][] b) { int[][] c = new int[n][n]; for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) { c[i][j] = (a[i][j] + b[i][j]) % m; } } return c; } // 输出矩阵 public static void printMatrix(int[][] a) { for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) { System.out.print(a[i][j] + " "); } System.out.println(); } } } ``` 解释: 1. 首先读入输入的n、k、m和矩阵A,同时初始化单位矩阵E。 2. 然后调用matrixPow函数求出A的k次矩阵res。 3. 最后将res和E相加得到结果ans,并输出。 4. matrixMul函数实现矩阵乘法,matrixPow函数实现矩阵快速,matrixAdd函数实现矩阵加法,printMatrix函数实现输出矩阵。

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值