Quelle  body.autodoc   Sprache: unbekannt

################# Outline of the documentation for AutoDoc #########################

@Chapter Introduction

@Chapter Examples

A few simple examples illustrating the use of the
package.  For more information see Chapter <Ref Chap="Chapter_AllFunctions"/>

@Section The 5-qubit code
@Section Hyperbolic codes from a file
@Section Randomly generated cyclic codes


@Chapter Description 
@ChapterTitle Description of the algorithm 

@Chapter AllFunctions
@ChapterTitle All Functions 

@Section DistanceFunctions 
@SectionLabel DistanceFunctions 
@SectionTitle Functions for computing the distance 

@Section IOFunctions 
@SectionLabel IOFunctions 
@SectionTitle Input/Output Functions 

@Section HelperFunctions
@SectionLabel HelperFunctions
@SectionTitle Helper Functions

@Chapter FileFormat
@ChapterTitle Extended MTX (MTXE) File Format
@ChapterLabel FileFormat

@Section General information

@Section Example MTXE files


@Chapter Introduction

The &GAP; package &QDistRnd; implements a
probabilistic algorithm for finding the distance of a $q$-ary quantum
low-density parity-check code linear over a finite field $F=\mathop{\rm GF}(q)$.
While there is no guarantee of the performance of the algorithm (the
existing bounds in the case of quantum LDPC codes are weak, see <Ref
Subsect="Subsection_TheAlgorithmDetails"/>), an empirical convergence
criterion is given to estimate the probability that a minimum weight
codeword has been found.  Versions for CSS and regular stabilizer
codes are given, see Section <Ref Sect="Section_DistanceFunctions"/>

In addition, a format for storing matrices associated with $q$-ary
quantum codes is introduced and implemented, see
Chapter <Ref Chap="Chapter_FileFormat"/> and
Sec. <Ref Sect="Section_IOFunctions"/>.  The
format is based on the well establised MaTrix market eXchange (MTX)
Coordinate format developed at NIST, and is designed for full backward
compatibility with this format.  Thus, the files are readable by any
software package which supports MTX.

The routines in the package are derived from the code originally written
by one of the authors (LPP).  A related Covering Set algorithm has a
provable performance for generic (non-LDPC) quantum codes based on
random matrices <Cite Key="Dumer-Kovalev-Pryadko-IEEE-2017"/>.
Implemented version is a variant of the random **information set**
(IS) algorithm based on random column permutations and Gauss'
elimination <Cite Key="Leon-1988"/> <Cite Key="Kruk-1989"/> <Cite
Key="Coffey-Goodman-1990"/>.

The &GAP; computer algebra system was chosen because of its excellent
support for linear algebra over finite fields.  Here we give a
reference implementation of the algorithm, with a focus on matrix
formats and generality, as opposed to performance.  Nevertheless, the
routines are sufficiently fast when dealing with codes of practically
important block lengths $n\lesssim 10^3$.

@Chapter Examples

@Chapter Description

@Section Elementary version
@SectionLabel SimpleVersion

@Subsection What it does?
@SubsectionLabel TheGoal


In the simplest possible terms, we are given a pair of matrices $P$
and $Q$ with orthogonal rows, $PQ^T=0$.  The matrices have entries in
a finite field $F=\mathop{\rm GF}(q)$, where $q$ is a power of a prime.  The goal
is to find the smallest weight of a non-zero vector $c$ over the same
field $F$, such that $c$ be orthogonal with the rows of $P$, $Pc^T=0$,
and linearly independent from the rows of $Q$.

@Subsection The algorithm
@SubsectionLabel TheProcedure

We first construct a generator matrix $G$ whose rows form
a basis of the $F$-linear space of all vectors orthogonal to the rows
of $P$.  At each step, a random permutation $S$ is generated and
applied to the columns of $G$.  Then, Gauss' elimination with back
substitution  renders the resulting matrix to the reduced
row echelon form, after which the inverse permutation $S^{-1}$ is
applied to the columns.  Rows of the resulting matrix $G_S$ that are
linearly independent from the rows of $Q$ are considered as candidates
for the minimum weight vectors.  Thus, after $N$ steps, we are getting
an upper bound on the distance which is improving with increasing $N$.

@Subsection Intuition
@SubsectionLabel Intuition


The intuition is that each row of $G_S$ is guaranteed to contain at
least `rank`$(G_S)-1$ zeros.  Thus, we are sampling mostly
lower-weight vectors from the linear space orthogonal to the rows of
$P$.  Further, it is easy to see that any vector obtained this way is
**irreducible** <Cite Key="Dumer-Kovalev-Pryadko-bnd-2015"/>, i.e., it
cannot be decomposed into a pair of zero-syndrome vectors with
non-overlapping supports.

Furthermore, the eventual convergence is guaranteed.  Indeed, if $c$
is a minimum-weight codeword of weight $d$, consider a permutation $S$
which places one position from its support into the 1st column, and
the remaining positions into the last $d-1$ columns.  Vector $c$ being
the lowest-weight non-trivial vector, no pivot column may be in the
block of last $d-1$ columns. This guarantees that vector $c$ is
obtained as the first row of $G_S$.
(This argument is adapted to degenerate quantum codes from
<Cite Key="Cuellar-etal-2020"/>).

@Subsection CSS version of the algorithm
@SubsectionLabel AlgorithmCSS

The described version of the algorithm is implemented in the function `DistRandCSS`
(<Ref Sect="Section_DistanceFunctions"/>).  It applies to the case of
Calderbank-Shor-Steane (CSS) codes, where the matrices $P=H_X$ and
$Q=H_Z$ are called the CSS generator matrices, and the computed
minimum weight is the distance $d_Z$ of the code.  The number of
columns $n$ is the block length of the code, and it encodes $k$
qudits, where $k=n-$`rank`$(H_X)-$`rank`$(H_Z)$.  To completely
characterize the code, we also need the distance $d_X$ which can be
obtained by calling the same function with the two matrices interchanged.
The conventional code distance $d$ is the minimum of $d_X$ and $d_Z$.
Parameters of such a $q$-ary CSS code are commonly denoted as
$[[n,k,(d_X,d_Z)]]_q$, or simply $[[n,k,d]]_q$ as for a general
$q$-ary stabilizer code.

@Subsection Generic version of the algorithm
@SubsectionLabel AlgorithmGeneric

CSS codes are a subclass of general $F$-linear stabilizer codes which
are specified by a single stabilizer generator matrix $H=(A|B)$
written in terms of two blocks of $n$ columns each.  The orthogonality
condition is given in a symplectic form, $$ A B^T-B A^T=0, $$ or,
equivalently, as orthogonality between the rows of $H$ and the
symplectic-dual matrix $\tilde H=(B|-A)$.  Non-trivial vectors in the code must be
orthogonal to the rows of $P=\tilde H$ and linearly independent from
the rows of $Q=H$.  The difference with the CSS version of the
algorithm is that we must minimize the **symplectic** weight of
$c=(a|b)$, given by the number of positions $i$, $1\le i\le n$, such
that either $a_i$ or $b_i$ (or both) be non-zero.

The parameters of such a code are denoted as $[[n,k,d]]_q$, where
$k=n-$`rank`$H$ is the number of encoded qudits, and $d$ is the
minimal symplectic weight of a non-trivial vector in the code.  It is
easy to check that a CSS code can also be represented in terms of a
single stabilizer generator matrix.  Namely, for a CSS code with
generators $H_X$ and $H_Z$, the stabilizer generator matrix has a
block-diagonal form, $H=$`diag`$(H_X,H_Z)$.

A version of the algorithm for general $F$-linear stabilizer codes is
implemented in the function `DistRandStab` (<Ref Sect="Section_DistanceFunctions"/>).

**Important Notice**: In general, here one could use most general
permutations of $2n$ columns, or restricted permutations of $n$
two-column blocks preserving the pair structure of the matrix.  While
the latter method would be much faster, there is no guarantee that
every vector would be found.  As a result, we decided to use general
permutations of $2n$ columns.

@Section Some more details 
@Subsection Quantum stabilizer codes

Representation of quantum codes in terms of linear spaces
is just a convenient map.  In the case $q=2$ (qubits), the details can
be found, e.g., in the book of Nielsen and Chuang, <Cite
Key="Nielsen-book"/>.  Further details on the theory of stabilizer
quantum error correcting codes based on qubits can be found in the
Caltech Ph.D. thesis of Daniel Gottesman <Cite Key="gottesman-thesis"/>
and in the definitive 1997 paper by Calderbank, Rains, Shor, and
Sloane <Cite Key="Calderbank-1997"/>.  Theory of stabilizer quantum
 codes based on qudits ($q$-state quantum systems)
was developed by Ashikhmin and Knill
<Cite Key="Ashikhmin-Knill-2001"/>
(prime fields with $q$ prime) and by Ketkar, Klappenecker, Kumar, & Sarvepalli <Cite
Key="Ketkar-Klappenecker-Kumar-Sarvepalli-2006"/> (extension fields
with $q$ a non-trivial power of a prime).

In the binary case (more generally, when $q$ is a prime), $F$-linear
codes coincide with **additive** codes.  The **linear** codes [e.g.,
over $\mathop{\rm GF}(4)$ in the binary case <Cite Key="Calderbank-1997"/>] is a
different construction which assumes an additional symmetry.  A brief
summary of $F$-linear quantum codes [where $F=\mathop{\rm GF}(q)$ with $q=p^m$,
$m>1$ a non-trivial power of a prime] can be found in the introduction
of Ref. <Cite Key="Zeng-Pryadko-hprod-2020"/>.  The construction is
equivalent to a more physical approach in terms of a lifted Pauli
group suggested by Gottesman <Cite
Key="Gottesman-prime-power-2014"/>.

@Subsection The algorithm
@SubsectionLabel TheAlgorithmDetails 

@LatexOnly \noindent\underline{%
Case of **classical linear codes** 
@LatexOnly }

The algorithm <Ref Subsect="Subsection_TheProcedure"/> is closely
related to the algorithm for finding minimum-weight codewords in a
classical linear code as presented by Leon <Cite Key="Leon-1988"/>,
and a related family of **information set** (IS) decoding algorithms <Cite
Key="Kruk-1989"/> <Cite Key="Coffey-Goodman-1990"/>.

Consider a classical linear $q$-ary code $[n,k,d]_q$ encoding $k$
symbols into $n$, specified by a generator matrix $G$ of rank $k$.
Using Gauss' algorithm and column permutations, the generator matrix
can be rendered into a **systematic form**, $G=(I|A)$, where the two
blocks are $I$, the size-$k$ identity matrix, and a $k$ by $n-k$
matrix $A$.  In such a representation, the first $k$ positions are
called the information set of the code (since the corresponding
symbols are transmitted directly) and the remaining $n-k$ symbols
provide the redundancy.  Any $k$ linearly-independent columns of $G$
can be chosen as the information set, which defines the systematic
form of $G$ up to a permutation of the rows of $A$.

The IS algorithm and the original performance bounds
<Cite Key="Leon-1988"/>
<Cite Key="Kruk-1989"/>
<Cite Key="Coffey-Goodman-1990"/>
are based on the observation that for a long random code a set of
$k+\Delta$ randomly selected columns, with $\Delta$ of order one, are
likely to contain an information set. ISs are (approximately) in
one-to-one correspondence with the column permutations, and a random
IS can thus be generated as a set of **pivot** columns in the Gauss'
algorithm after a random column permutation.  Thus, if there is a
codeword $c$ of weight $d$, the probability to find it among the rows
of reduced-row-echelon form $G_S$ after a column permutation $S$ can
be estimated as that for a randomly selected set of $k$ columns to hit
exactly one non-zero position in $c$.

The statistics of ISs is more complicated in other ensembles of random
codes, e.g., in linear **low-density parity-check** (LDPC) codes where
the check matrix $H$ (of rank $n-k$ and with rows orthogonal to those
of $G$) is additionally required to be sparse.  Nevertheless, a
provable bound can be obtained for a related **covering set** (CS)
algorithm where a randomly selected set of $s\ge k-1$ positions of a
putative codeword are set to be zero, and the remaining positions are
constructed with the help of linear algebra.  In this case, the
optimal choice <Cite Key="Dumer-Kovalev-Pryadko-IEEE-2017"/> is to
take $s\approx n(1-\theta)$, where $\theta $ is the erasure threshold
of the family of the codes under consideration.  Since $\theta\ge R$
(here $R=k/n$ is the code rate), here more zeros must be selected, and
the complexity would grow (assuming the distance $d$ remains the same,
which is usually **not** the case for LDPC codes).

Note however that rows of $G_S$ other than the last are not expected
to contain as many zeros (e.g., the first row is only guaranteed to
have $k-1$ zeros), so it is **possible** that the performance of the
IS algorithm on classical LDPC codes is actually closer to that on
random codes as estimated by Leon <Cite Key="Leon-1988"/>.

@LatexOnly \medskip\noindent\underline{%
Case of **quantum CSS codes**
@LatexOnly }

In the case of a random CSS code (with matrices $P$ and $Q$ selected
randomly, with the only requirement being the orthogonality between
the rows of $P$ and $Q$), the performance of the algorithm <Ref
Subsect="Subsection_TheProcedure"/> can be estimated as that of the CS
algorithm, in terms of the erasure threshold of a linear code with the
parity matrix $P$, see <Cite Key="Dumer-Kovalev-Pryadko-IEEE-2017"/>.

Unfortunately, such an estimate fails dramatically in the case of
**quantum LDPC codes**, where rows of $P$ and $Q$ have weights bounded
by some constant $w$.  This is a reasonable requirement since the
corresponding quantum operators (supported on $w$ qudits) have to
actually be measured frequently as a part of the operation of the
code, and it is reasonable to expect that the measurement accuracy
goes down (exponentially) quickly as $w$ is increased.  Then, the
linear code orthogonal to the rows of $P$ has the distance $\le w$
(the minimal weight of the rows of $Q$), and the corresponding erasure
threshold is exactly zero.  In other words, there is a finite
probability that a randomly selected $w$ symbols contain a vector
orthogonal to the rows of $P$ (and such a vector would likely have nothing to
do with non-trivial **quantum** codewords which must be linearly
independent from the rows of $Q$).

On the other hand, for every permutation $S$ in the algorithm <Ref
Subsect="Subsection_TheProcedure"/>, the matrix $G_S$ contains exactly
$k=n-$`rank`$(P)-$`rank`$(Q)$ rows orthogonal to rows of $P$ and
linearly independent from rows of $Q$ (with columns properly
permuted).  These vectors contain at least $s$ zeros, where
$[1-\theta_*(P,Q)] n\le s\le n-$`rank`$(Q)$, where $\theta_*(P,Q)$ is
the erasure threshold for $Z$-like codewords in the quantum CSS code
with $H_X=P$ and $H_Z=Q$.


@LatexOnly \medskip\noindent\underline{%
**What is it that we do not understand?**
@LatexOnly }

What missing is an understanding of the statistics of the ISs of
interest, namely, the ISs that overlap with a minimum-weight codeword
in one (or a few) positions.  

Second, we know that a given column permutation $S$ leads to the
unique information set, and that every information set can be obtained
by a suitably chosen column permutation.  However, there is no
guarantee that the resulting information sets have equal
probabilities.  In fact, it is easy to construct small matrices where
different information sets are obtained from different numbers of
column permutations (and thus have **different** probabilities).  It is
not clear whether some of the ISs may have vanishingly small
probabilities in the limit of large codes; in such a case the
algorithm may take an excessively long time to converge.

@Section Empirical estimate of the success probability
@SectionLabel Empirical

The probability to find a codeword after $N$ rounds of the algorithm
can be estimated empirically, by counting the number of times each
codeword of the minimum weight was discovered.  We **expect** the
probability $P(c)$ to discover a given codeword $c$ to depend only on
its (symplectic) weight `wgt`$(c)$, with the probability a
monotonously decreasing function of the weight.  If, after $N$ steps,
codewords $c_1$, $c_2$, $\ldots$ , $c_m$ of the same (minimal) weight $w$
are discovered $n_1$, $n_2$, $\ldots$ , $n_m$ times, respectively, we can
estimate the corresponding Poisson parameter as
$$ \lambda_w =\frac{1}{N m}\sum_{i=1}^m n_i. $$

Then, the probability that a codeword $c_0$ of the true minimal weight
$ d < w $ be **not** discovered after $N$ steps can be upper
bounded as (the inequalities tend to saturate and become equalities in
the limit of small $\lambda_w$) $$ P_{\rm fail} < (1-\lambda_w)^N
< e^{-N\lambda_w}=\exp\left(-m^{-1}\sum_{i=1}^m n_i\right)\equiv
\exp(-\langle n\rangle).  $$ Thus, the probability to fail is
decreasing as an exponent of the parameter $\langle n\rangle$, the
**average number of times a minimum-weight codeword has been found.**

The hypothesis about all $P(c_i)$ being equal to $\lambda_w$ is
testable, e.g., if one considers the distribution of the ratios
$x_i=n_i/N$, where $N=\sum_{i=1}^m n_i$ is the total number of
codewords found.  These quantities sum up to one and are distributed
according to multinomial distribution<Cite
Key="Steel-1953"/>.  Further, under our assumption
of all $P(c_i)$ being equal, we also expect the outcome probabilities
in the multinomial distribution to be all equal, $\pi_i=1/m$, $1\le
i\le m$.

This hypothesis can be tested using Pearson's $\chi^2$ test.  Namely,
in the limit where the total number of observations $N$ diverges, the
quantity $$ X^2=\sum_{i=1}^m \frac{(n_i-N \pi_i)^2}{ N\pi_i}=
N^{-1}\sum_{i=1}^m \frac{n_i^2}{\pi_i}-N
\stackrel{\pi_i=1/m}\to\frac{m}{N}\sum_{i=1}^m n_i^2-N, $$ is expected to be
distributed according to the $\chi^2_{m-1}$ distribution with $m-1$
parameters, see <Cite Key="Chernoff-Lehmann-1954"/> <Cite
Key="Cramer-book-1999"/>.

In practice, we can approximate with the $\chi^2_{m-1}$ distribution
as long as the total $N$ be large compared to the number $m$ of the
codewords found (i.e., the average $\langle n\rangle$ must be large,
which is the same condition as needed for confidence in the result.)

With `debug[4]` set (binary value 8) in `DistRandCSS` and
`DistRandStab` (<Ref Sect="Section_DistanceFunctions"/>), whenever more
than one minimum-weight vector is found, the quantity $X^2$ is
computed and output along with the average number of times $\langle
n\rangle$ a minimum-weight codeword has been found.  However, no
attempt is made to analyze the corresponding value or calculate the
likelihood of the null hypothesis that the codewords be equiprobable.

@Chapter FileFormat

@Section General information


The code supports reading matrices from an MTX file using `ReadMTXE`
and writing new MTX file using `WriteMTXE` functions.
Below a description of the format is given.

@Subsection Representation of field elements via integers

Every finite field is isomorphic to a Galois field $F=\mathop{\rm GF}(q)$, where
$q$ is a power of a prime, $q=p^m$.

* For a prime field, with $q=p$ a prime, $F$ is a prime field, isomorphic to the
  ring Z$(q)$ of integers modulo
  $q$.  In such a case, elements of the field are stored
  directly as integers.  After reading, these are taken modulo $p$,
  thus, e.g., with $p=7$, values $-1$, $6$, or $13$ are all equivalent.


* When $q=p^m$ with $m > 1$, $F$ is an extension field.  Non-zero
  elements of such a field can be represented as integer powers of a
  primitive element $\alpha$, a primitive root of unity in the field.
  Primitive element is a root of a primitive polynomial $f(x)$, that
  is, an irreducible polynomial with coefficients in the corresponding
  prime field $\mathop{\rm GF}(p)$, with the property that the smallest integer $n$
  such that $f(x)$ divides $x^n-1$ is $n=p^m-1$.  Alternatively, field
  elements can be represented as $p$-ary polynomials modulo the chosen
  primitive polynomial $f(x)$.
  
  Either definition requires that the primitive polynomial be
  specified.  By default, GAP uses Conway polynomials to represent
  field elements.  For details, as well as a large collection of
  Conway polynomials, please see the web page maintained by Frank
  Luebeck <Cite Key="ConwayPol-2022"/>.

In the actual file format, three different and mutually exclusive
storage options for elements of an extension field can in be used.
    
* First, assuming all elements needed are actually in the
  corresponding prime field $\mathop{\rm GF}(p)$, the integer values (mod $p$)
  directly correspond to the values of field elements.  This is the
  case, e.g., with $0,\pm1$ matrices which obey the orthogonality
  condition already over integers, and retain orthogonality over
  Z$(q)$ with any $q$ (what is more relevant here, the same matrix
  would work with any prime field).  (**this is currently not implemented**)
  
* Second option, and the only option currently implemented for
  extension fields, is to store the powers of the primitive element
  for each non-zero element in the field, and $-1$ for zero.
  
* Third option, is to take the coefficients of $p$-ary polynomial
  $a_0+a_1x+a_2x^2+\ldots +a_{m-1}x^{m-1}$ as
  digits of a $p$-ary integer $(a_{m-1}a_{m-2}\ldots a_2a_1)_p$. 
  (**this is currently not implemented**)
        
@Subsection Matrix storage format

There are two recommended storage formats.

1. For CSS matrices stored in
   separate files, the MTX header should use the `integer` type, with matrix
   elements stored in the usual order.

2. For general stabilizer codes, or to store both CSS matrices in a
   single file, the MTX header should use the `complex` type.  In this
   case the block matrix $(A,B)$ is stored as a complex matrix $A+iB$.

In both format versions, the number of columns specified in the file
coincides with the code length.

Two additional matrix format versions supported by `ReadMTXE` and
`WriteMTXE` are provided for compatibility.  Here, the columns $a_i$
and $b_i$ in the blocks $A$ and $B$ are listed individually, and are
either intercalated [the ordering $(a_1, b_1,
a_2,b_2,\ldots,a_n,b_n)$] or are separated into column blocks
$(a_1,\ldots,a_n,b_1,\ldots,b_n)$.  In both cases the number of
columns in the matrix stored is twice the code length.

The ordering of the columns is governed by a parameter
`pair`, optional in the function `ReadMTXE` and required in `WriteMTXE`.

* With `pair=0` the matrix elements are stored in the usual order.  In
  this case the MTX header should use the `integer` type.  This is the
  defailt storage format for stabilizer generator matrices of CSS
  codes, and also the internal matrix format for single-block matrices
  accepted by the function `DistRandCSS`, see Section <Ref
  Sect="Section_DistanceFunctions"/>.

* With `pair=1` the block matrix $(A,B)$ is stored with intercalated
  columns $(a_1,b_1,\ldots,a_n,b_n)$; the MTX header should use the
  `integer` type.  This is the internal matrix format for two-block
  matrices accepted by the functions `DistRandStab` (<Ref
  Sect="Section_DistanceFunctions"/>) and `WriteMTXE` (<Ref
  Sect="Section_IOFunctions"/>), and returned by `ReadMTXE` (<Ref
  Sect="Section_IOFunctions"/>).

* With `pair=2` the block matrix $(A,B)$ is stored with separated
  columns $(a_1,\ldots,a_n, b_1,\ldots,b_n)$; the MTX header should also
  use the `integer` type.

* With `pair=3` the block matrix $(A,B)$ is stored as a complex matrix
  $A+iB$, with columns $(a_1 + i b_1,\ldots,a_n + i b_n)$. In this
  case `type=complex`, since matrix elements are represented as
  complex integers.

By default, `pair=0` corresponds to `type=integer` and `pair=3`
corresponds to `type=complex`.  **It is strongly recommended** that
matrices intended for use by others should only use these two variants of
the MTXE format.

For efficiency reasons, the function `DistRandStab` (<Ref
Sect="Section_DistanceFunctions"/>) assumes the generator matrix with
intercalated columns.

@Subsection Explicit format of each line

The first line must have the following form:
@BeginCode MMXFormatLine1
%%MatrixMarket matrix coordinate `type` general
@EndCode
@InsertCode MMXFormatLine1
with `type` either `integer` or `complex`.

The second line is optional and specifies the field, the primitive
polynomial used (in the case of an extension field), and the storage
format of field elements. 

@InsertCode FieldHeaderA

Here the records should be separated by one or more spaces;
while `field`, `polynomial`, and `format` **should not contain any spaces.**
Any additional records in this line will be silently ignored.

The `field` option should specify a valid field, `GF(q)` or
`GF(p^m)`, where $q>1$ should be a power of the prime $p$. 

The `polynomial` should be a valid expanded monic polynomial with
integer coefficients, with a single independent variable `x`; it
should contain no spaces.  An error will be signaled if `polynomial`
is not a valid primitive polynomial of the `field`.  This argument is
optional; if not specified, one may assume that the Conway polynomial
should be used.

The optional `format` string should be "AdditiveInt" (the default for
prime fields), "PowerInt" (**currently** the default for extension
fields with $m>1$) or "VectorInt".

- `AdditiveInt` indicates that values listed
  are expected to be in the corresponding prime field and should be
  interpreted as integers mod $p$.  
  
- `PowerInt` indicates that field elements are
  represented as integers powers of the primitive element, root of
  the primitive polynomial, or $-1$ for the zero field element.  

- `VectorInt` corresponds to encoding coefficients of a degree-$(m-1)$
  $p$-ary polynomial representing field elements into a $p$-ary
  integer.  In this notation, any negative value will be taken mod
  $p$, thus $-1$ will be interpreted as $p-1$, the additive inverse of
  the field $1$.

The primitive polynomial must be written
explicitly as $x^m+a_{m-1}*x^{m-1}+\ldots+a_1*x+a_0$, where the
integer coefficients $a_i$ will be interpreted modulo `p`.  **The primitive polynomial
should not contain any spaces.**  

@BeginCode MMXFormatLine2 
% Field: GF(`q`) PrimitiveP(x): `polynomial` 
@EndCode 
@InsertCode MMXFormatLine2 

For example, with $q=5^2$, the Conway polynomial
$f_{5,2}(x)=x^2-x+2$, and the second line can read 
@BeginCode MMXFormatLine2A 
% Field: GF(25) PrimitiveP(x): x^2-x+2 Format: PowerInt
@EndCode
@InsertCode MMXFormatLine2A 

The following is an equivalent form of the
same polynomial and can also be used 
@BeginCode MMXFormatLine2B 
% Field: GF(25) PrimitiveP(x): x^2+4*x+2 
@EndCode 
@InsertCode MMXFormatLine2B 
The field may be left undefined; by default, it is $\mathop{\rm GF}(2)$, or it can
be specified by hand when reading the matrices. 
If the primitive polynomial is undefined, it will be assumed that the
Conway polynomial used internally by `GAP` should be used.

Next follows the comment section, with each line either empty or
starting with the `%` symbol:
@BeginCode MMXFormatLine3
% Example of the comment line
@EndCode
@InsertCode MMXFormatLine3

After the comment section, in agreement with MTX format, goes the line
giving the dimensions of the matrix and the number of non-zero elements:
@BeginCode MMXFormatLineA
rows     columns     `(number of non-zero elements)`
@EndCode
@InsertCode MMXFormatLineA

Then all non-zero elements are listed as three or four integers according to the `type`:

* `type=integer`:
@BeginCode  MMXFormatLineB
i     j     element[i,j]
@EndCode
@InsertCode MMXFormatLineB

* `type=complex`:       
@BeginCode  MMXFormatLineC
i     j     a[i,j]     b[i,j]
@EndCode
@InsertCode MMXFormatLineC

Notice that column and row numbers must start with 1, as prescribed in
the original MTX format.

@Subsection Matrix Storage Implementation Details

Neither the `WriteMTXE` nor `ReadMTXE` currently support the `Format:`
parameter.  Prime field elements are only stored "as is", i.e., as integers
to be taken modulo `p`, while extension field elements are only stored
in the `PowerInt` format, i.e., with the power of the primitive
element specified, or "-1" for zero.

The function `WriteMTXE` can only save field elements with the
primitive polynomial used internally by `GAP`, i.e., the Conway
polynomial.

The function `ReadMTXE` can read matrix elements specified (in the
case of an extension field) with any primitive polynomial as specified
in the file.

Given the field $\mathop{\rm GF}(p^m)$ and the primitive polynomial $p(x)$
specified in the file, the function `ReadMTXE` first checks that the
degree of $p(x)$ is indeed $m$ and that it is a primitive polynomial
in the corresponding prime field $\mathop{\rm GF}(p)$.  If either of these tests
fail, `ReadMTXE` produces an error.  Otherwise, it will attempt to
find the conversion coefficient $c$ such that $\alpha^c$ is a root of
$p(x)$, starting with $c=1$.  When found, the multiplicative inverse
$s$ such that $sc\equiv 1\bmod (q-1)$ will be used to convert the
elements being read, i.e., for any matrix element $y$ read, $y^s$ will
be used instead.

Notice that, unless the Conway polynomial was used (in which case
$c=s=1$, and the conversion is trivial), this search can be slow for
large fields, as all integer values in $[1,2,\ldots,q-2]$ will be
tested sequentially.  To help ensure that the correct polynomial is
used, it is recommended that orthogonality of matrices be checked.


@Section Example MTXE files

In this section we give two sample MTXE files storing the stabilizer generator
matrix of 5-qubit codes.

First, matrix (with one redundant linearly-dependent row) stored with
`type=integer` and `pair=1` (intercalated columns
$[a_1,b_1,a_2,b_2,\ldots]$) is presented.  Notice that the number of
columns is twice the actual length of the code. Even though the field
is specified explicitly, this matrix would work with any
prime field.

@BeginCode SampleFileA
%%MatrixMarket matrix coordinate integer general                
% Field: GF(7)
% 5-qubit code generator matrix / normal storage with intercalated cols
5 10 20
1 1 1
1 4 1
1 6 -1
1 7 -1
2 3 1
2 6 1
2 8 -1
2 9 -1
3 1 -1
3 5 1
3 8 1
3 10 -1
4 2 -1
4 3 -1
4 7 1
4 10 1
5 2 1
5 4 -1
5 5 -1
5 9 1
@EndCode
@InsertCode SampleFileA

This same matrix is stored in the file `matrices/n5k1A.mtx`.  This is
how the matrix can be read and distance calculated:
@BeginExample
        filedir:=DirectoriesPackageLibrary("QDistRnd","matrices");;
    lis:=ReadMTXE(Filename(filedir,"n5k1A.mtx" ));;
    Print("field ",lis[1],"\n");
#! field GF(7)
    dist:=DistRandStab(lis[3],100,0 : field:=lis[1]);
#! 3
@EndExample 


The same matrix can also be stored with `type=complex` and `pair=3`
(complex pairs $[a_1+i b_1,a_2+i b_2,\ldots]$).  In this format, the
number of columns equals the code length.  

@BeginCode SampleFileB
%%MatrixMarket matrix coordinate complex general
% works with any prime field
% 5-qubit code generator matrix / normal storage with intercalated cols
% [[5,1,3]]_p
4 5 16
1 1 1 0
1 2 0 1
1 3 0 -1
1 4 -1 0
2 2 1 0
2 3 0 1
2 4 0 -1
2 5 -1 0
3 1 -1 0
3 3 1 0
3 4 0 1
3 5 0 -1
4 1 0 -1
4 2 -1 0
4 4 1 0
4 5 0 1
@EndCode
@InsertCode SampleFileB

The matrix above is written in the file `matrices/n5k1.mtx`.  To
calculate the distance, we need to specify the field [unless we want
to use the default binary field].

@BeginExample
    lis:=ReadMTXE(Filename(filedir,"n5k1.mtx" ));;
    Print("field ",lis[1],"\n");
#! field GF(2)
    dist:=DistRandStab(lis[3],100,0,2 : field:=lis[1]);
#! 3
    q:=17;;
    lis:=ReadMTXE(Filename(filedir,"n5k1.mtx") : field:= GF(q));;
    Print("field ",lis[1],"\n");
#! field GF(17)
    dist:=DistRandStab(lis[3],100,0,2 : field:=lis[1]);
#! 3
@EndExample

Finally, the following is an example of a five-qudit code over
$\mathop{\rm GF}(2^3)$ constructed by the script `examples/cyclic.g`.
@BeginCode SampleFileC
%%MatrixMarket matrix coordinate complex general
% Field: GF(2^3) PrimitiveP(x): x^3+x+1
% code [[5,1,3]]_8
% cyclic w=4 x^6+Z(2^3)^4*x^5+Z(2^3)^4*x^3+Z(2)^0
% Powers of GF(8) primitive element and -1 for Zero are given
5 5 20
1 1 0 -1
1 2 -1 4
1 3 -1 4
1 4 0 -1
2 2 0 -1
2 3 -1 4
2 4 -1 4
2 5 0 -1
3 1 0 -1
3 3 0 -1
3 4 -1 4
3 5 -1 4
4 1 -1 4
4 2 0 -1
4 4 0 -1
4 5 -1 4
5 1 -1 4
5 2 -1 4
5 3 0 -1
5 5 0 -1
@EndCode
@InsertCode SampleFileC