1

I have a Galois Field GF(2^409) and irreducible polynomial f(x) = x^409 + x^15 + x^6 + x + 1 of which the coefficients can be only 1 or 0

If I have some element of this field a(x), how can I find reverse element a_1(x), such that 

a(x) * a_1(x) = 1 (mod f(x))

using the Extended Euklid Algorithm

  • I`m using polynomial basis

  • I know that I can find a_1(x) just using pow(a(x),2^409 - 2), but I`m using Python,and power operation takes too much time

  • I have troubles with defining the division operation under polynomials (without it I can`t use Extended Euklid Algorithm)

Patrick Haugh
  • 59,226
  • 13
  • 88
  • 96
Romanov Roman
  • 21
  • 4
  • Have you ever used Sage? It's a very powerful tool for these kind of problems. http://doc.sagemath.org/html/en/tutorial/index.html – Patrick Haugh Dec 13 '16 at 19:43

2 Answers2

2

One common approach to represent elements of GF(2k) would be using a big integer (long in Python 2, int in Python 3) and let each bit represent one coefficient. To perform polynomial division, you can essentially follow the steps from long division as you'd perform them with pen and paper.

Let's take some simpler example. Take GF(27) with modulus x7 + x + 1. Let's try to find the inverse of x6 + x3 + x. So you have to compute the GCD of these two polynomials using the extended Euclidean algorithm. And the first step is computing the quotient and remainder of one by the other. Namely the one with larger degree divided by the one with smaller.

x6 + x3 + x is 1x6 + 0x5 + 0x4 + 1x3 + 0x2 + 1x1 + 0x0 or 1 0 0 1 0 1 0 in short. The modulus is 1 0 0 0 0 0 1 1 by the same convention. So you divide

1 0 0 0 0 0 1 1 / 1 0 0 1 0 1 0 = 1 0 remainder 1 0 1 1 1
1 0 0 1 0 1 0
      1 0 1 1 1

So what did I do here? Dividing a by b I look for the highest set bit in each of these, i.e. for the degree of the polynomial. The difference between these degrees is a bit I want to set in the quotient. In this case a has degree 7 and b has degree 6. The difference is 1 so the quotient has to have a term x1. Now I write down b multiplied by that term, which in this case is simply the binary representation of b shifted left by the difference in degrees). Then I subtract that shifted value from the original a, but I do the subtraction in 2[x] so it's actually an XOR operation. The result of that is a new number, which I use as the value of a in the next iteration. I continue until the difference in degree becomes negative, i.e. I'm dividing a polynomial of smaller degree by one of larger. Then I'm done, and the last value of a is my remainder. In the division above, I was done after one step.

With this operation of performing a division with remainder, you should be able to figure out the Euclidean algorithm. The extended version will require a bit more work, but give it a try. In the end, you should be able to find that in the example given above, the inverse is x3 + x + 1 (computed by Sage). For comparison, the inverse of x6 + x3 + x in your original field GF(2409) would be the result of

"+".join("x^{}".format(i) for i in range(409,-1,-1) if (1 << i) & 71418312917235914488287388787154746126088900129923309868417397199063993100653429184865046255190862140902867842544110449930)

This number I computed in Sage:

x = polygen(ZZ)
m = x^ 409 + x^15 + x^6 + x + 1
F409 = GF(2^409, name="z", modulus=p)
z = F409.gen()
(1/(z^6+z^3+z)).polynomial().change_ring(ZZ)(2)
MvG
  • 57,380
  • 22
  • 148
  • 276
0

You can use SymPy and this awesome answer to work on finite fields.

You are on the ring Z/nZ with n = 2 so you want to use the GF class from the linked answer :

>>> %time field = GF(2, 409)
CPU times: user 1min 47s, sys: 6.64 ms, total: 1min 47s
Wall time: 1min 47s

As you can see the initialization of GF(2, 409) is quite long (~= 1-2 min)...
...but the computation of multiplication inverts is almost instant :

%%time

import numpy as np

# list of 410 coefficients regarding x^409 (left-most) up to x^0 (right-most)
coeffs = [0] * 410 
for idx in [0, 1, 6, 15, 409]: # exponents present in your polynome
    coeffs[409 - idx] = 1 # set to 1 to indicate existence

li = []
for i in range(4):
    '''
    To help SymPy make the polynome `x^409 + x^15 + x^6 + x + 1`,
    provide the `GF.inv()` method with its coefficients using `ZZ.map()`
    Actually, `inv()` will call `sympy.polys.galoistools.gf_gcdex()`
    '''
    coeffs = field.inv(ZZ.map(coeffs))
    li.append(np.array(coeffs))

CPU times: user 37.3 ms, sys: 2 µs, total: 37.3 ms
Wall time: 36 ms

The result returned by GF.inv() is the array of the invert's coefficients. Note that reapplying inversion again reverts to original polynome :

>>> print(np.equal(li[0], li[2]).all() and np.equal(li[1], li[3]).all())
True

You can also verify that p(x) * p(x)^(-1) = 1 :

>>> field.mul(li[0], li[1])
[1] # Ok !

Below is the stripped code from cited answer with only what you need for inversion :

import itertools
from sympy.polys.domains import ZZ
from sympy.polys.galoistools import (gf_irreducible_p, gf_gcdex)
from sympy.ntheory.primetest import isprime

class GF():
    def __init__(self, p, n=1):
        p, n = int(p), int(n)
        if not isprime(p):
            raise ValueError("p must be a prime number, not %s" % p)
        if n <= 0:
            raise ValueError("n must be a positive integer, not %s" % n)
        self.p = p
        self.n = n
        if n == 1:
            self.reducing = [1, 0]
        else:
            for c in itertools.product(range(p), repeat=n):
              poly = (1, *c)
              if gf_irreducible_p(poly, p, ZZ):
                  self.reducing = poly
                  break

    def inv(self, x):
        s, t, h = gf_gcdex(x, self.reducing, self.p, ZZ)
        return s


GF(2, 409).inv(ZZ.map(coeffs))
dilavasso
  • 11
  • 3