How to better plot elliptic curves over finite fields?

edit

I don't think it will be instructive, for learning purposes, to use a prime that big and plot it (assuming this is what you want to do) and I suspect the problem is the size of the field.

Here's the curve I use as an example, with cofactor 1, for "ecc" curves that are like NIST ones.

  E=EllipticCurve(GF(17),[3,5])
  E.order()
  E.is_supersingular()
  E.plot()

As you can see, it's a prime order group (so cofactor 1) and not supersingular, although it is obviously way too small for real world use. The plot however is nicer. Can't upload it to show you, but you can reproduce from this.

Let us understand first the called plot routine, asking for "help" and "full help". Using

E = EllipticCurve(GF(next_prime(20000)),[0,1])
E.plot?

we obtain also the information where the code is located:

Init docstring: Initialize self.  See help(type(self)) for accurate signature.
File:           /usr/lib/python3.8/site-packages/sage/schemes/elliptic_curves/ell_finite_field.py
Type:           method

Usually, E.plot?? would be enough to get the way the method works,

Source:   
    def plot(self, *args, **kwds):
        """
        !!! ... long doc string ... !!! PLEASE CHECK IT
        """
        R = self.base_ring()
        if not R.is_prime_field():
            raise NotImplementedError

        G = plot.Graphics()
        G += plot.points([P[0:2] for P in self.points() if not P.is_zero()], *args, **kwds)

        return G
File:      /usr/lib/python3.8/site-packages/sage/schemes/elliptic_curves/ell_finite_field.py

but in this case we need more, since the plot in the line G = plot.Graphics() is not the global plot, but the one in the lines...

import sage.plot.all as plot

class EllipticCurve_finite_field(EllipticCurve_field, HyperellipticCurve_finite_field):

at the beginning of the file. OK, we know the origin of the plot method, can also ask for its documentation:

sage: import sage.plot.all as xplot                                                                                                  
sage: ?xplot.points

and the last points method sends us to point2d, where we finally get the meaning of the options we can use, for instance

sage: xplot.point2d?

also contains examples using size and pointsize. Defaults are

sage: xplot.point2d.options                                                                                                          
{'alpha': 1,
 'aspect_ratio': 'automatic',
 'faceted': False,
 'legend_color': None,
 'legend_label': None,
 'marker': 'o',
 'markeredgecolor': None,
 'rgbcolor': (0, 0, 1),
 'size': 10}

So in order to get smaller points, one can try:

sage: E.plot(pointsize=1)

and note that the option pointsize=1 is passed to the plot through the line in the code of E.plot:

Source:   
    def plot(self, *args, **kwds):
        """etc
        """
        R = self.base_ring()
        if not R.is_prime_field():
            raise NotImplementedError

        G = plot.Graphics()
        G += plot.points([P[0:2] for P in self.points() if not P.is_zero()], *args, **kwds)

        return G

in the line with G += ... where we get the *args, **kwds from the E.plot caller, and insert them to the plot.points. To have access to the size of the graphics, we need to have access to G, so it is maybe a good idea to construct the G object with bare hands and use its G.show method with the right / needed figsize. See the example in xplot.Graphics?