/*
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 ge 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] ];
*/