//This is the ideal Fourier invariants.

ring rQ = 0,(q1,q2,q3,q4,q5,q6,q7),dp; 

ideal Invariants =
q5*q6-q4*q7,
q4*q6-q3*q7,
q4^2-q3*q5,
q2*q3-q1*q5,
q4^2*q5-q2*q7^2,
q4^3-q2*q6*q7,
q3*q4^2-q2*q6^2,
q3*q4^2-q1*q7^2,
q3^2*q4-q1*q6*q7,
q3^3-q1*q6^2;

// This is the inverse of the Fourier transform.

matrix ptoq[9][7] =
1/64,3/64,15/64,3/16,9/64,3/16,3/16,
3/32,-3/32,21/32,-3/8,-9/32,3/8,-3/8,
3/64,9/64,-3/64,-3/16,27/64,-3/16,-3/16,
3/32,-3/32,-3/32,3/8,-9/32,-3/8,3/8,
3/32,9/32,9/32,3/8,-9/32,-3/8,-3/8,
3/32,-3/32,9/32,-3/8,3/32,-3/8,3/8,
3/8,-3/8,-3/8,0,3/8,0,0,
3/32,9/32,-15/32,-3/8,-9/32,3/8,3/8,
3/32,-3/32,-15/32,3/8,3/32,3/8,-3/8;

// This is the ring of probability distributions. 

ring rP = 0,(p1,p2,p3,p4,p5,p6,p7,p8,p9),dp;

//This is the Fourier transform.

matrix qtop[7][9] =
1,1,1,1,1,1,1,1,1,
1,-1/3,1,-1/3,1,-1/3,-1/3,1,-1/3,
1,7/15,-1/15,-1/15,1/5,1/5,-1/15,-1/3,-1/3,
1,-1/3,-1/3,1/3,1/3,-1/3,0,-1/3,1/3,
1,-1/3,1,-1/3,-1/3,1/9,1/9,-1/3,1/9,
1,1/3,-1/3,-1/3,-1/3,-1/3,0,1/3,1/3,
1,-1/3,-1/3,1/3,-1/3,1/3,0,1/3,-1/3;


ideal Fourier = qtop*transpose(maxideal(1));

// This is the list of polynomial invariants.

map F = rQ, Fourier;
ideal PInvariants =  F(Invariants); 

// This is the polynomial parametrization.

ring r = 0,(c0,c1,e0,e1),dp;

ideal P = 
c0^3*e0^4+3*c0^3*e1^4+6*c0^2*c1*e0^3*e1+6*c0^2*c1*e0^2*e1^2+6*c0^2*c1*e0*e1^3+18*c0^2*c1*e1^4+6*c0*c1^2*e0^3*e1+24*c0*c1^2*e0^2*e1^2+42*c0*c1^2*e0*e1^3+36*c0*c1^2*e1^4+3*c1^3*e0^4+12*c1^3*e0^3*e1+18*c1^3*e0^2*e1^2+24*c1^3*e0*e1^3+51*c1^3*e1^4,
6*c0^3*e0^3*e1+6*c0^3*e0*e1^3+12*c0^3*e1^4+6*c0^2*c1*e0^3*e1+60*c0^2*c1*e0^2*e1^2+78*c0^2*c1*e0*e1^3+72*c0^2*c1*e1^4+6*c0*c1^2*e0^3*e1+132*c0*c1^2*e0^2*e1^2+366*c0*c1^2*e0*e1^3+144*c0*c1^2*e1^4+30*c1^3*e0^3*e1+144*c1^3*e0^2*e1^2+270*c1^3*e0*e1^3+204*c1^3*e1^4,
6*c0^3*e0^2*e1^2+6*c0^3*e1^4+3*c0^2*c1*e0^4+6*c0^2*c1*e0^3*e1+24*c0^2*c1*e0^2*e1^2+30*c0^2*c1*e0*e1^3+45*c0^2*c1*e1^4+3*c0*c1^2*e0^4+42*c0*c1^2*e0^3*e1+42*c0*c1^2*e0^2*e1^2+102*c0*c1^2*e0*e1^3+135*c0*c1^2*e1^4+6*c1^3*e0^4+24*c1^3*e0^3*e1+72*c1^3*e0^2*e1^2+84*c1^3*e0*e1^3+138*c1^3*e1^4,
6*c0^3*e0^2*e1^2+12*c0^3*e0*e1^3+6*c0^3*e1^4+12*c0^2*c1*e0^3*e1+30*c0^2*c1*e0^2*e1^2+120*c0^2*c1*e0*e1^3+54*c0^2*c1*e1^4+12*c0*c1^2*e0^3*e1+174*c0*c1^2*e0^2*e1^2+264*c0*c1^2*e0*e1^3+198*c0*c1^2*e1^4+24*c1^3*e0^3*e1+126*c1^3*e0^2*e1^2+324*c1^3*e0*e1^3+174*c1^3*e1^4,
6*c0^3*e0^3*e1+6*c0^3*e0*e1^3+12*c0^3*e1^4+6*c0^2*c1*e0^4+18*c0^2*c1*e0^3*e1+36*c0^2*c1*e0^2*e1^2+66*c0^2*c1*e0*e1^3+90*c0^2*c1*e1^4+6*c0*c1^2*e0^4+54*c0*c1^2*e0^3*e1+144*c0*c1^2*e0^2*e1^2+174*c0*c1^2*e0*e1^3+270*c0*c1^2*e1^4+12*c1^3*e0^4+66*c1^3*e0^3*e1+108*c1^3*e0^2*e1^2+186*c1^3*e0*e1^3+276*c1^3*e1^4,
12*c0^3*e0^2*e1^2+12*c0^3*e1^4+12*c0^2*c1*e0^3*e1+36*c0^2*c1*e0^2*e1^2+108*c0^2*c1*e0*e1^3+60*c0^2*c1*e1^4+12*c0*c1^2*e0^3*e1+144*c0*c1^2*e0^2*e1^2+324*c0*c1^2*e0*e1^3+168*c0*c1^2*e1^4+24*c1^3*e0^3*e1+144*c1^3*e0^2*e1^2+288*c1^3*e0*e1^3+192*c1^3*e1^4,
24*c0^3*e0^2*e1^2+48*c0^3*e0*e1^3+24*c0^3*e1^4+24*c0^2*c1*e0^3*e1+192*c0^2*c1*e0^2*e1^2+408*c0^2*c1*e0*e1^3+240*c0^2*c1*e1^4+96*c0*c1^2*e0^3*e1+552*c0*c1^2*e0^2*e1^2+1200*c0*c1^2*e0*e1^3+744*c0*c1^2*e1^4+72*c1^3*e0^3*e1+576*c1^3*e0^2*e1^2+1224*c1^3*e0*e1^3+720*c1^3*e1^4,
6*c0^3*e0^2*e1^2+12*c0^3*e0*e1^3+6*c0^3*e1^4+24*c0^2*c1*e0^3*e1+42*c0^2*c1*e0^2*e1^2+60*c0^2*c1*e0*e1^3+90*c0^2*c1*e1^4+18*c0*c1^2*e0^4+60*c0*c1^2*e0^3*e1+114*c0*c1^2*e0^2*e1^2+168*c0*c1^2*e0*e1^3+288*c0*c1^2*e1^4+6*c1^3*e0^4+60*c1^3*e0^3*e1+126*c1^3*e0^2*e1^2+192*c1^3*e0*e1^3+264*c1^3*e1^4,
24*c0^3*e0*e1^3+60*c0^2*c1*e0^2*e1^2+96*c0^2*c1*e0*e1^3+60*c0^2*c1*e1^4+36*c0*c1^2*e0^3*e1+132*c0*c1^2*e0^2*e1^2+276*c0*c1^2*e0*e1^3+204*c0*c1^2*e1^4+12*c1^3*e0^3*e1+144*c1^3*e0^2*e1^2+324*c1^3*e0*e1^3+168*c1^3*e1^4;

// This checks that the polynomial parametrization
// lies on the probability simplex.
// It requires suma.sing. Most likely, you should 
// change the directory where you saved this file. 

// If you do have this file, you should uncomment
// the following two lines.

// < "/home/lgp/singular/suma.sing";

// Suma(Substitute(1,P));

// This checks that the PInvariants vanish at
// the polynomial parametrization.

map Evaluate = rP, P;

// The following command takes a lot of space and time to
// finish for larger models.

// ideal Z = Evaluate(PInvariants);

setring rP;
ideal Z;
int i;
for (i=1; i<= size(PInvariants); i++)
{
  i;
  Z = PInvariants[i];
  setring r;
  Evaluate(Z);
  setring rP;
}