等比数列的二分求和

ACM

qsc视频中一道水题引发的一系列的知识

等比二分求和

一个很不错讲解原理的博客，我就直接搬运了
二分求解等比数列的和

二分求解等比数列的和

注意这个等比式的第一项为a, 不是1.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int mod = 1e9+7;
ll pow_mod(ll a, ll n){
    ll ans = 1;
    while(n){
        if(n&1){
            ans = ans*a%mod;
        }
        a = a*a;
        n>>=1;
    }
    return ans;
}
ll sum(ll a, ll n){
    if(n == 1) return a;
    ll half_sum = sum(a, n/2);
    ll ans;
    if(n&1){
        ll a1 = pow(a, n/2+1);
        ans = (half_sum+half_sum*a1%mod)%mod;
        ans = (ans+a1)%mod;
    }
    else{
        ll a1 = pow(a, n/2);
        ans = (half_sum+half_sum*a1%mod)%mod;
    }
    return ans;
}
ll a, n;
int main()
{
    while(~scanf("%I64d%I64d", &a, &n)){
        printf("%I64d\n", sum(a, n));
    }
    return 0;
}

二分求解等比矩阵的和

poj3233 Matrix Power Series
本题也可以采用构造的方式进行求解，但是我已经忘了怎么构造递推的矩阵了，我好菜啊。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
typedef long long ll;
const int maxn = 35;
int n, k, m;
struct Mat{
    ll mat[maxn][maxn];//以后一定用ll,不适时的进行取模就会导致错误。
    Mat(){
        memset(mat, 0, sizeof(mat));
   }
};
Mat mul(Mat a, Mat b, int n, int mod){
    Mat ans;
    for(int i=0; i<n; i++){
        for(int k=0; k<n; k++){
            for(int j=0; j<n; j++){
                ans.mat[i][j]+=(a.mat[i][k]*b.mat[k][j])%mod;
            }
        }
    }
    return ans;
}
Mat pow(Mat a, int b, int n, int mod){
    Mat ans;
    for(int i=0; i<n; i++) ans.mat[i][i] = 1;
    while(b){
        if(b&1){
            ans = mul(ans, a, n, mod);
        }
        a = mul(a, a, n, mod);
        b >>= 1;
    }
    return ans;
}
Mat add(Mat a, Mat b, int n, int mod){
    Mat ans;
    for(int i=0; i<n; i++){
        for(int j=0; j<n; j++){
            ans.mat[i][j] = (a.mat[i][j]+b.mat[i][j])%mod;
        }
    }
    return ans;
}
Mat sum(Mat a, int k, int n, int mod){
    if(k == 1) return a;
    Mat half_sum = sum(a, k/2, n, mod);
    Mat ans;
    if(k&1){
        Mat a1 = pow(a, k/2+1, n, mod);
        ans = add(half_sum, mul(half_sum, a1, n, mod), n, mod);
        ans = add(ans, a1, n, mod);
    }
    else{
        Mat a1 = pow(a, k/2, n, mod);
        ans = add(half_sum, mul(half_sum, a1, n, mod), n, mod);
    }
    return ans;
}
int main()
{
    while(~scanf("%d%d%d",&n, &k, &m)){
        Mat ma;
        for(int i=0; i<n; i++){
            for(int j=0; j<n; j++){
                scanf("%d", &ma.mat[i][j]);
            }
        }
        Mat ans = sum(ma, k, n, m);
        for(int i=0; i<n; i++){
            int j=0;
            for(j=0; j<n-1; j++){
                printf("%d ", ans.mat[i][j]);
            }
            printf("%d\n", ans.mat[i][j]);
        }
    }
    return 0;
}

题目

前面的都是为这做铺垫
poj1845 Subdiv

题意

题解

  1. 先用唯一分解定理分解该数
  2. 约数和定理

    第一个式子可以用二分求和来继续的求解，第二个化简的分式的形式可以用逆元来求解。
  3. 然后使用二分求等比数列的和(或者使用逆元来求解)

当然这一道题还可以用逆元来做，具体等以后的实现。
逆元的实现

若枚举5e7之间的所有的素数，那么结果就是T.
我觉的其中的缩小数字的搜索范围那一块处理的很奇妙，可以积累一下这个技巧。

1
2
3
4
5
6
7
8
9
10
11
12
13
for(int i=0; p[i]*p[i]<a&&i<cnt; i++){
        int tot=0;
        while(a%p[i] == 0){
            tot++;
            a /= p[i];
        }
        ans *= (sum(p[i], tot*b)%mod);
        ans %= mod;//注意及时的模数
    }
    if(a > 1){
        ans *= (sum(a, b)%mod);
        ans %= mod;//注意及时的模数
    }

AC代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
#include <cstdio>
#include <cstring>
using namespace std;
typedef long long ll;
const int mod = 9901;
int A, B;
const int maxn = 1e4+10;
typedef long long ll;
bool prime[maxn];
int p[maxn];
int cnt;
void get_prime(){
    for(int i=2; i<maxn; i++){
        if(prime[i]){
            p[cnt++] = i;
            for(int j=i+i; j<maxn; j+=i){
                prime[j] = false;
            }
        }
    }
}
ll pow_mod(int a, int n){
    ll ans = 1;
    a = a%mod;
    while(n){
        if(n&1){
            ans = ans*a%mod;
        }
        a = a*a%mod;
        n>>=1;
    }
    return ans;
}
int sum(int a, int n){
    if(n == 0) return 1;
    ll half_sum = sum(a, (n-1)/2);
    ll ans;
    if(n&1){
        ll a1 = pow_mod(a, (n+1)/2);
        ans = (half_sum+half_sum%mod*a1%mod)%mod;
    }
    else{
        ll a1 = pow_mod(a, (n+1)/2);
        ans = (half_sum+half_sum%mod*a1%mod)%mod;
        ans = (ans + pow_mod(a, n))%mod;
    }
    return ans;
}
int solve(int a, int b){
    ll ans = 1;
    int temp = a;
    for(int i=0; p[i]*p[i]<a&&i<cnt; i++){
        int tot=0;
        while(a%p[i] == 0){
            tot++;
            a /= p[i];
        }
        ans *= (sum(p[i], tot*b)%mod);
        ans %= mod;//注意及时的模数
    }
    if(a > 1){
        ans *= (sum(a, b)%mod);
        ans %= mod;//注意及时的模数
    }
    return ans;
}
int main()
{
    for(int i=0; i<maxn; i++) prime[i] = true;
    cnt = 0;
    get_prime();
    while(~scanf("%d%d", &A, &B)){
        printf("%d\n", solve(A, B));
    }
    return 0;
}

逆元AC代码

ll solve(ll a, ll b)

这个一开始没有使用ll类型害死人。

其中快速幂的时候使用了二分乘积。
因为mod·b的范围超过了int,a的范围在[0, mod·b]之间， 乘积的话会溢出long long ,而加法不会，所以要用乘积的快速幂。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
using namespace std;
typedef long long ll;
const int mod = 9901;//9901 is a prime number.
const int maxn = 1e5+10;
typedef long long ll;
ll a, b;
bool primes[maxn];
int cnt = 0;
int p[maxn];
void is_prime(){
    for(int i=2; i<=maxn; i++){
        if(primes[i]){
            p[cnt++] = i;
            for(int j=2*i; j<=maxn; j+=i){
                primes[j] = false;
            }
        }
    }
}
ll mul(ll a, ll b, ll mod){
    ll ans = 0;
    //a%=mod;
    while(b){
        if(b&1){
            ans = (ans+a)%mod;
        }
        b>>=1;
        a = (a+a)%mod;
    }
    return ans;
}
ll pow_mod(ll a, ll n, ll p){
    ll ans = 1;
    //a = a%p;
    while(n){
        if(n&1){
            ans = mul(ans, a, p);
        }
        a = mul(a, a, p);
        n>>=1;
    }
    return ans;
}
ll solve(ll a, ll b){
    ll ans = 1;
    for(int i=0; p[i]*p[i]<=a&&i<cnt; i++){
        int num = 0;
        if(a%p[i] == 0){
            while(a%p[i] == 0){
                num++;
                a/=p[i];
            }
            ll p_mod = (p[i]-1)*mod;
            ans *= ll(pow_mod(p[i], num*b+1, p_mod)+(p_mod-1))/(p[i]-1);
            ans %= mod;
        }
    }
    if(a>1){
        ll p_mod = (a-1)*mod;
        ans *= ll(pow_mod(a, b+1, p_mod)+(p_mod-1))/(a-1);
        ans %= mod;
    }
    return ans;
}
int main()
{
    for(int i=0; i<maxn; i++) primes[i] = true;
    is_prime();
    while(~scanf("%lld%lld", &a, &b)){
        printf("%lld\n", solve(a, b));
    }
    return 0;
}

未解决的问题

  • 等比矩阵的和的构造方法
  • 逆元来实现poj1845.(solved)

本文标题:等比数列的二分求和

文章作者:Babydragon

发布时间:2018-02-25, 15:02:53

最后更新:2018-03-23, 13:33:12

原始链接:http://baolintian.github.io/2018/02/25/等比数列的二分求和/

许可协议: "署名-非商用-相同方式共享 4.0" 转载请保留原文链接及作者。

文章目录
  1. 1. 等比二分求和
    1. 1.1. 二分求解等比数列的和
    2. 1.2. 二分求解等比矩阵的和
  2. 2. 题目
    1. 2.1. 题意
    2. 2.2. 题解
    3. 2.3. AC代码
    4. 2.4. 逆元AC代码
  3. 3. 未解决的问题
{{ live2d() }}