Square Roots in Fields of Prime Power Order

[TDQuering]

Greetings!

It appears the functionality with np.sqrt to compute square roots of elements is not properly implemented for extension fields of prime fields. I ran into this issue while trying to compute the square root of 6 in GF(11^2), which incorrectly returns a value of 1. I did a little more testing, and I found that the package computes the square root of all nonzero elements of GF(11^2) as 1. After some further investigation, it seems to be giving incorrect square root results consistently in fields GF(p^k) with odd characteristic p and even extension degree k.

It makes sense that there are no issues in characteristic 2, as the square root method for binary fields uses a specific implementation that appears correct:

galois/src/galois/_domains/_calculate.py

Lines 758 to 768 in c945356

	class sqrt_binary(_lookup.sqrt_ufunc):
	"""
	A ufunc dispatcher that provides the square root in binary extension fields.
	"""
	def implementation(self, a: Array) -> Array:
	"""
	Fact 3.42 from https://cacr.uwaterloo.ca/hac/about/chap3.pdf.
	"""
	return a ** (self.field.characteristic ** (self.field.degree - 1))

I believe the reason we have problems in the odd characteristic case is because the algorithms referenced here are designed for computing square roots modulo p, i.e. in the prime field, but not necessarily in extension fields:

galois/src/galois/_domains/_calculate.py

Lines 770 to 779 in c945356

	class sqrt(_lookup.sqrt_ufunc):
	"""
	A ufunc dispatcher that provides square root using NumPy array arithmetic.
	"""
	def implementation(self, a: Array) -> Array:
	"""
	Algorithm 3.34 from https://cacr.uwaterloo.ca/hac/about/chap3.pdf.
	Algorithm 3.36 from https://cacr.uwaterloo.ca/hac/about/chap3.pdf.
	"""

It's interesting that these algorithms appear (under admittedly limited testing) to work in extension fields of odd degree. Ideally, it's probably worth verifying that this is a genuine mathematical property and not just a fluke due to the small amount of testing I've done to research this issue. But it appears a different technique is definitely required to handle extension fields of even extension degree. I've found a paper on the problem of square roots in finite fields of even extension degree which looks promising, but I haven't had the chance to fully investigate it myself.

Also, I just wanted to say that I really appreciate those who have worked on this package - it's been very helpful to me. I hope this can be used to improve it even more.

[mhostetter]

I can reproduce this bug.

import galois

GF = galois.GF(11**2)
x = GF.squares
print(x)
y = np.sqrt(x)
print(y)
print(y**2)

[  0   1   2   3   4   5   6   7   8   9  10  13  14  15  16  20  26  28
  29  30  32  34  37  38  39  42  45  47  49  52  53  56  58  59  64  65
  67  68  73  74  76  79  80  83  85  87  90  93  94  95  98 100 102 103
 104 106 112 116 117 118 119]
[0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]
[0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]

[mhostetter]

Thank you for identifying this bug and doing as much debug as you already have.

From the textbook I referenced, it says

Perhaps I didn't "straightforwardly" extend the algorithms to extension fields? Maybe there's a bug in my Python code?

[TDQuering]

I've done some more reading and testing, and I believe I may have identified the issue. Luckily, I don't think it's as fundamental of a problem as I originally feared - the paper I link contends that the algorithms used here are not invalid in extension fields of even degree, just that they aren't "efficient" - I'm willing to wager that they're efficient enough for the purposes of this package, at least for the time being. But there was one other interesting thing in this paper that I believe identifies the bug in this code. In the notation of that paper, they consider fields GF(q) with q = p^k for some prime p. They give algorithms for the special cases where q is congruent to 3 modulo 4 or congruent to 5 modulo 8 (among others) which I believe line up with the algorithms implemented here. However, the implementation here applies those algorithms when p (rather than q) satisfies the relevant congruence.

For an odd extension degree k, this does not present a problem, as the congruences for p and q = p^k are equivalent. But they are not equivalent for an even extension degree. As an example: 11 is congruent to 3 modulo 4, but 11^2 = 169 is congruent to 1 modulo 4. This means that in these cases, galois relies on an algorithm that is not actually applicable, as it is checking the wrong congruence. If this theory is correct, there are some cases where the package should give the correct results for square roots in fields of even extension degree - this should occur when p and p^k are congruent to one another modulo 4 and 8. p = 17, k = 2 is such a case: both 17 and 17^2 are congruent to 1 modulo both 4 and 8. And indeed, galois gives correct square roots in this case:

import numpy as np
import galois

GF = galois.GF(17**2) 
squares = GF.squares
print(squares) 
roots = np.sqrt(squares)
print(roots) 
print(roots **2) 
print(np.all(squares == roots**2))

[  0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  19
  20  22  24  26  28  30  31  35  38  39  40  43  44  45  48  55  56  57
  59  60  61  66  67  69  70  71  73  76  78  79  80  86  87  89  93  95
  96  99 100 103 105 110 112 114 115 117 118 120 123 125 128 131 132 133
 134 138 139 140 141 142 143 146 152 154 160 163 164 165 166 167 168 172
 173 174 175 178 181 183 186 188 189 191 192 194 196 201 203 206 207 210
 211 213 217 219 220 226 227 228 230 233 235 236 237 239 240 245 246 247
 249 250 251 258 261 262 263 266 267 268 271 275 276 278 280 282 284 286
 287]
[  0   1   6 116   2  58  50  83   5   3  25 149  91   8 124   7   4  81
 148  61 135 109  86  46  17  42 129  79  64  98  26 102 147  38 123  93
  67 112 146  18  77  75  88 105  34 145 134  53  27  19 128  47 144 115
  73 100  56  43  59  95 143  71 108 122  28  90  69  20  62 118  39 142
 133 127 141  84  29  35 111  85  65  97  51 121 140  82  21  48 104 114
  44 139 132  92  80  30  54  78 107 126  22  40  57  87 138  99  60  31
 120  36 137 117  76  74  94 110  49 136 131  63  23  32  72  45 125  89
  66 103 152  41 119  70  52 101  24 113 151  68 150  55 130 106  96  37
  33]
[  0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  19
  20  22  24  26  28  30  31  35  38  39  40  43  44  45  48  55  56  57
  59  60  61  66  67  69  70  71  73  76  78  79  80  86  87  89  93  95
  96  99 100 103 105 110 112 114 115 117 118 120 123 125 128 131 132 133
 134 138 139 140 141 142 143 146 152 154 160 163 164 165 166 167 168 172
 173 174 175 178 181 183 186 188 189 191 192 194 196 201 203 206 207 210
 211 213 217 219 220 226 227 228 230 233 235 236 237 239 240 245 246 247
 249 250 251 258 261 262 263 266 267 268 271 275 276 278 280 282 284 286
 287]
True

As additional assurance, I wrote a quick rewrite of the code implemented in this library with the adjustments I believe are necessary. For this example, I also rewrote the function as a function which handles only one input, as I'm not experienced with wiring up numpy ufuncs. There's obviously no reason to do this in the actual implementation, I just did it to make the testing easier for myself. The code is as follows, and it indeed gives correct results for the case of GF(11^2) which presented the original issue.

import numpy as np
import galois

def galois_sqrt(a):
    if not a.is_square():
        raise ArithmeticError(
            f"{a} is not a square in {type(a).name}."
        )

    field = type(a)
    p = field.characteristic
    q = field.order

    # Swapped p for q here
    if q % 4 == 3:
        roots = a ** ((q + 1) // 4)

    # Swapped p for q here as well
    elif q % 8 == 5:
        d = a ** ((q - 1) // 4)
        roots = field(0)

        if d == 1: 
            roots = a ** ((q + 3) // 8)
        elif d == p - 1:
            roots = 2 * a * (4 * a) ** ((q - 5) // 8)

    else:
        # Find a non-square element `b`
        while True:
            b = field.Random(low=1)
            if not b.is_square():
                break

        # Write p - 1 = 2^s * t
        n = field.order - 1
        s = 0
        while n % 2 == 0:
            n >>= 1
            s += 1
        t = n
        assert q - 1 == 2**s * t

        roots = field(0)  # Roots will go here

        # Compute a root `r` if a is nonzero
        if a > 0:  # Where a has a reciprocal
            a_inv = np.reciprocal(a)
            c = b**t
            r = a ** ((t + 1) // 2)
            for i in range(1, s):
                d = (r**2 * a_inv) ** (2 ** (s - i - 1))
                if d == p - 1:
                    r *= c
                c = c**2
            roots = r  # Assign root if nonzero

    roots = field._view(np.minimum(roots, -roots))  # Return only the smaller root

    roots = field(roots)
    return roots

GF = galois.GF(11**2)
squares = GF.squares 
print(squares)
roots = GF([galois_sqrt(x) for x in squares])
print(roots)
print(roots**2)
print(np.all(squares == roots**2))

[  0   1   2   3   4   5   6   7   8   9  10  13  14  15  16  20  26  28
  29  30  32  34  37  38  39  42  45  47  49  52  53  56  58  59  64  65
  67  68  73  74  76  79  80  83  85  87  90  93  94  95  98 100 102 103
 104 106 112 116 117 118 119]
[ 0  1 20  5  2  4 56 38 29  3 47 54 15 32 55 40 50 65 24 42 21 16 33 46
 27 64 30 63 53 35 11 17 22 49 37 62 25 45 39 61 12 41 18 28 60 52 13 43
 48 31 59 44 58 23 34 19 36 51 14 26 57]
[  0   1   2   3   4   5   6   7   8   9  10  13  14  15  16  20  26  28
  29  30  32  34  37  38  39  42  45  47  49  52  53  56  58  59  64  65
  67  68  73  74  76  79  80  83  85  87  90  93  94  95  98 100 102 103
 104 106 112 116 117 118 119]
True

I tried out this same code on a few other finite fields of even extension degree, and it produced correct results there as well. So I believe the only change necessary may be a swap from p to self.field.order on lines 787 and 790 in the code linked below.

galois/src/galois/_domains/_calculate.py

Lines 775 to 799 in c945356

	def implementation(self, a: Array) -> Array:
	"""
	Algorithm 3.34 from https://cacr.uwaterloo.ca/hac/about/chap3.pdf.
	Algorithm 3.36 from https://cacr.uwaterloo.ca/hac/about/chap3.pdf.
	"""
	if not np.all(a.is_square()):
	raise ArithmeticError(
	f"Input array has elements that are non-squares in {self.field.name}.\n{a[~a.is_square()]}"
	)
	p = self.field.characteristic
	if p % 4 == 3:
	roots = a ** ((self.field.order + 1) // 4)
	elif p % 8 == 5:
	d = a ** ((self.field.order - 1) // 4)
	roots = self.field.Zeros(a.shape)
	idxs = np.where(d == 1)
	roots[idxs] = a[idxs] ** ((self.field.order + 3) // 8)
	idxs = np.where(d == p - 1)
	roots[idxs] = 2 * a[idxs] * (4 * a[idxs]) ** ((self.field.order - 5) // 8)

[mhostetter]

This was the first thing I tried, but apparently I didn't apply the idea correctly. I just re-tried, with your patch, and it does seem to work. Thanks for debugging this!

I can submit a PR soon, or you can if interested.

[TDQuering]

Done! See the PR here: #574