(* Helper functions *)

sorteddeglist[seppoly_, prime_] := 
 Sort[Map[Exponent[#, x] &, 
   Transpose[Drop[FactorList[seppoly, Modulus -> prime], 1]][[1]]]]

(* Computes the degree lists of the factorization of the polynomial mod p *)

dropstuff[list_, targetsum_] := 
 Select[list, Total[#[[1]]] == targetsum &]

(* Drops terms from a list that don't have the right total degree *)

discprimes[poly_] := 
 Transpose[FactorInteger[Abs[Discriminant[poly, x]]]][[1]]

(* Computes the primes dividing the discriminant *)




frobcount[poly_, numprimes_] := 
 dropstuff[
  Tally[Table[
    sorteddeglist[poly, Prime[n]], {n, 1, 
     numprimes + Length[discprimes[poly]]}]], Exponent[poly, x]]

(* Usage: Write the polynomial as a polynomial in x, and give the \
total number of primes you want to have the total counted for. The \
primes dividing the discriminant are automatically dropped from the \
list.  Count may be off if you don't pick a range that includes all \
of the primes dividing the discriminant of the polynomial, just \
subtract one or two from the range if needed. *)




frobcount[x^7 - 14 x^5 + 56 x^3 - 56 x - 22, 100]