Revision history [back]

This answer sort of cheats, but:

sage: p = 17
sage: quadratic_residues(p)
[0, 1, 2, 4, 8, 9, 13, 15, 16]

This is, after all, a well-known concept.

This works for your example too, even though 8 isn't prime

sage: quadratic_residues?
Signature:      quadratic_residues(n)
Docstring:
   Return a sorted list of all squares modulo the integer n in the
   range 0<= x < |n|.