标签: 快速幂

[BZOJ3209]花神的数论题 题解

[BZOJ3209]花神的数论题 题解

题目地址:洛谷:【P4317】花神的数论题 – 洛谷、BZOJ:Problem 3209. — 花神的数论题

题目描述

话说花神这天又来讲课了。课后照例有超级难的神题啦…… 我等蒟蒻又遭殃了。
花神的题目是这样的
设 sum(i) 表示 i 的二进制表示中 1 的个数。给出一个正整数 N ,花神要问你Π(sum(i)),也就是 sum(1)—sum(N) 的乘积。

输入输出格式

输入格式:
一个正整数 N。

输出格式:
一个数,答案模 10000007 的值。

输入输出样例

输入样例#1:

3

输出样例#1:

2

说明

对于样例一,1*1*2=2;
数据范围与约定
对于 100% 的数据,N≤10^15

题解

枚举二进制中含有多少个1,我们可以用数位DP处理出这种情况的方案数,然后快速幂算出结果,乘进答案中。

代码

// Code by KSkun, 2018/6
#include <cstdio>
#include <cctype>
#include <cstring>

#include <algorithm>

typedef long long LL;

inline char fgc() {
    static char buf[100000], *p1 = buf, *p2 = buf;
    return p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 100000, stdin), p1 == p2)
        ? EOF : *p1++;
}

inline LL readint() {
    register LL res = 0, neg = 1; register char c = fgc();
    for(; !isdigit(c); c = fgc()) if(c == '-') neg = -1;
    for(; isdigit(c); c = fgc()) res = (res << 1) + (res << 3) + c - '0';
    return res * neg;
}

const int MAXN = 60, MO = 10000007;

inline LL fpow(LL n, LL k) {
    LL res = 1;
    while(k) {
        if(k & 1) res = (res * n) % MO;
        n = n * n % MO;
        k >>= 1;
    }
    return res;
}

LL n, dp[MAXN][MAXN][MAXN][2], res[MAXN];
int m, num[MAXN];

inline LL dfs(int step, int cnt, int tot, bool lim) {
    if(!step) return cnt == tot;
    if(dp[step][cnt][tot][lim] != -1) return dp[step][cnt][tot][lim];
    LL res = 0; int limm = lim ? num[step] : 1;
    for(int i = 0; i <= limm; i++) {
        res += dfs(step - 1, cnt + (i == 1), tot, lim && i == limm);
    }
    return dp[step][cnt][tot][lim] = res;
}

int main() {
    memset(dp, -1, sizeof(dp));
    n = readint();
    while(n) {
        num[++m] = n & 1; n >>= 1;
    }
    for(int i = 1; i <= m; i++) {
        res[i] = dfs(m, 0, i, true);
    }
    LL ans = 1;
    for(int i = 1; i <= m; i++) {
        ans = ans * fpow(i, res[i]) % MO;
    }
    printf("%lld", ans);
    return 0;
}
[SDOI2015]序列统计 题解

[SDOI2015]序列统计 题解

题目地址:洛谷:【P3321】[SDOI2015]序列统计 – 洛谷、BZOJ:Problem 3992. — [SDOI2015]序列统计

题目描述

小C有一个集合S,里面的元素都是小于M的非负整数。他用程序编写了一个数列生成器,可以生成一个长度为N的数列,数列中的每个数都属于集合S。小C用这个生成器生成了许多这样的数列。但是小C有一个问题需要你的帮助:给定整数x,求所有可以生成出的,且满足数列中所有数的乘积mod M的值等于x的不同的数列的有多少个。小C认为,两个数列{Ai}和{Bi}不同,当且仅当至少存在一个整数i,满足Ai≠Bi。另外,小C认为这个问题的答案可能很大,因此他只需要你帮助他求出答案mod 1004535809的值就可以了。

输入输出格式

输入格式:
一行,四个整数,N、M、x、|S|,其中|S|为集合S中元素个数。第二行,|S|个整数,表示集合S中的所有元素。

输出格式:
一行,一个整数,表示你求出的种类数mod 1004535809的值。

输入输出样例

输入样例#1:

4 3 1 2
1 2

输出样例#1:

8

题解

参考资料:bzoj 3992: [SDOI2015]序列统计 (NTT+快速幂+DP) – CSDN博客
本题需要的数学姿势:数学笔记:数论(欧拉函数、阶、原根) | KSkun’s Blog
首先,计数自然是一个动态规划,我们考虑设计状态dp[i][j]为选了i个数且结果mod M的值是j,转移如下
dp[i][j] = \sum_{k \in S} dp[i-1][j'] \ (j'k \bmod M = j)
直观上感觉直接做是一个O(n|S|)的,好像只能得10分。
接下来我们要考虑一个问题,模数很特殊,并不是常见的1000000007,而是1004535809,它的原根为3,是常用的NTT模数,难道出题人在指引我们使用NTT?
想用NTT做,首先得搞出卷积,卷积下标是加减的事情没办法像上面那样乘吧。那么考虑利用原根把乘法转换为加法,dp[i][j]表示选到第i个数结果mod M的值为g^j,改变以后的转移方程就成这样了
dp[i][j] = \sum_{g^k \bmod M \in S} dp[i-1][j'] \ ((j' + k) \bmod (M - 1) = j)
或者这样?用f[i]表示g^i是否在集合中。
dp[i][j] = \sum_{k=1}^{M-1} dp[i-1][j-k] \cdot f[k]
咦,后面那部分现在可以卷积了,而且你会发现每一层的转移中f数组都是不变的,也就是说这是个线性变换,参照矩阵快速幂的操作,我们可以做多项式的快速幂,即dp[n] = dp[0] \times f^n。答案便是 dp[n][i] \ (g^i=x)
于是我们的思路就很明确了:求原根、算出f数组、多项式快速幂。
没完呢!这里有一个坑点,集合内的元素可以有0,0显然对答案没贡献,所以你要跳过它处理。

代码

// Code by KSkun, 2018/4
#include <cstdio>
#include <cmath>
#include <cstring>

#include <algorithm>

typedef long long LL;

inline char fgc() {
    static char buf[100000], *p1 = buf, *p2 = buf;
    return p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 100000, stdin), p1 == p2) ? EOF 
        : *p1++;
}

inline LL readint() {
    register LL res = 0, neg = 1;
    char c = fgc();
    while(c < '0' || c > '9') {
        if(c == '-') neg = -1;
        c = fgc();
    }
    while(c >= '0' && c <= '9') {
        res = res * 10 + c - '0';
        c = fgc();
    }
    return res * neg;
}

inline LL fpow(LL n, LL k, LL p) {
    n %= p;
    LL res = 1;
    while(k) {
        if(k & 1) res = res * n % p;
        n = n * n % p;
        k >>= 1;
    }
    return res;
}


const int MAXN = 1 << 15, G = 3, MO = 1004535809;

int n, m, x, s, t;
LL f[MAXN], g[MAXN];
bool vis[MAXN];

LL rev[MAXN], dft1[MAXN], dft2[MAXN], dftr[MAXN];

inline void calrev(int n, int len) {
    for(int i = 0; i < n; i++) {
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (len - 1));
    }
}

inline void dft(int n, LL *arr, int f) {
    for(int i = 0; i < n; i++) {
        if(i < rev[i]) std::swap(arr[i], arr[rev[i]]);
    }
    for(int i = 1; i < n; i <<= 1) {
        LL gn = fpow(G, (MO - 1) / (i << 1), MO);
        if(f == -1) gn = fpow(gn, MO - 2, MO);
        for(int j = 0; j < n; j += i << 1) {
            LL g = 1;
            for(int k = 0; k < i; k++) {
                LL x = arr[j + k], y = g * arr[j + k + i] % MO;
                arr[j + k] = (x + y) % MO;
                arr[j + k + i] = ((x - y) % MO + MO) % MO;
                g = g * gn % MO;
            }
        }
    }
    if(f == -1) {
        LL invn = fpow(n, MO - 2, MO);
        for(int i = 0; i < n; i++) {
            arr[i] = arr[i] * invn % MO;
        }
    }
}

inline void ntt(LL *a, LL *b) {
    memset(dft1, 0, sizeof(dft1));
    memset(dft2, 0, sizeof(dft2));
    memset(dftr, 0, sizeof(dftr));
    int n, len = 0;
    for(n = 1; n <= (m - 1) << 1; n <<= 1) len++;
    for(int i = 0; i < n; i++) {
        dft1[i] = a[i];
        dft2[i] = b[i];
    }
    calrev(n, len);
    dft(n, dft1, 1); dft(n, dft2, 1);
    for(int i = 0; i < n; i++) {
        dftr[i] = dft1[i] * dft2[i] % MO;
    }
    dft(n, dftr, -1);
    for(int i = 0; i < m - 1; i++) {
        a[i] = (dftr[i] + dftr[i + m - 1]) % MO;
    }
}

inline LL calg(LL n) {
    int fact[8005], tot = 0, x = n - 1;
    for(int i = 2; i * i <= x; i++) {
        if(x % i == 0) {
            fact[tot++] = (n - 1) / i;
            while(x % i == 0) x /= i;
        }
    }
    if(x > 1) fact[tot++] = (n - 1) / x;
    bool success;
    for(int i = 2;; i++) {
        success = true;
        for(int j = 0; j < tot; j++) {
            if(fpow(i, fact[j], n) == 1) {
                success = false;
                break;
            }
        }
        if(success) return i;
    }
}

int main() {
    n = readint(); m = readint(); x = readint(); s = readint();
    for(int i = 1; i <= s; i++) {
        t = readint();
        vis[t] = 1;
    }
    int xpow, G = calg(m), gn = 1;
    for(int i = 0; i < m - 1; i++) {
        if(vis[gn]) f[i] = 1;
        if(gn == x) xpow = i;
        gn = gn * G % m;
    }
    g[0] = 1;
    for(int i = n; i; i >>= 1) {
        if(i & 1) ntt(g, f);
        ntt(f, f);
    }
    printf("%lld", g[xpow]);
    return 0;
}
模运算性质、快速幂(取模)

模运算性质、快速幂(取模)

模运算的一些性质

同余公式

形如a≡b (mod d)的式子被称为同余公式,因为此式中a与b模d后的值相等,故被称为同余公式。相关性质将基于它展开。

它等价于(a-b)|d、a-(a/d)*d=b、a=b+nd (n∈Z)。

运算性质

  1. a≡a(mod d)
  2. 对称性 a≡b (mod d)→b≡a (mod d)
  3. 传递性 (a≡b (mod d),b≡c (mod d))→a≡c (mod d)
  4. 如果a≡x (mod d),b≡m (mod d),则
    a+b≡x+m (mod d)
    a-b≡x-m (mod d)
    a*b≡x*m (mod d )
    a/b≡x/m (mod d)
  5. a≡b (mod d)则a-b|d
  6. a≡b (mod d)则an≡bn (mod d)
  7. 如果ac≡bc (mod m),且c和m互质,则a≡b (mod m)

运算规则

(a + b)  mod  p = (a  mod  p + b  mod  p)  mod  p
(a – b)  mod  p = (a  mod  p – b  mod  p)  mod  p
(a * b)  mod  p = (a  mod  p * b  mod  p)  mod  p
ab  mod  p = ((a mod p)b)  mod  p
结合率: ((a + b)  mod  p + c)  mod  p = (a + (b + c)  mod  p)  mod  p
((a * b)  mod  p * c) mod  p = (a * (b * c)  mod  p)  mod  p
交换率: (a + b)  mod  p = (b+a)  mod  p
(a * b)  mod  p = (b * a)  mod  p
分配率: ((a + b) mod  p * c)  mod  p = ((a * c)  mod  p + (b * c)  mod  p)  mod  p
重要定理:若a≡b (mod  p),则对于任意的c,都有(a + c) ≡ (b + c) (mod p)
若a≡b (mod  p),则对于任意的c,都有(a * c) ≡ (b * c) (mod p)
若a≡b (mod  p),则对于任意的c,都有ac≡ bc (mod p)

快速幂(取模)

简介

快速幂是运用分治思想降低朴素算法复杂度的算法,是OI中大数据常用的优化小trick。

说明

假设我们现在正在解决求a^b的问题,其中b非常大,而且通常会遇到朴素算法运算中,中间得数就已经会爆int甚至爆long long的情况。

我们知道a^b=(a^(b/2))^2,因此我们可以将求a^b问题转化为求a^(b/2)问题,如此递归下去,总会遇到b/2=0的时候,此时递归就到头了。这样做的时间复杂度是O(logn)的(而朴素是O(n)的)。

我们还可以运用取模运算的性质(a * b)  mod  p = (a  mod  p * b  mod  p)  mod  p来对每一步的值取模缩小值的范围。此算法便是边幂边取模的算法。

如果b是一个单数怎么办?也不要紧:a^b=(a^(b/2))^2*a即可解决。

代码

下面是递归式快速幂。

int fpow(int x, int p) { // 求x^p 
    if(p == 0) return 1;
    int t = fpow(x, p / 2);
    if(p % 2 == 0) {
        return t * t;
    } else {
        return t * t * x;
    }
}

int fpow_m(int x, int p, int m) { // 求(x^p)%m
    if(p == 0) return 1;
    int t = fpow_m(x, p / 2, m) % m;
    if(p % 2 == 0) {
        return (t * t) % m;
    } else {
        return ((t * t) % m * (x % m)) % m; 
    }
}

下面是非递归式快速幂。

inline LL fpow(LL x, LL k) {
    LL t = 1;
    while(k) {
        if(k & 1) t = (t * x) % MO;
        x = (x * x) % MO;
        k >>= 1;
    }
    return t;
}

应用

应用不用说了,常用于数据大的NOIP提高组复赛T1之类的题目。

参考资料

模运算与同余公式的性质 – 根号下的麻辣烫 – CSDN博客

例题:2011 北大OpenJudge百练4010

描述

已知长度最大为200位的正整数n,请求出2011^n的后四位。

输入

第一行为一个正整数k,代表有k组数据,k<=200接下来的k行,

每行都有一个正整数n,n的位数<=200

输出

每一个n的结果为一个整数占一行,若不足4位,去除高位多余的0

样例输入

3
5
28
792

样例输出

1051
81
5521

分析

快速幂是肯定得用的,但是同时也考验了高精度幂运算的知识。

此处我用了一个小trick,打表知取2011的k次幂时,只要整数k后三位相等,他们运算后所得数后四位相等。故省去了写高精度的功夫。

代码

#include <cstdio>
#include <cstring>
#include <cstdlib> 

int fpow_m(int x, int p, int m) { // 求(x^p)%m
    if(p == 0) return 1;
    int t = fpow_m(x, p / 2, m) % m;
    if(p % 2 == 0) {
        return (t * t) % m;
    } else {
        return ((t * t) % m * (x % m)) % m; 
    }
}

int n, t;
char s[205];
char tmp[5];

int main() {
    scanf("%d", &n);
    for(int i = 0; i < n; i++) {
        memset(s, 0, sizeof s);
        memset(tmp, 0, sizeof tmp);
        scanf("%s", s);
        int len = strlen(s), cnt = 0;
        for(int i = 0; i < 3; i++) {
            if(len - 3 + i < 0) continue;
            tmp[cnt] = s[len - 3 + i];
            cnt++;
        }
        //printf("p: %s\n", tmp);
        printf("%d\n", fpow_m(2011, atoi(tmp), 10000));
    } 
    return 0;
}

参考资料

快速幂+快速幂经典例题 -  - CSDN博客