---Building some helpful functions

--Row reduced echelon form function
--Input: a matrix defined over any field
--Output: the row reduced echelon form
rref = M->(
    F := ring(M);
    R := F[vars(0..(numcols M)-1)];
    V := transpose vars R;
    N := sub(M,R);
    I := ideal(N*V);
    G := reverse flatten entries gens gb I;
    G = apply(G,g->(
	    c := (last(coefficients g))_(0,0);
	    1/sub(c,F)*g));
    T := gens R;
    matrix apply(G,g->apply(T,t->coefficient(t,g)))
    )

--Multiplication matrix for a ring element
--Input: a polynomial f and an ideal I (zero-dimensional!)
--Output: the multiplication matrix for f on the quotient S/I (provided I is zero-dimensional)
--The columns of the matrix represent multiplication of basis elements by f
multiplicationMatrix = (f,I)->(
    S := ring(I);
    if dim(S/I)>0 then(
	return("Error: not a zero dimensional ideal.")
	)else(
	B := flatten entries sub(basis (S/I),S);
	M := transpose matrix apply(B,b->(
		H := (f*b)%I;
		apply(B,c->coefficient(c,H))
		));
	K := coefficientRing(S);
	(B,map(K^(numcols M),K^(numcols M),M))
	))

--Division algorithm (from CLO p.65)
divisionAlgorithm=(f,L)->(
    --get the ring in which polynomials are defined
    S := ring f;
    --initialize quotient list and remainder to zero
    quot := new MutableList from apply(L,i->0);
    rem := 0;
    p := f;
    while p != 0 do(
	i := 0;
	divisionoccurred := false;
	while i<=(length(L)-1) and not divisionoccurred do(
	    diffi := (flatten exponents leadTerm(p))-(flatten exponents leadTerm(L_i));
	    if all(diffi,i->i>=0) then(
		nqt := sub(leadTerm(p)/leadTerm(L_i),S);
		quot#i = quot#i+nqt;
		p = p-nqt*L_i;
		divisionoccurred = true;
		)else(
		i = i+1;
		));
	if not divisionoccurred then(
	    rem = rem+leadTerm(p);
	    p = p-leadTerm(p);
	    );
	);
    (toList quot,rem)	
    )

--Computing S-polynomials
--Input: two polynomials f and g
--Output: the S-polynomial of f and g
sPolynomial=(f,g)->(
    S := ring f;
    lf := leadTerm f;
    lg := leadTerm g;
    LCM := lcm(lf,lg);
    m1 := sub(LCM/lf,S);
    m2 := sub(LCM/lg,S);
    m1*f-m2*g
    )

--Buchberger's algorithm
--Input: a list L of polynomials generating some ideal (monomial order should be determined in definition of polynomial ring)
--Output: a Grobner basis with respect to the monomial order
--Comment: Macaulay2 has a 'gb' command which is much more efficient than the following code.  What you find below is the most naive implementation!
--There is at least one glaring redundancy - see if you can spot it and make the code better.  Or write your own!
buchbergerAlgorithm=L->(
    GB := L;
    stop := false;
    while not stop do(
	rem := 0;
	i := 0;
	k := length GB;
	while (rem==0) and (i<(k-1)) do(
	    j := i+1;
	    while (rem==0) and (j<k) do(
		Sp := sPolynomial(GB_i,GB_j);
		rem = last divisionAlgorithm(Sp,GB);
		if rem != 0 then(
		    GB = append(GB,rem)
		    )else(
		    j = j+1
		    )
		);
	    if rem==0 then i = i+1
	    );
	if rem==0 then stop=true
	);
    GB
    )

--coupler curve for 4 bar mechanism with fixed vertices at (0,0) and (a,0)
--Input: x-coordinate of other fixed vertex and lengths for the five bars, as well
----as the ring you want the implicit function to be in.
--Output: equation for coupler curve of the corresponding mechanism. 
--(this may not have real points if the lengths were not chosen compatibly)
couplerCurve=(a,r1,r2,r3,r4,r5,R)->(
    S := QQ[gens R, (symbol x)_1, (symbol x)_2, (symbol y)_1, (symbol y)_2];
    l1 := x_1^2+y_1^2-r1^2;
    l2 := (x_2-a)^2+y_2^2-r2^2;
    l3 := (x_1-x_2)^2+(y_1-y_2)^2-r3^2;
    l4 := (sub(R_0,S)-x_1)^2+(sub(R_1,S)-y_1)^2-r4^2;
    l5 := (sub(R_0,S)-x_2)^2+(sub(R_1,S)-y_2)^2-r5^2;
    I := ideal(l1,l2,l3,l4,l5);
    J := eliminate(I,{x_1,x_2,y_1,y_2});
    sub(J_0,R)
    )

--Get implicit equation for trigonometric curve of the form r=cos(n/d*theta)
--Input: positive integers n and d, and the ring you want the implicit polynomial to be in
--Output: implicit equation in x and y for the trig curve r=cos(n/d*theta)
trigImplicitization=(n,d,R)->(
    S := QQ[gens R,symbol r,symbol s,symbol c,symbol a,symbol b,symbol i];
    eul := (a+b*i)^d-(c+i*s)^n;
    I := ideal(r-a,x^2+y^2-r^2,s^2+c^2-1,i^2+1,x-r*c,y-r*s,a^2+b^2-1,eul);
    J := eliminate(I,{r,s,c,a,b,i});
    use R;
    (sub(J,R))_0
    )

--Get real planar singular points of a planar curve
--This requires that the package NumericalAlgebraicGeometry is loaded
--Input: two variable polynomial f and an ideal I
--Output: the singular points of f that are not in V(I) (make I=ideal(0) to find all singular points)
needsPackage "NumericalAlgebraicGeometry"
getRealPlanarSingularPointsNAG=(f,I)->(
    R := ring(f);
    G := gens R;
    J := sub(I,R);
    F := flatten entries gens(saturate(ideal(f,diff(first G,f),diff(last G,f)),J));
    sols := solveSystem(F);
    RP := realPoints sols;
    apply(RP,p->flatten entries clean(10^(-6)_RR,matrix p))
    )

--Output the implicit equation and singular points of a trig curve to a text file
--Input: positive integers n,d, and a filename
--Output: a file located in filename containing the implicit equation of the curve r=cos(n/d*theta)
--along with the singular locus formatted well for copying and pasting into sage
getTrigSingularPointsNAG=(n,d,filename)->(
    R := QQ[x,y];
    J := ideal(x,y);
    f := trigImplicitization(n,d,R);
    sing := getRealPlanarSingularPointsNAG(f,J);
    filename<<toString(f)<<endl<<sing<<close
    )

--Compute the number of critical points of the distance function from a point P
--to a variety defined by an ideal I
--Input: an ideal I and a point P
--Output: the number of critical points of the distance function from a point P to a 
--variety defined by the ideal I (disregarding critical points in the singular locus)
--If P is a general point, this is the Euclidean distance degree.
euclideanDistanceDegree=(I,P)->(
    J := radical I;
    S := ring J;
    D := sum apply(length gens S,i->(S_i-(P_i)_S)^2);
    c := codim I;
    MS := minors(c,jacobian J);
    M := (jacobian J)|(jacobian ideal D);
    MN := minors(c+1,M);
    degree saturate(J+MN,MS)
    )

--A script to compute the affine Hilbert function
--Input: An ideal I and a degree d
--Output: The affine hilbert function dim (R/I)_(<=d)
affineHilbertFunction=(d,I)->(
    S := ring(I);
    afH := 0;
    for i from 0 to d do(
	afH = afH+numcols basis(i,S/I)
	);
    afH
    )

--A script to display affine Hilbert function up to a certain point
--Input: An ideal I and a degree d
--Output: A table displaying the values of the affine hilbert function up through degree d
affineHilbertFunctionTable=(d,I)->(
    S := ring(I);
    L := apply(d+1,i->numcols basis(i,S/I));
    A := accumulate(plus,prepend(0,L));
    tp := prepend("d",toList(0..d));
    bm := prepend("AffineHilbert(d)",A);
    netList {tp,bm}
    )