F(2^131)下有限域基本运算的实现(c++)_有限域运算编程

有限域实验

实验目标

实现$F_{2^{131}}$上的基本运算，并计算学号的逆元
不可约多项式为$f(x)=x^{131}+x^{13}+x^2+x+1$

算法原理

加法

在$F_{2^{131}}$下，对于定义在该域上的任意多项式$g(x)=\sum_{i=0}^{130}a_ix^i$，由于$a_i$只能取0或1，可将g(x)用长度为131的比特串来表示。
对于这样的多项式，加法和减法实际上是等价的，均为异或操作。(同次项系数相同则二者抵消为0，不同则保留为1的项)

取模

模运算建立在等式 $x^{131}+x^{13}+x^2+x+1=0 mod f(x)$ 的基础上
由上述加法法则可以得到 $x^{131}=x^{13}+x^2+x+1 mod f(x)$
记$r(x)=x^{13}+x^2+x+1,\ m=131$

假设x^j是一个高次项，那么可由下面的过程将其降次(图中省略了mod f(x)) 

$x^j=x^{j-m}*x^m=x^{j-m}* r(x) \\ \because deg(r)<m, deg(x^j)=j\\ \therefore deg(x^{j-m}*r(x))<j$

这里原来是latex语法，由于csdn的问题，显示有问题，实际文本为:

x^j=x^{j-m}*x^m=x^{j-m}* r(x) \\ \because deg(r)<m, deg(x^j)=j\\ \therefore deg(x^{j-m}*r(x))<j

经过这样的运算后，$x^j$从高次项变成了低次项。

那么对于需要取模的多项式，只要对其的每一个高次项进行这种运算，再将结果求和，就能得到取模后的结果

乘法

对于要进行乘法运算的两个多项式$a(x)=\sum_{i=0}^{130}a_ix^i,\ b(x)=\sum_{i=0}^{130}b_ix^i$
由分配律可知$a(x)*b(x)=\sum_{i=0}^{130}a(x)*b_i*x^i$
对于多项式a(x),将其乘以一个$x^i,$等价于将其比特串左移i位。
故算法可表示为

input:a(x),b(x)
output:c(x)=a(x)*b(x) mod f(x)
c=0
for i in (0,m):
    if b[i]==1:
        c ^= temp << i;
return c mod f(x)

平方

求逆

由拓展欧几里得算法,对任意多项式a(x)和不可约多项式f(x),以下简写为a,f

具体迭代过程为

x1=1;y1=0;
x2=0;y2=1;
k1=a;k2=f;
while deg(k1)!=0:
    j=deg(k1)-deg(k2)
    if j <0:
        swap(k1,k2)
        swap(x1,x2)
        j=-j;
    k1=k1+(k2<<j)
    x1=x1+(x2<<j)
return x1;

实验结果

Convert 16336163 to bits:
00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000111110010100100100001011
16337163's inverse is:
10010010100000111011101011101101100010101100010000110100100110110001100001101110000111010101101110010100111010110010110111100011101
The product of "16337163"and its inverse is:
00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001

具体代码实现

程序基于c++的bitset实现。
p为不可约多项式，r为p去掉最高次的余式

加法

bitset<M> add(const bitset<M> &a, const bitset<M> &b)
{
    bitset<M> ans = a ^ b;
    return ans;
}

取模

bitset<M> mod(const bitset<2 * M> &a, const bitset<M> &r)
{
    bitset<2 * M> x, y;
    bitset<2 * M> extend_r(r.to_string());
    x = a;
    y ^= bitset<2 * M>(x.to_string().substr(M)); //取比特串的低M位
    while (judge(x))                            //judge函数功能:判断多项式的高M位是否存在1
    {
        y.reset();
        y ^= bitset<2 * M>(x.to_string().substr(M));
        for (int i = M; i < 2 * M; i++)
            if (x[i])
                y ^= extend_r << (i - M);
        x = y;
    }
    return bitset<M>(y.to_string().substr(M));
}

乘法

bitset<M> multiply(const bitset<M> &a, const bitset<M> &b, const bitset<M> &r)
{
    bitset<2 * M> x;
    bitset<2 * M> temp(a.to_string());
    for (int i = 0; i < M; i++)
    {
        if (b[i])
            x ^= temp << i;
    }
    return mod(x, r);
}

平方

bitset<M> square(const bitset<M> &a, const bitset<M> &r)
{
    bitset<2 * M> x;
    for (int i = 0; i < M; i++)
        x[i * 2] = a[i];
    return mod(x, r);
}

求逆

bitset<M> inverse(const bitset<M> &a, const bitset<M + 1> &p)
{
    bitset<2 * M> b, c, u, v, temp;
    bitset<M> r(p.to_string().substr(1));
    int j;
    b[0] = 1;
    u = bitset<2 * M>(a.to_string());
    v = bitset<2 * M>(p.to_string());
    while (degree(u))
    {
        j = degree(u) - degree(v);
        if (j < 0)
        {
            j = -j;
            temp = u;
            u = v;
            v = temp;
            temp = b;
            b = c;
            c = temp;
        }
        u = u ^ (v << j);
        b = b ^ (c << j);
    }
    return mod(b, r);
}

code.cpp

#include <iostream>
#include <string.h>
#include <bitset>
#include <vector>
#include <string>
using namespace std;
#define M 131
int degree(const bitset<2 * M> &a)
{
    for (int i = 2 * M - 1; i >= 0; i--)
        if (a[i] || i == 0)
            return i;
}
bool judge(const bitset<2 * M> &a)
{
    bitset<M> temp;
    bitset<2 * M> t(temp.set().to_string() + temp.to_string());
    t &= a;
    return t.any();
}
bitset<M> add(const bitset<M> &a, const bitset<M> &b)
{
    bitset<M> ans = a ^ b;
    return ans;
}
bitset<M> mod(const bitset<2 * M> &a, const bitset<M> &r)
{
    bitset<2 * M> x, y;
    bitset<2 * M> extend_r(r.to_string());
    x = a;
    y ^= bitset<2 * M>(x.to_string().substr(M));
    while (judge(x))
    {
        y.reset();
        y ^= bitset<2 * M>(x.to_string().substr(M));
        for (int i = M; i < 2 * M; i++)
            if (x[i])
                y ^= extend_r << (i - M);
        x = y;
    }
    return bitset<M>(y.to_string().substr(M));
}
bitset<M> multiply(const bitset<M> &a, const bitset<M> &b, const bitset<M> &r)
{
    bitset<2 * M> x;
    bitset<2 * M> temp(a.to_string());
    for (int i = 0; i < M; i++)
    {
        if (b[i])
            x ^= temp << i;
    }
    return mod(x, r);
}
bitset<M> square(const bitset<M> &a, const bitset<M> &r)
{
    bitset<2 * M> x;
    for (int i = 0; i < M; i++)
        x[i * 2] = a[i];
    return mod(x, r);
}
bitset<M> inverse(const bitset<M> &a, const bitset<M + 1> &p)
{
    bitset<2 * M> b, c, u, v, temp;
    bitset<M> r(p.to_string().substr(1));
    int j;
    b[0] = 1;
    u = bitset<2 * M>(a.to_string());
    v = bitset<2 * M>(p.to_string());
    while (degree(u))
    {
        j = degree(u) - degree(v);
        if (j < 0)
        {
            j = -j;
            temp = u;
            u = v;
            v = temp;
            temp = b;
            b = c;
            c = temp;
        }
        u = u ^ (v << j);
        b = b ^ (c << j);
    }
    return mod(b, r);
}
int main()
{
    bitset<M + 1> p("111");
    p[131] = 1;
    p[13] = 1;
    bitset<M> r("111");
    r[13] = 1;
    bitset<M> id(16337163);
    cout << "Convert 16336163 to bits:" << endl
         << id << endl
         << "16337163's inverse is:" << endl
         << inverse(id, p) << endl
         << "The product of \"16337163\"and its inverse is:" << endl
         << multiply(inverse(id, p), id, r) << endl;
}