/* This file contains a rather straightforward implementation of the algorithms developed in my Ph.D. thesis "Deterministic equation solving over finite fields", Universiteit Leiden 2006. Implementation in Magma by Christiaan van de Woestijne, Spring 2006. E-mail of author: c.vandewoestijne@tugraz.at. Features: - Tonelli-Shanks algorithm for taking roots in finite fields - deterministic version of T-S - n'th power generator of a finite field - writing finite field elements as sums of n'th powers - solving diagonal equations over finite fields whose number of variables exceeds the (true) degree and the Main Feature: - all algorithms here are deterministic and efficient (except those, like TonelliShanks and FindRationalPointProbabilistically, that were provided for purposes of comparison) At the end of the file there are some examples (commented out, uncomment to include these). MORE DETAILS ON THE INDIVIDUAL FUNCTIONS Thesis CHAPTER 2: ConvertPower := function( a, n ) Let d = gcd(n,#F-1), where F is the finite field containing a. We want b in F such that b^n = a^d. By the cyclic structure of F*, it's easy to see that b is in the group generated by a. The function returns b and d. ExtractPower := function( a, d ) Extract the highest power of d out of a, where a and d are integers. We also return the quotient of a by the highest power of d that divides it. PowerOrder := function( a, l ) Compute w such that l^w is the multiplicative order of the finite field element a. NOTICE: this generates an infinite loop if the order of a is not a power of l! If successful and if the order of a is nontrivial, the function also returns the last nontrivial power of a (of l-power exponent). FindRationalPointProbabilistically := function( a, b, n ) Find one nonzero solution to the equation a1 x1^n + ... + a_v x_v^n = b where the a_i and b are elements of a finite field F, with a_i nonzero. Use a probabilistic method: guess x1 up to x_{v-1}, then try to compute x_v by taking an n'th root. Might fail, of course; if v and n are reasonable with respect to #F, we do 2d tries, where d=gcd(n,#F-1). Thesis CHAPTER 3: TonelliShanks := function( a, n ) Compute an n'th root of the finite field element a (if it has such a root), using the algorithm of Tonelli and Shanks. DeterministicTonelliShanks := function( a, g, n ) Compute an n'th root of the finite field element a, using the algorithm of Tonelli and Shanks in my deterministic variant. The multiplicative n-part of the root will be contained in the n-part of the group generated by g. This may fail, even if a is an n'th power, if the group generated by g does not contain an element of sufficiently large order. Here 'sufficiently large' means that g generates a sufficiently large portion of all l-Sylow subgroups of F*, for all primes dividing gcd( n, #F-1 ). If g is a primitive root of F, then the algorithm succeeds if and only if a is an n'th power in F. SelectiveRoot := function( a, n : verbose := false ) RecursiveSelectiveRoot := function( a, n : verbose := false ) From a list a of (at least) n+1 elements of a finite field, select two such that the quotient of the two is an n'th power inside the group generated by all the n+1 elements, and compute the root. Output: i,j,b, with 1 le i lt j le n+1, such that b^n eq a[i]/a[j]. If the option verbose is set to true, quite some diagnostic output is given. The first function is an iterative implementation, the second follows the naive recursive formulation given in the thesis. Thesis CHAPTER 4: MultiplicativePrimitiveElement := function( a, E ) Let M be the (finite) field generated by the a[i] over E. This function computes one generator for M over E that is a multiplicative combination of the a[i]. PowerGenerator := function( F, E, n ) Computes an element beta in F such that E(beta^n) is equal to the subfield of F generated by all n'th powers in F. E(beta^n) is equal to F whenever #F > gcd(n,#F-1)^2. Makes use of MultiplicativePrimitiveElement. E and F should be finite fields with E \subseteq F, and n should be a positive integer dividing #E-1. Thesis CHAPTER 5: SumOfIntegerPowers := function( a, n ) Reduce the size of a by writing it as a sum of integer n'th powers, plus some small part that is not usefully written as such a sum. More details: we return a sequence of integers r and an integer a' s.t. (1) a = sum [ rr^n : rr in r ] + a' (2) a' le n^n ExpandOnRationalBase := function( a, t, s ) Expand the integer a on the rational base t/s. See my thesis for more comments on what this does. If the input is right, then afterwards, we have a = c[1] + c[2]*(t/s) + ... + c[d]*(t/s)^(d-1) where d=#c a should be a nonnegative integer, and t and s should be integers with 0 < s < t. PowerSumRepresentation := function( a, n : verbose := false ) Represent the nonzero finite field element a as a sum of n'th powers in the same field, using at most n terms. Diagnostic output is given if the verbose option is set to true. ZeroPowerSum := function( F, n ) Find a set of nonzero elements of F such that their n'th powers sum to zero. Note that this problem is always solvable over the prime field already. This implementation just uses PowerSumRepresentation to give a sequence of n'th powers summing to -1. In my thesis, I give an asymptotically faster way for doing this, but it's a bit involved to implement, and this works also. Thesis CHAPTER 6: ZeroOfDiagonalForm := function( a, n : verbose := false ) Find a nontrivial zero to the diagonal form given by sum [ a[i]*x[i]^n : i in [1..#a] ] over the finite field F = Universe( a ). We should have #a gt n. Uses the trapezium method as defined in my thesis. RepresentationByDiagonalForm := function( a, b, n : verbose := false ) Find a representation of the element b by the diagonal form given by sum [ a[i]*x[i]^n : i in [1..#a] ] over the finite field F = Universe( a ). We should have #a gt n, and b in F. Uses the trapezium method as defined in my thesis. If #F < n^2, the only algorithm is just enumerating all solutions; this was not yet implemented. */ ConvertPower := function( a, n ) // Let d = gcd(n,#F-1), where F is the finite field containing a. // We want b in F such that b^n = a^d. By the cyclic structure of F*, // it's easy to see that b is in the group generated by a. F := Parent( a ); // finite field, hopefully if Type( F ) ne FldFin then error "Expecting a finite field element"; end if; d,x,y := XGCD( n, #F-1 ); if x le 0 then x +:= #F-1; y -:= n; end if; // Get x positive (if x eq 0, we replace it by q-1 to avoid problems with // a eq 0 return a^x, d ; end function; _FFConvertPower := function( F, a, n ) // Same as ConvertPower, but low level; assume the finite field is given. d,x,y := XGCD( n, #F-1 ); if x le 0 then x +:= #F-1; y -:= n; end if; // Get x positive (if x eq 0, we replace it by q-1 to avoid problems with // a eq 0 return a^x, d ; end function; ExtractPower := function( a, d ) // Extract the highest power of d out of a, where a and d are integers. // Yes, I know that Bernstein developed an essentially linear method for // this, but I don't feel like implementing that here. // We also return the quotient of a and the highest power of d that divides it. r := 0; while a mod d eq 0 do a := a div d; r := r + 1; end while; return r, a; end function; PowerOrder := function( a, l ) // Compute w such that l^w is the multiplicative order of the finite field // element a. // NOTICE: this generates an infinite loop if the order of a is not a power // of l! // If successful and if the order of a is nontrivial, the function also returns // the last nontrivial power of a. if a eq 0 then error "Order not defined for zero"; end if; w := 0; while a ne 1 do aprev := a; a := a^l; w := w+1; end while; if w gt 0 then return w, aprev; else return w, _; end if; end function; FindRationalPointProbabilistically := function( a, b, n ) // Find one nonzero solution to the equation // a1 x1^n + ... + a_v x_v^n = b // where the a_i and b are elements of a finite field F, with a_i nonzero. // Use a probabilistic method: guess x1 up to x_{v-1}, then try to compute // x_v by taking an n'th root. Might fail, of course; if v and n are // reasonable with respect to #F, we do 2d tries, where d=gcd(n,#F-1). if #a lt 1 then error "Expecting a nonempty list of coefficients"; end if; F := Parent( b ); if Type( F ) ne FldFin then F := Parent( a[1] ); // all elements of a list are in the same ring if Type( F ) ne FldFin then error "Expecting a finite field element"; end if; if not IsCoercible( F, b ) then error "Elements should belong to one finite field"; end if; else if not IsCoercible( F, a[1] ) then error "Elements should belong to one finite field"; end if; end if; if n le 0 then error "Expecting a positive integer for the degree"; end if; d := Gcd( n, #F-1 ); if d eq 1 then // trivial case: equation is linear return true, [ ConvertPower( a[1]/b, n ) ] cat [ 0 : i in [1..#a-1] ]; end if; if #F le (d-1)^2 then error "Too high degree to expect a solution"; else quality := 2*( Log( d-1 ) + Log( 2 ) ) / ( Log( #F ) - 2*Log( d-1 ) ); if #a lt quality then error "Too few variables to expect a solution in reasonable time"; end if; end if; for j in [1..2*d] do // Do 2d tries x := [ Random( F ) : i in [1..#a-1] ]; // guess first #a-1 coordinates testvalue := ( b - &+[ a[i] * x[i]^d : i in [1..#a-1] ] ) / a[#a]; // print "Testing:", x, d; bo, rt := IsPower( testvalue, d ); // try to compute last coordinate if bo then zeroq := &+[ a[i] * x[i]^d : i in [1..#a-1] ] + a[#a] * rt^d - b; if zeroq ne 0 then error "Some error occurred!"; end if; return true, [ ConvertPower( z, n ) : z in x cat [ rt ] ]; end if; end for; return false, _; end function; _PrimePowerTonelliShanks := function( F, a, l, f ) // Compute an l^f'th root of the finite field (F) element a, using the // algorithm of Tonelli and Shanks. Low level function, not meant to be // called by the end user. Is used by TonelliShanks. if not IsPrime( l ) or f lt 0 then error "Expecting an integer prime power exponent"; end if; if a eq 0 then // Root of 0 return a; end if; q := #F; r, u := ExtractPower( q-1, l ); if r lt f then // l^f doesn't divide q-1 error "The exponent should divide the group order"; end if; A := a^u; B := 1; g, x, y := XGCD( u, l^f ); if PowerOrder( A, l ) + f gt r then // Order of A is too large, no root exists return false, _; end if; h := F!1; // Find a non-l'th power by guessing while h^((q-1) div l) eq 1 do repeat h := Random( F ); until h ne 0; // Nonzero, please end while; g := h^u; // Remove all factors other than l from order gpowers := [ g ]; for i in [1..r-1] do gpowers[i+1] := gpowers[i]^l; end for; zeta := gpowers[r]; // really g^(l^(r-1)); zeta is an l'th root of unity zetapowers := [ zeta^i : i in [1..l-1] ]; while A ne 1 do // We reduce the order of A by at least one factor l // in every iteration. Notice: A*B^(l^f) is a loop // invariant, so when A=1, we are done. w, Apr := PowerOrder( A, l ); // Apr = A^(l^(w-1)), A^(l^w) = 1 z := Position( zetapowers, Apr ); // zeta^z = Apr A := A * g^( (l^w-z)*l^(r-w) ); // use l^w-z to avoid inversion B := B * g^( z*l^(r-w-f) ); // These computations can be done much more efficiently if one stores // powers of g with l-power exponents end while; return true, B^x * a^y; // Convert root of A into root of a end function; TonelliShanks := function( a, n ) // Compute an n'th root of the finite field element a, using the // algorithm of Tonelli and Shanks. if n le 0 or not IsIntegral( n ) then error "Expecting a positive integer exponent"; end if; F := Parent( a ); if Type( F ) ne FldFin then error "Expecting a finite field element"; end if; if a eq 0 or a eq 1 then // Root of 0 or 1 return a; end if; d := GCD( #F-1, n ); // We may reduce the exponent n to d if d eq 1 then // Trivial root return a; end if; dfac := Factorisation( d ); rootie := a; for pp in dfac do // Compute root prime power by prime power bo, rootie := _PrimePowerTonelliShanks( F, rootie, pp[1], pp[2] ); if not bo then return false, _; end if; end for; return true, _FFConvertPower( F, rootie, n ); // Restore exponent n end function; _PrimePowerDeterministicTonelliShanks := function( F, a, g, l, f ) // Compute an l^f'th root of the nonzero finite field element a, using the // algorithm of Tonelli and Shanks in my deterministic variant. This // may fail, even if a is an l^f'th power, if the group generated by g // does not contain an element of sufficiently large order. // We assume that a and g are elements of the finite field F, or at least // coercible into it. // Low level function, called by DeterministicTonelliShanks, not meant // to be called directly by the user. if not IsPrime( l ) or f lt 0 then error "Expecting an integer prime power exponent"; end if; q := #F; r, u := ExtractPower( q-1, l ); if r lt f then // l^f doesn't divide q-1 error "The exponent should divide the group order"; end if; A := a^u; B := 1; _, x, y := XGCD( u, l^f ); if A ne 1 then // Test if order of a is prime to l g := g^u; // Make sure order of g is a power of l ("ell") s := PowerOrder( g, l ); // print "In PrimePower: g has order", l, "^", s, "and a has order", // l, "^", PowerOrder( A, l ); if PowerOrder( A, l ) + f gt s then // Order of A is too large, no root exists (at least not inside the group // generated by g) return false, _; end if; gpowers := [ g ]; for i in [1..s-1] do gpowers[i+1] := gpowers[i]^l; end for; zeta := gpowers[s]; // really g^(l^(s-1)); zeta is an l'th root of unity zetapowers := [ zeta^i : i in [1..l-1] ]; while A ne 1 do w, Apr := PowerOrder( A, l ); // Apr = A^(l^(w-1)), A^(l^w) = 1 z := Position( zetapowers, Apr ); // zeta^z = Apr A := A * g^( (l^w-z)*l^(s-w) ); // use l^w-z instead of -z to avoid inversion B := B * g^( z*l^(s-w-f) ); // These computations can be done much more efficiently if one stores // powers of g with l-power exponents end while; end if; // End of nontrivial stuff that has to be done if A ne 1 return true, B^x * a^y; // Convert root of A into root of a end function; DeterministicTonelliShanks := function( a, g, n ) // Compute an n'th root of the finite field element a, using the // algorithm of Tonelli and Shanks in my deterministic variant. This // may fail, even if a is an n'th power, if the group generated by g // does not contain an element of sufficiently large order. Here // 'sufficiently large' means that g generates a sufficiently large // portion of all l-Sylow subgroups of F*, for all primes dividing // gcd( n, #F-1 ). If g is a primitive root of F, then the algorithm // succeeds if and only if a is an n'th power in F. if n le 0 or not IsIntegral( n ) then error "Expecting a positive integer exponent"; end if; F := Parent( a ); // Try to find an ambient finite field if Type( F ) ne FldFin then if Type( Parent( g ) ) eq FldFin then if IsCoercible( a, Parent(g) ) then F := Parent( g ); a := F!a; else error "Input should be elements of one finite field"; end if; else error "Expecting a finite field element"; end if; end if; if a eq 0 or a eq 1 then // Root of 0 or 1 return a; end if; d := GCD( #F-1, n ); // We may reduce the exponent n to d if d eq 1 then // Trivial root return a; end if; dfac := Factorisation( d ); rootie := a; for pp in dfac do // Compute root prime power by prime power bo, rootie := _PrimePowerDeterministicTonelliShanks( F, rootie, g, pp[1], pp[2] ); if not bo then return false, _; // Group generated by g doesn't contain a root end if; end for; return true, _FFConvertPower( F, rootie, n ); // Restore exponent n end function; forward RecursiveSelectiveRoot; RecursiveSelectiveRoot := function( a, n : verbose := false ) // From a list of (at least) n+1 elements of a finite field, select two // such that the quotient of the two is an n'th power inside the group // generated by all the n+1 elements, and compute the root. // Output: i,j,b, with 1 le i lt j le n+1, such that b^n eq a[i]/a[j]. F := Universe( a ); if Type( F ) ne FldFin then error "Expecting a list of finite field elements"; end if; if n le 0 or not IsIntegral( n ) then error "Expecting a positive integral exponent"; end if; q := #F; d := Gcd( n, q-1 ); // Can reduce the exponent n to d if #a le d then error "Length of list should exceed the (true) exponent"; end if; if 0 in a[ 1..d+1 ] then error "Expecting nonzero finite field elements"; end if; if d eq 1 then // Base case of recursion: trivial root return 1,2,_FFConvertPower( F, a[1]/a[2], n ); // Restore exponent n end if; pd := PrimeDivisors( d ); ell := pd[#pd]; // Largest prime divisor of d r, u := ExtractPower( q-1, ell ); f, _ := ExtractPower( d, ell ); if verbose then print "Considering prime power", ell, "^", f; end if; if verbose then print "Sylow group has depth", r, "and index", u; end if; A := [ a[i]^u : i in [1..d+1] ]; Aord := [ PowerOrder( A[i], ell ) : i in [1..d+1] ]; if verbose then print "The ell-ified elements are", A; end if; if verbose then print "Their (logarithmic) orders are", Aord; end if; maxord,maxloc := Maximum( Aord ); // A[maxloc] has maximal order ell^maxord if maxord le f then // Group gen. by the a[i] has no more factors ell than n? B := A; else B := [ A[i]^(ell^(maxord-f)) : i in [1..d+1] ]; end if; // At this moment, the following holds: B[i][2] eq B[j][2] if and only if the // group generated by the a[i] contains an ell^f'th root of a[i]/a[j]. newd := d div ell^f; if verbose then print "B is", B; end if; Sort( ~B, ~perm ); // Easiest way to find equals in B if verbose then print "Sorted B:", B, "with permutation", perm; end if; // Semantics of the permutations used with sorting are as follows: // the element that was at position i before sorting is at position i^perm // after. // In particular: if I find that after sorting the elements at positions // i and i+1 are equal, then the original positions of those elements are // i^Inverse(perm) and (i+1)^Inverse(perm). equals := 1; equalstart := 1; i := 1; while equals le newd and i lt d+1 do i := i+1; if B[i] eq B[i-1] then equals := equals+1; // Sequence of equals in B continues else equals := 1; equalstart := i; // Sequence breaks off; start from 1 again end if; end while; if equals le newd then error "Couldn't find enough equals"; else if verbose then print "Found", newd+1, "equals in sorted B starting from position", equalstart; end if; end if; if verbose then print "Continuing with a", [ i^perm : i in [equalstart..equalstart+newd] ]; end if; newa := a[ [ i^perm : i in [equalstart..equalstart+newd] ] ]; // Find the a[i] that correspond to the sequence of equals in B isub, jsub, rtsub := RecursiveSelectiveRoot( newa, newd ); // Recursive call ilift := ( equalstart + isub - 1 )^perm; jlift := ( equalstart + jsub - 1 )^perm; if verbose then print "Recursive call gives", isub, jsub, rtsub; end if; if verbose then print "Original indices are", ilift, "and", jlift; end if; bo, rt := _PrimePowerDeterministicTonelliShanks( F, rtsub, a[maxloc], ell, f ); if not bo then error "Oops! Could not compute root"; end if; if verbose then print "We are about to say that a[", ilift, "]/a[", jlift, "] = ", a[ilift]/a[jlift], "is equal to", [_FFConvertPower( F, rt, n )][1], "^", n, "=", rt^d; end if; if ilift lt jlift then return ilift, jlift, _FFConvertPower( F, rt, n ); else return jlift, ilift, _FFConvertPower( F, 1/rt, n ); end if; end function; FindEquals := function( S, newd ) equals := 1; equalstart := 1; i := 1; while equals le newd and i lt #S do i := i+1; if S[i] eq S[i-1] then equals := equals+1; // Sequence of equals in B continues else equals := 1; equalstart := i; // Sequence breaks off; start from 1 again end if; end while; if equals le newd then error "Couldn't find enough equals"; end if; return equalstart; end function; SelectiveRoot := function( a, n : verbose := false ) // From a list of (at least) n+1 nonzero elements of a finite field, select two // such that the quotient of the two is an n'th power inside the group // generated by all the n+1 elements, and compute the root. // Output: i,j,b, with 1 le i lt j le n+1, such that b^n eq a[i]/a[j]. F := Universe( a ); if Type( F ) ne FldFin then error "Expecting a list of finite field elements"; end if; if n le 0 or not IsIntegral( n ) then error "Expecting a positive integral exponent"; end if; q := #F; d := Gcd( n, q-1 ); // Can reduce the exponent n to d if d ne n then if verbose then print "It is enough to take a", d, "th root and restore to exponent", n, "at the end"; end if; end if; if #a le d then error "Length of list should exceed the (true) exponent"; end if; if 0 in a[ 1..d+1 ] then print "The list is", a[1..d+1]; error "Expecting nonzero finite field elements"; end if; if d eq 1 then // Trivial root return 1,2,_FFConvertPower( F, a[1]/a[2], n ); // Restore exponent n end if; dfac := Reverse( Factorisation( d ) ); // It is profitable in this // function to consider the largest // prime factors of d first tempd := d; generator := F!1; indexset := [1..d+1]; for pp in dfac do // Loop over prime power divisors of d ell := pp[1]; f := pp[2]; r, u := ExtractPower( q-1, ell ); if verbose then print "Considering prime power", ell, "^", f; end if; if verbose then print "Sylow group has depth", r, "and index", u; end if; A := [ a[i]^u : i in indexset ]; Aord := [ PowerOrder( x, ell ) : x in A ]; if verbose then print "The ell-ified elements are", A; end if; if verbose then print "Their (logarithmic) orders are", Aord; end if; maxord,maxloc := Maximum( Aord ); // A[maxloc] has maximal order ell^maxord generator := generator * A[maxloc]; // Add power of ell to order of generator if maxord le f then // Group gen. by the a[i] has no more factors ell than n? B := A; else B := [ x^(ell^(maxord-f)) : x in A ]; end if; // At this moment, the following holds: B[i][2] eq B[j][2] if and only if the // group generated by the A[i] contains an ell^f'th root of A[i]/A[j]. // Equivalently, if the group generated by the a[k], with k running over // indexset, contains an ell^f'th root of the corresponding quotient of a's. tempd := tempd div ell^f; // Prime ell is ready: remove it from d if verbose then print "B is", B; end if; Sort( ~B, ~perm ); // Easiest way to find equals in B if verbose then print "Sorted B:", B, "with permutation", perm; end if; // Semantics of the permutations used with sorting are as follows: // the element that was at position i before sorting is at position i^perm // after. // In particular: if I find that after sorting the elements at positions // i and i+1 are equal, then the original positions of those elements are // i^Inverse(perm) and (i+1)^Inverse(perm). equals := 1; equalstart := 1; i := 1; while equals le tempd and i lt #B do i := i+1; if B[i] eq B[i-1] then equals := equals+1; // Sequence of equals in B continues else equals := 1; equalstart := i; // Sequence breaks off; start from 1 again end if; end while; if equals le tempd then error "Couldn't find enough equals"; else if verbose then print "Found", tempd+1, "equals in sorted B starting from position", equalstart; end if; end if; if verbose then print "Continuing with indices", Sort( [ indexset[ i^perm ] : i in [equalstart..equalstart+tempd]] ); end if; indexset := [ indexset[ (equalstart+i)^perm ] : i in [0..tempd] ]; Sort( ~indexset ); // This is complicated: we start the loop considering the full list // a[1]..a[d+1] of a's. In every loop we restrict consideration to // a sublist of these, viz., those a[k] for which k is in indexset. // We construct the list B, with equal length as indexset, then we // select among the entries of B a sublist all whose elements are // equal, and whose length is the length of indexset for the next // loop. To be able to select these B-entries, we sort B and keep // the permutation that sorted it. In the sorted list B, the entries // B[equalstart] up to (including) B[equalstart+tempd] are equal; // therefore, before sorting, the entries B[equalstart^perm] up to // B[(equalstart+tempd)^perm] were equal; and using indexset, these // B-entries correspond to the entries a[ indexset[equalstart^perm] ] // up to (incl.) a[ indexset[(equalstart+tempd)^perm] ]. // This last assertion is used to obtain the new indexset. // Then, finally, we have to SORT the indexset, because the internal // list sorting routine of Magma is NOT conservative (it's quicksort). // Thus if it sorts a list, then equal elements may end up in a different // order. This is not very important, but at the end we want to give two // indices i and j with i gcd(n,#F-1)^2. Makes use of MultiplicativePrimitiveElement. // E and F should be finite fields with E \subseteq F, and n should be a // positive integer dividing #E-1. q := #E; e := Degree( F, E ); if n le 0 then error "Exponent should be a positive integer"; end if; n2 := Gcd( n, q^e-1 ); // print "Real exponent is", n2; if q^e le n2^2 then // Small field F? // print "The field is small."; t := (q^e-1) div n2; // Number of distinct n'th powers in F gammas := []; gammapowers := []; for x in F do // This is a rather stupid way of computing all // distinct n'th powers in F... if x^n2 notin gammapowers then Append( ~gammas, x ); Append( ~gammapowers, x^n2 ); end if; if #gammas ge t then // This already makes it less stupid! break; end if; end for; d, pos := Maximum( [ _FFDegree( gp, E ) : gp in gammapowers ] ); return _FFConvertPower( F, gammas[pos], n ), d; // gammapowers is a cyclic subgroup of F*, so a generator of this subgroup // certainly generates the subfield generated by all n'th powers in F. // We return an n'th root of such a generator, as well as the degree of the // subfield in question. // Converting the power is necessary as the computation was done with // n2'th powers, where n2 = gcd( n, #F-1 ). end if; beta := Generator( F, E ); // Important thingie if #E le n2 then // Small field E? // print "The base field is small."; gammas := []; // Get n elements in F whose n'th powers gammapowers := []; // are distinct for x in F do if x^n2 notin gammapowers then Append( ~gammas, x ); Append( ~gammapowers, x^n2 ); end if; if #gammas ge n2 then break; end if; end for; // We adjoin the n'th powers in gammapowers to E, so that E has at // least n+1 elements (including zero); next we want n+1 distinct // elements in E (the enlarged E), so why don't choose the elements // just computed plus zero. xx := MultiplicativePrimitiveElement( gammapowers cat [ ( beta + gamma )^n2 : gamma in gammas cat [0] ], E ); // Compute a generator of F over E that is a multiplicative // combination of // 1) the n'th powers of the elements in gammas // 2) the n'th powers of (beta + c_i) where c_i runs over // the n elements in gammas, plus zero // See my thesis for an argument why these n'th powers together // generate F over E. bases := gammas cat [ beta + gamma : gamma in gammas ] cat [ beta ]; // print bases, " as generators gives exponents", xx; alpha := &*[ bases[i]^xx[i] : i in [1..#bases] ]; // Compute the generator as a multiplicative combination // of the bases. else // print "Neither the field nor the base field are small."; c := []; // Find n+1 distinct elements in E for x in E do Append( ~c, x ); if #c ge n2+1 then break; end if; end for; bases := [ beta + cc : cc in c ]; xx := MultiplicativePrimitiveElement( [ bb^n2 : bb in bases ], E ); // print c, bases, "as generators gives exponents", xx; // Compute a generator of F over E that is a multiplicative // combination of the n'th powers of (beta + c_i), where // c_i runs through a list of n+1 distinct elements of E // computed above. alpha := &*[ bases[i]^xx[i] : i in [1..#bases] ]; // Compute the generator as a multiplicative combination // of the bases. end if; return _FFConvertPower( F, alpha, n ), e; // Converting the power is necessary as the computation was done with // n2'th powers, where n2 = gcd( n, #F-1 ). end function; SumOfIntegerPowers := function( a, n ) // Reduce the size of a by writing it as a sum of integer n'th powers, // plus some small part that is not usefully written as such a sum. // More details: we return a sequence of integers r and an integer a' s.t. // (1) a = sum [ rr^n : rr in r ] + a' // (2) a' le n^n r := []; while a gt n^n do if n eq 1 then // Stupid!! Magma's Iroot function does not allow rt := a; // taking roots of exponent 1. Why is this ?! else rt := Iroot( a, n ); end if; Append( ~r, rt ); a -:= rt^n; end while; return r, a; end function; ExpandOnRationalBase := function( a, t, s ) // Expand the integer a on the rational base t/s. // See my thesis for more comments on what this does. // If the input is right, then afterwards, we have // a = c[1] + c[2]*(t/s) + ... + c[d]*(t/s)^(d-1) where d=#c // a should be a nonnegative integer, and t and s should be integers with // 0 < s < t. c := []; while a gt 0 do q,r := Quotrem( a, t ); // a = qt + r, with 0 le r le t-1 Append( ~c, r ); // r is the next digit a := q*s; // decrease a (note that s; if not( a in PowerField ) then // a is not a sum of n'th powers in F return false, _; end if; coeffs := ElementToSequence( PowerField!a ); if verbose then print "The coefficients of", a, "relative to the power basis are", coeffs; end if; // so a = sum [ coeffs[i]*(alpha^n2)^(i-1) : i in [1..deg] ] if #Fp gt n2+1 then // Big base field, apply some coefficient reduction // Precomputations for coefficient reduction. Only useful if #Fp > n+1. ix,jx,beta := SelectiveRoot( [ Fp!i : i in [ 1..n2+1 ] ], n2 ); beta := 1/beta; if verbose then print "In F", Characteristic( Fp ), "we have", jx, "/", ix, "=", beta, "^", n2; end if; // We know now that jx/ix = beta^n in Fp. powerlist := []; for i in [1..deg] do sumpow, rest := SumOfIntegerPowers( Integers()!coeffs[i], n2 ); expan := ExpandOnRationalBase( rest, jx, ix ); if verbose then print "The expansion of", coeffs[i], "minus the powers of", sumpow, "on the base", jx, "/", ix, "is", expan; end if; // So, coeffs[i] = sum [ sp^n : sp in sumpow ] + // sum [ expan[i]*(jx/ix)^(i-1) : i in [1..#expan] ]. powerlist cat:= [ sp*alpha^(i-1) : sp in sumpow ]; powerlist cat:= &cat[ [ alpha^(i-1)*beta^(j-1) : k in [1..expan[j]] ] : j in [1..#expan] ]; // So, for each term expan[j]*(jx/ix)^(j-1), we add expan[j] terms // alpha^(i-1)*beta^(j-1) to our list of powers. end for; else // Small base field, take coefficients as such powerlist := []; for i in [1..deg] do powerlist cat:= [ alpha^(i-1) : k in [1..Integers()!coeffs[i]] ]; end for; end if; // Now, we have a = &+[ x^n2 : x in powerlist ] !! if verbose then print a, "is equal to the sum of the", n2, "th powers of", powerlist; print "We thus found a sum of", #powerlist, "powers. Reducing..."; end if; while #powerlist gt n2 do // Reduce length of sum powerlist := _FFPowerSumReduce( F, powerlist, n2, verbose ); if verbose then print "Number of terms dropped to", #powerlist; end if; end while; truepowerlist := [ _FFConvertPower( F, pw, n ) : pw in powerlist ]; return true, truepowerlist; end function; ZeroPowerSum := function( F, n ) // Find a set of nonzero elements of F such that their n'th powers sum to // zero. // Note that this problem is always solvable over the prime field already. // In my thesis, I give an asymptotically faster way for doing this, but // it's a bit involved to implement, and this works also. _, psp := PowerSumRepresentation( PrimeField(F)!(-1), n ); return [ F!1 ] cat psp; end function; ZeroOfDiagonalForm := function( a, n : verbose := false ) // Find a nontrivial zero to the diagonal form given by // sum [ a[i]*x[i]^n : i in [1..#a] ] // over the finite field F = Universe( a ). // We should have #a gt n. // Uses the trapezium method as defined in my thesis. if #a le n then error "The number of variables should exceed the degree."; end if; F := Universe( a ); if Type( F ) ne FldFin then error "Expecting a sequence of finite field elements."; end if; // Here as well, we could reduce to the case where n divides #F-1 // Didn't implement that yet if 0 in a[ 1..n+1 ] then error "Expecting nonzero finite field elements"; end if; zs := ZeroPowerSum( F, n ); // Sequence whose n'th powers sum to zero if verbose then print "Found zero sequence", zs; end if; zspartsums := [ zs[2]^n ]; // Make sure no prefix of zs works as well for i in [2..#zs-1] do // Actually, not prefix, but indices 2..k zspartsums[i] := zspartsums[i-1] + zs[i+1]^n; end for; zeroset := [ k : k in [1..#zs-1] | zspartsums[k] eq 0 ]; if zeroset ne [] then k := Maximum( zeroset ); zs := zs[k+2..#zs]; if verbose then print "Restricting zero sequence to", zs; end if; end if; zerolengths := []; x := IdentityMatrix( F, n+1 ); S := []; for i in [1..n+1] do // I'd rather use indices 0..n, but that's not allowed zerolengths[i] := #zs; x[i][i] := zs[1]; S[i] := a[i]*zs[1]^n; end for; while Minimum( zerolengths ) gt 1 do // No sequence exhausted yet? if verbose then print "The sums are", S; print "The variables are", x; print "The progress is", zerolengths; end if; k, l, beta := SelectiveRoot( S, n ); // So k := FiniteField( 5, 12 ); Fp := PrimeField(F); // F_5 Fbla := sub< F | alpha^7813 >; // 2*7813 = (5^12-1) div (5^6-1) // Fbla is isomorphic to F, but its generator Fbla.1 is not primitive. // For example, Fbla.1^2 does not generate F over Fp. PowerGenerator( Fbla, Fp, 2 ); // Answer: Fbla.1 + 1. Nice, isn't it? // This thing is almost primitive, it generates a subgroup of index 2 in F*. Fbla!alpha; // Some complicated expression. As the generator is not primitive, the elements // of Fbla are represented on the power basis generated by Fbla.1, instead of // by powers of the generator, which is done for alpha. PowerSumRepresentation( Fbla!alpha, 2 ); // Write alpha (considered as an element of Fbla) as a sum of squares in Fbla. // If you go through the whole example with 5 replaced by 3, some differences // arise. Of course, the 7813 doesn't work; in fact, (3^12-1)/(3^6-1) is also // even, so we can replace 7813 by 365 to get the same effect. // Furthermore, when computing a representation of elements as sums of squares // in F, we see that the base field (F_3) has become so small that reduction // of coefficients is not necessary --- in fact, it doesn't work anymore. // Other nice example: pri := NextPrime( 12345678 ); F := FiniteField( pri, 6 ); Factorisation( #F-1 ); PowerSumRepresentation( alpha, 2 : verbose := true ); PowerSumRepresentation( alpha, 3 : verbose := true ); // and the same for exponent 5, 7, 13, 37, 43, 59, and some bigger primes. // Note that the order of 2 in this field is tiny: it's 5x37x333667! // So in particular, 2 is a square, a cube, a seventh power, and so on. // That's very handy for coefficient reduction! [ ZeroPowerSum( F, i ) : i in [1..10] ]; */