Dear Michael,

Here are a few remarks, which are just to continue our conversation.

The following MAGMA program is meant to illustrate our idea to use
densities to determine (with some unknown probability), whether or not
a mod-p representation is surjective.  The first intrinsic computes an
approximate density for the number of p with a given characteristic
polynomial.  The second computes, in a stupid way, how many elements
of GF(p) have a given charpoly.   There's some magma code at the bottom
that illustrates usage of this code.   The examples aren't interesting
because they are semistable.

Do you know an example of a pair (E,p), where E is an elliptic curve,
E[p] is not surjective, and there is no curve in the isogeny class of
E that has a rational p-torsion point?  There must be interesting
papers about examples of irreducible reps like that.  I vaguely recall
that Iman Chen might have written one of them.

// --- START ---------------------------------------------------------

intrinsic DensityWithCharpoly(E::CrvEll, p::RngIntElt,
                  f::RngUPolElt, ell_max::RngIntElt) -> FldRatElt 
{The density of primes ell <= ell_max such that the charpoly
of Frob_ell on E[p] is f.}

   R := PolynomialRing(GF(p));
   f := R!f;
   a := -Coefficient(f,1);
   det := Coefficient(f,0);
   i := 0;
   j := 0;
   ell := 2;
   for a_ell in TracesOfFrobenius(E,ell_max) do
      if a_ell eq a and (ell mod p) eq det then
         i := i + 1;
      end if;
      ell := NextPrime(ell);
      j := j+1;
   end for;
   return i/j;
end intrinsic;

intrinsic DensityIfSurjective(p::RngIntElt, f::RngUPolElt) 
  -> FldRatElt
{The density of elements z in GL(2,GF(p)) such that the
characteristic polynomial of z is f.}

   R := PolynomialRing(GF(p));
   f := R!f;
   G := GL(2,p);
   a := [z : z in G | CharacteristicPolynomial(z) eq f];
   return #a / #G;

end intrinsic;

// --- END ---------------------------------------------------------


> >> Attach("surj.m");
> E := EC("37A");
> DensityWithCharpoly(E,5,(x-1)^2, 100);
1/25
> DensityWithCharpoly(E,5,(x-1)^2, 100)*1.0;
0.0399999999999999999999999999994
> DensityWithCharpoly(E,5,(x-1)^2, 1000)*1.0;
0.065476190476190476190476190475
> DensityWithCharpoly(E,5,(x-1)^2, 10000)*1.0;
0.0480065093572009764035801464597
> DensityWithCharpoly(E,5,(x-1)^2, 100000)*1.0;
0.050458715596330275229357798164
> E := EC("11A");
> DensityWithCharpoly(E,5,(x-1)^2, 10000)*1.0;
0.248169243287225386493083807973
> DensityWithCharpoly(E,5,(x-1)^2, 50000)*1.0;
0.248003117085525034093122930058
> DensityIfSurjective(3,(x-1)^2)*1.0;
3/16
> DensityIfSurjective(3,(x-1)^2)*1.0;
0.187500000000000000000000000000
> E := EC("37A");
> DensityWithCharpoly(E,3,(x-1)^2, 10000)*1.0;
213/1229
> DensityWithCharpoly(E,3,(x-1)^2, 10000)*1.0;
0.173311635475996745321399511795
> E := EC("37B");
> DensityWithCharpoly(E,3,(x-1)^2, 10000)*1.0;
0.496338486574450772986167615947
>