###### Sage code to visualize the proof that a^p = a (mod p) if p is a prime #####
# if you draw a^p necklaces then they can be placed in groups
# where the necklaces are the same if they are equal up to rotation 
# the number of elements in a group will be p 
# UNLESS the necklace is a solid color (because then it is in a group by itself).
# this shows:
# a^p = p*(number of groups of size p) + (number of solid colored necklaces)
# therefore a^p = a (mod p) because there are "a" solid colored necklaces
#############################################################################
# you can copy and paste this code in
# https://sagecell.sagemath.org
# to see it work
#
# Sage may refuse to draw the picture if a and p are too big.
# but try to explain what goes wrong if we instead use a number of beads
# not equal to a prime number.  For example:
# 3^4 = 1 (mod 4) so it is not the case in general that a^n = a (mod n)
# the code should work if you execute arrange_necklaces(3,4) as well.

def necklace(x,y,w):
    """
    Draw a single necklace of radius 1/2 at the point (x,y) using the colors
    specified by the word w

    INPUT:

    - ``x, y`` -- coordinates of the center of the necklace
    - ``w`` -- a list of positive integers (e.g. [1,1,2,3])

    OUPUT:

    - a graphics object
    """
    return circle((x,y),1/2,color='black',axes=False)\
        +sum(circle((x+cos(2*pi*r/len(w))/2,y+sin(2*pi*r/len(w))/2)\
        ,1/8,fill=True,color=clr[w[r]]) for r in range(len(w)))

# clr is a list of colors of beads to use in the necklaces
clr = [colors['white'], colors['black'], colors['red'], colors['green'],\
        colors['blue'],colors['orange'], colors['yellow'], colors['cyan'],\
        colors['magenta'], colors['brown'], colors['grey']]


def mcs(w):
    """
    return the minimal cyclic shift of w

    INPUT:

    - ``w`` -- a word or a list

    OUPUT:

    - a list
    """
    if len(uniq(w))==1:
        return [100]*len(w)
    return min(list(w[i:]+w[:i]) for i in range(len(w)))

def arrange_necklaces(a,n):
    """
    Arrange all necklaces of length ``n`` with ``a`` colors of beads
    in a square and draw them

    INPUT:

    - ``a`` and ``n`` -- integers

    OUTPUT:

    - a graphics object
    """
    d = n*ceil(sqrt(a^n)/n)                                                                                                                                      
    wds = map(lambda x: x[1], sorted([[mcs(w),w] for w in  Words(a,n)]))
    return sum(necklace(2*int(mod(r,d)),2*int(r//d),wds[r])\
        for r in range(len(wds)))

arrange_necklaces(2,5)