FORM version 2.3 Dec 10 1997

    f dx,dy;
    s n,x,y,intx,inty,a,b,c,d,e,f,g,h,l,m;
    s f00,f01,f10,f02,f20,f11,f03,f30,f12,f21;
    s A,B,C,D,E,F,G,H,I,L,dmy;
    *
    .global

    *
    g F00= 1;
    g F01= dy*y*(1-x-y);
    g F10= dx*x*(1-x-y);
    g F02= dy*dy*y^2*(1-x-y)^2;
    g F20= dx*dx*x^2*(1-x-y)^2;
    g F11= dx*dy*x*y*(1-x-y)^2;
    g F03= dy*dy*dy*y^3*(1-x-y)^3;
    g F30= dx*dx*dx*x^3*(1-x-y)^3;
    g F12= dx*dy*dy*x*y^2*(1-x-y)^3;
    g F21= dx*dx*dy*x^2*y*(1-x-y)^3;
    nwrite statistics;
    *
    repeat;
     id dx*x^n?= n*x^(n-1);
     id dy*y^n?= n*y^(n-1);
    endrepeat;
    *
    print;
    *
    .store

   F00 =
      1;

   F01 =
      1 - x - 2*y;

   F10 =
      1 - 2*x - y;

   F02 =
      2 + 12*x*y - 4*x + 2*x^2 - 12*y + 12*y^2;

   F20 =
      2 + 12*x*y - 12*x + 12*x^2 - 4*y + 2*y^2;

   F11 =
      1 + 8*x*y - 4*x + 3*x^2 - 4*y + 3*y^2;

   F03 =
      6 + 144*x*y - 180*x*y^2 - 18*x - 72*x^2*y + 18*x^2 - 6*x^3 - 72*y + 180*
      y^2 - 120*y^3;

   F30 =
      6 + 144*x*y - 72*x*y^2 - 72*x - 180*x^2*y + 180*x^2 - 120*x^3 - 18*y + 
      18*y^2 - 6*y^3;

   F12 =
      2 + 72*x*y - 72*x*y^2 - 12*x - 54*x^2*y + 18*x^2 - 8*x^3 - 18*y + 36*y^2
       - 20*y^3;

   F21 =
      2 + 72*x*y - 54*x*y^2 - 18*x - 72*x^2*y + 36*x^2 - 20*x^3 - 12*y + 18*
      y^2 - 8*y^3;


    *
    g I1= intx*inty*F01*F10;
    g I2= intx*inty*F01*F02;
    g I3= intx*inty*F01*F20;
    g I4= intx*inty*F01*F11;
    g I5= intx*inty*F01*F03;
    g I6= intx*inty*F01*F30;
    g I7= intx*inty*F01*F12;
    g I8= intx*inty*F01*F21;
    *
    g I9= intx*inty*F10*F10;
    g I10= intx*inty*F10*F02;
    g I11= intx*inty*F10*F20;
    g I12= intx*inty*F10*F11;
    g I13= intx*inty*F10*F03;
    g I14= intx*inty*F10*F30;
    g I15= intx*inty*F10*F12;
    g I16= intx*inty*F10*F21;
    *
    g I17= intx*inty*F02*F10;
    g I18= intx*inty*F02*F02;
    g I19= intx*inty*F02*F20;
    g I20= intx*inty*F02*F11;
    g I21= intx*inty*F02*F03;
    g I22= intx*inty*F02*F30;
    g I23= intx*inty*F02*F12;
    g I24= intx*inty*F02*F21;
    *
    g I25= intx*inty*F20*F10;
    g I26= intx*inty*F20*F02;
    g I27= intx*inty*F20*F20;
    g I28= intx*inty*F20*F11;
    g I29= intx*inty*F20*F03;
    g I30= intx*inty*F20*F30;
    g I31= intx*inty*F20*F12;
    g I32= intx*inty*F20*F21;
    *
    g I33= intx*inty*F11*F10;
    g I34= intx*inty*F11*F02;
    g I35= intx*inty*F11*F20;
    g I36= intx*inty*F11*F11;
    g I37= intx*inty*F11*F03;
    g I38= intx*inty*F11*F30;
    g I39= intx*inty*F11*F12;
    g I40= intx*inty*F11*F21;
    *
    g I41= intx*inty*F03*F10;
    g I42= intx*inty*F03*F02;
    g I43= intx*inty*F03*F20;
    g I44= intx*inty*F03*F11;
    g I45= intx*inty*F03*F03;
    g I46= intx*inty*F03*F30;
    g I47= intx*inty*F03*F12;
    g I48= intx*inty*F03*F21;
    *
    g I49= intx*inty*F30*F10;
    g I50= intx*inty*F30*F02;
    g I51= intx*inty*F30*F20;
    g I52= intx*inty*F30*F11;
    g I53= intx*inty*F30*F03;
    g I54= intx*inty*F30*F30;
    g I55= intx*inty*F30*F12;
    g I56= intx*inty*F30*F21;
    *
    g I57= intx*inty*F12*F10;
    g I58= intx*inty*F12*F02;
    g I59= intx*inty*F12*F20;
    g I60= intx*inty*F12*F11;
    g I61= intx*inty*F12*F03;
    g I62= intx*inty*F12*F30;
    g I63= intx*inty*F12*F12;
    g I64= intx*inty*F12*F21;
    *
    g I65= intx*inty*F21*F10;
    g I66= intx*inty*F21*F02;
    g I67= intx*inty*F21*F20;
    g I68= intx*inty*F21*F11;
    g I69= intx*inty*F21*F03;
    g I70= intx*inty*F21*F30;
    g I71= intx*inty*F21*F12;
    g I72= intx*inty*F21*F21;
    nwrite statistics;
    *
    id inty*y^n?= (1-x)^(n+1)/(n+1);
    id intx*x^n?= 1/(n+1);
    *
    print;
    *
    .store

   I1 =
      1/24;

   I2 = 0;

   I3 = 0;

   I4 = 0;

   I5 = 0;

   I6 = 0;

   I7 = 0;

   I8 = 0;

   I9 =
      1/12;

   I10 = 0;

   I11 = 0;

   I12 = 0;

   I13 = 0;

   I14 = 0;

   I15 = 0;

   I16 = 0;

   I17 = 0;

   I18 =
      2/15;

   I19 =
      1/45;

   I20 =
      1/30;

   I21 = 0;

   I22 = 0;

   I23 = 0;

   I24 = 0;

   I25 = 0;

   I26 =
      1/45;

   I27 =
      2/15;

   I28 =
      1/30;

   I29 = 0;

   I30 = 0;

   I31 = 0;

   I32 = 0;

   I33 = 0;

   I34 =
      1/30;

   I35 =
      1/30;

   I36 =
      1/45;

   I37 = 0;

   I38 = 0;

   I39 = 0;

   I40 = 0;

   I41 = 0;

   I42 = 0;

   I43 = 0;

   I44 = 0;

   I45 =
      9/14;

   I46 =
      9/280;

   I47 =
      3/28;

   I48 =
      3/70;

   I49 = 0;

   I50 = 0;

   I51 = 0;

   I52 = 0;

   I53 =
      9/280;

   I54 =
      9/14;

   I55 =
      3/70;

   I56 =
      3/28;

   I57 = 0;

   I58 = 0;

   I59 = 0;

   I60 = 0;

   I61 =
      3/28;

   I62 =
      3/70;

   I63 =
      3/70;

   I64 =
      9/280;

   I65 = 0;

   I66 = 0;

   I67 = 0;

   I68 = 0;

   I69 =
      3/70;

   I70 =
      3/28;

   I71 =
      9/280;

   I72 =
      3/70;


    *
    g eq1= intx*inty*F00*(x-a*F00-b*F01-c*F10);
    g eq2= intx*inty*F01*(x-a*F00-b*F01-c*F10);
    g eq3= intx*inty*F10*(x-a*F00-b*F01-c*F10);
    g eq4= intx*inty*F00*(y-d*F00-e*F01-f*F10);
    g eq5= intx*inty*F01*(y-d*F00-e*F01-f*F10);
    g eq6= intx*inty*F10*(y-d*F00-e*F01-f*F10);
    nwrite statistics;
    *
    id inty*y^n?= (1-x)^(n+1)/(n+1);
    id intx*x^n?= 1/(n+1);
    *
    id a= 1/3;
    id b= -1/2*c;
    id c= -2/3;
    *
    id d= 1/3;
    id f= -1/2*e;
    id e= -2/3;
    *
    b a,b,c;
    print;
    *
    .store

   eq1 = 0;

   eq2 = 0;

   eq3 = 0;

   eq4 = 0;

   eq5 = 0;

   eq6 = 0;


    *
    g eq8= intx*inty*F00*(x^2-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11);
    g eq9= intx*inty*F01*(x^2-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11);
    g eq10= intx*inty*F10*(x^2-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11);
    g eq11= intx*inty*F02*(x^2-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11);
    g eq12= intx*inty*F20*(x^2-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11);
    g eq13= intx*inty*F11*(x^2-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11);
    nwrite statistics;
    *
    id inty*y^n?= (1-x)^(n+1)/(n+1);
    id intx*x^n?= 1/(n+1);
    *
    id a= 1/6;
    id b= -1/2*c;
    id c= -8/15;
    id f= 30*(- 2/15*d - 1/45*e);
    id d= 1/3*e;
    id e= 3/20;
    *
    print;
    *
    .store

   eq8 = 0;

   eq9 = 0;

   eq10 = 0;

   eq11 = 0;

   eq12 = 0;

   eq13 = 0;


    *
    g eq14= intx*inty*F00*(y^2-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11);
    g eq15= intx*inty*F01*(y^2-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11);
    g eq16= intx*inty*F10*(y^2-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11);
    g eq17= intx*inty*F02*(y^2-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11);
    g eq18= intx*inty*F20*(y^2-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11);
    g eq19= intx*inty*F11*(y^2-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11);
    nwrite statistics;
    *
    id inty*y^n?= (1-x)^(n+1)/(n+1);
    id intx*x^n?= 1/(n+1);
    *
    id a= 1/6;
    id b= -2*c;
    id f= -3/2*(d+e);
    id c= 4/15;
    id e= 1/3*d;
    id d= 3/20;
    *
    print;
    *
    .store

   eq14 = 0;

   eq15 = 0;

   eq16 = 0;

   eq17 = 0;

   eq18 = 0;

   eq19 = 0;


    *
    g eq20= intx*inty*F00*(x*y-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11);
    g eq21= intx*inty*F01*(x*y-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11);
    g eq22= intx*inty*F10*(x*y-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11);
    g eq23= intx*inty*F02*(x*y-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11);
    g eq24= intx*inty*F20*(x*y-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11);
    g eq25= intx*inty*F11*(x*y-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11);
    nwrite statistics;
    *
    id inty*y^n?= (1-x)^(n+1)/(n+1);
    id intx*x^n?= 1/(n+1);
    *
    id a= 1/12;
    id c= -(1/5+2*b);
    id b= -1/15;
    id f= -30*(2/15*d+1/45*e);
    id e= d;
    id d= -27/360;
    *
    print;
    *
    .store

   eq20 = 0;

   eq21 = 0;

   eq22 = 0;

   eq23 = 0;

   eq24 = 0;

   eq25 = 0;


    *
    g R1= dx*F02;
    g R2= dy*F02;
    g R3= dx*F20;
    g R4= dy*F20;
    g R5= dx*F11;
    g R6= dy*F11;
    g R7= dx*dx*F20;
    g R8= dx*dy*F20;
    g R9= dy*dy*F20;
    g R10= dx*dx*F02;
    g R11= dx*dy*F02;
    g R12= dy*dy*F02;
    g R13= dx*dx*F11;
    g R14= dx*dy*F11;
    g R15= dy*dy*F11;
    g R16= dx*F03;
    g R17= dy*F03;
    g R18= dx*dx*F03;
    g R19= dx*dy*F03;
    g R20= dy*dy*F03;
    g R21= dx*F30;
    g R22= dy*F30;
    g R23= dx*dx*F30;
    g R24= dx*dy*F30;
    g R25= dy*dy*F30;
    g R26= dx*F12;
    g R27= dy*F12;
    g R28= dx*dx*F12;
    g R29= dx*dy*F12;
    g R30= dy*dy*F12;
    g R31= dx*F21;
    g R32= dy*F21;
    g R33= dx*dx*F21;
    g R34= dx*dy*F21;
    g R35= dy*dy*F21;
    nwrite statistics;
    *
    repeat;
     id dx*x^n?= n*x^(n-1);
     id dy*y^n?= n*y^(n-1);
    endrepeat;
    *
    id x^2= a*f00+b*f01+c*f10+d*f02+e*f20+f*f11;
    id a= 1/6;
    id b= -1/2*c;
    id c= -8/15;
    id f= 30*(- 2/15*d - 1/45*e);
    id d= 1/3*e;
    id e= 3/20;
    *
    id y^2= a*f00+b*f01+c*f10+d*f02+e*f20+f*f11;
    id a= 1/6;
    id b= -2*c;
    id f= -3/2*(d+e);
    id c= 4/15;
    id e= 1/3*d;
    id d= 3/20;
    *
    id x*y= a*f00+b*f01+c*f10+d*f02+e*f20+f*f11;
    id a= 1/12;
    id c= -(1/5+2*b);
    id b= -1/15;
    id f= -30*(2/15*d+1/45*e);
    id e= d;
    id d= -27/360;
    *
    id x= 1/3+1/3*f01-2/3*f10;
    id y= 1/3-2/3*f01+1/3*f10;
    *
    print;
    *
    .store

   R1 =
      4/3 - 20/3*f01 + 4/3*f10;

   R2 =
       - 12*f01;

   R3 =
       - 12*f10;

   R4 =
      4/3 + 4/3*f01 - 20/3*f10;

   R5 =
      2/3 - 10/3*f01 - 4/3*f10;

   R6 =
      2/3 - 4/3*f01 - 10/3*f10;

   R7 =
      24;

   R8 =
      12;

   R9 =
      4;

   R10 =
      4;

   R11 =
      12;

   R12 =
      24;

   R13 =
      6;

   R14 =
      8;

   R15 =
      6;

   R16 =
      42 - 45*f00 + 84/5*f01 - 24/5*f10 - 171/10*f02 - 9/10*f20 + 9*f11;

   R17 =
      96 - 102*f00 + 24/5*f01 - 48/5*f10 - 153/5*f02 - 9/5*f20 + 18/5*f11;

   R18 =
       - 24 + 84*f01 - 24*f10;

   R19 =
       - 24 + 192*f01 - 24*f10;

   R20 =
      360*f01;

   R21 =
      96 - 102*f00 - 48/5*f01 + 24/5*f10 - 9/5*f02 - 153/5*f20 + 18/5*f11;

   R22 =
      42 - 45*f00 - 24/5*f01 + 84/5*f10 - 9/10*f02 - 171/10*f20 + 9*f11;

   R23 =
      360*f10;

   R24 =
       - 24 - 24*f01 + 192*f10;

   R25 =
       - 24 - 24*f01 + 84*f10;

   R26 =
      24 - 25*f00 + 16/5*f01 + 4/5*f10 - 39/10*f02 + 9/10*f20 - 9*f11;

   R27 =
      30 - 31*f00 + 16/5*f01 - 8/5*f10 - 9/10*f02 - 3/10*f20 - 81/5*f11;

   R28 =
       - 16 + 56*f01 - 4*f10;

   R29 =
       - 12 + 60*f01 + 24*f10;

   R30 =
       - 16 + 32*f01 + 56*f10;

   R31 =
      30 - 31*f00 - 8/5*f01 + 16/5*f10 - 3/10*f02 - 9/10*f20 - 81/5*f11;

   R32 =
      24 - 25*f00 + 4/5*f01 + 16/5*f10 + 9/10*f02 - 39/10*f20 - 9*f11;

   R33 =
       - 16 + 56*f01 + 32*f10;

   R34 =
       - 12 + 24*f01 + 60*f10;

   R35 =
       - 16 - 4*f01 + 56*f10;


    *
    g P1= F01*F01;
    g P2= F01*F10;
    g P3= F01*F02;
    g P4= F01*F20;
    g P5= F01*F11;
    *
    g P6= F10*F01;
    g P7= F10*F10;
    g P8= F10*F02;
    g P9= F10*F20;
    g P10= F10*F11;
    nwrite statistics;
    *
    id x^3= a*f00+b*f01+c*f10+d*f02+e*f20+f*f11+g*f03+h*f30+l*f12+m*f21;
    id a= 1/10;
    id b= -1/2*c;
    id d= -15/2*(e/45+f/30);
    id g= 14/9*( - 9/280*h - 3/28*l - 3/70*m);
    id f= -72/36*e;
    id c= -16/40;
    id e= 27/140;
    id h= 80/3*( - 1/40*l - 1/40*m);
    id l= -200/9*3/100*m;
    id m= 16*3/560;
    *
    id y^3= a*f00+b*f01+c*f10+d*f02+e*f20+f*f11+g*f03+h*f30+l*f12+m*f21;
    id a= 1/10;
    id b= -2*c;
    id d= 45*( - 2/15*e - 1/30*f);
    id g= 70/3*( - 3/28*h - 9/280*l - 3/70*m);
    id c= 8/40;
    id e= -6/36*f;
    id h= -40/9*( + 3/80*l + 3/40*m);
    id f= -27/70;
    id l= -40/3*9/80*m;
    id m= -32/3*3/560;
    *
    id x^2*y= a*f00+b*f01+c*f10+d*f02+e*f20+f*f11+g*f03+h*f30+l*f12+m*f21;
    id a= 1/30;
    id c= 24*( - 1/360 - 1/12*b);
    id d= 15/2*( - 1/45*e - 1/30*f);
    id g= 14/9*( - 9/280*h - 3/28*l - 3/70*m);
    id b= 0;
    id f= 36*(1/630-7/54*e);
    id h= 800/513*( - 3/80*l - 21/200*m);
    id e= -27/1260;
    id l= -570/13*43/2280*m;
    id m= -416/3/1680;
    *
    id x*y^2= a*f00+b*f01+c*f10+d*f02+e*f20+f*f11+g*f03+h*f30+l*f12+m*f21;
    id a= 1/30;
    id c= -24*(1/180+b/12);
    id d= 45*( - 2/15*e - 1/30*f);
    id g= 14/9*( - 9/280*h - 3/28*l - 3/70*m);
    id b= -8/120;
    id f= -6*( 1/630 + 7/9*e);
    id h= 800/513*( - 3/80*l - 21/200*m);
    id e= -27/756;
    id l= -2280/43*13/570*m;
    id m= 344/3/1680;
    *
    id x^2= a*f00+b*f01+c*f10+d*f02+e*f20+f*f11;
    id a= 1/6;
    id b= -1/2*c;
    id c= -8/15;
    id f= 30*(- 2/15*d - 1/45*e);
    id d= 1/3*e;
    id e= 3/20;
    *
    id y^2= a*f00+b*f01+c*f10+d*f02+e*f20+f*f11;
    id a= 1/6;
    id b= -2*c;
    id f= -3/2*(d+e);
    id c= 4/15;
    id e= 1/3*d;
    id d= 3/20;
    *
    id x*y= a*f00+b*f01+c*f10+d*f02+e*f20+f*f11;
    id a= 1/12;
    id c= -(1/5+2*b);
    id b= -1/15;
    id f= -30*(2/15*d+1/45*e);
    id e= d;
    id d= -27/360;
    *
    id x= 1/3+1/3*f01-2/3*f10;
    id y= 1/3-2/3*f01+1/3*f10;
    *
    print;
    *
    .store

   P1 =
       - 1 + 7/6*f00 - 2/15*f01 + 4/15*f10 + 7/20*f02 + 1/20*f20 - 1/10*f11;

   P2 =
       - 1 + 13/12*f00 + 2/15*f01 + 2/15*f10 + 1/40*f02 + 1/40*f20 + 11/20*f11
      ;

   P3 =
       - 16/3 + 16/3*f00 + 8/15*f01 - 4/35*f02 + 16/35*f11 + 13/63*f03 - 4/63*
      f12 + 4/63*f21;

   P4 =
       - 16/3 + 16/3*f00 + 4/15*f10 + 8/35*f20 + 8/35*f11 + 2/63*f30 + 2/63*
      f12 + 25/63*f21;

   P5 =
       - 8/3 + 8/3*f00 + 1/15*f01 + 2/15*f10 + 2/35*f02 + 4/35*f20 - 4/35*f11
       + 1/126*f03 + 1/63*f30 + 17/63*f12 - 1/18*f21;

   P6 =
       - 1 + 13/12*f00 + 2/15*f01 + 2/15*f10 + 1/40*f02 + 1/40*f20 + 11/20*f11
      ;

   P7 =
       - 1 + 7/6*f00 + 4/15*f01 - 2/15*f10 + 1/20*f02 + 7/20*f20 - 1/10*f11;

   P8 =
       - 16/3 + 16/3*f00 + 4/15*f01 + 8/35*f02 + 8/35*f11 + 2/63*f03 + 25/63*
      f12 + 2/63*f21;

   P9 =
       - 16/3 + 16/3*f00 + 8/15*f10 - 4/35*f20 + 16/35*f11 + 13/63*f30 + 4/63*
      f12 - 4/63*f21;

   P10 =
       - 8/3 + 8/3*f00 + 2/15*f01 + 1/15*f10 + 4/35*f02 + 2/35*f20 - 4/35*f11
       + 1/63*f03 + 1/126*f30 - 1/18*f12 + 17/63*f21;


    *
    g neq1= intx*inty*F00*(x^3-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11-
                           g*F03-h*F30-l*F12-m*F21);
    g neq2= intx*inty*F01*(x^3-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11-
                           g*F03-h*F30-l*F12-m*F21);
    g neq3= intx*inty*F10*(x^3-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11-
                           g*F03-h*F30-l*F12-m*F21);
    g neq4= intx*inty*F02*(x^3-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11-
                           g*F03-h*F30-l*F12-m*F21);
    g neq5= intx*inty*F20*(x^3-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11-
                           g*F03-h*F30-l*F12-m*F21);
    g neq6= intx*inty*F11*(x^3-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11-
                           g*F03-h*F30-l*F12-m*F21);
    g neq7= intx*inty*F03*(x^3-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11-
                           g*F03-h*F30-l*F12-m*F21);
    g neq8= intx*inty*F30*(x^3-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11-
                           g*F03-h*F30-l*F12-m*F21);
    g neq9= intx*inty*F12*(x^3-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11-
                           g*F03-h*F30-l*F12-m*F21);
    g neq10= intx*inty*F21*(x^3-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11-
                           g*F03-h*F30-l*F12-m*F21);
    nwrite statistics;
    *
    id inty*y^n?= (1-x)^(n+1)/(n+1);
    id intx*x^n?= 1/(n+1);
    *
    id a= 1/10;
    id b= -1/2*c;
    id d= -15/2*(e/45+f/30);
    id g= 14/9*( - 9/280*h - 3/28*l - 3/70*m);
    id f= -72/36*e;
    id c= -16/40;
    id e= 27/140;
    id h= 80/3*( - 1/40*l - 1/40*m);
    id l= -200/9*3/100*m;
    id m= 16*3/560;
    *
    print;
    *
    .store

   neq1 = 0;

   neq2 = 0;

   neq3 = 0;

   neq4 = 0;

   neq5 = 0;

   neq6 = 0;

   neq7 = 0;

   neq8 = 0;

   neq9 = 0;

   neq10 = 0;


    *
    g neq11= intx*inty*F00*(y^3-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11-
                           g*F03-h*F30-l*F12-m*F21);
    g neq12= intx*inty*F01*(y^3-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11-
                           g*F03-h*F30-l*F12-m*F21);
    g neq13= intx*inty*F10*(y^3-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11-
                           g*F03-h*F30-l*F12-m*F21);
    g neq14= intx*inty*F02*(y^3-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11-
                           g*F03-h*F30-l*F12-m*F21);
    g neq15= intx*inty*F20*(y^3-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11-
                           g*F03-h*F30-l*F12-m*F21);
    g neq16= intx*inty*F11*(y^3-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11-
                           g*F03-h*F30-l*F12-m*F21);
    g neq17= intx*inty*F03*(y^3-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11-
                           g*F03-h*F30-l*F12-m*F21);
    g neq18= intx*inty*F30*(y^3-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11-
                           g*F03-h*F30-l*F12-m*F21);
    g neq19= intx*inty*F12*(y^3-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11-
                           g*F03-h*F30-l*F12-m*F21);
    g neq20= intx*inty*F21*(y^3-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11-
                           g*F03-h*F30-l*F12-m*F21);
    nwrite statistics;
    *
    id inty*y^n?= (1-x)^(n+1)/(n+1);
    id intx*x^n?= 1/(n+1);
    *
    id a= 1/10;
    id b= -2*c;
    id d= 45*( - 2/15*e - 1/30*f);
    id g= 70/3*( - 3/28*h - 9/280*l - 3/70*m);
    id c= 8/40;
    id e= -6/36*f;
    id h= -40/9*( + 3/80*l + 3/40*m);
    id f= -27/70;
    id l= -40/3*9/80*m;
    id m= -32/3*3/560;
    *
    print;
    *
    .store

   neq11 = 0;

   neq12 = 0;

   neq13 = 0;

   neq14 = 0;

   neq15 = 0;

   neq16 = 0;

   neq17 = 0;

   neq18 = 0;

   neq19 = 0;

   neq20 = 0;


    *
    g neq21= intx*inty*F00*(x^2*y-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11-
                           g*F03-h*F30-l*F12-m*F21);
    g neq22= intx*inty*F01*(x^2*y-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11-
                           g*F03-h*F30-l*F12-m*F21);
    g neq23= intx*inty*F10*(x^2*y-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11-
                           g*F03-h*F30-l*F12-m*F21);
    g neq24= intx*inty*F02*(x^2*y-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11-
                           g*F03-h*F30-l*F12-m*F21);
    g neq25= intx*inty*F20*(x^2*y-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11-
                           g*F03-h*F30-l*F12-m*F21);
    g neq26= intx*inty*F11*(x^2*y-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11-
                           g*F03-h*F30-l*F12-m*F21);
    g neq27= intx*inty*F03*(x^2*y-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11-
                           g*F03-h*F30-l*F12-m*F21);
    g neq28= intx*inty*F30*(x^2*y-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11-
                           g*F03-h*F30-l*F12-m*F21);
    g neq29= intx*inty*F12*(x^2*y-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11-
                           g*F03-h*F30-l*F12-m*F21);
    g neq30= intx*inty*F21*(x^2*y-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11-
                           g*F03-h*F30-l*F12-m*F21);
    nwrite statistics;
    *
    id inty*y^n?= (1-x)^(n+1)/(n+1);
    id intx*x^n?= 1/(n+1);
    *
    id a= 1/30;
    id c= 24*( - 1/360 - 1/12*b);
    id d= 15/2*( - 1/45*e - 1/30*f);
    id g= 14/9*( - 9/280*h - 3/28*l - 3/70*m);
    id b= 0;
    id f= 36*(1/630-7/54*e);
    id h= 800/513*( - 3/80*l - 21/200*m);
    id e= -27/1260;
    id l= -570/13*43/2280*m;
    id m= -416/3/1680;
    *
    print;
    *
    .store

   neq21 = 0;

   neq22 = 0;

   neq23 = 0;

   neq24 = 0;

   neq25 = 0;

   neq26 = 0;

   neq27 = 0;

   neq28 = 0;

   neq29 = 0;

   neq30 = 0;


    *
    g neq31= intx*inty*F00*(x*y^2-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11-
                           g*F03-h*F30-l*F12-m*F21);
    g neq32= intx*inty*F01*(x*y^2-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11-
                           g*F03-h*F30-l*F12-m*F21);
    g neq33= intx*inty*F10*(x*y^2-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11-
                           g*F03-h*F30-l*F12-m*F21);
    g neq34= intx*inty*F02*(x*y^2-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11-
                           g*F03-h*F30-l*F12-m*F21);
    g neq35= intx*inty*F20*(x*y^2-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11-
                           g*F03-h*F30-l*F12-m*F21);
    g neq36= intx*inty*F11*(x*y^2-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11-
                           g*F03-h*F30-l*F12-m*F21);
    g neq37= intx*inty*F03*(x*y^2-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11-
                           g*F03-h*F30-l*F12-m*F21);
    g neq38= intx*inty*F30*(x*y^2-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11-
                           g*F03-h*F30-l*F12-m*F21);
    g neq39= intx*inty*F12*(x*y^2-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11-
                           g*F03-h*F30-l*F12-m*F21);
    g neq40= intx*inty*F21*(x*y^2-a*F00-b*F01-c*F10-d*F02-e*F20-f*F11-
                           g*F03-h*F30-l*F12-m*F21);
    nwrite statistics;
    *
    id inty*y^n?= (1-x)^(n+1)/(n+1);
    id intx*x^n?= 1/(n+1);
    *
    id a= 1/30;
    id c= -24*(1/180+b/12);
    id d= 45*( - 2/15*e - 1/30*f);
    id g= 14/9*( - 9/280*h - 3/28*l - 3/70*m);
    id b= -8/120;
    id f= -6*( 1/630 + 7/9*e);
    id h= 800/513*( - 3/80*l - 21/200*m);
    id e= -27/756;
    id l= -2280/43*13/570*m;
    id m= 344/3/1680;
    *
    print;
    *
    .store

   neq31 = 0;

   neq32 = 0;

   neq33 = 0;

   neq34 = 0;

   neq35 = 0;

   neq36 = 0;

   neq37 = 0;

   neq38 = 0;

   neq39 = 0;

   neq40 = 0;


    *
    g Quad= A*x^2+B*y^2+C*x*y+D*x+E*y+F;
    g Cub= A*x^3+B*y^3+C*x^2*y+D*x*y^2+E*x^2+F*y^2+G*x*y+H*x+I*y+L;
    
    
    nwrite statistics;
    *
    id x^3= a*f00+b*f01+c*f10+d*f02+e*f20+f*f11+g*f03+h*f30+l*f12+m*f21;
    id a= 1/10;
    id b= -1/2*c;
    id d= -15/2*(e/45+f/30);
    id g= 14/9*( - 9/280*h - 3/28*l - 3/70*m);
    id f= -72/36*e;
    id c= -16/40;
    id e= 27/140;
    id h= 80/3*( - 1/40*l - 1/40*m);
    id l= -200/9*3/100*m;
    id m= 16*3/560;
    *
    id y^3= a*f00+b*f01+c*f10+d*f02+e*f20+f*f11+g*f03+h*f30+l*f12+m*f21;
    id a= 1/10;
    id b= -2*c;
    id d= 45*( - 2/15*e - 1/30*f);
    id g= 70/3*( - 3/28*h - 9/280*l - 3/70*m);
    id c= 8/40;
    id e= -6/36*f;
    id h= -40/9*( + 3/80*l + 3/40*m);
    id f= -27/70;
    id l= -40/3*9/80*m;
    id m= -32/3*3/560;
    *
    id x^2*y= a*f00+b*f01+c*f10+d*f02+e*f20+f*f11+g*f03+h*f30+l*f12+m*f21;
    id a= 1/30;
    id c= 24*( - 1/360 - 1/12*b);
    id d= 15/2*( - 1/45*e - 1/30*f);
    id g= 14/9*( - 9/280*h - 3/28*l - 3/70*m);
    id b= 0;
    id f= 36*(1/630-7/54*e);
    id h= 800/513*( - 3/80*l - 21/200*m);
    id e= -27/1260;
    id l= -570/13*43/2280*m;
    id m= -416/3/1680;
    *
    id x*y^2= a*f00+b*f01+c*f10+d*f02+e*f20+f*f11+g*f03+h*f30+l*f12+m*f21;
    id a= 1/30;
    id c= -24*(1/180+b/12);
    id d= 45*( - 2/15*e - 1/30*f);
    id g= 14/9*( - 9/280*h - 3/28*l - 3/70*m);
    id b= -8/120;
    id f= -6*( 1/630 + 7/9*e);
    id h= 800/513*( - 3/80*l - 21/200*m);
    id e= -27/756;
    id l= -2280/43*13/570*m;
    id m= 344/3/1680;
    *
    id x^2= a*f00+b*f01+c*f10+d*f02+e*f20+f*f11;
    id a= 1/6;
    id b= -1/2*c;
    id c= -8/15;
    id f= 30*(- 2/15*d - 1/45*e);
    id d= 1/3*e;
    id e= 3/20;
    *
    id y^2= a*f00+b*f01+c*f10+d*f02+e*f20+f*f11;
    id a= 1/6;
    id b= -2*c;
    id f= -3/2*(d+e);
    id c= 4/15;
    id e= 1/3*d;
    id d= 3/20;
    *
    id x*y= a*f00+b*f01+c*f10+d*f02+e*f20+f*f11;
    id a= 1/12;
    id c= -(1/5+2*b);
    id b= -1/15;
    id f= -30*(2/15*d+1/45*e);
    id e= d;
    id d= -27/360;
    *
    id x= 1/3*f00+1/3*f01-2/3*f10;
    id y= 1/3*f00-2/3*f01+1/3*f10;
    *
    id F= f00;
    *
    id f00^2= 1;
    id f00*dmy?= dmy;
    *
    b f00,f01,f10,f02,f20,f11,f03,f30,f12,f21,f03,f30,f12,f21;
    print;
    *
    .end

   Quad =
       + f00 * ( 1 )

       + f01 * ( 4/15*A - 8/15*B - 1/15*C + 1/3*D - 2/3*E )

       + f10 * (  - 8/15*A + 4/15*B - 1/15*C - 2/3*D + 1/3*E )

       + f02 * ( 1/20*A + 3/20*B - 3/40*C )

       + f20 * ( 3/20*A + 1/20*B - 3/40*C )

       + f11 * (  - 3/10*A - 3/10*B + 7/20*C )

       + 1/6*A + 1/6*B + 1/12*C + 1/3*D + 1/3*E;

   Cub =
       + f01 * (  - 8/15 + 1/5*A - 2/5*B - 1/15*D + 4/15*E - 1/15*G + 1/3*H - 
         2/3*I )

       + f10 * ( 4/15 - 2/5*A + 1/5*B - 1/15*C - 8/15*E - 1/15*G - 2/3*H + 1/3
         *I )

       + f02 * ( 3/20 + 9/140*A + 27/140*B - 1/28*C - 3/140*D + 1/20*E - 3/40*
         G )

       + f20 * ( 1/20 + 27/140*A + 9/140*B - 3/140*C - 1/28*D + 3/20*E - 3/40*
         G )

       + f11 * (  - 3/10 - 27/70*A - 27/70*B + 11/70*C + 11/70*D - 3/10*E + 7/
         20*G )

       + f03 * ( 1/210*A - 2/105*B - 2/315*C + 1/105*D )

       + f30 * (  - 2/105*A + 1/210*B + 1/105*C - 2/315*D )

       + f12 * (  - 2/35*A + 3/35*B + 43/630*C - 26/315*D )

       + f21 * ( 3/35*A - 2/35*B - 26/315*C + 43/630*D )

       + 1/6 + 1/10*A + 1/10*B + 1/30*C + 1/30*D + 1/6*E + 1/12*G + 1/3*H + 1/
         3*I + L;