//This is the ideal Fourier invariants.

ring rQ = 0, (q1,q2,q3,q4,q5,q6,q7,q8,q9,q10), dp;
ideal Invariants = 
q9^2-q8*q10,
q7*q9-q6*q10,
q6*q9-q5*q10,
q5*q9-q4*q10,
q3*q9-q2*q10,
q7*q8-q5*q10,
q6*q8-q4*q10,
q5*q8-q4*q9,
q3*q8-q2*q9,
q2*q8-q1*q10,
q6^2-q5*q7,
q5*q6-q4*q7,
q3*q6-q2*q7,
q5^2-q4*q6,
q3*q5-q2*q6,
q2*q5-q1*q7,
q3*q4-q1*q7,
q2*q4-q1*q6,
q5*q7^2-q3*q10^2,
q4*q7^2-q2*q10^2,
q4*q6*q7-q2*q9*q10,
q4*q5*q7-q1*q10^2,
q4^2*q7-q1*q9*q10,
q4^2*q6-q1*q8*q10,
q4^2*q5-q1*q8*q9,
q4^3-q1*q8^2,
q2^3-q1*q3^2;

// This is the inverse of the Fourier transform.

matrix ptoq[15][10] =
1/256,9/256,3/128,21/256,9/64,69/256,9/128,9/128,9/64,21/128,
9/256,9/256,-9/128,117/256,9/64,-99/256,-27/128,45/128,9/64,-63/128,
9/256,9/256,-9/128,45/256,-27/64,117/256,-27/128,9/128,-27/64,45/128,
9/128,-27/128,9/64,45/128,-9/32,-63/128,27/64,9/64,-9/32,9/64,
3/256,27/256,9/128,-9/256,-9/64,39/256,27/128,-9/128,-9/64,-21/128,
9/128,9/128,-9/64,-27/128,9/32,45/128,-27/64,-27/64,9/32,9/64,
3/128,-9/128,3/64,-9/128,9/32,-45/128,9/64,-9/64,9/32,-9/64,
3/128,27/128,9/64,15/128,9/32,15/128,-9/64,-9/64,-9/32,-21/64,
9/128,9/128,-9/64,45/128,-9/32,-27/128,9/64,-27/64,9/32,9/64,
9/64,9/64,-9/32,9/64,9/16,-63/64,9/32,-9/32,-9/16,27/32,
9/64,-27/64,9/32,9/64,-9/16,45/64,-9/32,-9/32,9/16,-9/32,
9/64,9/64,-9/32,-27/64,-9/16,45/64,9/32,9/32,9/16,-27/32,
9/64,-27/64,9/32,-27/64,9/16,9/64,-9/32,9/32,-9/16,9/32,
3/128,27/128,9/64,-21/128,-9/32,-69/128,-9/64,9/64,9/32,21/64,
9/128,9/128,-9/64,-63/128,9/32,9/128,9/64,27/64,-9/32,-9/64;

// This is the ring of probability distributions. 

ring rP = 0,(p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12,p13,p14,p15),dp;

//This is the Fourier transform.

matrix qtop[10][15] =
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
1,1/9,1/9,-1/3,1,1/9,-1/3,1,1/9,1/9,-1/3,1/9,-1/3,1,1/9,
1,-1/3,-1/3,1/3,1,-1/3,1/3,1,-1/3,-1/3,1/3,-1/3,1/3,1,-1/3,
1,13/21,5/21,5/21,-1/7,-1/7,-1/7,5/21,5/21,1/21,1/21,-1/7,-1/7,-1/3,-1/3,
1,1/9,-1/3,-1/9,-1/3,1/9,1/3,1/3,-1/9,1/9,-1/9,-1/9,1/9,-1/3,1/9,
1,-11/69,13/69,-7/69,13/69,5/69,-5/23,5/69,-1/23,-7/69,5/69,5/69,1/69,-1/3,1/69,
1,-1/3,-1/3,1/3,1,-1/3,1/3,-1/3,1/9,1/9,-1/9,1/9,-1/9,-1/3,1/9,
1,5/9,1/9,1/9,-1/3,-1/3,-1/3,-1/3,-1/3,-1/9,-1/9,1/9,1/9,1/3,1/3,
1,1/9,-1/3,-1/9,-1/3,1/9,1/3,-1/3,1/9,-1/9,1/9,1/9,-1/9,1/3,-1/9,
1,-1/3,5/21,1/21,-1/3,1/21,-1/7,-1/3,1/21,1/7,-1/21,-1/7,1/21,1/3,-1/21;

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,f0,f1),dp;

ideal P = 
c0^3*f0^5+3*c0^3*f1^5+6*c0^2*c1*f0^4*f1+3*c0^2*c1*f0^3*f1^2+3*c0^2*c1*f0^2*f1^3+6*c0^2*c1*f0*f1^4+18*c0^2*c1*f1^5+6*c0*c1^2*f0^4*f1+21*c0*c1^2*f0^3*f1^2+3*c0*c1^2*f0^2*f1^3+42*c0*c1^2*f0*f1^4+36*c0*c1^2*f1^5+3*c1^3*f0^5+12*c1^3*f0^4*f1+12*c1^3*f0^3*f1^2+6*c1^3*f0^2*f1^3+24*c1^3*f0*f1^4+51*c1^3*f1^5,
9*c0^3*f0^4*f1+9*c0^3*f0*f1^4+18*c0^3*f1^5+63*c0^2*c1*f0^3*f1^2+63*c0^2*c1*f0^2*f1^3+90*c0^2*c1*f0*f1^4+108*c0^2*c1*f1^5+63*c0*c1^2*f0^3*f1^2+333*c0*c1^2*f0^2*f1^3+360*c0*c1^2*f0*f1^4+216*c0*c1^2*f1^5+27*c1^3*f0^4*f1+126*c1^3*f0^3*f1^2+216*c1^3*f0^2*f1^3+297*c1^3*f0*f1^4+306*c1^3*f1^5,
9*c0^3*f0^3*f1^2+9*c0^3*f0^2*f1^3+18*c0^3*f1^5+9*c0^2*c1*f0^4*f1+18*c0^2*c1*f0^3*f1^2+126*c0^2*c1*f0^2*f1^3+63*c0^2*c1*f0*f1^4+108*c0^2*c1*f1^5+9*c0*c1^2*f0^4*f1+126*c0*c1^2*f0^3*f1^2+180*c0*c1^2*f0^2*f1^3+441*c0*c1^2*f0*f1^4+216*c0*c1^2*f1^5+18*c1^3*f0^4*f1+99*c1^3*f0^3*f1^2+297*c1^3*f0^2*f1^3+252*c1^3*f0*f1^4+306*c1^3*f1^5,
18*c0^3*f0^3*f1^2+36*c0^3*f0*f1^4+18*c0^3*f1^5+36*c0^2*c1*f0^3*f1^2+198*c0^2*c1*f0^2*f1^3+306*c0^2*c1*f0*f1^4+108*c0^2*c1*f1^5+36*c0*c1^2*f0^3*f1^2+630*c0*c1^2*f0^2*f1^3+1062*c0*c1^2*f0*f1^4+216*c0*c1^2*f1^5+126*c1^3*f0^3*f1^2+540*c1^3*f0^2*f1^3+972*c1^3*f0*f1^4+306*c1^3*f1^5,
3*c0^3*f0^3*f1^2+3*c0^3*f0^2*f1^3+6*c0^3*f1^5+3*c0^2*c1*f0^5+6*c0^2*c1*f0^4*f1+18*c0^2*c1*f0^3*f1^2+6*c0^2*c1*f0^2*f1^3+30*c0^2*c1*f0*f1^4+45*c0^2*c1*f1^5+3*c0*c1^2*f0^5+42*c0*c1^2*f0^4*f1+36*c0*c1^2*f0^3*f1^2+6*c0*c1^2*f0^2*f1^3+102*c0*c1^2*f0*f1^4+135*c0*c1^2*f1^5+6*c1^3*f0^5+24*c1^3*f0^4*f1+51*c1^3*f0^3*f1^2+21*c1^3*f0^2*f1^3+84*c1^3*f0*f1^4+138*c1^3*f1^5,
36*c0^3*f0^2*f1^3+18*c0^3*f0*f1^4+18*c0^3*f1^5+18*c0^2*c1*f0^4*f1+54*c0^2*c1*f0^3*f1^2+162*c0^2*c1*f0^2*f1^3+252*c0^2*c1*f0*f1^4+162*c0^2*c1*f1^5+18*c0*c1^2*f0^4*f1+270*c0*c1^2*f0^3*f1^2+486*c0*c1^2*f0^2*f1^3+576*c0*c1^2*f0*f1^4+594*c0*c1^2*f1^5+36*c1^3*f0^4*f1+180*c1^3*f0^3*f1^2+540*c1^3*f0^2*f1^3+666*c1^3*f0*f1^4+522*c1^3*f1^5,
6*c0^3*f0^2*f1^3+18*c0^3*f0*f1^4+18*c0^2*c1*f0^3*f1^2+36*c0^2*c1*f0^2*f1^3+144*c0^2*c1*f0*f1^4+18*c0^2*c1*f1^5+18*c0*c1^2*f0^3*f1^2+252*c0*c1^2*f0^2*f1^3+252*c0*c1^2*f0*f1^4+126*c0*c1^2*f1^5+36*c1^3*f0^3*f1^2+162*c1^3*f0^2*f1^3+378*c1^3*f0*f1^4+72*c1^3*f1^5,
6*c0^3*f0^4*f1+6*c0^3*f0*f1^4+12*c0^3*f1^5+6*c0^2*c1*f0^5+18*c0^2*c1*f0^4*f1+30*c0^2*c1*f0^3*f1^2+6*c0^2*c1*f0^2*f1^3+66*c0^2*c1*f0*f1^4+90*c0^2*c1*f1^5+6*c0*c1^2*f0^5+54*c0*c1^2*f0^4*f1+102*c0*c1^2*f0^3*f1^2+42*c0*c1^2*f0^2*f1^3+174*c0*c1^2*f0*f1^4+270*c0*c1^2*f1^5+12*c1^3*f0^5+66*c1^3*f0^4*f1+84*c1^3*f0^3*f1^2+24*c1^3*f0^2*f1^3+186*c1^3*f0*f1^4+276*c1^3*f1^5,
18*c0^3*f0^3*f1^2+18*c0^3*f0^2*f1^3+36*c0^3*f1^5+18*c0^2*c1*f0^4*f1+72*c0^2*c1*f0^3*f1^2+144*c0^2*c1*f0^2*f1^3+234*c0^2*c1*f0*f1^4+180*c0^2*c1*f1^5+18*c0*c1^2*f0^4*f1+180*c0*c1^2*f0^3*f1^2+576*c0*c1^2*f0^2*f1^3+666*c0*c1^2*f0*f1^4+504*c0*c1^2*f1^5+36*c1^3*f0^4*f1+234*c1^3*f0^3*f1^2+486*c1^3*f0^2*f1^3+612*c1^3*f0*f1^4+576*c1^3*f1^5,
36*c0^3*f0^3*f1^2+72*c0^3*f0*f1^4+36*c0^3*f1^5+36*c0^2*c1*f0^4*f1+108*c0^2*c1*f0^3*f1^2+360*c0^2*c1*f0^2*f1^3+432*c0^2*c1*f0*f1^4+360*c0^2*c1*f1^5+36*c0*c1^2*f0^4*f1+432*c0*c1^2*f0^3*f1^2+1116*c0*c1^2*f0^2*f1^3+1188*c0*c1^2*f0*f1^4+1116*c0*c1^2*f1^5+72*c1^3*f0^4*f1+432*c1^3*f0^3*f1^2+972*c1^3*f0^2*f1^3+1332*c1^3*f0*f1^4+1080*c1^3*f1^5,
72*c0^3*f0^2*f1^3+36*c0^3*f0*f1^4+36*c0^3*f1^5+72*c0^2*c1*f0^3*f1^2+360*c0^2*c1*f0^2*f1^3+684*c0^2*c1*f0*f1^4+180*c0^2*c1*f1^5+180*c0*c1^2*f0^3*f1^2+1116*c0*c1^2*f0^2*f1^3+2088*c0*c1^2*f0*f1^4+504*c0*c1^2*f1^5+180*c1^3*f0^3*f1^2+1188*c1^3*f0^2*f1^3+1944*c1^3*f0*f1^4+576*c1^3*f1^5,
72*c0^3*f0^2*f1^3+36*c0^3*f0*f1^4+36*c0^3*f1^5+180*c0^2*c1*f0^3*f1^2+324*c0^2*c1*f0^2*f1^3+432*c0^2*c1*f0*f1^4+360*c0^2*c1*f1^5+108*c0*c1^2*f0^4*f1+396*c0*c1^2*f0^3*f1^2+972*c0*c1^2*f0^2*f1^3+1296*c0*c1^2*f0*f1^4+1116*c0*c1^2*f1^5+36*c1^3*f0^4*f1+432*c1^3*f0^3*f1^2+1080*c1^3*f0^2*f1^3+1260*c1^3*f0*f1^4+1080*c1^3*f1^5,
36*c0^3*f0^2*f1^3+108*c0^3*f0*f1^4+36*c0^2*c1*f0^3*f1^2+432*c0^2*c1*f0^2*f1^3+648*c0^2*c1*f0*f1^4+180*c0^2*c1*f1^5+252*c0*c1^2*f0^3*f1^2+1080*c0*c1^2*f0^2*f1^3+1944*c0*c1^2*f0*f1^4+612*c0*c1^2*f1^5+144*c1^3*f0^3*f1^2+1188*c1^3*f0^2*f1^3+2052*c1^3*f0*f1^4+504*c1^3*f1^5,
6*c0^3*f0^3*f1^2+12*c0^3*f0*f1^4+6*c0^3*f1^5+24*c0^2*c1*f0^4*f1+30*c0^2*c1*f0^3*f1^2+12*c0^2*c1*f0^2*f1^3+60*c0^2*c1*f0*f1^4+90*c0^2*c1*f1^5+18*c0*c1^2*f0^5+60*c0*c1^2*f0^4*f1+84*c0*c1^2*f0^3*f1^2+30*c0*c1^2*f0^2*f1^3+168*c0*c1^2*f0*f1^4+288*c0*c1^2*f1^5+6*c1^3*f0^5+60*c1^3*f0^4*f1+96*c1^3*f0^3*f1^2+30*c1^3*f0^2*f1^3+192*c1^3*f0*f1^4+264*c1^3*f1^5,
18*c0^3*f0^2*f1^3+54*c0^3*f0*f1^4+72*c0^2*c1*f0^3*f1^2+198*c0^2*c1*f0^2*f1^3+198*c0^2*c1*f0*f1^4+180*c0^2*c1*f1^5+54*c0*c1^2*f0^4*f1+234*c0*c1^2*f0^3*f1^2+468*c0*c1^2*f0^2*f1^3+576*c0*c1^2*f0*f1^4+612*c0*c1^2*f1^5+18*c1^3*f0^4*f1+198*c1^3*f0^3*f1^2+540*c1^3*f0^2*f1^3+684*c1^3*f0*f1^4+504*c1^3*f1^5;

// 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;
}