ARTICLE DETAIL

资讯详情

深耕网站建设与运营推广的一线实战洞察。

有限域不可约多项式与本原多项式:判定、手算与Python枚举

有限域不可约多项式与本原多项式:判定、手算与Python枚举 不可约多项式和本原多项式这两个词第一次撞见基本都在有限域、纠错编码或者 LFSR 的教材里公式一摆、符号一多人就开始飘。我第一次做 GF(2^8) 上的乘法器时随手从网上抄了个八次多项式当模仿真跑起来乘法逆元全错寄存器值满屏乱飞查了大半天才发现那个多项式是可约的——它连一个有限域都撑不起来谈何求逆。从那以后我养成一个习惯凡是涉及不可约多项式、本原多项式的选型一律自己算一遍、跑一遍代码再往工程里填绝不信手里的表。这篇内容就是把这件事说透。我会先把不可约和本原这两个概念钉死再用大量手算例子把判定过程走一遍包括 GF(2) 和 GF(3) 上的一堆具体多项式然后给出一套完整的 Python 实现让你不查表也能枚举任意次数的不可约多项式和本原多项式。如果你在做 CRC、LFSR、Reed-Solomon、AES 的 S 盒、随机数发生器或者只是被教材上的定理绕晕了这篇应该能帮你省下几个小时。全文计算过程我都写全了你可以拿纸笔跟着对一遍。1. 先把概念钉死不可约多项式到底不可约在哪1.1 从没有实根这个直觉说起不可约多项式的定义课本上一句话在域 F 上一个次数大于等于 1 的多项式 f(x)如果它不能写成两个次数都小于 deg(f) 且大于 0 的多项式的乘积就叫在 F 上不可约。这个定义里有两个关键词最容易被忽略一个是在哪个域上另一个是次数都小于 deg(f) 且大于 0。这两个点直接决定了后面所有的坑。先说在哪个域上。x^21 在实数域上不可约因为没有任何实数 x 让 x^210判别式是负的。但放到复数域上x^21(xi)(x-i)就变成可约的了。同样的道理x^2-2 在有理数域上不可约但在实数域上可约成 (x-√2)(x√2)。所以不可约永远是一个相对概念离开域的定义谈不可约是没有意义的。我们后面讨论的域主要是有限域 GF(p) 和 GF(p^n)其中最常用的是 GF(2)。再说次数都小于 deg(f) 且大于 0。这里排除了乘一个常数的情况。比如在 GF(2) 上x^2x x·(x1)两个因子都是 1 次所以 x^2x 可约。但 2x^22x 这种写法在 GF(2) 里就等于 0根本不用讨论。换成 GF(3)3x^2 这种常数倍数不算因式分解因为常数不影响能不能再分常数是单位元不是真正的因子。我见过最常见的误判就是有人把在 GF(2) 上没有根直接等同于不可约。这个等价关系只在次数不超过 3 的时候成立。到了 4 次及以上就开始失效因为 4 次多项式可以拆成两个 2 次的乘积而这两个 2 次因子各自都没有一次因子所以整个 4 次多项式在 GF(2) 上照样没有根却明明是可约的。最经典的例子就是 f(x) x^4 x^2 1。代 x0 进去得到 1代 x1 进去得到 1111两个取值都不是 0所以在 GF(2) 上没有根。但它实际上等于 (x^2x1)^2展开一下x^4 x^2 1中间项 2x^3、2x^2 在 GF(2) 里全被 2 消掉2 次项剩一个 x^2。这就是一个标准的 4 次可约、却无根的例子。类似的还有 x^41在 GF(2) 上它等于 (x1)^4因为 (x1)^2 x^21再平方就是 x^41。这两个例子建议你记牢面试和排错都用得上。提示判断无根只是排除了一次因子可约还包含两个高次因子相乘的情况。次数大于 3 时无根是最低门槛不是充分条件。1.2 本原多项式是不可约多项式里的优选生不可约多项式只是分不开的多项式而本原多项式是在此基础上再挑一层出来它不仅不可约而且它的某个根还是扩域乘法群的生成元。具体说设 f(x) 是 GF(q) 上次数为 n 的不可约多项式α 是它的一个根那么 α 一定落在扩域 GF(q^n) 里。GF(q^n) 去掉零元之后是一个乘法群这个群的阶是 q^n - 1。如果 α 的阶恰好等于 q^n - 1也就是说 α 的幂能遍历所有非零元素那么 α 就是一个本原元f(x) 就叫做 GF(q) 上的 n 次本原多项式。反过来说不可约多项式的根 α 的阶一定是 q^n - 1 的因子但不一定等于它。只有当这个阶刚好顶满 q^n - 1 时才升级为本原多项式。所以本原多项式一定不可约不可约多项式的却不一定本原。这个方向千万不能记反。举一个具体的例子。在 GF(2) 上x^4x^3x^2x1 是 4 次不可约的但它的根 α 满足 α^51阶只有 5而 2^4-1155 只是 15 的因子。所以它是不可约但非本原。而 x^4x1 和 x^4x^31 这两个多项式的根阶都是 15它们都是本原的。这个区别在工程上非常要命。用 LFSR 做伪随机序列时特征多项式必须是本原的才能保证序列周期达到最大长度 2^n-1如果只保证不可约周期可能只有 2^n-1 的某个因子序列会提前重复随机性和统计特性直接崩掉。反过来做有限域乘法的时候只要不可约就够了因为乘法逆元的存在性只依赖域结构不依赖是不是本原。所以选型之前先问自己一句我要的是能构成域还是要有满周期这个问题的答案决定了你得验到哪一层。1.3 一张表看清两者关系与计数规律不可约多项式和本原多项式的个数都有闭式公式这在做枚举和交叉验证时特别有用。GF(q) 上首一 n 次不可约多项式的个数是N_q(n) (1/n) · Σ_{d|n} μ(d) · q^(n/d)其中 μ 是莫比乌斯函数取值规则是μ(1)1μ(素数)-1μ(素数平方)0μ(两个不同素数之积)1。GF(q) 上 n 次本原多项式的个数是P_q(n) φ(q^n - 1) / nφ 是欧拉函数。注意本原多项式的个数只跟 q^n - 1 的因式分解有关跟 n 的因数关系不大这一点很多人第一次看会觉得别扭。把 q2 代进去算一张常用表n不可约个数 N_2(n)本原个数 P_2(n)本原占比12150%211100%322100%43267%566100%69667%71818100%8301653%9564886%10996061%1118617695%1233514443%164080204850%表里 n2、3、5、7 这几个本原占比 100% 的次数都是因为 2^n-1 是素数。2^2-13、2^3-17、2^5-131、2^7-1127 都是素数而 q^n-1 为素数时除 1 以外的任何元素阶都是 q^n-1所有不可约多项式的根都自动是生成元。这就是为什么小次数的效果看起来很整齐一旦到了 n82^8-12553×5×17 有多个因子情况立刻复杂起来。拿 n4 手算一遍确认公式不是骗人的。n4 的因子是 1、2、4μ(1)1、μ(2)-1、μ(4)0所以N_2(4) (1/4) × [1×2^4 (-1)×2^2 0×2^1] (1/4) × [16 - 4] 3确实只有 3 个 4 次不可约多项式。本原的个数2^4-115φ(15)φ(3)×φ(5)2×48P_2(4)8/42正好 2 个。而 GF(2) 上 4 次首一多项式总共只有 2^416 个可约的有 13 个不可约 3 个本原 2 个。这个比例在选型时挺有用的你以为随手一抓就不可约其实命中率只有两成不到。再看 n8N_2(8) (1/8)[2^8 - 2^4] (256-16)/8 30。P_2(8)2553×5×17φ(255)255×(2/3)×(4/5)×(16/17)128128/816。30 个不可约里只有 16 个本原另外 14 个的根阶分别是 85、51、17、15 这些因子。所以如果你从网上随便抓一个八次不可约多项式当 LFSR 的特征多项式有超过一半的概率拿到的是非本原的周期直接腰斩甚至只有原来的十几分之一。2. GF(p)上判定不可约的三条实用路线2.1 低次2 次、3 次直接试根就够次数是 2 或 3 的时候判定逻辑最简单把域里所有元素代进去只要有一个取值为 0就说明有一次因子可约如果全都不为 0就不可约。道理很直白。一个 2 次多项式如果能分解只能是1 次 × 1 次那就必须有一个一次因子 (x-a)也就意味着 f(a)0。3 次同理能分解的话必然是 1 次乘 2 次同样要有一个根。所以 2 次、3 次不存在两个高次因子相乘的退化情况试根法在这里是完全充分的。在 GF(2) 上只有 0 和 1 两个元素要试。所有 2 次首一多项式有 4 个x^2、x^21、x^2x、x^2x1。逐个代x^2代 0 得 0有根可约。x^21代 1 得 110有根可约等于 (x1)^2。x^2x代 0 得 0有根可约等于 x(x1)。x^2x1代 0 得 1代 1 得 1111都非 0不可约。所以 GF(2) 上只有 x^2x1 一个 2 次不可约多项式跟公式 N_2(2)(1/2)(4-2)1 完全对得上。3 次在 GF(2) 上有 8 个首一多项式其中能分解的占 6 个剩下 2 个不可约x^3x1 和 x^3x^21。这两个互为反序多项式系数顺序倒过来在 GF(2) 上这种配对很常见后面枚举时你会发现大量这种情况。换到 GF(3) 上做一遍体会一下 p2 的感觉。GF(3) 的元素是 0、1、2。所有 2 次首一多项式形如 x^2bxcb、c 各有 3 种取法共 9 个。逐个筛多项式x0x1x2判定x^2011可约x^21122不可约x^22200可约x^2x020可约x^2x1101可约x^2x2212不可约x^22x002可约x^22x1110可约x^22x2221不可约三个不可约x^21、x^2x2、x^22x2。跟公式 N_3(2)(1/2)(9-3)3 一致。这里有个小技巧值得说在 GF(p) 上用试根法时与其老老实实代 p 个值不如先算一下判别式。对 2 次多项式 x^2bxc它不可约当且仅当 b^2-4c 不是 GF(p) 中的平方数且 p 为奇素数。在 GF(3) 里0^20、1^21、2^21平方数集合是 {0,1}。x^21 的判别式是 0-4-4≡2 mod 32 不在 {0,1} 里所以不可约x^2x2 的判别式是 1-8-7≡2 mod 3同样不在平方数集合里不可约x^22x2 的判别式是 4-8-4≡2 mod 3也不可约。三次都能秒判不用逐个代入。这个判据只对 2 次有效3 次以上就不成立了。注意试根法在次数 ≥4 时只是必要条件。真正做工程的时候我通常先用试根法快速筛掉一大半可约的剩下的再上 Rabin 或者直接跑代码。2.2 4 次以上先看阶再用分圆多项式拆结构次数到 4 以上试根法就不够用了得换个思路。一个非常有效的观察角度是分圆多项式。赛道上有个常用结论在 GF(q) 上x^m - 1 的因式分解跟 m 的素因子以及 q 在模 m 下的乘法阶密切相关。具体说x^m-1 分解出的每个不可约因子其次数都等于 ord_m(q)也就是 q 模 m 的乘法阶m 与 q 互素时。所有不可约因子次数相同个数就是 φ(m)/ord_m(q)。拿 GF(2) 上的 x^5-1 举例。在 GF(2) 里减法和加法一样x^5-1 x^51。因式分解出来是x^5 1 (x1)(x^4x^3x^2x1)m5q2ord_5(2) 是让 2^k ≡ 1 (mod 5) 的最小 k。2^12、2^24、2^38≡3、2^416≡1所以 ord_5(2)4。φ(5)44/41说明除了那个一次因子 (x1) 之外剩下的是一个 4 次不可约多项式。这就直接证明了 x^4x^3x^2x1 在 GF(2) 上不可约完全不需要试根。这个方法特别好用因为它把判定不可约变成了算一个模意义下的乘法阶。再举 x^71 的例子。m7ord_7(2)2^12、2^24、2^38≡1所以 ord_7(2)3。φ(7)66/32说明除了 (x1) 之外剩下 6 次的部分 x^6x^5x^4x^3x^2x1 会拆成两个 3 次不可约多项式。实际拆出来是x^6x^5x^4x^3x^2x1 (x^3x1)(x^3x^21)正好是前面提到的那两个 3 次不可约多项式完全吻合。再看 9 次的情况这个例子后面讲本原的时候还要用。x^91 在 GF(2) 上分解涉及的是 Φ_9(x) x^6x^31 这个九次分圆的一部分。ord_9(2)2^12、2^24、2^38≡-1、2^416≡7、2^514≡5、2^610≡1所以 ord_9(2)6。φ(9)66/61说明 Φ_9 本身就是一个 6 次不可约多项式。也就是 x^6x^31 不可约。顺手再算一个 21 次的为后面 6 次多项式的阶分布做准备。ord_21(2)2^38、2^664≡1 (mod 21)中间 2^12、2^24 都不等于 1所以 ord_21(2)6。φ(21)1212/62说明 Φ_21 会拆成两个 6 次不可约多项式。这两只的根阶都是 21。所以 6 次不可约多项式一共有三组来源两组来自 Φ_21阶 21一个来自 Φ_9阶 9再加上本原的那 6 个阶 63正好凑出 N_2(6)2169。这个分解思路在做枚举验证时非常有用可以拿它当标准答案去对代码输出。2.3 Rabin 判定法一条公式管到底附手算全过程分圆法靠的是恰好能对上 x^m-1这个条件遇到不对应的多项式就不好使了。要一个普适的判定得请出 Rabin 不可约性测试。定理是这么说的设 f 是 GF(q) 上次数为 n 的多项式f 不可约当且仅当下面两条同时成立一是 x^(q^n) ≡ x (mod f) 二是对 n 的每个素因子 d都有 gcd(x^(q^(n/d)) - x, f) 1。第一条保证 f 的根都在 GF(q^n) 里第二条保证根不会掉进任何一个真子域 GF(q^(n/d)) 里。两条合起来说明 f 的所有根都是真正的 n 次代数元也就是 f 不可约。先用它验一个已知答案f x^4x^3x^2x1n4q2。n4 的素因子只有 d2n/d2第二步需要算 gcd(x^(2^2) - x, f) gcd(x^4x, f)。先算第一步。记 f 的关系式x^4 x^3x^2x1因为 x^4x^3x^2x10移项即可GF(2) 上加减一样。x^4 x^3x^2x1 x^5 x·x^4 x^4x^3x^2x (x^3x^2x1)x^3x^2x 1x^5 居然等于 1这个结果很关键。继续 x^8 x^5·x^3 x^3 x^16 x^15·x (x^5)^3·x 1·x x所以 x^16 ≡ x (mod f)第一条满足。再算 gcd。x^4x x(x^31) x(x1)(x^2x1)它的所有不可约因子是 x、x1、x^2x1。f 是 4 次如果它跟 x^4x 有公因子公因子只能是 x^2x1因为 f 没有一次因子前面试根已经确认。而 (x^2x1)^2 x^4x^21 ≠ f所以 x^2x1 不整除 f两者互素gcd1第二条也满足。结论f 不可约。跟分圆法的判断一致两条路殊途同归。再用它验一个可约的例子看看它是怎么抓住漏洞的。取 f x^4x^21 (x^2x1)^2n 还是 4d2。第一步f 的关系式是 x^4 x^21。 x^4 x^21 x^8 (x^4)^2 (x^21)^2 x^41 (x^21)1 x^2 x^16 (x^8)^2 (x^2)^2 x^4 x^21我们需要 x^16 ≡ x也就是 x^21 ≡ x显然不成立两者差 x^2x1非零。所以第一条就直接挂了根本不用算第二步。这个对比很有意思本原与否会影响 x 的阶但可约与否会直接让 x^(q^n) ≡ x 这个整体条件失效。注意 x^16 ≡ x 这一条本质上是要求f 的所有根都落在 GF(q^n) 里对可约多项式来说只要有一个因子的次数不整除 n这个条件就崩。而 gcd 那一条是专门用来排除根的次数是 n 的真因子这种情况的。举个需要靠第二条才能抓住的例子。取 f (x^2x1)(x^2x1)不行刚用过。换 f x^4x^3x^2x1 的兄弟四次的另一个可约情形f (x^2x1)(x^21)。展开x^4x^2x^3xx^21 x^4x^31……等一下x^2x^20x 剩一个所以是 x^4x^3x1。验一下 x111110有根是可约的这个例子太容易被试根法抓住。要构造一个第一条满足、第二条不满足的最典型的是 f 本身不可约但次数不是 n……那不可能。真正的场景是 f 不可约且次数为真因子。比如取 f x^2x12 次不可约如果我们错把它当成 4 次来测n 会算错实际它满足 x^4 ≡ x (mod f)因为 x^31x^4x。而此时 n4 的 d2n/d2要算 gcd(x^4x, f)。x^4x mod f xx 0gcd 就是 f 本身不等于 1第二条直接判死。这就是第二条的作用它把根的阶太小、落在了真子域里这种情况筛掉。手算 Rabin 的时候有两个提速技巧我自己常用。一个是别硬算 x^(q^n)一路用平方递推x^2 → 平方得 x^4 → 平方得 x^8 → 平方得 x^16每一步做完立刻用 f 化简。另一个是 gcd 那一步别真的去展开 x^(q^(n/d)) 这个多项式直接在商环里算出它的化简结果再求 gcd因为 gcd 只跟化简结果有关。次数一高这两点能省掉大量时间。3. 本原多项式判定把阶算明白就赢了一半3.1 阶的三个性质记住就不会错判定本原的核心就是算一个数根 α 的乘法阶。围绕它有三条性质必须刻在脑子里。第一条阶一定整除 q^n - 1。这不是巧合是拉格朗日定理的直接推论GF(q^n) 的非零元素构成一个阶为 q^n-1 的乘法群群中任何元素的阶都整除群的阶。所以算阶的时候不用一个个往上试只需要在 q^n-1 的因子集合里找。第二条不可约多项式所有根的阶都相同。这一点非常省事。f 的 n 个根是彼此共轭的α、α^q、α^(q^2)、…、α^(q^(n-1))共轭元素的阶必然相等。所以你随便取哪个根算阶结果都一样也正因为如此我们才能说这个多项式的阶而不用特指是哪个根。第三条本原等价于阶等于 q^n - 1。不是能整除且最大是严格的相等一个字都不能差。由此可以得到一个很实用的判定流程先确认不可约然后算出 x 在商环 GF(q)[x]/(f) 中的阶如果等于 q^n-1 就是本原。之所以能拿 x 代替 α 来算是因为 α 就是 x 在商环里的像两者阶完全一致。有个细节我踩过坑判断 x^k 是不是 1 的时候别忘了 k 要最小的那个。有人算到 x^151 就宣布本原完全没检查 x^3、x^5 是不是也等于 1。万一是 x^51那阶就是 5 而不是 15结论直接反了。正确的做法是先确认 x^(q^n-1)1再对 q^n-1 的每个素因子 p验证 x^((q^n-1)/p) ≠ 1。全部通过阶才等于 q^n-1。3.2 手算阶的完整流程以 x^4x^31 为例光说流程太干直接上完整的幂表。取 f x^4x^31在 GF(2) 上n4目标阶是 2^4-115。关系式x^4 x^31。逐次往上推kx^k 的化简结果说明011x2x^23x^34x^31用关系式5x^3x1x·x^4 x^4x6x^3x^2x17x^2x18x^3x^2x9x^2110x^3x11x^3x^2112x113x^2x14x^3x^2151回到单位元x^15 确实等于 1而且从表里能直接看出来1 到 14 次幂没有任何一个等于 1。所以阶就是 15等于 2^4-1f 是本原多项式。再看一个反面例子f x^4x^3x^2x1。前面已经算过 x^4 x^3x^2x1一步就得到 x^5 1。所以阶是 5而 2^4-1155≠15不是本原的。这个多项式在 LFSR 里用起来周期只有 5序列短得可怜。顺手把另一个本原的 4 次多项式 x^4x1 的幂表也贴出来方便你对照kx^k 化简结果011x2x^23x^34x15x^2x6x^3x^27x^3x18x^219x^3x10x^2x111x^3x^2x12x^3x^2x113x^3x^2114x^31151两张表对比着看你会发现 x^7 和 x^7 就不一样了在 x^4x1 里 x^7x^3x1在 x^4x^31 里 x^7x^2x1。这说明两个本原多项式虽然都撑起 GF(16) 的乘法群但元素的编号也就是离散对数表是不同的。所以工程里换多项式就等价于换了一套对数表S 盒、CRC 表全得重算。这也是为什么我在项目里从不轻易改特征多项式。三个阶段性的检查点总结一下先看 x^(q^n-1) 是否等于 1再对 q^n-1 的每个素因子 p 检查 x^((q^n-1)/p) 是否不等于 1两条都过阶就是满的。以 15 为例153×5素因子是 3 和 5所以要额外验 x^5≠1 和 x^3≠1。15/35、15/53都在表里能查到x^5x^3x1≠1、x^3x^3≠1通过。3.3 6 次不可约多项式的分歧63、21、9 三种阶前面用分圆法算过GF(2) 上 6 次不可约多项式有 9 个其中 6 个本原阶 632 个阶为 211 个阶为 9。这个分布是理解不可约不等于本原最好的教材。2^6-1633^2×7。63 的因子有 1、3、7、9、21、63。一个 6 次不可约多项式的根阶必须是 63 的因子。理论上可能的阶是 3、7、9、21、63但受次数必须等于 ord_m(2) 整除性的约束实际能出现的就是 63、21、9 三种。阶为 9 的那一个是 x^6x^31也就是 Φ_9(x)。验证一下它的阶设 α 是根α^6 α^31那么 α^9 α^3·α^6 α^3(α^31) α^6α^3 (α^31)α^3 1。确实 α^91。而且 α^3≠1否则 α^6 α^31 11... 等等如果 α^31那 α^6α^31 111 1 ≠ 0矛盾α^1≠1所以阶就是 9。阶为 21 的那两个来自 Φ_21 的分解它们的幂次只会落在 21 个值上永远到不了 63。这两个多项式具体是什么不重要重要的是这个事实同样是 6 次不可约多项式用它们做 LFSR周期会从 63 掉到 21 或 9掉了三分之二还多。我在实际项目里见过更惨的有人用了一个 16 次的多项式做流密码的驱动序列2^16-1655353×5×17×257阶的因子很多。他挑的那个多项式的阶只有 255周期短了两个数量级测试时看着随机性还行实际在长时间抓包下序列重复得非常明显。这类问题在短测试里根本暴露不出来只能靠提前把阶算清楚。所以我的建议是任何人让你用一个特征多项式先问一句话它的阶是多少答不上来的就用代码跑一遍再谈。3.4 三项式、五项式怎么选以及 8 次为什么没有三项式选本原多项式的时候工程上有两个偏好项数越少越好最高次项系数固定为 1中间只有少量非零系数。原因是硬件上每个非零系数对应一个异或门或者一个抽头项数越少LFSR 的反馈逻辑越省资源软件上计算也越快。系数只有三项的形如 x^n x^k 1叫三项式是最理想的形态。但三项式不存在万能的情况。一个基本事实是如果 n 是偶数任何 x^n x^k 1 只要 k 也是偶数它就是完全平方必然可约。比如 x^8x^41(x^4x^21)^2、x^8x^21(x^4x1)^2、x^8x^61(x^4x^31)^2全是这个套路。所以在 GF(2) 上要找不可约三项式k 必须是奇数。更麻烦的是有些次数上根本不存在不可约三项式。8 次就是最著名的例子13 次、16 次、19 次这些也在列。这一点你自己跑一遍枚举就能确认代码在后面第 4 节。所以做 8 位 LFSR 的时候市面上流传的抽头组合没有一个三项式的全得用五项式比如 x^8x^4x^3x^21、x^8x^6x^5x^41 这类。这里我要提醒一句网上流传的本原多项式表质量参差不齐我遇到过至少两次抄错系数的情况。有的是印刷错误有的是把不可约当成本原列进去了。所以我的习惯是候选多项式一律先用代码验一遍不可约和本原两个性质验过了再往 RTL 里写。花五分钟验证比在板子上调三天强得多。还有一个容易被忽略的取舍项数少不等于性能好。在软件实现里非零项的分布位置会影响缓存和指令流水有时候一个五项式的实际跑分反而比三项式更稳。这个只能实测别光看公式。4. 从零写代码验证不查表也能自己算4.1 GF(2) 多项式运算的位运算实现前面全在纸上推现在把它工程化。GF(2) 上的多项式有个天然的好处系数只有 0 和 1可以用一个整数的二进制位直接表示。比如 0b1011 表示 x^3x1bit i 对应 x^i 的系数。加法和减法在 GF(2) 上都是异或乘法和取模也都是位运算写起来非常短。先把基础函数搭起来# GF(2) 多项式用整数位掩码表示bit i x^i 的系数 # 例如 0b1011 表示 x^3 x 1 def deg(a): return a.bit_length() - 1 def poly_mod(a, m): 在 GF(2) 上求 a mod m加减统一用异或 dm deg(m) while a and deg(a) dm: a ^ m (deg(a) - dm) return a def poly_mul(a, b, mNone): 多项式乘法给了 m 就在商环里做 r 0 while b: if b 1: r ^ a b 1 a 1 return r if m is None else poly_mod(r, m) def poly_gcd(a, b): 欧几里得算法求最大公因式 while b: a, b b, poly_mod(a, b) return apoly_mod用的是长除法只要被除式的次数不低于除式次数就把除式左移相应位数再异或回去。这里能直接用异或代替减法是因为 GF(2) 上 1-10 和 110 完全一样。这个特性也解释了为什么代码这么短。poly_mul是标准的俄罗斯农民乘法乘数逐位右移被乘数逐位左移遇到 1 就累加。注意这里我用的是边乘边不进模最后统一取模对小次数完全够用。如果要做 32 次以上的多项式建议改成边乘边取模否则中间结果会膨胀得很大。再补一个快速幂Rabin 判定和阶的计算都要用def poly_powmod(base, e, m): base^e mod m快速幂 r 1 base poly_mod(base, m) while e: if e 1: r poly_mul(r, base, m) base poly_mul(base, base, m) e 1 return r4.2 枚举所有 n 次不可约多项式并与公式对照有了基础运算不可约判定可以先写个朴素版把所有次数不超过 n/2 的首一多项式都拿来试除一遍。逻辑上无脑但胜在正确性一目了然。def is_irreducible(f): 朴素试除法用所有次数 deg(f)/2 的首一多项式试除 n deg(f) if n 0: return False for d in range(1, n // 2 1): for g in range(1 d, 1 (d 1)): # d 次首一多项式全集 if poly_mod(f, g) 0: return False return Truerange(1 d, 1 (d1))这个写法值得说一下。所有次数恰好为 d 的首一多项式最高位一定是 1所以最低的值是 1d最高是 1(d1)-1正好是这么个左闭右开区间。这个技巧在后面枚举时反复用到。跑一遍枚举把所有 n 次不可约多项式列出来for n in range(1, 13): irr [f for f in range(1 n, 1 (n 1)) if is_irreducible(f)] print(n, len(irr), [bin(f) for f in irr])n4 时会输出 3 个0b10011 (x^4x1)、0b11011 (x^4x^31)、0b11111 (x^4x^3x^2x1)。跟我前面手算的完全一致。n6 时会输出 9 个。n8 时会输出 30 个。这些数字跟公式 N_2(n) 算出来的就一一对上了这也是验证代码有没有写错的最好方式。顺便说一句性能。这个朴素版本枚举 12 次多项式时每个候选要试除大约 Σ(2^d) ≈ 2^(n/21) 个除式总计算量对 n≤16 完全能接受跑起来就是几秒钟。真要上 20 次以上就得换成 Rabin 或者 Berlekamp别硬扛。4.3 用 Rabin 和阶来判本原性不可约判完了接着判本原。先把 Rabin 实现出来它比试除法快得多def prime_factors(n): fs, d set(), 2 while d * d n: while n % d 0: fs.add(d) n // d d 1 if n 1: fs.add(n) return fs def rabin_irreducible(f): Rabin 不可约判定仅针对 GF(2) n deg(f) if n 1: return False # 条件一x^(2^n) ≡ x (mod f) if poly_powmod(2, 1 n, f) ! poly_mod(2, f): return False # 条件二对每个素因子 dgcd(x^(2^(n/d)) x, f) 1 for d in prime_factors(n): h poly_powmod(2, 1 (n // d), f) ^ 2 # 异或 2 就是减 x if poly_gcd(poly_mod(h, f), f) ! 1: return False return True这里有个实现细节要注意条件二里的减 x在 GF(2) 上就是异或 2因为 x 的位掩码是 0b10。写别的域的时候要老老实实用减法GF(2) 才能这么省事。再看阶的计算这是判本原的核心def poly_order(f): 求 x 在 GF(2)[x]/(f) 里的阶返回 None 表示 f 不可约性有问题 n deg(f) cur, k 1, 0 while True: cur poly_mul(cur, 2, f) # 每次乘一个 x k 1 if cur 1: return k if k (1 n): return None def is_primitive(f): n deg(f) return is_irreducible(f) and poly_order(f) (1 n) - 1最后把两者串起来把前面那张表跑出来for n in range(1, 13): irr [f for f in range(1 n, 1 (n 1)) if is_irreducible(f)] pri [f for f in irr if poly_order(f) (1 n) - 1] print(fn{n:2d} 不可约{len(irr):3d} 本原{len(pri):3d})输出会是 1:2/1、2:1/1、3:2/2、4:3/2、5:6/6、6:9/6、7:18/18、8:30/16、9:56/48、10:99/60、11:186/176、12:335/144跟第 1.3 节那张表一字不差。到这一步你已经有了一个完全可以自证的验证工具以后不管谁给你一个多项式三十行代码就能验完再也不用求人。有个小优化值得提一句。poly_order是老老实实一次乘一个 x对 n16 要循环 65535 次还行n32 就是 40 亿次跑不动了。真要做大次数应该改成先算 x^(2^n-1) 确认是 1再对 2^n-1 的每个素因子 p 验证 x^((2^n-1)/p)≠1这样只需 O(log n × 素因子个数) 次快速幂。这是把前面 3.1 节的理论直接用上了。4.4 用现成库交叉验证别只信自己的代码自己写实现最大的风险是错得一致——逻辑本身有 bug但输出内部自洽你反而更放心。所以我一定会拿现成库对一遍。用 sympy 验不可约性最省事from sympy import Poly, GF, symbols x symbols(x) f Poly(x**8 x**4 x**3 x**2 1, x, domainGF(2)) print(f.is_irreducible) # True用 galois 库可以直接把不可约和本原的列表一起拉出来import galois print(galois.irreducible_polys(2, 8, reverseTrue)) print(galois.primitive_polys(2, 8, reverseTrue))如果你有 Sage 环境那更简单一行搞定R.x PolynomialRing(GF(2)) f x^8 x^4 x^3 x^2 1 print(f.is_irreducible(), f.is_primitive())我习惯的做法是用自己写的代码生成一份列表用库再生成一份两边取集合做差。只要差集为空就说明实现没问题。这个方法帮我抓出过两次 bug一次是poly_mod里忘了处理 a 归零的情况导致死循环一次是阶的循环上限设成了 2^n 而不是 2^n-1边界上多算了一轮。这种错误光靠看代码很难发现必须交叉验证。提示不同库的版本 API 可能略有差异尤其是 galois 的返回类型和参数名。以上代码以常见版本为准跑之前建议先确认一下自己环境里的函数签名。5. 踩坑记录与速查表5.1 五个最容易翻车的点第一个把无根当成不可约。前面反复强调过了次数 ≥4 时这个等价关系不成立。典型的反例 x^4x^21(x^2x1)^2。我现在看新人代码的习惯动作就是搜一遍有没有只判根就返回不可约的实现十个里面能抓到三个。第二个阶算错。最常见的两种错法一种是查到 x^k1 就收工没确认 k 是不是最小的另一种是循环边界写错把 q^n-1 写成 q^n 或者 q^n-2。这两种错误都不会报异常只会让结论静悄悄地错掉。我的做法是在代码里加一条断言算出来的阶必须整除 q^n-1不整除就说明实现有问题。第三个特征多项式选成了不可约但非本原的。这个坑最隐蔽因为仿真跑起来一切都正常只是周期比预期短。判断办法很简单如果 LFSR 的周期不是 2^n-1先怀疑多项式的阶。第四个混淆 p 和 n。GF(p^n) 里 p 是特征是模的素数n 是扩张次数是多项式的次数。这两个量在很多公式里同时出现抄公式的时候特别容易串位。我一般会在草稿上把每个符号的含义写一遍再代入。第五个跨域使用结论。GF(2) 上的表不能直接搬到 GF(3) 或者 GF(2^8) 上用。比如判别式是不是平方数这个判据只对奇素数 p 的 2 次多项式成立在 GF(2) 上根本不适用因为 2 是偶素数判别式公式退化。换域就得重算没有捷径。5.2 症状—原因—处理速查表现象可能原因处理方式有限域乘法求逆失败或结果异常模多项式可约商环不是域用 Rabin 或试除确认不可约LFSR 周期远小于 2^n-1特征多项式不可约但非本原计算 x 的阶换成本原多项式代码判不可约但试除能拆开只做了试根漏了高次因子补上 d≤n/2 的完整试除阶算出来不整除 q^n-1循环边界或关系式化简写错加整除断言重查 poly_mod同一多项式在不同资料里结论冲突一份资料把不可约当本原了自己跑一遍代码以计算结果为准8 次找不到三项式该次数确实不存在不可约三项式改用五项式如 x^8x^4x^3x^21换了特征多项式后 CRC/S 盒全错换多项式等于换了整个对数表所有派生表全部重新生成这张表基本覆盖了我自己遇到过的所有情况。最后一条尤其值得强调模多项式和它派生出来的表是强绑定的改一个系数整张表就全变了不存在改一点点、影响一点点的情况。5.3 我自己用的验证流程跑过这么多项目之后我形成了一个固定动作几乎每次都按这个顺序走一遍。第一步拿到候选多项式先转成位掩码跑is_irreducible。这一步过滤掉最基础的错误。第二步跑poly_order拿到 x 的阶跟 2^n-1 比。相等才算过。第三步用库交叉验证一遍确认自己的代码没写歪。第四步如果这个多项式要用在硬件上再检查一下非零项的位置和数量评估异或门的开销。第五步如果它要派生 S 盒、CRC 表或者对数表生成完之后再做一次一致性自测比如验证任意元素乘它的逆等于 1。这五步加起来可能要多花二十分钟但相比在板子上调三天的成本完全值。我印象最深的一次是做一个 16 位 CRC掐着时间上线多项式是从一份老文档里直接抄的。结果 CRC 的检错能力明显不达标突发错误测试通过率异常最后发现那个多项式压根就是可约的。从那以后我再也没跳过第一步。最后分享一个小习惯我会把每次项目里用到的本原多项式、它的阶、验证时间记在一个自己的小本子上同时把对应的枚举代码段一起存下来。用的时候直接拿出来对一遍比重新推一遍快得多也比查网上的表放心得多。次数多了之后你会发现常用的就那么几十个比如 4 次用 x^4x1、8 次用 x^8x^4x^3x^21、16 次用 x^16x^12x^3x1 这类手熟了之后看一眼系数就知道是不是三项式、奇偶性对不对直觉会变得很准。
返回列表