Factorizing the polynomial $x^6-2x^3+8$ over finite fields This is a very specific question. This equation came up when attempting to compute certain subgroups of groups of Lie type. The polynomial $x^6-2x^3+8$ splits over $\mathbb{F}_q$ if and only if a particular subgroup ($\mathrm{SL}_2(8).3$) embeds in a particular group of Lie type ($E_6(q)$).
More specifically, let $q$ be a prime power congruent to $1,2,4\bmod 7$, and not a power of $2$ or $3$. We ask for solutions to the equation $x^6-2x^3+8=0$ over $\mathbb{F}_q$. If $q\equiv 2\bmod 3$ then there is always a solution, but the polynomial does not split. Thus we may in addition assume that $q\equiv 1\bmod 3$. This polynomial has a solution if and only if $1+\sqrt{-7}$ is a cube in $\mathbb{F}_q$.
If it is of help, notice that $1+\sqrt{-7}$ is the sum $2(\zeta+\zeta^2+\zeta^4)$, where $\zeta$ is a primitive $7$th root of unity.
Ideally, one would like a statement of the form 'This equation has a solution if and only if $q$ is congruent to one of [...] modulo $n$'.
 A: A comment asks for emperical data. I do not have a solution but in a search for inspiration I computed the first bunch of prime powers with $1 + \sqrt{-7}$ a cube (under the conditions $q \equiv 1 \mod 3$)

The $q \equiv 1 \mod 7$
169, 211, 421, 673, 841, 967, 1051, 1093, 2017, 2143, 2269, 2311, 2647, 2731


The $q \equiv 2 \mod 7$
121, 457, 751, 961, 1129, 1549, 2053, 2221, 2557, 2809, 3019, 3187, 3271


The $q \equiv 4 \mod 7$
67, 361, 487, 529, 613, 739, 1117, 1201, 1453, 1663, 1747, 1999, 2671, 2713

The Magma is below for those who wish to play (apologies for ad-hoc awfulness):
cubes_in_class_mod_7 := function(c)
    our_cubes := [* *];
    for p in [2..10000] do
        if IsPrimePower(p) then
            if (p mod 2) eq 1 then
                if ((p mod 7) eq c) and ((p mod 3) eq 1) then
                    Fp := FiniteField(p);
                    F<x> := PolynomialRing(Fp);
                    a := Roots(x^2 + 7)[1][1];
                    is_cube := Roots(x^3 - (1 + a));
                    b := #is_cube;
                    if b eq 3 then
                       our_cubes := Append(our_cubes, p);
                    end if;
                end if;
            end if;
        end if;
    end for;

    return our_cubes;
end function;

cubes_in_class_mod_7(1);
cubes_in_class_mod_7(2);
cubes_in_class_mod_7(4);

A: I am a little surprised to see that my program for finding the positive primes represented by an indefinite binary quadratic form is, well, self contained. It accepts a triple of integers that are the coefficients of the original form, then applies Gauss-Lagrange reduction. Once the form is in the correct shape, it calculates large portions of the J. H. Conway Topograph. Conway's method is well suited to find just positive terms. Here is the C++ program. Note that the discriminant $b^2 - 4ac$ which must be positive but not a square, for this program must be smaller than $2^{31}-1 = 2147483647.$  I remember changing the IntSqrt function because an infinite loop occurred when finding the square root of a number too close to $2^{31}.$ From what I see, it is best to take the coefficients $|a|, |b|, |c| < 10000.$ It would not be much of a problem to call in the GMP libraries (for C++) and change all integer variables to mpz_class.
#include <iostream>
#include <stdlib.h>
#include <fstream>
#include <sstream>
#include <list>
#include <set>
#include <math.h>
#include <iomanip>
#include <string>
#include <algorithm>
#include <iterator>
#include <strstream>


using namespace std;




//  this file is named     Conway_Positive_Primes.cc 
//  main writing   June 2014
  
 
// to compile
//        g++  -o     Conway_Positive_Primes     Conway_Positive_Primes.cc  -lm




void insert_primitive_reps(unsigned int a, unsigned int h, unsigned int b, unsigned int M,  set<unsigned int>  *setPtr)
{
 // cout << setw(12) << a  << setw(12) << h  << setw(12) << b << "   insert_primitive_reps"    << endl;
  if ( a <= M )
  {
    (*setPtr).insert(a);
    if ( b <= M )
    {
      (*setPtr).insert(b);
      if ( a <= M - b && h <= M - a - b)
      {
        if( a <= M - a - h ) insert_primitive_reps(a, h + 2 * a, a + b + h, M, setPtr);
        if( b <= M - b - h ) insert_primitive_reps(a + b + h, h + 2 * b,b, M, setPtr);
         // comment: once a+b+h <= M, min(2a+h, 2b+h) <= M
      }  // if a + b + h
    } // if  b
  } // if a

} // end insert_primitive_rep



int IntSqrt(int i)
{
  // cerr << "called intsqrt  with   " << i << endl;
  if ( i <= 0 ) return 0;
  else if ( i <= 3) return 1;
  else if ( i >= 2147395600) return 46340;
  else
  {
    float r;
    r = 1.0 * i;
    r = sqrt(r);
    int root = (int) ceil(r);
    while ( root * root   <= i ) ++root;
    while ( root * root   > i ) --root;
    return  root ;
  }
}


string stringify(int x)
 {
   ostringstream o;
   o << x  ;
   return o.str();
 }


string Factored(unsigned int i)
{
  string fac;
  fac = " = ";
  int p = 2;
 unsigned int temp = i;
  if (temp < 0 )
  {
    temp *= -1;
    fac += " -1 * ";
  }

  if ( 1 == temp) fac += " 1 ";
  if ( temp > 1)
  {
    int primefac = 0;
    while( temp > 1 && p * p <= temp)
    {
      if (temp % p == 0)
      {
        ++primefac;
        if (primefac > 1) fac += " * ";
         fac += stringify( p) ;
        temp /= p;
        int exponent = 1;
        while (temp % p == 0)
        {
          temp /= p;
          ++exponent;
        } // while p is fac
        if ( exponent > 1)
        {
          fac += "^" ;
          fac += stringify( exponent) ;
        }
      }  // if p is factor
      ++p;
    } // while p
    if (temp > 1 && primefac >= 1) fac += " * ";
    if (temp > 1 ) fac += stringify( temp)  ;
  } // temp > 1
  return fac;
} // Factored


int PrimeQ(int i)
{
  if ( i < 0 ) i *= -1;
  if ( i <= 3) return 1;
  else
  {
    int boo = 1;
    int j = 2;

    while (boo && j * j   <= i )
    {
      if ( i %  j  == 0) boo = 0;
      ++j;
    }
    return boo;
  }
}


int main(int argc, char *argv[])
{
  if ( argc != 6) cout << "Usage: ./Conway_Positive_Primes A B C      Bound    Modulus " << endl;
  else {
 


  int a,b,c, discr;
     int N,M;
    a = atoi(argv[1]);
        b = atoi(argv[2]);
    c = atoi(argv[3]);
   M = atoi(argv[4]);
     N = atoi(argv[5]);
 cout << setw(12) << atoi(argv[1])  << setw(12) << atoi(argv[2])  << setw(12) << atoi(argv[3]) << "   original form " << endl    << endl;
      int d = b * b - 4 * a * c ;
      int droot = IntSqrt(d) ;
      if ( d <= 0) cout << "nonpositive discriminant  " << d << endl << endl;
       if (d == droot * droot) cout << "square discriminant  " << d << endl << endl;
 
      if (d > 0 && d != droot * droot)
      {



      int aa,bb,cc;
      while ( a <= 0 || c >= 0 || b <= abs(a+c) )
      {
        int delta, cAbs;
       cAbs = c;
        if (cAbs < 0) cAbs *= -1;

       delta =   (b + droot)/( 2 * cAbs)  ;
  if (c < 0) delta *= -1;
     aa = c ; bb = 2 * c * delta - b ; cc =  c * delta * delta - b * delta + a ;
       a = aa; b = bb; c = cc;
       } // while not reduced with a > 0


  cout << setw(12) << a  << setw(12) << b  << setw(12) << c << "   Lagrange-Gauss reduced " << endl    << endl;


        int a_old = a;
        int b_old = b;
        int c_old = c;

        int goon = 1;
       

set<unsigned int>  S;





        while (goon )
        {
         
      
         int newval = a+b+c;

          if ( newval > 0 )
          {
       //      cout << setw(65) << b + 2 * a << endl;
             insert_primitive_reps(a, b + 2 * a,newval ,M,  &S); // note ampersand
             b+= 2 * c ;
            
             a = newval;

           } // newval > 0
         else if ( newval < 0 )
          {
     //        cout << setw(5) << -1 * ( b + 2 * c) << endl;
 
             b+= 2 * a ;
             c = newval;
 
           } // newval < 0

  
         goon = (a != a_old)  || (b != b_old)  || (c != c_old)  ;

        } // while goon

     
  cout << endl << endl << " Represented (positive) primes up to  " << M << endl << endl;
 set<unsigned int> mods ;

 int rount = 0;
 set<unsigned int>::iterator iterU;
 for(iterU = S.begin() ;   iterU != S.end() ; ++iterU)
    {
      
      unsigned int p = *iterU;
     if ( p > 1 && PrimeQ(p) )
      {
       cout << setw(6) << p ;
       ++rount;
       if (rount % 10 == 0 ) cout << endl;
     // cout << setw(12) << p << setw(12) << p %  N << endl ;
      mods.insert( p % N);
      

      } // if prime
    }

  int count = 0;
  cout << endl << endl;
  cout << "=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=   " << endl;
  cout<< " these are the collection of remainders when dividing by   " << N << endl << endl;
 for(iterU = mods.begin() ;   iterU != mods.end() ; ++iterU )
  {
    int    u = *iterU ;
    ++count;
    {
      cout << setw(7) << u  ;
      if (0 == count % 10) cout << endl;
    }
  } // for iterU

  cout << endl << endl;

   


  cout << endl << endl << " Represented (positive) primes up to  " << M <<  "  and value mod    " << N << endl << endl;
 cout << setw(12) << atoi(argv[1])  << setw(12) << atoi(argv[2])  << setw(12) << atoi(argv[3]) << "   original form " << endl    << endl;
       } // not square


    } // else argc
    return 0 ;
}





 
//  g++  -o     Conway_Positive_Primes     Conway_Positive_Primes.cc  -lm

