restart
-- macaulay2 can be used as a calculator

3/5 + 7/11
4*17*31
4^3
96!
sqrt 2
2^(4/3)

-- more interestingly, macaulay2 is a computer algebra system that
-- was created to support research in commutative algebra and 
-- algebraic geometry. 

-- with this purpose, it has been designed allow the construction and
-- manipulation of many sorts of algebraic objects.  

------------
-- RINGS ---
------------

-- to start the first one does is declare the  ambient ring where 
-- the magic happens

-- examples:
R=ZZ --ring of integers
R=ZZ[x] --polynomial ring over the integers
R=(ZZ/32003)[x,y] --polynomial ring over finite field
R=GF(2,3)[z] --polynomial ring over finite field finite field of order 2^3
R=QQ[w,x,y,z] --polynomial ring over the rationals

-- By convention, all calculations that follow are assumed to occur
-- in the last ring declared, which is taken as the ambient ring
(x+y+z)^4
factor oo
1+b
-- here b hasn't been identified as a variable, so the sum doesn't make sense

-------------
-- IDEALS ---
-------------

-- also the subsequent ideals are defined in the current ambient ring

I = ideal (w*x-y^2,z^3)
J = ideal (w^3,x*y-z^2)

-- macaulay2 can do many operations on ideals
I+J -- sum
I*J -- product
I^3 -- powers of an ideal
isPrime I --test for primeness
isPrimary I^2 -- test for being primary
radical I -- radical
J:ideal(w^2) -- colon ideals
IcapJ=intersect(I,J) -- intersections
isSubset(I,J) -- test ideal containment
isSubset(I,I+J)
primaryDecomposition IcapJ -- primary decompositions

--------------------
-- QUOTIENT RINGS --
--------------------

-- a quotient ring is defined as one would expect
S=R/I

-- since S was the last ring defined, it is replaced as the current
-- ambient ring: e.g. any polynomial we write is taken modulo I
-- example:
{x, w*x, z^2, z^3}

-- macaulay2 can obtain the generators of a ring 
-- i.e. as algebra over field:
gens S
-- this are the same as the generators of the unquotiented ring
-- however, they are subject to the relations in I
gens R

-- any further object defined now will be considered over the 
-- quotient ring S=R/I:
i=ideal(x-y)

-------------
-- MODULES --
-------------

-- let' start anew
restart
R=QQ[x,y,z,w]
I=ideal(x*y^2-z^3,y^3-x^2*z, x*y*z*w)

-- we can have macaulay2 view R/I as an R-module, thus as an 
-- R-algebra:
N=R^1/I
generators N -- obviously N is generated as an R-algebra by 1..
relations N -- and its elements are subject to the relations in I

-- an advanced feature of M2 is computation of free resolutions:
FR = res N 
FR.dd -- along with all their decorations...

-- here's how a free module over R looks like
M=R^3 
generators M -- this should be seen as R*e_1 + R*e_2 + R*e_3

-- of course, there are more general types of modules besides 
-- free modules. 
-- the first is quotients of a free module, which can be presented
-- as the cokernel of a module homomorphism
A=matrix{{x,y,z},{y^2,z^2,w^2}}
M=coker A
generators M
relations M
presentation M
D=res M
D.dd

-- the second is (finitely generated) submodules of a free module,
-- which can be presented as the image of a homomorphism:
B=matrix{{x,y,z,w},{x^2,y^2,z^2,w^2}}
N=ker B 
generators N
relations N
presentation N
E=res N
E.dd

-- incorporating both classes of modules, we have subquotients,
-- the main data structure for representing arbitrary modules
-- in macaulay2. these are defined as (finitely generated) 
-- submodules of quotient modules
M
Mod=ideal(x,y)*M
Mod==x*M+y*M
generators Mod
relations Mod
presentation Mod
F=res Mod
F.dd

--------------------------------------
-- EXAMPLES FROM THE EXERCISES -------
--------------------------------------

--SHEET1

-- we can use macaulay2 to find relations among generators of
-- invariant rings. In sheet 1, for instance, we saw that the 
-- invariant ring of QQ[x,y] under the ZZ/2 action (x,y)->(-x,-y)
-- is QQ[x^2,x*y,y^2], which can be retrieved as a quotient ring 
-- as follows:
restart
R=QQ[x,y]
S=QQ[u,v,w]
f=map(R,S,{x^2,x*y,y^2})
I=ker f
FR=res I
FR.dd

--SHEET2

-- remember we proved in sheet2 that the preimage in R of the 
-- localization of an ideal I at U equals the sum of the ideals
-- J:f^infty as f ranges over f. let's apply this for 
-- U={1,f,f^2,...}

restart
R=QQ[x,y,z]
I=ideal(x^4-y^5, x^3-y^7)
-- to compute the "saturation" we compute successive quotients until
-- no newer ideals appear
If1=I:x
If2=I:x^2
If3=I:x^3
If4=I:x^4
-- let's compare the ideals so far
I==If1
If1==If2
If2==If3
If3==If4 -- this says I:x^infty=I:x^3
-- So I:x^3 is the preimage of the localization of I at U={1,x,x^2,...}
-- In this case, we can check that by alternatively working in the 
-- larger polynomial ring
R1=QQ[x,y,z,t]
I=ideal(x^4-y^5, x^3-y^7)
J=I+ideal(1-x*t)
Jp=eliminate(t,J) -- this is the preimage
Jp==If3 -- OOPS, n order to compare Jp to the ideal above, we have to move
-- it to the polynomial ring R=QQ[x,y,z]
sub(Jp,R)==If3 

-- M2 can be used to compute primary decompositions. for our 
-- preimage we have
primaryDecomposition If3 / print

-- comparing with the original ideal, we see that the preimage
-- of the localization got rid of the component (x^3,y^5) of I. 
primaryDecomposition I /print
-- all polynomials in this component give functions that vanish 
-- on x=0, which we wanted to exclude

-- one can also have macaulay2 compute Hom-modules, here are the 
-- solutions to the last exercise of sheet 2:
restart
R=ZZ
M=R^1/ideal(3)
N=R^1/ideal(12)
H=Hom(M,N)
minimalPresentation H

M=R^1/ideal(63)
N=R^1/ideal(28)
H=Hom(M,N)
minimalPresentation H

M=R^1/ideal(3)
N=R^1/ideal(7)
H=Hom(M,N)
minimalPresentation H 
 
R=QQ[x]
M=R^1/ideal(x^7)
N=R^1/ideal(x^12)
H=Hom(M,N)
minimalPresentation H
H=Hom(N,M) -- here any element of R/(x^7) can be chosen for the image of 1
minimalPresentation H

-- here's the code to compute tensor products
restart
R=QQ[x]
S=QQ[x]
M=R**S
describe M
gens M 
-- macaulay2 doesn't know it should represent the generators of
-- this tensor product with different symbols, although it handles
-- the variables as different:
M_0==x
M_1==x
random({2,0},M)
random({0,2},M)
random({1,1},M)

-- however, we can force M2 to use different symbols:
N=tensor(R,S,Variables=>{a,b})
describe N
gens N
random({2,2},N)

--SHEET3

-- as mentioned before, M2 can compute free resolutions.
-- here's how the first exercise in sheet3 would look like:
restart 
R=QQ[x,y]
m=ideal(x,y)
k=R^1/m
FR=res k
FR.dd

-- let's close with more examples of primary decomposition

--EXAMPLE0

-- let's try to understand what colon ideals mean, and what 
-- their role in primary decomposition is

restart
R=QQ[x,y,z]
I=ideal(x*z-y^2, x^3-y*z)
-- the two polynomials below are a gift from heaven. 
-- let's compute the colon ideals
J1=I:(x^2*y-z^2)
J2=I:(x)
-- we compare the result with a primary decomposition:
primaryDecomposition I / print
-- observe that it consists of the colon ideals we just 
-- computed. this means that the effect of coloning is 
-- "getting rid" of certain ideals in the intersection
-- representing I

--QUIZ: Let V,W be algebraic sets in k^n, where k is 
--algebraically closed. prove that 
--     I(V):I(W)=I(V-W)
-- hint: first prove that for general I,J ideals in k[x1..xn]
--    I:J is contained in I(Z(I) - Z(J))
-- you are free to assume that if x is not in Z(I), then 
-- there is some f in I with f(x) different from 0


--EXAMPLE1
restart
R=QQ[a..d]
I = ideal(a*b + c*d, a*c + b*d, a*d + b*c)
ass(I)
minimalPrimes I
-- we can check that I is radical as follows
I== intersect(ass(I))

-- we can check that the associated primes are of the form 
-- P=I:x for some x. for instance
I:(c^2 - d^2)*d
I:(a*d^2 + b*d^2 - c*d^2 + c^2*d)

--EXAMPLE2

R=QQ[a..j]
A=genericMatrix(R,2,5) 
Iall = minors(2,A)
-- we know that the ideal defined by the 2by2 minors is prime:
-- indeed, this is the ideal of polynomials vanishing on the 
-- 2 by 5 matrices of rank one
isPrime Iall

-- in a typical application, we would be given an ideal generated 
-- by a subset of these generators, say minors spanned by adjacent
-- columns. 

-- For instance, we would like to know which of polynomial
-- functions that vanish on the set of rank-1 matrices can be 
-- written in terms of the adjacent minors

Iadj = ideal (a*d - c*b, c*f - d*e, e*h - g*f, g*j - i*h)
isPrime Iadj -- not surprisingly, the ideal will not be prime
radical Iadj == Iadj

-- the answer to our question can be deducted from 
-- a primary decomposition
L = primaryDecomposition Iadj / print

-- note that the first component is the ideal generated by all 
-- the minors of A. 

-- from the 4th component, we see that a necessary condition 
-- is that each term of the polynomial function has variables
-- h,g,d,c

-- A fun fact is that the number of associated primes of the ideal
-- of adjacent 2by2 minors of a 2 by n matrix is the nth fibonacci
-- number f(n), defined by the recurrence f(n) = f(n-1) + f(n-2)
-- with the initial conditions f(1)=f(2)=1