算法竞赛中的数学模板

最新推荐文章于 2025-03-02 12:33:32 发布

原创最新推荐文章于 2025-03-02 12:33:32 发布 · 371 阅读

3 ·

CC 4.0 BY-SA版权

Mathematics 专栏收录该内容

0 篇文章

订阅专栏

博客介绍了多种数论相关内容，包括快速幂取模、最大公约数与最小公倍数、扩展欧几里得、求解线性同余方程、逆元求解等。还提及费马小定理、中国剩余定理等重要定理，以及斐波那契数列、Catalan数等，同时介绍了素数判断、筛法求素数和分解质因数的方法。

快速幂取模（不取模就把取模的位置删除）：

#include <bits/stdc++.h>

using namespace std;

typedef long long ll;
ll mod = 1e9 + 7;
//fast modular exponentiation O(logn)
ll qpow(ll a, ll b)
{
    ll re = 1;
    while(b)
    {
        if(b & 1)
            re = (re * a) % mod;
        b >>= 1;
        a = (a * a) % mod;
    }
    return re % mod;
}

int main()
{
    ll a, b;
    cin >> a >> b;
    cout << qpow(a, b) << endl;
    return 0;
}

最大公约数与最小公倍数：

#include <bits/stdc++.h>
using namespace std;
//greatest common divisor O(logb)
//gcd recursion
int gcd(int a, int b)
{
    return b ? gcd(b, a%b) : a;
}

//gcd iteration
int igcd(int a, int b)
{
    int tmp;
    if(a == 0 || b == 0) return a + b;
    while(b)
        tmp = b, b = a%b, a = tmp;
    return tmp;
}

//gcd optimization
int optgcd(int a, int b)
{
    int i, j;
    if(a == 0 || b == 0) return a + b;
    for(i = 0; (a&1) == 0; ++i) a >>= 1;
    for(j = 0; (b&1) == 0; ++j) b >>= 1;
    if(j < i) i = j;
    while(1)
    {
        if(a < b) a ^= b, b ^= a, a ^= b;
        if((a -= b) == 0) return b << i;
        while((a&1) == 0) a >>= 1;
    }
}

int main()
{
    int a, b;
    while(cin >> a >> b)
    {
        cout << igcd(a, b) << endl;
        cout << gcd(a, b) << endl;
        cout << optgcd(a, b) << endl;
        //lowest common multiple
        cout << a * b / gcd(a, b) << endl;
    }
    return 0;
}

扩展欧几里得（求解贝祖等式）：

#include <bits/stdc++.h>

using namespace std;
//extended Euclidean algorithm O(lnn)
//ax + by = gcd(a, b)
//first code 
int exgcd(int a, int b, int& x, int& y)
{
    int ret, tmp;
    if(!b)
    {
        x = 1, y = 0;
        return a;
    }
    ret = exgcd(b, a%b, x, y);
    tmp = x, x = y, y = tmp - a / b * y;
    return ret;
}

//second code
void egcd(int a, int b, int& gcd, int& x, int& y)
{
    if(!b)
        x = 1, y = 0, gcd = a;
    else
    {
        egcd(b, a%b, gcd, y, x);
        y -= x * (a / b);
    }
}

int main()
{
    int a = 84, b = 36, x, y, gcd;
    gcd = exgcd(a, b, x, y);
    cout << gcd << " " << x << " " << y << endl;
    egcd(a, b, gcd, x, y);
    cout << gcd << " " << x << " " << y << endl;
    return 0;
}

求解线性同余方程：

#include <bits/stdc++.h>

using namespace std;
//judged O(lnn)
//extended greatest common divisor
void exgcd(int a, int b, int& gcd, int& x, int& y)
{
    if(!b)
        x = 1, y = 0, gcd = a;
    else
    {
        exgcd(b, a%b, gcd, y, x);
        y -= a / b * x;
    }
}

//Linear Congruence Theorem
//it is a specific solution of the Linear Congruence Theorem
bool lct(int a, int b, int c, int& gcd, int& x, int& y)
{
    exgcd(a, b, gcd, x, y);
    if(c%gcd)
        return false;
    int t = c / gcd;
    x *= t;
    y *= t;
    return true;
}

int main()
{
    int a, b, c, gcd, x, y;
    //sample input: 17 50 3 sample output: 1 9 -3
    cin >> a >> b >> c;
    bool flag = lct(a, b, c, gcd, x, y);
    if(flag)
        cout << gcd << " " << x << " " << y << endl;
    else
        cout << "no answer" << endl;
    return 0;
}

逆元及求解a/b%mod：

费马小定理(Fermat's little theorem)是数论中的一个重要定理；一个质数p，而整数a不是p的倍数，则有a^（p-1）≡1（mod p），

GCD(a, p) = 1。

扩展欧几里得ax+by=gcd(a,b)，当gcd(a,b)=1的时候，ax+by=1,根据逆元的定义可知x的解就是a关于模b的逆元，y的解就是b关于模a的逆元，当gcd(a,b)=1时成立。

#include <bits/stdc++.h>

using namespace std;
typedef long long ll;
//Fermat's little theorem get inverse element O(log2n)
ll qpow(ll a, ll b, ll mod)
{
    ll re = 1;
    while(b)
    {
        if(b & 1)
            re = (re * a) % mod;
        b >>= 1;
        a = (a * a) % mod;
    }
    return re % mod;
}

ll fermatInv(ll a, ll b)
{
    return qpow(a, b - 2, b);
}

//extended greatest common divisor get inverse element O(lnn)
void exgcd(ll a, ll b, ll& gcd, ll& x, ll& y)
{
    if(!b)
        x = 1, y = 0, gcd = a;
    else
    {
        exgcd(b, a%b, gcd, y, x);
        y -= a / b * x;
    }
}

ll exgcdInv(ll a, ll b)
{
    ll gcd, x, y;
    exgcd(a, b, gcd, x, y);
    return gcd == 1 ? (x + b) % b : -1;
}

int main()
{
    ll a, b, mod;
    //sample input: 1000 53 9973
    //sample output: 8844 7922
    cin >> a >> b >> mod;
    ll exi = exgcdInv(b, mod);
    cout << "exgcdInv: " << exi << endl;
    cout << ((a % mod) * (exi % mod)) % mod << endl;
    ll fi = fermatInv(b, mod);
    cout << "fermatInv: " << fi << endl;
    cout << ((a % mod) * (fi % mod)) % mod << endl;
    return 0;
}

逆元表：给定一个质数p，一个整数n小于p，求1到n的模p逆元表。

#include <bits/stdc++.h>

using namespace std;
//judged O(n)
typedef long long ll;
const int maxn = 1e5 + 10;
ll inv[maxn];

//get inverse table
void invs(ll inv[], ll n, ll p)
{
    inv[1] = 1;
    for(ll i = 2; i <= n; ++i)
        inv[i] = (p - p/i) * inv[p % i] % p;
}

int main()
{
    ll n, p;
    cin >> n >> p;
    invs(inv, n, p);
    for(int i = 1; i <= n; ++i)
        cout << inv[i] << " ";
    cout << endl;
    return 0;
}

中国剩余定理（Chinese Remainder Theorem）：

#include <bits/stdc++.h>

using namespace std;
//judged
const int maxn = 1e5 + 10;
int mod[maxn], a[maxn], n;
void exgcd(int a, int b, int& x, int& y)
{
    if(!b) x = 1, y = 0;
    else
    {
        exgcd(b, a%b, y, x);
        y -= a / b * x;
    }
}

//Chinese Remainder Theorem
int crt(int n)
{
    int ans = 0, m = 1, x, y, tmp;
    for(int i = 0; i < n; ++i)
        m *= mod[i];
    for(int i = 0; i < n; ++i)
    {
        tmp = m / mod[i];
        exgcd(tmp, mod[i], x, y);
        x = (x % mod[i] + mod[i]) % mod[i];
        ans =(ans + x * a[i] * tmp) % m;
    }
    return ans;
}

int main()
{
    while(cin >> n)
    {
        for(int i = 0; i < n; ++i)
            cin >> a[i] >> mod[i];
        cout << crt(n) << endl;
    }
    return 0;
}

扩展中国剩余定理（extended Chinese Remainder Theorem）（解决模数m不能互质的情况）：

#include <bits/stdc++.h>

using namespace std;
//judged 
typedef long long ll;
const ll maxn = 1e6 + 10;

ll m[maxn], a[maxn], n;

ll gcd(ll a, ll b)
{
    return b ? gcd(b, a % b) : a;
}

void exgcd(ll a, ll b, ll& x, ll& y)
{
    if(!b) x = 1, y = 0;
    else
    {
        exgcd(b, a % b, y, x);
        y -= a / b * x;
    }
}

ll inv(ll a, ll b)
{
    ll x, y;
    exgcd(a, b, x, y);
    return (x + b) % b;
}

//extended Chinese Remainder Theorem
ll excrt(int n)
{
    ll ans, A, M, mg1, mg2, mg, g;
    ans = a[0], M = m[0];
    for(int i = 1; i < n; ++i)
    {
        g = gcd(m[i], M);
        mg1 = m[i] / g, mg2 = M / g, mg = m[i] * M / g;
        ans = (ans - a[i]) / g * inv(mg1, mg2) % mg2 * m[i] + a[i];
        M = m[i] * M / g;
    }
    return (ans % M + M) % M;
}

int main()
{
    ll n;
    while(cin >> n && n >= 2)
    {
        for(int i = 0; i < n; ++i)
            cin >> m[i] >> a[i];
/*        ll A, mg1, mg2, mg, g;
        for(int i = 1; i < n; ++i)
        {
            g = gcd(m[i], m[i -1]);
            cout << "gcd m1 m2: " << g << endl;
            a[i] = (a[i - 1] - a[i]) / g * inv(m[i] / g, m[i - 1] / g) % (m[i - 1] / g) * m[i] + a[i];
            m[i] = m[i] * m[i - 1] / g;
        }
        cout << (a[n - 1] + m[n - 1]) % m[n - 1] << endl;
        */
        cout << excrt(n) << endl;
    }
    return 0;
}

斐波那契数列（Fabonacci Sequence）：

#include <iostream>
#include <cstdlib>
#include <cmath>

using namespace std;
typedef long long ll;
const double sqrtFive = sqrt(5.0);

ll fab(int n)
{
   return (pow((sqrtFive + 1.0) / 2.0, n) - pow((1.0 - sqrtFive) / 2.0, n)) * (sqrtFive / 5.0);
}

ll f(int n)
{
    return n == 1 || n == 2 ? 1 : f(n - 1) + f(n - 2);
}

ll ff(int n)
{
    ll a[100]={0, 1, 1};
    if(n <= 2) return a[n];
    for(int i = 2; i <= n; ++i)
        a[i] = a[i - 1] + a[i - 2];
    return a[n];
}

int main()
{
    cout << sqrtFive << endl;
    int n;
    while(cin >> n)
    {
        cout << fab(n) << endl;
        cout << f(n) << endl;
    }
    return 0;
}

Catalan数（括号化、出栈次序、凸多边形三角划分、给定节点组成二叉搜索树、n对括号正确匹配数目）：

#include <iostream>
#include <cstring>

using namespace std;

typedef long long ll;

//f(n){(i=0--n-1){f(i)*f(n-i-1)}}
ll a[100];
ll f(int n)
{
    a[0] = 1;
    for(int i = 1; i <= n; ++i)
    {
        for(int j = 0; j < i; ++j)
            a[i] += a[j] * a[i - j - 1];
    }
    return a[n];
}

//f(n){f(n-1)*(4*n-2)/(n+1)}
ll b[100];
ll ff(int n)
{
    b[0] = 1;
    for(int i = 1; i <= n; ++i)
        b[i] = b[i - 1] * (4 * n - 2) / (n + 1);
    return b[n];
}

//f(n){C(n, 2n)/(n+1)}
ll fff(int n)
{
    ll up = 1, down = 1;
    for(int i = 0; i < n; ++i)
    {
        up *= (2 * n - i);
        down *= (n - i);
        if(up % down == 0) {up /= down; down = 1;}
    }
    return up / (n + 1);
}

//f(n){C(n, 2n)-C(n-1, 2n)}
ll ffff(int n)
{
    ll up = 1, down = 1;
    for(int i = 0; i < n; ++i)
    {
        up *= (2 * n - i);
        down *= ( n - i);
        if(up % down == 0) {up /= down; down = 1;}
    }
    return up - up * n / (n + 1);
}

int main()
{
    int n;
    while(cin >> n)
    {
        memset(a, 0, sizeof(a));
        cout << fff(n) << endl;
        cout << ffff(n) << endl;
    }
    return 0;
}

判断单个素数：

#include <bits/stdc++.h>
using namespace std;
//judged O(sqrt(n))
bool isPrime(int n)
{
    if(n == 2) return true;
    if(n < 2 || n % 2 == 0) return false;
    int len = sqrt(n);
    for(int i = 3; i <= len; i += 2)
        if(n % i == 0)
            return false;
    return true;
}

int main()
{
    for(int i = 0; i < 105; ++i)
        if(isPrime(i))
            cout << i << " ";
    cout << endl;
    return 0;
}

埃拉托色尼（Eratosthenes）筛法求素数：

#include <bits/stdc++.h>
using namespace std;

const int maxn = 1e6 + 10;
int vis[maxn];
//judged O(nloglogn)
void EPrimes(int n)
{
    //合数n一定有不超过根号n的质因子，因此外层循环只需要循环到根号n
    int len = sqrt(n);
    for(int i = 2; i <= len; ++i)
        if(!vis[i])
            for(int j = i; j <= n/i; ++j)
                vis[i * j] = 1;
}

int main()
{
    double dur;
    clock_t st, ed;
    st = clock();
    int n = 1000000;
    EPrimes(n);
    int cnt = 0;
    for(int i = 2; i <= n; ++i)
        if(vis[i] == 0)
            cnt++;
    cout << cnt << endl;
    ed = clock();
    dur = (double)(ed - st);
    cout << "Use Time: " << (dur / CLOCKS_PER_SEC) << endl;
    return 0;
}

线性筛法求素数：

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
using namespace std;

const int maxn = 1e5 + 10;
int q[maxn], num = 0;
bool a[maxn];
//judged O(n)
void linearPrimes(int n)
{
    memset(a, 1, sizeof(a));
    a[0] = false, a[1] = false;
    for(int i = 2; i <= n; ++i)
    {
        if(a[i]) q[++num] = i;
        for(int j = 1; j <= num && i * q[j] <= n; ++j)
        {
            a[i * q[j]] = false;
            if(i % q[j] == 0) break;
        }
    }
}

int main()
{
    double dur;
    clock_t st, ed;
    st = clock();
    int n = 100000;
    linearPrimes(n);
    cout << num << endl;
    ed = clock();
    dur = (double)(ed - st);
    cout << "Use Time: " << (dur / CLOCKS_PER_SEC) << endl;
    return 0;
}

分解质因数（Eratosthenes筛法和试除法）：

#include <bits/stdc++.h>

using namespace std;

const int maxn = 1e5 + 10;
bool prime[maxn];
int factor[(int)sqrt(maxn)], fcnt; 
//judged O(sqrt(n)) Eratosthenes and "try to divide"
void divide(int n)
{
    memset(prime, 1, sizeof(prime)); fcnt = 0;
    int num = n, len = sqrt(n);
    for(int i = 2; i <= n; ++i)
    {
        if(prime[i])
        {
            while(num % i == 0)
            {
                factor[fcnt++] = i;
                num /= i;
                if(num == 1) return ;
            }
            for(int j = i; j <= len; j += 1)
               prime[j * i] = 0;
        }
    }
}
int main()
{
    double dur;
    clock_t st, ed;
    st = clock();
    int n;
    cin >> n;
	divide(n);
    for(int i = 0; i < fcnt; ++i)
        cout << factor[i] << " ";
    cout << endl;
    ed = clock();
    dur = (double)(ed - st);
    cout << "Use Time: " << (dur / CLOCKS_PER_SEC) << endl;
    return 0;
}

分解质因数（线性筛法和试除法）：

#include <bits/stdc++.h>

using namespace std;

const int maxn = 1e7 + 10;
bool prime[maxn];
int factor[(int)sqrt(maxn)], fcnt;
vector<int> vec;

void divide(int n)
{
    memset(prime, 1, sizeof(prime));
    vec.clear();
    fcnt = 0;
    int num = n;
    for(int i = 2; i <= n; ++i)
    {
        if(prime[i])
        {
            factor[fcnt++] = i;
            while(num % i == 0)
            {
                vec.push_back(i);
                num /= i;
                if(num == 1) return ;
            }
        }
        for(int j = 0; j < fcnt; ++j)
        {
            prime[factor[j] * i] = 0;
            if(i % factor[j] == 0) break;
        }
    }
}

int main()
{
    double dur;
    clock_t st, ed;
    st = clock();
    int n;
    cin >> n;
    divide(n);
    for(int i = 0; i < vec.size(); ++i)
        cout << vec[i] << " ";
    cout << endl;
    ed = clock();
    dur = (double)(ed - st);
    cout << "Use Time: " << (dur / CLOCKS_PER_SEC) << endl;
    return 0;
}