read `gpolyutil.txt`; 
read `definitive_basis.txt`;

bigC:=table([seq(i=simplify(c[i]*delta^2),i=0..3)]):
index_fn:=proc(idxa,idxb)
  global x,y, bigC;
  local this_c, idx_param;
  this_c:=expand(bigC[idxa-1]);
  if idxb=1 then subs(seq(y[idx_param]=0,idx_param=1..3),this_c)
    else if idxb=2 then coeff(coeff(this_c,y[1]),y[2])
      else if idxb=3 then coeff(coeff(this_c,y[1]),y[3]) 
        else if idxb=4 then coeff(coeff(this_c,y[2]),y[3])
        fi
      fi
    fi
  fi;
  simplify(%);
end:

Csplit:=linalg[matrix](4,4,index_fn):

y1y2coeff:=proc(arga::integer,argb::integer)
  global F,x,y;
  local mult_table,i,j,G,total,loop_param;

  G:=[F(x[1]),F(x[2]),F(x[3])];

  if nargs=2 then
    mult_table   :=table([1=[1,2],G[3]=[3,4]]):
  else if nargs=3 then
      mult_table  :=table([1=[1,1,2],G[3]=[1,3,4],G[1]*G[2]=[2,2,2],G[1]*G[3]=[2,3,3],
                        G[2]*G[3]=[2,4,4]]):
      else if nargs=4 then
        mult_table :=table([1=[1,1,1,2],G[3]=[1,1,3,4],G[1]*G[2]=[1,2,2,2],G[1]*G[3]=[1,2,3,3],
                      G[2]*G[3]=[1,2,4,4],G[1]*G[2]*G[3]=[2,2,3,4],
                      G[1]*G[3]^2=[3,3,3,4],G[2]*G[3]^2=[3,4,4,4]]):
        else ERROR(`only 2,3, or 4 arguments allowed`);
      fi
    fi
  fi;
  total:=0:
  for loop_param in indices(mult_table) do
    combinat[permute](mult_table[op(loop_param)]):
    total:=total+op(loop_param)*sum('product(Csplit[args['j']+1,%[i]['j']],'j'=1..nargs)','i'=1..nops(%));
  od:

  total;
end:

# expand(simplify((c[0]*delta^2)*(c[2]*delta^2))):
# subs(seq(y[i]^2=F(x[i]),i=1..3),%):
# coeff(coeff(%,y[1]),y[2]):
# simplify(%-y1y2coeff(0,2));