# *****************************************************************************
#                            Gauss Sums
#    1. Computing quadratic Gauss sums; checking Gauss Theorem
#    2. Plotting Gauss Sums and their argument g(c)/p^.5 ................
#
#    3. Is Gauss Sum a class function ord(c)=ord(d) => g(c)=g(d)? ....
#    4. Cubic Gauss Sums: g(c) for ord(c)=3; special case: SG-primes ..
#
# *****************************************************************************

# ****************************************************************
#          Computing Quadratic Gauss Sums
# ****************************************************************
p=7 # Pick a prime (for now)
G=DirichletGroup(p) # Group of Dirichlet characters
c=G[1]              # select a generator for characters
L=c^ZZ((p-1)/2)    # This is the only quadratic character: Legendre symbol
print("Legendre symbol:", L)
print( "L^2=", L^2)
print( (" "))
GS=L.gauss_sum()/p^(1/2); GS
P=(real_part(GS), imag_part(GS)); P

a = plot([],figsize=(3,3),title='Gauss Sums for p='+str(p),frame=True,axes_labels=['$x$-axis','$y$-axis'])
a+=circle((0.0,0.0),1,color='blue')
a+=point(P,color='red', size=20)
show(a)

print("Numeric values")
z=CC(GS); ord=2; print("order=",ord, " Arg(g(c))=", z)
print( "z^4=", z^4)
k=2*pi.n()/z.argument(); print("2 Pi/Arg = ", k, " ", floor(k))

# *********************************************************************************************
#                   Plotting the argument of Gauss Sums (unit complex numbers): g(c)/p^(1/2)
# *********************************************************************************************
P=Primes()
for k in range(1,30):
    p=P.unrank(k); # print "\nPrime p=", p
    G=DirichletGroup(p) # Group of Dirichlet characters
    c=G[1]              # select a generator
    dv=divisors(p-1); ndv=len(dv)   # possible orders m|p-1
    print( "\n**************************************************************")
    print( "Prime p=", p, "      p-1=", factor(p-1), " has ", ndv, " divisors:"), dv
    gs=[c.gauss_sum() for j in range(ndv)] # initialize the list of Gauss Sums
    arggs=[CC(0) for j in range(ndv)]      # initialize the list of arguments of Gauss sums ("normalizing" to modulus 1)
    a = plot([],figsize=(3,3),title='Gauss Sums for p='+str(p),frame=True,axes_labels=['$x$-axis','$y$-axis'])
    for j in range(ndv):
        # print "j=",j
        m=dv[j]                                     # selecting a divisor m of p-1
        # print ("Order of character:",m)
        c1=c^ZZ((p-1)/m)                            # constructing the character of order m
        gs[j]=c1.gauss_sum()
#        print ("Ord(c)=",m, " Gauss sum g(c):", gs[j])  # the Gauss Sum g(c) - too ugly to show :)
        if m==1:
            arggs[j]=gs[j]
        else:
            arggs[j]=gs[j]/p^(1/2)                     # argument of Gauss sum (unit complex number)
#        print("Normalized:", arggs[j])
#    print ("gs:", gs)
    a+=list_plot(gs,color='blue', size=20)
    for j in range(ndv):                            # plotting the orders of characters as labels for the respective Gauss Sums
        z=CC(gs[j])+0.1
        a+=text(str(dv[j]), z, color='red')
#    show(a) # For each prime        -              SKIP for now ...
    a = plot([],figsize=(4,4),title='Argument of Gauss Sums for p='+str(p),frame=True,axes_labels=['Red label = order of character',' '])
#    print "arg:", arggs
    a+=list_plot(arggs,color='red', size=20, xmin=-1.1,xmax=1.1,ymin=-1.1,ymax=1.1)
    a+=circle((0.0,0.0),1,color='blue')
    for j in range(ndv):                        # adding the orders of characters as labels for the respective Gauss Sums
        z=CC(arggs[j])
        a+=text(str(dv[j]), z*CC(1.1), color='red')     # shifting the label radially
    show(a) # For each prime

# Questions / variations
# 1. Test quadratic GS - formula
# 2. Test cubic GS
# 3. Test SG-primes
# 4. Test divisibility of GSs