数论

jiazhaopeng

2019-12-08 19:14:17

Personal

本文主要是数论的基础部分,高级部分已给出链接:

下面进入正题:

初等数论

部分摘自秦岳学长的课件《初等数论》。

π(pi^qi)约数的和=(1+p1+p1^2...)*(...
//设S=1+p+p^2+...+p^q, pS-S=p^(q+1)-1, S=(p^(q+1)-1)/(p-1)
①gcd(a,b)=gcd(b,a-b)
②gcd(2a,2b)=2gcd(a,b)
③gcd(2a,b)=gcd(a,b) (b是奇数)
用途:高精GCD
n个整数间的裴蜀定理

设a1,a2,a3......an为n个整数,d是它们的最大公约数,
那么存在整数x1......xn使得x1*a1+x2*a2+...xn*an=d。

特别来说,如果a1...an互质(不是两两互质),
那么存在整数x1......xn使得x1*a1+x2*a2+...xn*an=1。
//摘自:百度百科:https://baike.baidu.com/item/%E8%A3%B4%E8%9C%80%E5%AE%9A%E7%90%86/5186593?fr=aladdin

例题:P4571 [JSOI2009]瓶子和燃料

这里给出求n个数选k个数,求最大化的gcd的方法:

把n个数的因数都存到一个数组里面,然后排序,从大往小扫,如果某个数出现了至少k次,则说明这个因数在至少k个数中出现过,这个因数即为答案。

Code:

for (register int i = 1; i <= n; ++i) {
    read(num);
    for (register int j = 1; j * j <= num; ++j) {
        if (num % j == 0) {
            yin[++tot] = j;
            if (j * j != num) yin[++tot] = num / j;
        }
    }
}
sort(yin + 1, yin + 1 + tot);
int ct = 1;
for (register int i = tot - 1; i > 0; --i) {
    if (yin[i] == yin[i + 1])   ct++;
    else    ct = 1;
    if (ct == k) {
        printf("%d\n", yin[i]);
        return 0;
    }
}
O(n)计算1~n的逆元
p%i+[p/i]*i = p 
p%i+[p/i]*i ≡ 0(mod p) 
[p/i]*i ≡ -p%i (mod p) 
[p/i] ≡ -p%i*i^(-1) (mod p)
[p/i] * (p%i)^(-1) ≡ -i^(-1) (mod p) 
i^-1 ≡ -[p/i]*(p%i)^(-1) (mod p) 
在数论中的积性函数:
对于正整数n的一个算术函数 f(n),若f(1)=1,且当a,b互质时f(ab)=f(a)f(b),在数论上就称它为积性函数。
若对于某积性函数 f(n) ,就算a, b不互质,也有f(ab)=f(a)f(b),则称它为完全积性的。
φ(n) -欧拉函数,计算与n互质的正整数之数目
μ(n) -莫比乌斯函数,关于非平方数的质因子数目
gcd(n,k) -最大公因子,当k固定的情况
d(n) -n的正因子数目
σ(n) -n的所有正因子之和
σk(n) - 因子函数,n的所有正因子的k次幂之和,当中k可为任何复数。
1(n) -不变的函数,定义为 1(n) = 1 (完全积性)
Id(n) -单位函数,定义为 Id(n) = n(完全积性)
Idk(n) -幂函数,对于任何复数、实数k,定义为Idk(n) = n^k (完全积性)
ε(n) -定义为:若n = 1,ε(n)=1;若 n > 1,ε(n)=0。别称为“对于狄利克雷卷积的乘法单位”(完全积性)
λ(n) -刘维尔函数,关于能整除n的质因子的数目
γ(n),定义为γ(n)=(-1)^ω(n),在此加性函数ω(n)是不同能整除n的质数的数目
另外,所有狄利克雷特征均是完全积性的

较小数的 \phi

φ(n)=πφ(pi^qi)
φ(x=p^q)=p^(q-1)*(p-1)=x*(1-1/p)(p为质数)

阶乘的 \phi

\phi(n) = n * \Pi \frac{p_i - 1}{p_i}

注意:与质因数的指数 c_i 无关!

方法:预处理所有质数,然后对于每一个数的阶乘都包含小于该数的所有质因子,不包含大于该数的所有质因子。根据此预处理处上式的右边部分即可。

my code

x^k mod m = x^kmod m, k < \phi(m) x^k mod m = x^{k mod \phi(m) + \phi(m)}modm, k >= \phi(m)

CRT

(插播)博客:中国剩余定理 CRT

x = a_1(mod m_1) x = a_2(mod m_2) ... x = a_n(mod m_n)

设M为所有模数之积,Mi为除mi以外的所有模数之积,ki为Mi在mi意义下的逆元(即除mi以外的所有模数之积在mi意义下的逆元)。有:

answer = ∑(a_ik_iM_i)mod M

或者这么说:设 P 为模数之积(总模数); p_i 为模数(分模数),则:

Ans = \sum{a_i * (\frac{P}{p_i}) * inv(\frac{P}{p_i}(mod~p_i))} ~ ~ ~ ~ ~ ~(mod~ ~P)

可以这么想:

因为我们要求 a_i 在其它的所有方程中都为0,因此要乘一个 \frac{P}{p_i},只要乘了它,就能搞定其它的所有方程。但是,在第 i 个方程这里,就多乘了一个 \frac{P}{p_i},因此我们要i 这里把它除掉。

x≡A (mod N)  且  x≡B (mod M)

x*=(A*M*ni(M,N)+B*N*ni(N,M))
ni(a,b)表示a在模b意义下的乘法逆元
考虑全部解x=x*+δ,必有δ≡0 (mod N or M)
通解:x=x*+k×lcm(N,M)    (k=...-2,-1,0,1,2...)

EXCRT

直接看题解,讲得很清楚。题解 P4777 【【模板】扩展中国剩余定理(EXCRT)】

//x = b[i](mod p[i])
long long gcd(long long x, long long y)
inline long long lcm(long long x, long long y)
void exgcd(long long aa, long long bb, long long &x, long long &y)
int main() {
    for (register int i = 2; i <= n; ++i) {
        g = exgcd(P, p[i], k1, k2);
        k1 = k1 * (((bb - b[i]) % p[i] + p[i]) % p[i]) / g % p[i]);
        bb -= P * k1;
        P = lcm(P, p[i]);
        bb = (bb % P + P) % P;
    }
   //the answer is bb
}

屠龙勇士

EXCRT 的板子题。需要龟速乘。

实际上EXCRT 真的不难,自己稍微推一推就好了。

关键代码:

namespace jzp2 {
    inline ll slowmul(ll x, ll k, ll P)
    inline ll get_gcd(ll a, ll b)
    inline ll exgcd(ll a, ll b, ll &x, ll &y)
    inline PLL calc(ll a, ll b, ll c) {
        ll g = get_gcd(a, b);
        if (c % g)  return MP(-RND, -RND);
        a /= g; b /= g; c /= g;
        ll x = -RND, y = -RND;
        exgcd(a, b, x, y);
        return MP(slowmul(c, x, b), b);
    }
    inline void sol() {
        ll nwx = 0, nwp = 1;
        for (int i = 1; i <= n; ++i) {
            multiset<ll>::iterator it = st.upper_bound(a[i]);
            ll t;
            if (it == st.begin())   t = *st.begin();
            else    --it, t = *it;
            st.erase(it);
            PLL pr = calc(t, p[i], a[i]);
            if (pr.first == -RND && pr.second == -RND) { puts("-1"); return ; }
            ll newx = pr.first, newp = pr.second;
            pr = calc(nwp, newp, newx - nwx);
            if (pr.first == -RND && pr.second == -RND) { puts("-1"); return ; }
            ll tmp = nwp;
            nwp = nwp / get_gcd(nwp, newp) * newp;
            nwx = ((slowmul(pr.first % nwp, tmp % nwp, nwp) + nwx % nwp) % nwp + nwp) % nwp;
            st.insert(s[i]);
        }
        printf("%lld\n", nwx);
    }
}

阶和原根

阶和原根

阶:

给定 a 和模数 m,满足:

a^t =1(modm)

的最小的正整数 tam 的阶,记作ord_m(a)

ord_m(a) 一定为 \phi(m) 的约数。

原根

给定模数 m,满足:

a^t(modm)

在t∈\left[1, \phi(m) \right ]两两不同的 am 的原根。

  1. 若a为m的原根,则 ord_m(a) = \phi(m)(其实这个才是定义)

  2. 不是所有数都有原根,只有 2, 4, P^a, 2P^a 有原根(P 为奇素数),并且原根个数为 \phi(\phi(m))

  3. 如果数 m 有原根,则其最小原根的规模为 O(m^{\frac{1}{4}})

求一个原根:

先将\phi(p)质因数分解,\phi(p) = P_1^{k_1} * P_2^{k_2} * ... * P_i^{k_i},然后从 2 开始从小到大枚举 a如果对于 \phi(p) 的所有质因数 P_i,都有 a^{\frac{m}{P_i}} \not= 1 \bmod m,那么 am 的原根。

(证明:如果不是原根,那么其阶一定不为 \phi(m),这是定义;并且阶是 \phi(m) 的因数,这是性质。那么肯定存在一个 p_i | \phi(m),这个因数被“漏掉了”,也就是说,不要这个因数也可以同余于一,因此 a^{\frac{\phi(m)}{p_i}} \equiv 1 \pmod m

复杂度:小于 O(m^{\frac{1}{4}}log^2m)

求所有原根:

先求出一个原根 g,则所有原根为 {g^t}(gcd(t, \phi(m)) = 1),因此原根个数为 \phi(\phi(m)) 个。

原根有类似对数(ln)的性质,可以看作模意义下的对数。两个数在模意义下相乘可以看作其对应的原根幂的质数之和。因此可以把乘除变成加减。

例题:P3321 [SDOI2015]序列统计

BSGS

问题:已知 y, z, P,保证 gcd(y, P) = 1。求使得 y^x = z(mod P)成立的最小正整数 x。(求模意义下的对数)

解:

如果把 x 看作 am-bm 自己选定),其中b<=m,a<=\frac{P}{m},那么问题就转化为 y^{am} = zy^b(modP)成立的 a,m,b。我们把等式右边的所有可能取值算出来,存到哈希表(或map)里面。然后枚举左边的 a。如果算出的左边答案在哈希表里面有,就可以得出答案,否则继续枚举 a

复杂度:O(max(m, \frac{P}{m}))。当 m\sqrt{P} 时最优,为 O(\sqrt{P})。用map还会多一个log。

Code:
//求最小非负整数x
namespace Bsgs{
    map<ll, ll> mp;
    inline void sol() {
        mp.clear();
        y = ((y % P) + P) % P; z = ((z % P) + P) % P;
        if (y == 0) {
            if (z == 0) return puts("1"), void();
            failed(); return ;
        }
        if (y == 1) {
            if (z == 1) return puts("0"), void();
            failed(); return ;
        }
        if (z == 1) {
            puts("0"); return ;
        }

        ll m = sqrt(P);
        for (register ll i = 0, t = z; i < m; ++i, t = t * y % P) 
            mp[t] = i;
        for (register ll i = 1, T = quickpow(y, m), t = T; i <= P / m + 10; ++i, t = t * T % P)
            if (mp.count(t)) 
                return printf("%lld\n", m * i - mp[t]), void();
        failed();
    }
}

注意!!

练习题:

P4884 多少个1?

P2485 [SDOI2011]计算器

exBSGS

还是那个问题:已知 y, z, P,求使得 y^x = z(mod P)成立的最小正整数 x

然而这回不保证 gcd(y, P) = 1 了。

我们可以感性地进行特判:设gcd(y, P) = g,如果 g | z 成立,且z还不为 1,那么该方程就无解。(毕竟y^x在模P意义下都为g的倍数)。

如果g | z,那么 g 同时是y, z, P的约数,我们可以全除以 g,就成:

\frac{y}{d}y^{x-1} = \frac{z}{d}(mod\frac{P}{d})

然后我们把 \dfrac{y}{d} 移到右边:

y^{x-1} = inv(\frac{y}{d})\frac{z}{d}(mod\frac{P}{d})

就可以递归解决了,一直到 yP 互质。

然而有可能 y^{x-1} 中的 x-1 为负数了,那样我们就不好办了。(虽然 y^{x-1} 会出现 1 的情况,但是 y 并不与 P 互质。)不过由于递归层数不超过 \log P,我们可以在一开始就把所有的 x=0,1,...,\log P 特判一下。

Code:
int exbsgs(int y, int z, int p) {//y^x = z(mod p)
    if (z == 1) return 0;
    int g = gcd(y, p);
    if (g == 1) return bsgs(y, z, p);
    if (z % g)  return failed(), -1;
    int invv = get_inv(y / g, p / g);
    int res = exbsgs(y, z / g * invv % p, p / g);
    if (res == -1)  return -1;
    return res + 1;
}

...

y %= P; z %= P;
if (y == 0) {
    if (z == 0) {
        puts("1");
        continue;
    }
    failed();
}
for (i = 0 -> 100)  check(i);
int res = exbsgs(y, z, P);
if (res != -1) {
    printf("%lld\n", res);
}

注意:

在保证正确性的前提下,要尽可能地加特判(主要是对0和1的特判)。

除法分块(数论分块)

例题:1257: [CQOI2007]余数之和

//暴力算法
for (register int i = 1; i <= n; ++i)   ans += k % i;

我们发现k % i = k - (k / i) * i;其中k的和容易求,为n × k,重点在于求sigma(k / i * i)

对于k/i相同的情况,实际上就是等差数列求和。我们有发现k/i的值只可能有sqrt(k)个。证明如下:

i <= sqrt(k) :k/i 最多有sqrt(k)个。因为i最多有sqrt(k)个。
i > sqrt(k) : k/i 最多有sqrt(k)个。因为k/i的结果最多有sqrt(k)个。

因此,我们对每个答案等差数列求一遍和。找每个答案的最开始的点和最后的点代码如下:

for (int i = 1; i <= n; ++i) {
            //i为最小的点
    int j = n / i;//j为答案
   int k = n / j;//k为最大的点
   k = min(k, n);
   //x ∈ i ... k : n / x == j

   ...
   i = k;
}
//    22/1  22
//    22/2  11
//    22/3  7
//    22/4  5
//    22/5  4
      22/6  3    22/x = y  x是 22/k = y最小的k
                 22/y = x' x'是 22/k = y最大的k
      22/7  3
      22/8  2
      ...
      22/11 2

这道题的代码(简化版):

for (register int i = 1; i <= n && i <= k; ++i) {
    int res = k / i;
    int lst = k / res;
    lst = min(lst, n);
    ans -= 1ll * res * (i + lst) * (lst - i + 1) / 2;
    i = lst;
}

注意!!

除法分块中,一定要保证i<=n,否则res = n / i出现0时,lst = n / i会RE!!

线性推逆元/阶乘逆元

注意:

一、推阶乘逆元的时候先要求的是jie[up] 的逆元,而不是 up 的逆元!!

二、注意推阶乘逆元的上届:不能超过模数 P,否则阶乘逆元将会推出一堆0!!(毕竟0无逆元)

一些数论神题

P2150 [NOI2015]寿司晚宴

n <= 20 的暴力这里就不说了 ,就是给简单的状压。

n <= 100,我们发现这 n-1 个数的质因数不会很多,因此我们把每个数看作一个质因数集合,要求选的两个集合无交集。

我们依然选择 状压。

设定 f[i][s1][s2] 为:考虑前 i 个数,其中甲(指小G)选的集合为 s1,乙(指小W)选的集合为 s2,只有当s1和s2的交集为空时状态合法。

然后考虑推 i时第i个集合的贡献,就有:

f[i][s1 | state_i][s2]+=f[i-1][s1][s2] f[i][s1][s2 | state_i]+=f[i - 1][s1][s2]

然后发现第一维可以滚动数组优化,或者像背包那样优化。

这种方法只适用于 state 不会很大的时候。

n <= 500,因为每个数的质因数集合中超过 \sqrt{n} 的最多只有一个,因此我们可以将那些大质因数相同的一块考虑,强制他们不被同一个人选。

具体实现:

先将数按大质因数排序(无 则为-1),这样大质因数就成段分布了。

定义 f[s1][s2] 为甲(指小G)选的集合为 s1(可为空),乙(指小W)选的集合为 s2 的方案数(可为空),f1[s1][s2] 为强制甲选当前段的大质因数,...的方案数,

每段开头:把 $f[][]$ 复制(赋值)给 $f1[][],f2[][]$。 每段结束:汇总 $f[][]$:$f[s1][s2] = f1[s1][s2] + f2[s1][s2] - f[s1][s2]$,其中减去 $f[s1][s2]$ 是因为两人都空集会算两遍。 答案:$\sum_{i, j}f[i][j]~ * ~ [i ~ and ~ j = 0]

关键代码:

inline void dp() {
    f[0][0] = 1;
    for (register int i = 1; i < n; ++i) {
        if (i == 1 || nd[i].mx != nd[i - 1].mx || nd[i].mx == -1) {
            memcpy(f1, f, sizeof(f));
            memcpy(f2, f, sizeof(f));
        }
        register int st = nd[i].state;
        for (register int s1 = All; ~s1; --s1) {
            for (register int s2 = All; ~s2; --s2) {
                if (s1 & s2)    continue;
                if ((st & s2) == 0) f1[s1 | st][s2] = (f1[s1 | st][s2] + f1[s1][s2]) % P;
                if ((st & s1) == 0) f2[s1][s2 | st] = (f2[s1][s2 | st] + f2[s1][s2]) % P;
            }
        }
        if (i == n - 1 || nd[i].mx != nd[i + 1].mx || nd[i].mx == -1) {
            for (register int s1 = All; ~s1; --s1) {
                for (register int s2 = All; ~s2; --s2) {
                    if (s1 & s2)    continue;
                    f[s1][s2] = (f1[s1][s2] + f2[s1][s2] - f[s1][s2] + P) % P;
                }
            }
        }
    }
}

类似题目:T121257 简单的数学题 (团队题目)