ACM模板 数论

ACM模板 数论

一、素数筛法

1.Eratosthenes筛法（埃氏筛法）

生成1~n的素数表，生成了tot个
复杂度$O(n\log log\left( n\right))$

void getPrimes(int n, int vis[], int prime[], int &tot)
{
    int m = sqrt(n+0.5);
    memset(vis, 0, sizeof(int)*(n+1));
    for(int i = 2; i <= m; i++) if(!vis[i])
        for(int j = i*i; j <= n; j+=i) vis[j] = true;
    tot = 0;
    for(int i=2;i<=n;i++) if(!vis[i]) prime[tot++]= i;
}

2.线性筛法

const int maxp = 1e6;
int prime[maxp] = {0},cnt = 0;
bool notPri[maxp] = {1, 1};
void linearGetPrimes()
{
    for(int i = 2; i < maxp; i++)
    {
        if(!notPri[i]) prime[cnt++] = i;
        for(int j = 0; j < cnt && (long long)i*prime[j] < maxp; j++)
        {
            notPri[i*prime[j]] = 1;
            if(!(i%prime[j])) break;
        }
    }
}

二、分解质因数

1.$\sqrt{n}扫质因数$

n是待分解的数，tot是质因数种类数，d是质因数，fre是每种质因数的个数。

void getFactors(int n,int fre[],int d[],int& tot)
{
    tot = 0;
    int m = sqrt(n+0.5);
    for(int i=2;i<=m;i++)
    {
        if(n%i==0)
        {
            d[tot] = i;
            fre[tot]=0;
            while(n%i==0) n/=i,fre[tot]++;
            tot++;
        }
    }
    if(n>1)
        d[tot] = n,
        fre[tot++] = 1;
}

2.素筛扫质因数

三、欧几里得、扩展欧几里得

1.GCD

int gcd(int a, int b)
{
    if(!b) return a;
    return gcd(b,a%b);
}

2.扩展欧几里得

• 求a和b的最大公约数,
• 求整数x和y，使得$ax+by=gcd(a,b)$

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

得到x和y后，可继续得到$ax+by=c$的整数解，下面计算了x的最小正值。
若无解，返回false

int d;
bool solve(int a, int b,int c, int& x,int& y)
{
    exgcd(a, b, d, x, y);
    if(!d || c%d) return false;
    x *= c/d;
    y *= c/d;
    a /= d;
    b /= d;
    if(a<0) a=-a;
    if(b<0) b=-b;
    ((x %= b) += b) %= b;//此时x+mb，y-ma都是解，m为自然数  
    return true;
}

四、常用运算取模大全

1.加减乘取模

$(a+b)\bmod n = (a\bmod n + b\bmod n)\bmod n$
$(a-b)\bmod n = (a\bmod n - b\%n +n)\bmod n$
$(a \times b)\bmod n = ((a\bmod n) \times (b\bmod n))\bmod n$
注意溢出，如当n为int时，下面代码可用：

    int mul_m(int a, int b, int n)
    {
        a %= n; b %= n;
        return (int)((long long)a * b % n);
    }

若扔可能溢出，用快速乘取模

2.大整数取模

大整数存在字符数组n里

    char n[10000];
    int m;
    scanf("%s%d", n, &m);
    int len = strlen(n);
    int ans = 0 ;
    for(int i = 0; i < len; i++)
        ans = (int)((long long)ans*10 + n[i] - 48) %m;
    printf("%d\n", ans);

3.快速幂取模

求$a^n \bmod{}m$

long long pow_m(LL a, LL n, int m)
{
    a %= m;
    if(!n) return 1;
    LL x = pow_m(a, n/2, m);
    LL ans = x * x %m;
    if(n&1) (ans *= a) %= m;
    return ans;
}

4.快速乘取模

LL mul_m(LL a, LL b , int mod)
{
    LL ans = 0;
    a %= mod;
    b %= mod;
    while (b)
    {
        if (b&1) ans = (ans + a ) % mod;
        b >>= 1;
        (a <<= 1) %= mod;
    }
    return ans;
}

5.等比数列求和取模

求$1+p+p^2...+p^n\pmod{mod}$

LL sumg_m(LL p,LL n,LL mod)
{
    if(p==0)return 0;
    if(n==0)return 1;
    if(n&1) return mul_m(
                         pow_m(p,n/2+1,mod)+1,
                         sumg_m(p,n/2,mod),
                         mod);
    else return (
                    mul_m(
                          pow_m(p,n/2+1,mod)+1,
                          sumg_m(p,n/2-1,mod),
                          mod)
                    +pow_m(p,n/2,mod)
                )%mod;
}

五、求模逆元

需要用到
- 扩展欧几里得算法

$ax\equiv 1\pmod{m}$，求$x$

int a,b,c,x,y,d;
void modinv(int a, int m, int& x)
{
    solve(a,m,1,x,y);
}

六、因子和

1.求某个数字的因子和

因子和即该数字所有因子的和。

int n,sum[maxn];
void sigma(int n,int sum[])
{
    memset(sum,0,sizeof(int)*(n+1));
    for(int i=1;i<=n;i++)
        for(int j=i;j<=n;j+=i) sum[j]+=i;
    //for(int i=1;i<=n;i++) sum[i]-=i;//因数不含自身时用
}

2.求幂的因子和的模

求$\sigma(a^b)\pmod{mod}$
需要用到

• 分解质因数
• 快速乘取模
• 等比数列求和取模

要注意不能a%=mod

LL sumFactors(LL a,LL b,LL mod)
{
    getFactors(a,fre,d,tot);
    LL ans = 1;
    for(int i=0;i<tot;i++)
    {
        LL p = d[i];
        LL n = b*fre[i];
        ans = mul_m(ans, sumg_m(p,n,mod), mod);
    }
    if(a==0) ans = 0;
    return ans;
}

七、同余方程组 中国剩余定理

1.同余方程组（中国剩余定理）不互质版

已知同余方程组$x\equiv r_i \pmod{m_i}$(其中$m_1,m_2...m_n$不保证互质)的$m_i$和$r_i$，求x.
先看前两个方程，$x \bmod m_1 = r_1$，即

$$x = r_1 + m_1k_1 ......(1)$$
$x \bmod m_2 = r_2$，即 $$x = r_2 + m_2k_2......(2)$$
由上面两式可得， $r_1 + m_1k_1 = r_2 + m_2k_2$，
移项，得 $$m_1k_1 - m_2k_2 = r_2 - r_1......(3)$$
令 $a=m_1,b=-m_2,c=r_2-r_1$，
得 $$ak_1+bk_2=c......(4)$$
(此时要判断是否有解)
用扩展欧几里得解出 $k_1$后带入(1)，得到满足前两个同余方程的一个特解 $x_2$.
在(4)中, $k_1+p·b/gcd(m_1,m_2)$(p∈Z)也是解，则 $k_1+p·m_2/gcd(m_1,m_2)$也是解，
将 $k_1+p·m_2/gcd(m_1,m_2)$带入(1)中，得
$$x = r_1 + m_1k_1+p·m_1m_2/gcd(m_1,m_2)$$ ，
即
$$x = x_0+p·lcm(m_1,m_2)......(5)$$
也就是
$$x \equiv x_0 \pmod{lcm(m_1,m_2)}·......(6)$$
至此，我们求出了满足方程(1)和(2)的解，它就是(6)的解
将(6)和 $x \bmod m_3 = r_3$联立，会发现这和解(1)(2)是同一个问题！于是就可以愉快的往下递推了！！！
下面代码求了x的 最小正值

int excrt(int m[], int r[], int n)
{
     int M = m[0], R = r[0], x, y;
     for (int i = 1; i < n; i++)
     {
         if(!solve(M, m[i], r[i] - R, x, y)) return -1;//无解输出-1
         R += x * M;
         M /= d;
         M *= m[i];
         (R += M) %= M;
     }
     return R;
 }

2.同余方程组（中国剩余定理）互质版

$x\equiv r_i\pmod{m_i}$，其中$m_1,m_2...m_n$互质。
需要用到扩展欧几里得算法
下面返回了最小的正数x
设M为$m_i$的乘积，则x+kM也是解

int crt(int n,int r[],int m[])
{
    int lcm = 1;
    for(int i = 0; i < n; i++) lcm *= m[i];
    int ans = 0;
    int M,t,tmp;
    for(int i = 0; i < n; i++)
    {
        M = lcm/m[i];
        if(!solve(M, m[i], 1, t, tmp)) return -1;
        ans+= r[i]*t*M;
    }
    (ans += lcm) %= lcm;
    return ans;
}

八、欧拉函数

对正整数n，欧拉函数是小于等于n的数中与n互质的数的数目。
公式$\varphi(x) = x(1- \frac{1}{p_1}) (1- \frac{1}{p_2}) \cdots (1- \frac{1}{p_{n-1}})(1- \frac{1}{p_n})$，其中$p_1,p_2, \cdots p_n$为$x$的所有质因数。注意$\varphi(1) = 1$.

1. 求欧拉函数

此为实现2的基础。

int euler(int n)
{
     int res = n, a = n;
     for(ll i = 2; i*i <= a; i++)
     {
         if(0 == a%i)
         {
             res = res/i*(i-1);
             while(a%i==0) a /= i;
         }
     }
     if(a > 1) res=res/a*(a-1);
     return res;
}

2. 求欧拉函数表

const int maxn = 1000001;  
int euler[maxn];  
void Init(){   
     euler[1]=1;  
     for(int i=2;i<maxn;i++) euler[i]=i;  
     for(int i=2;i<maxn;i++)  
        if(euler[i]==i) 
           for(int j=i;j<maxn;j+=i)  
              euler[j]=euler[j]/i*(i-1);//先进行除法是为了防止中间数据的溢出   
}