read(`morscrpt.txt`):
n:=4;xx:=[seq(x.i,i=1..n)];
s2fL:=[seq(x.1*x.i,i=1..n),seq(x.2*x.i,i=2..n)];mon(2,omit(xx,x1,x2));
s2f:=[op(s2fL),op(")];
q1:=s2fL[1]+sum('a.(i-2)*s2fL[i]','i'=3..nops(s2fL));
q2:=s2fL[2]+sum('b.(i-2)*s2fL[i]','i'=3..nops(s2fL));

vars0:=[seq(a.i,i=1..nops(s2fL)-2),seq(b.i,i=1..nops(s2fL)-2)];
f11:=quo(q1,x1,x1,'r');f12:=si(r/x2);
f21:=quo(q2,x1,x1,'r');f22:=si(r/x2);
q3:=mtaylor(f11*f22-f12*f21,xx,3);qs:=[q1,q2,q3]:

factor(origin([x3,x4],qs));
subs(x2=1,x1=-b3,");
incid_line:="[1];

M:=matrix([seq([seq(coefmon(q.i,s2f[k]),k=1..nops(s2f))],i=1..3)]):
rows_no_pivo(M):
pivo_M:="[2];rows_no_pivo_M:=""[1];M:="""[3]:
 eqs(M,rows_no_pivo_M,vars0):
 	rows_with_eqn_M:="[1]:"":
 	position_eqn_M:="[2]:"":
 	eqs_M:="[3]:"":
 	eqs_M_1:="[4]:

 prep(M,position_eqn_M,eqs_M, eqs_M_1, s2f):
 eqs_M:="[1]:s2f:=""[2]:M:=eval("""[3]):
 B:=sol(eqs_M, eqs_M_1);

for i in position_eqn_M do print(i,s2f[i[2]])od:
(factor(subs(B,qs)));
(origin(vars0,qs));

collect(rem(-b3*x2+x1+a2*x4+a1*x3,(x2+b1*x3+b2*x4),x2),x3);

i_0:=1;relexcn(eqs_M,eqs_M_1,'c',i_0 );
subs(seq(c[i]=c.i,i=2..nops(eqs_M)),");
rel_1:="[1]:
eqsM:=""[2]:
eqsM1:="""[3]:
exc1:=op(solve({eqsM[1]},{eqsM1[1]}));
vars1:=[eqsM1[1]]:
for i in vars0 do if member(i,{seq(eqsM1[ii],
ii=1..nops(eqsM))})=false then
[op("),i]fi:
od:
vars.1:=[op("),seq(c.i,i=2..nops(eqsM))];

 M:=subs(rel_1,copy(M)):
 local_triv(M,rows_with_eqn_M,eqsM[1]);
 rows_no_pivo(M):
 pivo_M:="[2];rows_no_pivo_M:=""[1];M:="""[3]:
 M:=si(M):s2f:=subs(rel_1,copy(s2f));
 qs:=[seq(sum('M[i,j]*s2f[j]','j'=1..coldim(M)),i=1..3)]:

omit(vars1,vars0);
rel_1:=collect(rel_1,"):

#save B,incid_line,vars0,eqsM1, s2f,s2fL,vars1,rel_1,qs,exc1,`vlast1.m`;

#read(`vlast1.m`);
subs(rel_1,incid_line);
incid_line1:=c3;
s3f:=fixorder(mon(3,xx),xx):


 mu:=matrix([ seq(seq(
 [seq(coefmon(qs[i]*x.j,s3f[k]),k=1..nops(s3f))],j=1..4),i=1..3)
 ]):

 rows_no_pivo(mu):
 pivo_mu:="[2];rows_no_pivo_mu:=""[1];mu:="""[3]:
 eqs(mu,rows_no_pivo_mu,vars1);
rows_with_eqn_mu:="[1];"":
position_eqn_mu:="[2]:"":
eqs_mu:="[3]:"":
eqs_mu_1:="[4]:

# save position_eqn_mu,s3f,incid_line1,
 #mu,rows_no_pivo_mu,rows_with_eqn_mu,eqs_mu,eqs_mu_1,`vlast2.m`;
#read(`vlast2.m`);

Y0:=sol(eqs_mu,eqs_mu_1):
eqs_mu_1:=[c3, c4, c5, a1, a2];#c4 replaces b3,a1 replaces b4
Y0:=sol(eqs_mu,eqs_mu_1):
mu0:=convert(row(mu,10),list);
factor(subs(Y0,mu0)):omit(",0);
mdc(");seq(""[i]/",i=1..nops(""));
Y2:=sol({op(Y0),""});Y1:=sol({op(Y0),b1=0,b2=0});
{}:for i in Y1 do if rhs(i)<>lhs(i)then {op("),i}fi od:Y1:=" ;
{}:for i in Y2 do if rhs(i)<>lhs(i)then {op("),i}fi od:Y2:=" ;

prep(mu,position_eqn_mu,eqs_mu, eqs_mu_1, s3f):
eqs_mu:="[1]:s3f:=""[2]:mu:=eval("""[3]):

eqs_mun:=[op(eqs_mu),b1,b2];eqs_mun_1:=[op(eqs_mu_1),b1,b2];
sol("",");

Y1 minus";
 position_eqn_mu;

[]:for i in position_eqn_mu do [op("),[si(mu[op(i)]),s3f[i[2]]]]od:
(origin(vars1,[op("),op(qs)]));


 i_0:=6;
 relexcn(eqs_mun,eqs_mun_1,'d',i_0 );
 subs(seq(d[i]=d.i,i=2..nops(eqs_mun)),");
 rel_2:="[1]:
 eqsmu:=""[2]:
 eqsmu1:="""[3]:
 exc2:=op(solve({eqsmu[1]},{eqsmu1[1]}));
 	vars2:=[eqsmu1[1]]:
 	for i in vars1 do if member(i,{seq(eqsmu1[ii],
 			ii=1..nops(eqsmu))})
 			=false then
 				[op("),i]fi:
 	od:
 	vars.2:=[op("),seq(d.i,i=2..nops(eqsmu))];

mun:=subs(rel_2,copy(mu)):
mu0:=factor(convert(row(mun,rows_with_eqn_mu[1]),list));

local_triv(mun,rows_with_eqn_mu,eqsmu[1]);
mu0:=factor(convert(row(mun,rows_with_eqn_mu[1]),list));

s3f:=subs(rel_2,s3f);

g:=mtaylor(sum('mu0[j]*s3f[j]','j'=1..nops(s3f)),xx,4);
 factor(subs(Y1,qs));
 newplane:=mdc(");newpt:=sol([seq(""[i]/",i=1..3)],xx);


# save Y1,Y2,Y0,mun,vars2,exc2,rel_2,eqsmu1,eqsmu,newpt,
 #newplane,g,s3f,`vlast3.m`;
#read(`vlast3.m`);

subs(rel_2,incid_line1);
incid_line2:=d2;

subs(rel_2,[seq(collect(lhs(Y2[i])-rhs(Y2[i]),b1),i=1..nops(Y2))]):
sature(",b1);

Y21:=sol(
[d4, d2, d3, d5, d6, b3*b1-b4],[d4, d2, d3, d5, d6, b4]);

subs(Y21,g);
qs:=subs(rel_2,qs):

 factor(subs(Y21,qs));
lastplane:= mdc(");
seq(""[i]/",i=1..3);
lastline:=seq("[i],i=1..2);
lastpt :=sol([""],xx);

o:=sol([ (lastline),x1,x2],xx);
si(subs(o,g));si(subs(lastpt,grad(g,xx)));

eqsn:=[b4-b1*b3, d2,d3, d4, d5, d6]:
eqsn_1:=[b4, d2,d3, d4, d5, d6]:


i_0:=3;
relexcn(eqsn,eqsn_1,'e',i_0 ):
subs(seq(e[i]=e.i,i=2..nops(eqsn)),");
rel_3:="[1]:
eqslast:=""[2]:
eqslast1:="""[3]:
exc3:=op(solve({eqslast[1]},{eqslast1[1]}));
	vars3:=[eqslast1[1]]:
	for i in vars2 do if member(i,{seq(eqslast1[ii],
			ii=1..nops(eqslast))})
			=false then
				[op("),i]fi:
	od:
	vars3:=[op("),seq(e.i,i=2..nops(eqslast))];
subs(rel_3, subs(rel_2,exc1));factor(lhs(")-rhs("));exc3;exc1;
subs(rel_3,exc2);

mu0:=factor(convert(row(mun,rows_with_eqn_mu[1]),list)):

mun0:=matrix([subs(rel_3,copy(mu0))]);

local_triv(mun0,[1],eqslast[1]);
eval(mun0);

qs:=ssubs(rel_3,qs):indets(qs)minus indets(vars3);

g:=mtaylor(sum('mun0[1,j]*s3f[j]','j'=1..nops(s3f)),xx,4);

#save vars3,exc3,rel_3,eqslast1,eqslast,g,qs,`vlast4.m`;
#read(`vlast4.m`);

mtaylor(collect(g,b1),xx,4);