#########################################
# Congruence arithmetic - basics;
# see also http://userpages.umbc.edu/~rcampbel/Computers/SAGE/numbthy.html
# and http://www.math.ucla.edu/~jimc/mathnet_d/sage/reference/sage/rings/finite_rings/integer_mod_ring.html
# or http://doc.sagemath.org/html/en/reference/rings_standard/sage/rings/finite_rings/integer_mod_ring.html
#########################################
# Part I

# Defining the ring Z/NZ 
R=Zmod(20); R # Z/20Z
R.order()
R.list()

# Generators
R.gens()
R.gen(0) # Why not all Z/nZ*?

# Specifying elements
x=R(5); print "Converting the integer 5 into 5 mod 20=", x

a=Mod(7,20); a # a modulo n as element of Z/nZ: Mod(a, n)

# Addition
a+ 19 # 7+15 mod 20=2

# Multiplicative group for R=Z/20Z
U=R.list_of_elements_of_multiplicative_group() # The group of units U=(Z/nZ)*
print "Z/20Z*= ", U, " has ", len(U), " elements"
U[0] # is always 1 
# Accesing individual units, e.g. the 3rd unit
U[2]
G=R.unit_gens(); G                  # multiplicative generators, i.e. Primitive roots
print "No. of generators (rank/dimension):", len(G)
m=1
for k in range(len(G)):
    o=multiplicative_order(G[k])
    print "gen.", k, "is ", G[k], " additive order is ", G[k].order(), " mult. order:", o
    m=m*o
# print m, len(G), len(U)
if m==len(U): print "No of elements", len(U), " equals prod. of orders of gen.", m
# How to access the group structure?

# The multiplicative orbit of generators for (Z/20Z)*
for i in range(len(G)):
    g=G[i]; [g^k for k in range(g.multiplicative_order())]

# Multiplication in Z/nZ
R=IntegerModRing(20) # also denoted Zmod(20) or Integers(20)
if R==Integers(20): print "IntegerModRing(n) == Integers(n)"
# It is a list of 20 elements
R.list()
# Defining an element
R(5) # converting integer 5 to 5 mod 20; NOT R[5] ...
# ... but labeling the list allows to access it:
L=R.list(); print L[0], L[7] # L[20] is not defined!

# Multiplication mod 20
a=21 % 20; b=7 % 20; c=a*b # multiplication mod 20
print a,b,c
if c==R(7): print "TRUE" # equal as elements in Z/20Z

# inverse of n (mod m): n.inverse mod(m)
# power a
# n (mod m): power mod(a, n, m)

# square root of a (mod n) =Mod(a,n).sqrt()
Mod(1,7).sqrt()
Mod(6,7).sqrt() # 6, i.e. -1, is not a quadratic residue mod 7 ... see below

# **************************************************************************************
#           Part II:        More Number Theory - functions & multiplication
# From http://wiki.sagemath.org/quickref?action=AttachFile&do=get&target=quickref-nt.pdf
#

# Ring Z/nZ = Zmod(n) = IntegerModRing(n)

# Euler’s φ(n) function: 
euler_phi(6) # equals the number of elements of multiplicative group: |Z/nZ|*
Zmod(6).list_of_elements_of_multiplicative_group()

# Kronecker symbol kronecker_symbol(a,b)
kronecker_symbol(2,7) # QR
Mod(2,7).sqrt() # sqrt is 3
print "Is 3^2=2? Yes: ", Mod(2,7).sqrt()^2 # checking 3^2=2

kronecker_symbol(6,7) # not a QR
print "has no square root: ", Mod(6,7).sqrt()

# Quadratic residues: 
quadratic_residues(7)

# Quadratic non-residues: 
[k for k in range(7) if kronecker(k,7)==-1]

# primitive root modulo n
g=primitive_root(7); g
g.multiplicative_order() 
# To make g an element of Z/7Z

g1=mod(g,7); g1
o=g1.multiplicative_order(); o
# so ... order of a (mod n) =Mod(a,n).multiplicative order()
# checking that it generates Z/7Z* ...
[3^k % 7 for k in range(o)]

# discrete log: log(Mod(6,7), Mod(3,7))
b=Mod(3,7); a=Mod(6,7); x=log(Mod(6,7), Mod(3,7))
print "a=",a," b=",b, " x=",x, "b^x=a i.e. x=log_b(a): ", b^x

# Chinese remainder theorem: x = crt(a,b,m,n)
#    finds x with x ≡ a (mod m) and x ≡ b (mod n)
a=1; m=2; b=3; n=5
print a,m,b,n
x=crt(a,b,m,n); x 

# Checking
x % m # is a=1
x % n # is b=3