Gabidulin Codes

[arpitdm]

A Linear Gabidulin Code Gab[n,k] over F_{q^m} of length n <= m and dimension k <= n is the set of all words, formed by the operator evaluation of a q-degree restricted skew polynomial (with frobenius endomorphism and a finite field) f(x) \in S_{F_{q^m}}[x].

i.e. Gab[n,k] = {(f(g_0), f(g_1),..., f(g_{n-1})) = f(g): deg(f) < k}

where g_0,...,g_{n-1} are fixed elements belonging to F_{q^m} and are linearly independent over F_{q}

This ticket proposes a new class for Gabidulin Codes along with encoders and decoders for it.

[arpitdm]

A few points:

1. I've not added the tests and documentation for now since the code design is still in flux.
2. The file compiles successfully and one can create the code and codewords in vector form, i.e. F_{q^m}ⁿ.
3. One problem right now is that, although the dual code and the parity check matrix that are constructed are correct (verified by trying G*H^T = 0, as well as checking linear independence of h_0, h_1,..., h_{n-1}), I am using a random element from the right kernel of the coefficient matrix (Equation 2.28 of Antonia Wachter-Zeh's PhD thesis). And so dual of the dual code does not give me the original code. solve_right returns the trivial (all-zeros) solution.
4. Should there be a def generator_matrix method within the code class itself?

---

New commits:

​ddaef41

Initial construction of the Gabidulin Code class. PolynomialEvaluationEncoder and GeneratorMatrixEncoder also included.

[dlucas]

Hello,

To answer your question (4.), for now, yes, I guess you need to have such a method. Later on, when the abstract class to manage rank-metric codes will be implemented we can replicate the design used in AbstractLinearCode, but for now, copy-pasting generator_matrix from AbstractLinearCode and put it in GabidulinCode is fine.

Some remarks:

1. You don't have to use backslashes to jump lines when you're in a parenthesis/square bracket block. For instance, you can write:

   % (self.length(), self.dimension(), \
   self.minimum_distance(), self.base_field())
   

   as

   % (self.length(), self.dimension(),
   self.minimum_distance(), self.base_field())
   

2. Lines 19-34 are more or less copy-pasted from relative_finite_field_extension.py, and are not necessary. Take lines 19-20 for instance: it's a carbon copy of lines 116-117 of relative_finite_field_extension.py (except you replaced must be by has to be in the error message)... Which basically means you run the same sanity check twice: first time on lines 19-20 and second time line 49 when calling the constructing a RelativeFiniteFieldExtension. I think you should just remove all those tests and let RelativeFiniteFieldExtension do its job on sanitizing.

3. In the same vein: you copy-pasted a lot of getter from relative_finite_field_extension.py, and I'm not sure it makes sense. The user might be a bit puzzled by absolute_field_power which returns something of little relevance wrt. the Gabidulin code. Furthermore, if the user desperately needs to access this value, they can still access it by the RelativeFiniteFieldExtension returned by relative_finite_field_extension (this one we should keep).

4. There's some kind of naming conflict: evaluation_points and linearly_independent_elements are the same (see line 68), but while the user has to write linearly_independent_elements=stuff at construction time, they have to write my_code.evaluation_points() to access them from the constructed code. I think you should choose one name and stick to it. I vote for evaluation_points for consistency with grs.py.

5. Remember to not call methods in for loops except it's absolutely necessary. You do this on lines 96-97, 180-184 and 212. E.g. line 182, as return value of p.coefficients() never changes throughout the iterations, I would store it in a local variable and use this local variable instead.

Otherwise it seems good.

Best,

David

[jsrn]

Noticed the following: In one place, you should be computing sigma^t(g) for some power t and some element g, where sigma is the Frobenius automorphism. You implement this as pow(sigma(g), t), but that's not what it means. Rather, it means sigma( sigma( ... sigma(g)...), where there is t sigmas. In another place, you compute simply g^(q^t) instead of calling sigma, which is also bad style; it should again be sigma called t times on g.

[jsrn]

Also, avoid "log" if you can:

• s = log(q,p) should be s = relative_field.degree(). Similar with sm.
• m should be sm / m. Check if the division goes well.

[git]

Branch pushed to git repo; I updated commit sha1. Last 10 new commits:

​de61ca4	Importing ring_element, integral_domain and principal_ideal_domain were throwing deprecation warnings. Updated the import statements according to Tickets 19167 and 20011.
​9d21b54	module that implements class 'Side' with Left and Right as the two instances.
​9de7bc1	Merge branch 'u/arpitdm/skew_polynomials' of git://trac.sagemath.org/sage into skew_polynomials
​9c526ad	Fixed doctests caused by deprecation warnings, inability to access private variables, AttributeError in accessing of parent, and .
​a06c2bf	Merge branch 'u/arpitdm/skew_polynomials' of git://trac.sagemath.org/sage into skew_polynomials
​18c7982	Fixed bug with integer coercion
​c6183bc	Merge branch 'develop' into temp3
​dd5c575	Editing declarations of cython functions so as to make them compatible.
​60b7869	Merge branch 'u/arpitdm/skew_polynomials' of git://trac.sagemath.org/sage into gabidulin_codes
​40a3d37	Refactored code based on comments 3, 4 and 5. Added support for vector_to_matrix and matrix_to_vector representations of codewords. Added methods for rank_weight of a codeword and rank_distance between two codewords.

[arpitdm]

Oh! I thought I only added the gabidulin.py file to my commit. Is there a way I can remove the skew polynomial files from this ticket?

Also, @jsrn, you mention I computed g^qt in the generator matrix of the code. According to Equation 2.27 of the PhD thesis, there is no sigma in the formation of the elements of the generator matrix.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​799d871

Added support for in the PolynomialEvaluationEncoder. This includes methods for fast lagrange-type interpolation of skew polynomials, minimum subspace polynomials and multi point evaluation of skew polynomials.

[arpitdm]

unencode_nocheck method in the Polynomial Evaluation Encoder does not work right now because it calls the right division of one skew polynomial by another which fails because of the integer coercion error in #13215. Apart from that, as far as I am able to check, the four new methods are implemented correctly and they do not cause any errors during build or runtime.

[jsrn]

Replying to arpitdm:

Oh! I thought I only added the gabidulin.py file to my commit. Is there a way I can remove the skew polynomial files from this ticket?

I guess you merged in the current tip of #13215 and that's what is showing up as a long list of commits? That's perfectly ok (#13215 is a dependency of this ticket).

Also, @jsrn, you mention I computed g^qt in the generator matrix of the code. According to Equation 2.27 of the PhD thesis, there is no sigma in the formation of the elements of the generator matrix.

Well, due to sigma always being the Frobenius, then sigma^t(g) is, by definition, g^(q^t). Wachter-Zeh just chose to write things using explicit q-powering notation everywhere, instead of writing sigma. But it is more mathematically elegant to write sigma, and it fits better with the rest of the code. (say for instance that we choose to make a parameter that you can construct Gabidulin codes with sigma be a power of Frobenius, e.g. sigma(g) = g^(q^2). These are sometimes, somewhat ridiculously, called "Generalized Gabidulin codes").

[git]

Branch pushed to git repo; I updated commit sha1. Last 10 new commits:

​9064e0e	Fix division bug in interpolation
​6618340	Fix divide-by-two potential problem in interpolation
​0635ec9	Merge branch '13215_skew_polynomials' into 21131_interpolation_skew_poly
​e67411f	merge mvp_mpe ticket
​392047b	Merge branch 'u/arpitdm/gabidulin_codes' of git://trac.sagemath.org/sage into gabidulin_codes
​d97d9c2	removed commented code
​4e3e143	updated encode function
​16d2d3d	fixed generator matrix method
​3c1261f	changed code constructor method to add linearized polynomial ring as mandatory argument
​e1bf28f	added Gao decoder

[jsrn]

Hi Arpit,

Just saw your Gao decoder, and I have a few initial comments:

• r_ (perhaps R is a better name) is an interpolation polynomial: construct it using your new interpolation_polynomial method instead of what you're manually doing right now (which is much slower).

• After you have performed the division and checked that rem is zero (use if rem.is_zero()), then you need to check that quo indeed corresponds to the message of a codeword close to r. See how the Gao decoder for Reed-Solomon codes does this with a helper function, covering both decode_to_message and decode_to_codeword in an efficient way.

• You _partial_xgcd would be more readable in a compact form, IMO. Say you call r_current -> r_c and r_previous -> r_p, then you can replace the loop body with:

    (quo, r_c), r_p = r_c.right_quo_rem(r_p), r_c
    u_c, u_p = u_p - q*u_c, u_c
    v_c, v_p = v_p - q*v_c, v_c

Best, Johan

[jsrn]

Another comment: The Gabidulin constructor shouldn't take the skew polynomial ring (NOT linearized polynomial ring) as an argument, just like Reed-Solomon and Reed-Muller codes don't take it a polynomial ring as argument. Indeed: the code itself (being a sub-vector space of GF(q^m)^n) is invariant of the skew polynomial ring used for encoding the code from a message space. By this logic, the polynomial encoder objects for Reed-Solomon and Reed-Muller codes can optionally take a polynomial ring, but if the user doesn't supply one, one will be created automatically.

[arpitdm]

Replying to jsrn:

Another comment: The Gabidulin constructor shouldn't take the skew polynomial ring (NOT linearized polynomial ring) as an argument, just like Reed-Solomon and Reed-Muller codes don't take it a polynomial ring as argument. Indeed: the code itself (being a sub-vector space of GF(q^m)^n) is invariant of the skew polynomial ring used for encoding the code from a message space. By this logic, the polynomial encoder objects for Reed-Solomon and Reed-Muller codes can optionally take a polynomial ring, but if the user doesn't supply one, one will be created automatically.

The reason I put it there is because although the code is a subspace, the skew polynomial ring takes an additional argument, the base ring automorphism. This is not required in a GRS code. It could be Frobenius or a power of the Frobenius, right? And so, a given polynomial could belong to different rings. In the earlier commits I had assumed that the automorphism is the Frobenius.

[jsrn]

The reason I put it there is because although the code is a subspace, the skew polynomial ring takes an additional argument, the base ring automorphism. This is not required in a GRS code. It could be Frobenius or a power of the Frobenius, right? And so, a given polynomial could belong to different rings. In the earlier commits I had assumed that the automorphism is the Frobenius.

True. The Gabidulin code constructor could therefore take as argument /which/ power of the Frobenius (integer between 1 and m-1, both inclusive). That feels more natural for the user, and is much less work to specify.

Best, Johan

[caruso]

Actually, thanks to the recent PhD thesis of Robert, Gabidulin codes make sense (and seem to have applications) for extensions of infinite fields (e.g. number fields) as well. Moreover I strongly believe (although I have not checked it carefully) and Wechter-Zeh's algorithms extend readily to the general case.

For this reason, I would say that it makes sense to allow more general morphisms that powers of Frobenius. Concretely I propose to replace the current linearized_polynomial_ring by twisting_homomorphism (and possibly make relative_field optional and let the default be the fixed ring under the twisting homomorphism).

[jsrn]

Replying to caruso:

Actually, thanks to the recent PhD thesis of Robert, Gabidulin codes make sense (and seem to have applications) for extensions of infinite fields (e.g. number fields) as well.

That's true and a good point. However, codes over infinite fields needs to be handled specially in many ways, and it would require special casing and some additional testing to do make sure things work. Also, our current plan is to let GabidulinCode inherit from AbstractLinearCode, and that class has many things which I know break over infinite fields.

Moreover I strongly believe (although I have not checked it carefully) and Wechter-Zeh's algorithms extend readily to the general case.

I'm also pretty sure it "just works". But it will have terrible performance in the Euclidean algorithm because of unmitigated coefficient growth.

For this reason, I would say that it makes sense to allow more general morphisms that powers of Frobenius. Concretely I propose to replace the current linearized_polynomial_ring by twisting_homomorphism (and possibly make relative_field optional and let the default be the fixed ring under the twisting homomorphism).

For the sake of getting this GSoC finished in time, I would down-vote officially supporting infinite fields. But of course, we can still think about the best interface that will scale well in the future. I think your suggestion wrt. these arguments is good if it is practical in the current Sage. For instance, how would I create a Frobenius automorphism from a relative field to an extension field, and such that the Gabidulin code constructor can query the fixed field? The following works poorly:

sage: q = 3^2
sage: F.<a> = GF(q^2, 'a')
sage: sigma = F.hom([a^q])
sage: sigma
Ring endomorphism of Finite Field in a of size 3^4
  Defn: a |--> a^3 + a^2 + 2*a
sage: sigma.fixed_<tab gives nothing>

Best, Johan

[caruso]

OK for not supporting infinite fields for now. (And I agree that there are issues with growing of coefficients. By the way, do you know if there exists an analogue of subresultants for skew polynomials?)

To answer your last question:

sage: q = 5^2
sage: k.<a> = GF(q^2)
sage: Frob = k.frobenius_endomorphism(2)
sage: Frob.fixed_field()
(Finite Field in a_fixed of size 5^2, Ring morphism:
   From: Finite Field in a_fixed of size 5^2
   To:   Finite Field in a of size 5^4
   Defn: a_fixed |--> 4*a^3 + 4*a^2 + 4*a + 3)

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​de37ab8

fixed Gao decoder class

[jsrn]

@arpitdm: Remember to add a message after you push. Your changes to the decoder looks good. I'm somewhat baffled why you changed GabidulinCode to inherit from AbstractLinearCode since we have been discussing for some time now that you should create a class AbstractLinearRankMetricCode.

@caruso: Thanks for the snippet, I hadn't noticed that but of course it's pretty obvious. I'm not sure, though, that I would find this behaviour the most intuitive for Gabidulin codes, especially if relative_field was automatically determined from the fixed field of twisting_homomorphism. Fir instance, an unsuspecting user might well do something like:

sage: F = GF(5^20,'a')
sage: evals = [ a^i for i in range(20) ]
sage: C = GabidulinCode(evals, k=5) # twisting_homomorphism taken to be absolute Frob.
sage: <fun experiments with C.>

sage: C2 = GabidulinCode(evals, k=5, twisting_homomorphism=F.frobenius_endomorphism(2))
<BOOM: Some error probably stating that the evaluation points are not linearly independent over the subfield>

A carefully phrased exception at the end would help the user find the bug in an interactive session, but I think this is the kind of surprising behaviour that one keeps on running into every time one uses the functionality.

I would prefer a session looking something like this:

sage: F = GF(5^20,'a')
sage: evals = [ a^i for i in range(20) ]
# unless specified, subfield = prime field, frobenius power = 1
sage: C = GabidulinCode(evals, k=5)
sage: C
[20,5,16] Gabidulin code over Finite Field in a of size 5^20 with subfield of size 5

sage: C2 = GabidulinCode(evals, k=5, frobenius_power=2) #subfield is still prime field
sage: C2
[20,5,16] Gabidulin code with Frobenius power 2 over Finite Field in a of size 5^20 with subfield of size 5

sage: C3 = GabidulinCode([a^i for i in range(10)], k=5, frobenius_power=2, subfield=GF(5^2,'b'))
sage: C3
[10,5,6] Gabidulin code over Finite Field in a of size 5^20 with subfield of size 5^2.

For your other question: I'm not aware of subresultants for skew polynomials, no. I've recently started looking at computations in skew polynomials in the more general setting of derivations, and over infinite fields. All algorithms I know of perform pretty terribly compared to the case without coefficient growth: the usual approaches of homomorphic imaging doesn't seem to work. Such rings seem more widely known as Ore polynomials. We should add a note about this in the doc of #13215.

[arpitdm]

Replying to jsrn:

@arpitdm: Remember to add a message after you push. Your changes to the decoder looks good. I'm somewhat baffled why you changed GabidulinCode to inherit from AbstractLinearCode since we have been discussing for some time now that you should create a class AbstractLinearRankMetricCode.

I was about to when discussion started on #13215. Anyway, in the _decode_to_code_and_message, I needed to use the connected_encoder method which is available in the AbstractLinearCode? class. Since ARMC is not yet opened, I inherited this temporarily from that. I will update to ARMC once it is possible.

[jsrn]

Replying to arpitdm:

I was about to when discussion started on #13215. Anyway, in the _decode_to_code_and_message, I needed to use the connected_encoder method which is available in the AbstractLinearCode? class. Since ARMC is not yet opened, I inherited this temporarily from that. I will update to ARMC once it is possible.

It's a method on Decoder - not AbstractLinearCode. But I can see how it could be practical to temporarily sidestep ARMC until it's there.

[git]

Branch pushed to git repo; I updated commit sha1. Last 10 new commits:

​b9f8239	merging latest changes
​34bebbc	created top level functions for interpolation and fraction field skew polynomial ring. added documentation and indirect doctests.
​b523806	changed names of private methods
​e435f56	Merge branch 'u/arpitdm/mvp_mpe' of git://trac.sagemath.org/sage into mvp_mpe
​d403fae	fixes to documentation
​7e4b3a3	added a very basic constructor, a couple of basic getter methods and deleted a couple of methods from ALC class. added some rank distance, rank weight, to_matrix_representation and from_matrix_representation.
​0657293	Merge branch 'u/arpitdm/abstract_linear_rank_metric_code_class' of git://trac.sagemath.org/sage into abstract_rank_metric_code
​08aef71	merging with updates from other tickets
​da374e9	small fix to docs
​7278453	added random_element method. refactored some code. added documentation and tests.

[arpitdm]

I've added refactored and scrubbed the code to make it cleaner where ever possible. I've created proper inheritances for the classes, and added some methods. And I've added documentation and tests for nearly every method and class. There are two issues that remain. One is def parity_evaluation_points where the problem is to find a concrete, non-trivial solution for the equation Ax = 0 where A is the coefficient matrix. I am not sure how to do that in Sage.

And the second is, finding the right name for the getter method on line 321 (currently it is left as def m).

[jsrn]

Replying to arpitdm:

I've added refactored and scrubbed the code to make it cleaner where ever possible. I've created proper inheritances for the classes, and added some methods. And I've added documentation and tests for nearly every method and class.

I've looked through the code in overview, and it looks good! Some comments on the structure:

• Constructor of Gabidulin code: you require both subfield and twist map. I've thought a bit more about it, and I think that the subfield should always be the fixed field of the twist map: if it is smaller than the fixed field, then the code is not MRD, and if it is larger, then the code is not linear. Therefore, I suggest that both twist map and subfield is optional: if twist map is set and subfield is not, then subfield is the fixed field of the twist map; throw an exception if the twist map does not have the fixed_field method. If subfield is given but twist map is not, then default to Frobenius of that field extension. If both twist map and subfield is set, then throw an exception if the twist map has a method fixed_field and it does not return the given subfield.

• You don't have an example using field_extension.

• Your parity_check_matrix is strange. Why not simply return self.dual_code().generator_matrix()?

• What's the point of _vector_space and def vector_space?

• You shouldn't override generator_matrix and random_element.

• What's the point of GabidulinCode.message_space?

• Doc: There's no such thing a as a q-degree of a skew polynomial. The doc should mention that over finite fields, the skew polynomials are very similar to linearized polynomials (isomorphic when the twist map is power-by-q).

There are two issues that remain. One is def parity_evaluation_points where the problem is to find a concrete, non-trivial solution for the equation Ax = 0 where A is the coefficient matrix. I am not sure how to do that in Sage.

Just proceed exactly as you have done, until you get solution_space. Then a non-zero element is solution_space.basis()[0], assuming that solution_space is not empty.

And the second is, finding the right name for the getter method on line 321 (currently it is left as def m).

That should be a method on AbstractLinearRankMetricCode. It could be called field_extension_degree.

Best, Johan

[arpitdm]

Replying to jsrn:

• Constructor of Gabidulin code: you require both subfield and twist map. I've thought a bit more about it, and I think that the subfield should always be the fixed field of the twist map: if it is smaller than the fixed field, then the code is not MRD, and if it is larger, then the code is not linear. Therefore, I suggest that both twist map and subfield is optional: if twist map is set and subfield is not, then subfield is the fixed field of the twist map; throw an exception if the twist map does not have the fixed_field method. If subfield is given but twist map is not, then default to Frobenius of that field extension. If both twist map and subfield is set, then throw an exception if the twist map has a method fixed_field and it does not return the given subfield.

What do I do if both subfield and twist_map are not given?

Also, should I remove the dependency on #21088?

[jsrn]

Replying to arpitdm:

What do I do if both subfield and twist_map are not given?

Then take twist_map to be absolute frobenius, i.e. frobenius wrt. the prime field, and subfield to be the prime field.

Also, should I remove the dependency on #21088?

Oh, yeah I guess you can, since finite fields are already formally supported without #21088.