Array Arithmetic¶
After creating a FieldArray subclass and one or two arrays, nearly any arithmetic operation can be
performed using normal NumPy ufuncs or Python operators.
In the sections below, the finite field \(\mathrm{GF}(3^5)\) and arrays \(x\) and \(y\) are used.
In [1]: GF = galois.GF(3**5)
In [2]: x = GF([184, 25, 157, 31]); x
Out[2]: GF([184, 25, 157, 31], order=3^5)
In [3]: y = GF([179, 9, 139, 27]); y
Out[3]: GF([179, 9, 139, 27], order=3^5)
Ufuncs¶
NumPy ufuncs are universal functions that operate on scalars. Unary ufuncs operate on a single scalar and binary ufuncs operate on two scalars. NumPy extends the scalar operation of ufuncs to operate on arrays in various ways. This extensibility enables NumPy broadcasting.
Expand any section for more details.
Addition: x + y == np.add(x, y)
In [4]: x + y
Out[4]: GF([ 81, 7, 215, 58], order=3^5)
In [5]: np.add(x, y)
Out[5]: GF([ 81, 7, 215, 58], order=3^5)
Additive inverse: -x == np.negative(x)
In [6]: -x
Out[6]: GF([ 98, 14, 206, 62], order=3^5)
In [7]: np.negative(x)
Out[7]: GF([ 98, 14, 206, 62], order=3^5)
Any array added to its additive inverse results in zero.
In [8]: x + np.negative(x)
Out[8]: GF([0, 0, 0, 0], order=3^5)
Subtraction: x - y == np.subtract(x, y)
In [9]: x - y
Out[9]: GF([17, 16, 18, 4], order=3^5)
In [10]: np.subtract(x, y)
Out[10]: GF([17, 16, 18, 4], order=3^5)
Multiplication: x * y == np.multiply(x, y)
In [11]: x * y
Out[11]: GF([ 41, 225, 106, 123], order=3^5)
In [12]: np.multiply(x, y)
Out[12]: GF([ 41, 225, 106, 123], order=3^5)
Scalar multiplication: x * z == np.multiply(x, z)
Scalar multiplication is essentially repeated addition. It is the “multiplication” of finite field elements and integers. The integer value indicates how many additions of the field element to sum.
In [13]: x * 4
Out[13]: GF([184, 25, 157, 31], order=3^5)
In [14]: np.multiply(x, 4)
Out[14]: GF([184, 25, 157, 31], order=3^5)
In [15]: x + x + x + x
Out[15]: GF([184, 25, 157, 31], order=3^5)
In finite fields \(\mathrm{GF}(p^m)\), the characteristic \(p\) is the smallest value when multiplied by any non-zero field element that results in \(0\).
In [16]: p = GF.characteristic; p
Out[16]: 3
In [17]: x * p
Out[17]: GF([0, 0, 0, 0], order=3^5)
Multiplicative inverse: y ** -1 == np.reciprocal(y)
In [18]: y ** -1
Out[18]: GF([ 71, 217, 213, 235], order=3^5)
In [19]: GF(1) / y
Out[19]: GF([ 71, 217, 213, 235], order=3^5)
In [20]: np.reciprocal(y)
Out[20]: GF([ 71, 217, 213, 235], order=3^5)
Any array multiplied by its multiplicative inverse results in one.
In [21]: y * np.reciprocal(y)
Out[21]: GF([1, 1, 1, 1], order=3^5)
Division: x / y == x // y == np.divide(x, y)
In [22]: x / y
Out[22]: GF([237, 56, 122, 126], order=3^5)
In [23]: x // y
Out[23]: GF([237, 56, 122, 126], order=3^5)
In [24]: np.divide(x, y)
Out[24]: GF([237, 56, 122, 126], order=3^5)
Remainder: x % y == np.remainder(x, y)
In [25]: x % y
Out[25]: GF([0, 0, 0, 0], order=3^5)
In [26]: np.remainder(x, y)
Out[26]: GF([0, 0, 0, 0], order=3^5)
Divmod: divmod(x, y) == np.divmod(x, y)
In [27]: x / y, x % y
Out[27]: (GF([237, 56, 122, 126], order=3^5), GF([0, 0, 0, 0], order=3^5))
In [28]: divmod(x, y)
Out[28]: (GF([237, 56, 122, 126], order=3^5), GF([0, 0, 0, 0], order=3^5))
In [29]: np.divmod(x, y)
Out[29]: (GF([237, 56, 122, 126], order=3^5), GF([0, 0, 0, 0], order=3^5))
In [30]: q, r = divmod(x, y)
In [31]: q*y + r == x
Out[31]: array([ True, True, True, True])
Exponentiation: x ** z == np.power(x, z)
In [32]: x ** 3
Out[32]: GF([175, 76, 218, 192], order=3^5)
In [33]: np.power(x, 3)
Out[33]: GF([175, 76, 218, 192], order=3^5)
In [34]: x * x * x
Out[34]: GF([175, 76, 218, 192], order=3^5)
Square root: np.sqrt(x)
# Ensure the elements of x have square roots
In [35]: x.is_quadratic_residue()
Out[35]: array([ True, True, True, True])
In [36]: z = np.sqrt(x); z
Out[36]: GF([102, 109, 14, 111], order=3^5)
In [37]: z ** 2 == x
Out[37]: array([ True, True, True, True])
Logarithm: np.log(x)
In [38]: z = np.log(y); z
Out[38]: array([60, 2, 59, 3])
In [39]: α = GF.primitive_element; α
Out[39]: GF(3, order=3^5)
In [40]: α ** z == y
Out[40]: array([ True, True, True, True])
Ufunc methods¶
FieldArray instances support NumPy ufunc methods. Ufunc methods allow
a user to apply a NumPy ufunc in a unique way across the target array. All arithmetic ufuncs are supported.
Expand any section for more details.
reduce()
The reduce methods reduce the input array’s dimension by one, applying the ufunc across one axis.
In [41]: np.add.reduce(x)
Out[41]: GF(7, order=3^5)
In [42]: x[0] + x[1] + x[2] + x[3]
Out[42]: GF(7, order=3^5)
In [43]: np.multiply.reduce(x)
Out[43]: GF(105, order=3^5)
In [44]: x[0] * x[1] * x[2] * x[3]
Out[44]: GF(105, order=3^5)
accumulate()
The accumulate methods accumulate the result of the ufunc across a specified axis.
In [45]: np.add.accumulate(x)
Out[45]: GF([184, 173, 57, 7], order=3^5)
In [46]: GF([x[0], x[0] + x[1], x[0] + x[1] + x[2], x[0] + x[1] + x[2] + x[3]])
Out[46]: GF([184, 173, 57, 7], order=3^5)
In [47]: np.multiply.accumulate(x)
Out[47]: GF([184, 9, 211, 105], order=3^5)
In [48]: GF([x[0], x[0] * x[1], x[0] * x[1] * x[2], x[0] * x[1] * x[2] * x[3]])
Out[48]: GF([184, 9, 211, 105], order=3^5)
reduceat()
The reduceat methods reduces the input array’s dimension by one, applying the ufunc across one axis
in-between certain indices.
In [49]: np.add.reduceat(x, [0, 3])
Out[49]: GF([57, 31], order=3^5)
In [50]: GF([x[0] + x[1] + x[2], x[3]])
Out[50]: GF([57, 31], order=3^5)
In [51]: np.multiply.reduceat(x, [0, 3])
Out[51]: GF([211, 31], order=3^5)
In [52]: GF([x[0] * x[1] * x[2], x[3]])
Out[52]: GF([211, 31], order=3^5)
outer()
The outer methods applies the ufunc to each pair of inputs.
In [53]: np.add.outer(x, y)
Out[53]:
GF([[ 81, 166, 80, 211],
[165, 7, 155, 52],
[ 54, 139, 215, 103],
[198, 40, 89, 58]], order=3^5)
In [54]: np.multiply.outer(x, y)
Out[54]:
GF([[ 41, 192, 93, 97],
[ 91, 225, 193, 196],
[ 80, 211, 106, 145],
[149, 41, 129, 123]], order=3^5)
at()
The at methods performs the ufunc in-place at the specified indices.
In [55]: z = x.copy()
# Negate indices 0 and 1 in-place
In [56]: np.negative.at(x, [0, 1]); x
Out[56]: GF([ 98, 14, 157, 31], order=3^5)
In [57]: z[0:1] *= -1; z
Out[57]: GF([ 98, 25, 157, 31], order=3^5)
Advanced arithmetic¶
Convolution: np.convolve(x, y)
In [58]: np.convolve(x, y)
Out[58]: GF([ 79, 80, 5, 79, 167, 166, 123], order=3^5)
FFT: np.fft.fft(x)
The Discrete Fourier Transform (DFT) of size \(n\) over the finite field \(\mathrm{GF}(p^m)\) exists when there exists a primitive \(n\)-th root of unity. This occurs when \(n\ |\ p^m - 1\).
In [59]: GF = galois.GF(7**5)
In [60]: n = 6
# n divides p^m - 1
In [61]: (GF.order - 1) % n
Out[61]: 0
In [62]: x = GF.Random(n); x
Out[62]: GF([1646, 6326, 1752, 7345, 5189, 330], order=7^5)
In [63]: X = np.fft.fft(x); X
Out[63]: GF([ 139, 10300, 11803, 11743, 14247, 11050], order=7^5)
In [64]: np.fft.ifft(X)
Out[64]: GF([1646, 6326, 1752, 7345, 5189, 330], order=7^5)
See also ntt() and primitive_root_of_unity.
Inverse FFT: np.fft.ifft(X)
The inverse Discrete Fourier Transform (DFT) of size \(n\) over the finite field \(\mathrm{GF}(p^m)\) exists when there exists a primitive \(n\)-th root of unity. This occurs when \(n\ |\ p^m - 1\).
In [65]: GF = galois.GF(7**5)
In [66]: n = 6
# n divides p^m - 1
In [67]: (GF.order - 1) % n
Out[67]: 0
In [68]: x = GF.Random(n); x
Out[68]: GF([ 7908, 12487, 1602, 4632, 15761, 12631], order=7^5)
In [69]: X = np.fft.fft(x); X
Out[69]: GF([15744, 8433, 10265, 12377, 7082, 13770], order=7^5)
In [70]: np.fft.ifft(X)
Out[70]: GF([ 7908, 12487, 1602, 4632, 15761, 12631], order=7^5)
See also ntt() and primitive_root_of_unity.
Linear algebra¶
Linear algebra on FieldArray arrays/matrices is supported through both native NumPy linear algebra functions
in numpy.linalg and additional linear algebra routines not included in NumPy.
Expand any section for more details.
Dot product: np.dot(a, b)
In [71]: GF = galois.GF(31)
In [72]: a = GF([29, 0, 2, 21]); a
Out[72]: GF([29, 0, 2, 21], order=31)
In [73]: b = GF([23, 5, 15, 12]); b
Out[73]: GF([23, 5, 15, 12], order=31)
In [74]: np.dot(a, b)
Out[74]: GF(19, order=31)
Vector dot product: np.vdot(a, b)
In [75]: GF = galois.GF(31)
In [76]: a = GF([29, 0, 2, 21]); a
Out[76]: GF([29, 0, 2, 21], order=31)
In [77]: b = GF([23, 5, 15, 12]); b
Out[77]: GF([23, 5, 15, 12], order=31)
In [78]: np.vdot(a, b)
Out[78]: GF(19, order=31)
Inner product: np.inner(a, b)
In [79]: GF = galois.GF(31)
In [80]: a = GF([29, 0, 2, 21]); a
Out[80]: GF([29, 0, 2, 21], order=31)
In [81]: b = GF([23, 5, 15, 12]); b
Out[81]: GF([23, 5, 15, 12], order=31)
In [82]: np.inner(a, b)
Out[82]: GF(19, order=31)
Outer product: np.outer(a, b)
In [83]: GF = galois.GF(31)
In [84]: a = GF([29, 0, 2, 21]); a
Out[84]: GF([29, 0, 2, 21], order=31)
In [85]: b = GF([23, 5, 15, 12]); b
Out[85]: GF([23, 5, 15, 12], order=31)
In [86]: np.outer(a, b)
Out[86]:
GF([[16, 21, 1, 7],
[ 0, 0, 0, 0],
[15, 10, 30, 24],
[18, 12, 5, 4]], order=31)
Matrix multiplication: A @ B == np.matmul(A, B)
In [87]: GF = galois.GF(31)
In [88]: A = GF([[17, 25, 18, 8], [7, 9, 21, 15], [6, 16, 6, 30]]); A
Out[88]:
GF([[17, 25, 18, 8],
[ 7, 9, 21, 15],
[ 6, 16, 6, 30]], order=31)
In [89]: B = GF([[8, 18], [22, 0], [7, 8], [20, 24]]); B
Out[89]:
GF([[ 8, 18],
[22, 0],
[ 7, 8],
[20, 24]], order=31)
In [90]: A @ B
Out[90]:
GF([[11, 22],
[19, 3],
[19, 8]], order=31)
In [91]: np.matmul(A, B)
Out[91]:
GF([[11, 22],
[19, 3],
[19, 8]], order=31)
Matrix exponentiation: np.linalg.matrix_power(A, z)
In [92]: GF = galois.GF(31)
In [93]: A = GF([[14, 1, 5], [3, 23, 6], [24, 27, 4]]); A
Out[93]:
GF([[14, 1, 5],
[ 3, 23, 6],
[24, 27, 4]], order=31)
In [94]: np.linalg.matrix_power(A, 3)
Out[94]:
GF([[ 1, 16, 4],
[11, 9, 9],
[ 8, 24, 29]], order=31)
In [95]: A @ A @ A
Out[95]:
GF([[ 1, 16, 4],
[11, 9, 9],
[ 8, 24, 29]], order=31)
Matrix determinant: np.linalg.det(A)
In [96]: GF = galois.GF(31)
In [97]: A = GF([[23, 11, 3, 3], [13, 6, 16, 4], [12, 10, 5, 3], [17, 23, 15, 28]]); A
Out[97]:
GF([[23, 11, 3, 3],
[13, 6, 16, 4],
[12, 10, 5, 3],
[17, 23, 15, 28]], order=31)
In [98]: np.linalg.det(A)
Out[98]: GF(0, order=31)
Matrix rank: np.linalg.matrix_rank(A, z)
In [99]: GF = galois.GF(31)
In [100]: A = GF([[23, 11, 3, 3], [13, 6, 16, 4], [12, 10, 5, 3], [17, 23, 15, 28]]); A
Out[100]:
GF([[23, 11, 3, 3],
[13, 6, 16, 4],
[12, 10, 5, 3],
[17, 23, 15, 28]], order=31)
In [101]: np.linalg.matrix_rank(A)
Out[101]: 3
In [102]: A.row_reduce()
Out[102]:
GF([[ 1, 0, 0, 11],
[ 0, 1, 0, 25],
[ 0, 0, 1, 11],
[ 0, 0, 0, 0]], order=31)
Matrix trace: np.trace(A)
In [103]: GF = galois.GF(31)
In [104]: A = GF([[23, 11, 3, 3], [13, 6, 16, 4], [12, 10, 5, 3], [17, 23, 15, 28]]); A
Out[104]:
GF([[23, 11, 3, 3],
[13, 6, 16, 4],
[12, 10, 5, 3],
[17, 23, 15, 28]], order=31)
In [105]: np.trace(A)
Out[105]: GF(0, order=31)
In [106]: A[0,0] + A[1,1] + A[2,2] + A[3,3]
Out[106]: GF(0, order=31)
Solve a system of equations: np.linalg.solve(A, b)
In [107]: GF = galois.GF(31)
In [108]: A = GF([[14, 21, 14, 28], [24, 22, 23, 23], [16, 30, 26, 18], [4, 23, 18, 3]]); A
Out[108]:
GF([[14, 21, 14, 28],
[24, 22, 23, 23],
[16, 30, 26, 18],
[ 4, 23, 18, 3]], order=31)
In [109]: b = GF([15, 11, 6, 29]); b
Out[109]: GF([15, 11, 6, 29], order=31)
In [110]: x = np.linalg.solve(A, b)
In [111]: np.array_equal(A @ x, b)
Out[111]: True
Matrix inverse: np.linalg.inv(A)
In [112]: GF = galois.GF(31)
In [113]: A = GF([[14, 21, 14, 28], [24, 22, 23, 23], [16, 30, 26, 18], [4, 23, 18, 3]]); A
Out[113]:
GF([[14, 21, 14, 28],
[24, 22, 23, 23],
[16, 30, 26, 18],
[ 4, 23, 18, 3]], order=31)
In [114]: A_inv = np.linalg.inv(A); A_inv
Out[114]:
GF([[27, 17, 9, 8],
[20, 21, 12, 4],
[30, 10, 23, 22],
[13, 25, 6, 13]], order=31)
In [115]: A @ A_inv
Out[115]:
GF([[1, 0, 0, 0],
[0, 1, 0, 0],
[0, 0, 1, 0],
[0, 0, 0, 1]], order=31)
Additional linear algebra¶
Below are additional linear algebra routines provided for FieldArray arrays/matrices that are
not included in NumPy.
Row space: A.row_space()
In [116]: GF = galois.GF(31)
In [117]: A = GF([[23, 11, 3, 3], [13, 6, 16, 4], [12, 10, 5, 3], [17, 23, 15, 28]]); A
Out[117]:
GF([[23, 11, 3, 3],
[13, 6, 16, 4],
[12, 10, 5, 3],
[17, 23, 15, 28]], order=31)
In [118]: A.row_space()
Out[118]:
GF([[ 1, 0, 0, 11],
[ 0, 1, 0, 25],
[ 0, 0, 1, 11]], order=31)
See row_space() for more details.
Column space: A.column_space()
In [119]: GF = galois.GF(31)
In [120]: A = GF([[23, 11, 3, 3], [13, 6, 16, 4], [12, 10, 5, 3], [17, 23, 15, 28]]); A
Out[120]:
GF([[23, 11, 3, 3],
[13, 6, 16, 4],
[12, 10, 5, 3],
[17, 23, 15, 28]], order=31)
In [121]: A.column_space()
Out[121]:
GF([[ 1, 0, 0, 0],
[ 0, 1, 0, 11],
[ 0, 0, 1, 5]], order=31)
See column_space() for more details.
Left null space: A.left_null_space()
In [122]: GF = galois.GF(31)
In [123]: A = GF([[23, 11, 3, 3], [13, 6, 16, 4], [12, 10, 5, 3], [17, 23, 15, 28]]); A
Out[123]:
GF([[23, 11, 3, 3],
[13, 6, 16, 4],
[12, 10, 5, 3],
[17, 23, 15, 28]], order=31)
In [124]: A.left_null_space()
Out[124]: GF([[ 0, 1, 23, 14]], order=31)
See left_null_space() for more details.
Null space: A.null_space()
In [125]: GF = galois.GF(31)
In [126]: A = GF([[23, 11, 3, 3], [13, 6, 16, 4], [12, 10, 5, 3], [17, 23, 15, 28]]); A
Out[126]:
GF([[23, 11, 3, 3],
[13, 6, 16, 4],
[12, 10, 5, 3],
[17, 23, 15, 28]], order=31)
In [127]: A.null_space()
Out[127]: GF([[ 1, 22, 1, 14]], order=31)
See null_space() for more details.
Gaussian elimination: A.row_reduce()
In [128]: GF = galois.GF(31)
In [129]: A = GF([[23, 11, 3, 3], [13, 6, 16, 4], [12, 10, 5, 3], [17, 23, 15, 28]]); A
Out[129]:
GF([[23, 11, 3, 3],
[13, 6, 16, 4],
[12, 10, 5, 3],
[17, 23, 15, 28]], order=31)
In [130]: A.row_reduce()
Out[130]:
GF([[ 1, 0, 0, 11],
[ 0, 1, 0, 25],
[ 0, 0, 1, 11],
[ 0, 0, 0, 0]], order=31)
See row_reduce() for more details.
LU decomposition: A.lu_decompose()
In [131]: GF = galois.GF(31)
In [132]: A = GF([[4, 1, 24], [7, 6, 1], [11, 20, 2]]); A
Out[132]:
GF([[ 4, 1, 24],
[ 7, 6, 1],
[11, 20, 2]], order=31)
In [133]: L, U = A.lu_decompose()
In [134]: L
Out[134]:
GF([[ 1, 0, 0],
[25, 1, 0],
[26, 15, 1]], order=31)
In [135]: U
Out[135]:
GF([[ 4, 1, 24],
[ 0, 12, 21],
[ 0, 0, 24]], order=31)
In [136]: np.array_equal(L @ U, A)
Out[136]: True
See lu_decompose() for more details.
PLU decomposition: A.plu_decompose()
In [137]: GF = galois.GF(31)
In [138]: A = GF([[15, 4, 11], [7, 6, 1], [11, 20, 2]]); A
Out[138]:
GF([[15, 4, 11],
[ 7, 6, 1],
[11, 20, 2]], order=31)
In [139]: P, L, U = A.plu_decompose()
In [140]: P
Out[140]:
GF([[1, 0, 0],
[0, 0, 1],
[0, 1, 0]], order=31)
In [141]: L
Out[141]:
GF([[ 1, 0, 0],
[ 9, 1, 0],
[17, 0, 1]], order=31)
In [142]: U
Out[142]:
GF([[15, 4, 11],
[ 0, 15, 27],
[ 0, 0, 0]], order=31)
In [143]: np.array_equal(P @ L @ U, A)
Out[143]: True
See plu_decompose() for more details.
Last update:
May 17, 2022