Pensieve header: Finding the YB element for NOE-1t.

SetDirectory["C:\\drorbn\\AcademicPensieve\\Projects\\OneCo-1606"];
<< NOE-1t.m

AddRule[ff_,rule_]:=(
Print["As ",e0,", adding ",rule];
done=False; EC=EC⋃{rule}
);
MF[ϕ_,v_]:=Module[{t=ϕ,t1}, If[Simplify[t]===0, 1,
While[(t1=Simplify[t/.v->0])===0, t=D[t,v]]; t1
]];

RunMess:=(
done=False; While[!done, done=True;
E0=DeleteCases[Simplify[E0//.EC], 0] // SortBy[LeafCount];
Print["Length[E0]==", Length[E0], "; Length[EC]==", Length[EC], "; {gn,hn}==", {gn,hn}];
For[k=1, k≤Length@E0, ++k,
e1=Factor[e0=E0[[k]]];
If[Head[e1]=!=Times, e2=e1,
E0[[k]]=e2=Select[e1, !FreeQ[#, (Subscript[f, _])[_,_]|(Subscript[g, _])[_]|(Subscript[h, _])[]]&] ];
If[e2==1,Print["Panic at ",e0,"! No solutions."]; Break[]];
If[!FreeQ[e2,(Subscript[f, _])[_,_]]∧(FreeQ[e2,x]∨FreeQ[e2,y]∨FreeQ[e2,z]),
{ff}=Cases[e2, (Subscript[f, _])[_,_], {0,∞}, 1];
{{sol}}=Solve[e2==0, ff];
rule=((ff/.{x-> x_, y->y_,z->z_}) -> (ff/.sol)) /. Rule->RuleDelayed;
AddRule[ff,rule]; Break[]
];
If[!FreeQ[e2,(Subscript[g, _])[_]]∧(FreeQ[e2,y|z]∨FreeQ[e2,x|z]∨FreeQ[e2,z|y]),
{gg}=Cases[e2, (Subscript[g, _])[_], ∞, 1];
{{sol}}=Solve[e2==0, gg];
rule=((gg/.{x-> x_, y->y_,z->z_}) -> (gg/.sol)) /. Rule->RuleDelayed;
AddRule[gg,rule]; Break[]
];
If[Head[e2]===Plus,
s=List@@Collect[e2, (Subscript[f, _])[_,_], Factor]; s1=Select[s, FreeQ[(Subscript[f, _])[_,_]]];
sxy=Cases[s,a_.*(Subscript[f, _])[x,y]]; sxz=Cases[s,a_.*(Subscript[f, _])[x,z]]; syz=Cases[s,a_.*(Subscript[f, _])[y,z]];
Which[
sxy=={}∧sxz=!={}∧syz=!={}, (
{ff}=Cases[sxz, a_.*(Subscript[f, k_])[x,z]:>(Subscript[f, k])[x,z],{1},1];
mf=MF[First@sxz /. (Subscript[f, _])[x,z]->1,x]; mf*=MF[First@syz /. (Subscript[f, _])[y,z]->1,y];
s1=Plus@@Simplify[s1/mf]; sxz=Plus@@Simplify[sxz/mf]; syz=Plus@@Simplify[syz/mf];
If[FreeQ[sxz,y]∧FreeQ[syz,x]∧FreeQ[s1,x|y]∧Simplify[(sxz/.x->y)+syz ==0],
{{sol}}=Solve[sxz==(Subscript[g, ++gn])[z], ff];
rule=((ff/.{x-> x_, y->y_,z->z_}) -> (ff/.sol)) /. Rule->RuleDelayed;
AddRule[ff,rule]; Break[]
]
),
syz=={}∧sxy=!={}∧sxz=!={}, (
{ff}=Cases[sxy, a_.*(Subscript[f, k_])[x,y]:>(Subscript[f, k])[x,y],{1},1];
mf=MF[First@sxy /. (Subscript[f, _])[x,y]->1,y]; mf*=MF[First@sxz /. (Subscript[f, _])[x,z]->1,z];
s1=Plus@@Simplify[s1/mf]; sxy=Plus@@Simplify[sxy/mf]; sxz=Plus@@Simplify[sxz/mf];
If[FreeQ[sxy,z]∧FreeQ[sxz,y]∧FreeQ[s1,y|z]∧Simplify[(sxz/.z->y)+sxy ==0],
{{sol}}=Solve[sxy==(Subscript[g, ++gn])[x], ff];
rule=((ff/.{x-> x_, y->y_,z->z_}) -> (ff/.sol)) /. Rule->RuleDelayed;
AddRule[ff,rule]; Break[]
]
),
sxy=!={}∧sxz=!={}∧syz=!={}, (
kk=Union@Cases[e2, a_.*(Subscript[f, k_])[x,y]:>k,∞];
If[Length[kk]==1,
{kk}=kk;
{{sol}}=Solve[e2==0, (Subscript[f, kk])[x,y]];
sol=(Subscript[f, kk])[x,y]/.sol;
e3=D[sol,z]//Factor;
If[FreeQ[e3,(Subscript[f, _])[_,_]],
If[Head[e3]===Times, e3=Select[e3, !FreeQ[#, Derivative[0,1][(Subscript[f, _])][_,_]|(Subscript[g, _])[_]|(Subscript[h, _])[]]&] ];
s=Collect[e3, Derivative[0,1][(Subscript[f, _])][_,_], Factor];
s1=Select[s, FreeQ[Derivative[0,1][(Subscript[f, _])][_,_]]];
pxz=Coefficient[s,Derivative[0,1][(Subscript[f, kk])][x,z]];
pyz=Coefficient[s,Derivative[0,1][(Subscript[f, kk])][y,z]];
mf=MF[pxz,x]; mf*=MF[pyz,y];
{s1,pxz,pyz}=Simplify[{s1,pxz,pyz}/mf];
If[FreeQ[pxz,y]∧FreeQ[pyz,x]∧FreeQ[s1,x|y]∧Simplify[(pyz/.y->x)+pxz ==0],
rule=((Subscript[f, kk])[x_,z_]->(Subscript[g, ++gn])[z]/pxz+(Subscript[g, ++gn])[x]) /. Rule->RuleDelayed;
AddRule[(Subscript[f, kk])[x,z],rule]; Break[]
]
]
]
)
];
If[FreeQ[e2, (Subscript[f, _])[_,_]]∧!FreeQ[e2,(Subscript[g, _])[_]],
s=List@@Collect[e2, (Subscript[g, _])[_], Factor]; s1=Select[s, FreeQ[(Subscript[g, _])[_]]];
sx=Cases[s,a_.*(Subscript[g, _])[x]]; sy=Cases[s,a_.*(Subscript[g, _])[y]]; sz=Cases[s,a_.*(Subscript[g, _])[z]];
Which[
FreeQ[e2,x]∧sy=!={}∧sz=!={}, (
{gg}=Cases[sy, a_.*(Subscript[g, k_])[y]:>(Subscript[g, k])[y],{1},1];
mf=MF[First@sy /. (Subscript[g, _])[y]->1,y]; mf*=MF[First@sz /. (Subscript[g, _])[z]->1,z];
s1=Plus@@Simplify[s1/mf]; sy=Plus@@Simplify[sy/mf]; sz=Plus@@Simplify[sz/mf];
If[FreeQ[sx,y]∧FreeQ[sz,y]∧FreeQ[s1,y|z]∧Simplify[(sz/.z->y)+sy ==0],
{{sol}}=Solve[sy==(Subscript[h, ++hn])[], gg];
rule=((gg/.{x-> x_, y->y_,z->z_}) -> (gg/.sol)) /. Rule->RuleDelayed;
AddRule[gg,rule]; Break[]
]
),
FreeQ[e2,z]∧sy=!={}∧sx=!={}, (
{gg}=Cases[sy, a_.*(Subscript[g, k_])[y]:>(Subscript[g, k])[y],{1},1];
mf=MF[First@sy /. (Subscript[g, _])[y]->1,y]; mf*=MF[First@sx /. (Subscript[g, _])[x]->1,x];
s1=Plus@@Simplify[s1/mf]; sy=Plus@@Simplify[sy/mf]; sx=Plus@@Simplify[sx/mf];
If[FreeQ[sz,y]∧FreeQ[sx,y]∧FreeQ[s1,y|x]∧Simplify[(sx/.x->y)+sy ==0],
{{sol}}=Solve[sy==(Subscript[h, ++hn])[], gg];
rule=((gg/.{x-> x_, y->y_,z->z_}) -> (gg/.sol)) /. Rule->RuleDelayed;
AddRule[gg,rule]; Break[]
]
)
] (* Which *)
] (* If *)
] (* If *)
] (* For *)
]; (* While *)
E0=Union[DeleteCases[Simplify[E0//.EC], 0]] // SortBy[LeafCount];
)

Rp[i_,j_]:=𝔼[1, (Subscript[b, i])(Subscript[c, j]),(Subscript[u, i])(Subscript[w, j]),(Subscript[c, i]) (Subscript[f, 1])[(Subscript[t, i]),(Subscript[t, j])]+(Subscript[c, j]) (Subscript[f, 2])[(Subscript[t, i]),(Subscript[t, j])]+(Subscript[c, i]^2) (Subscript[f, 3])[(Subscript[t, i]),(Subscript[t, j])]+(Subscript[c, i]) (Subscript[c, j]) (Subscript[f, 4])[(Subscript[t, i]),(Subscript[t, j])]+(Subscript[c, j]^2) (Subscript[f, 5])[(Subscript[t, i]),(Subscript[t, j])]+(Subscript[u, i]) (Subscript[w, i]) (Subscript[f, 6])[(Subscript[t, i]),(Subscript[t, j])]+(Subscript[u, i]) (Subscript[w, j]) (Subscript[f, 7])[(Subscript[t, i]),(Subscript[t, j])]+(Subscript[u, j]) (Subscript[w, i]) (Subscript[f, 8])[(Subscript[t, i]),(Subscript[t, j])]+(Subscript[u, j]) (Subscript[w, j]) (Subscript[f, 9])[(Subscript[t, i]),(Subscript[t, j])]+(Subscript[c, i]) (Subscript[u, i]) (Subscript[w, i]) (Subscript[f, 10])[(Subscript[t, i]),(Subscript[t, j])]+(Subscript[c, i]) (Subscript[u, i]) (Subscript[w, j]) (Subscript[f, 11])[(Subscript[t, i]),(Subscript[t, j])]+(Subscript[c, i]) (Subscript[u, j]) (Subscript[w, i]) (Subscript[f, 12])[(Subscript[t, i]),(Subscript[t, j])]+(Subscript[c, i]) (Subscript[u, j]) (Subscript[w, j]) (Subscript[f, 13])[(Subscript[t, i]),(Subscript[t, j])]+(Subscript[c, j]) (Subscript[u, i]) (Subscript[w, i]) (Subscript[f, 14])[(Subscript[t, i]),(Subscript[t, j])]+(Subscript[c, j]) (Subscript[u, i]) (Subscript[w, j]) (Subscript[f, 15])[(Subscript[t, i]),(Subscript[t, j])]+(Subscript[c, j]) (Subscript[u, j]) (Subscript[w, i]) (Subscript[f, 16])[(Subscript[t, i]),(Subscript[t, j])]+(Subscript[c, j]) (Subscript[u, j]) (Subscript[w, j]) (Subscript[f, 17])[(Subscript[t, i]),(Subscript[t, j])]+(Subscript[u, i]^2) (Subscript[w, i]^2) (Subscript[f, 18])[(Subscript[t, i]),(Subscript[t, j])]+(Subscript[u, i]^2) (Subscript[w, i]) (Subscript[w, j]) (Subscript[f, 19])[(Subscript[t, i]),(Subscript[t, j])]+(Subscript[u, i]^2) (Subscript[w, j]^2) (Subscript[f, 20])[(Subscript[t, i]),(Subscript[t, j])]+(Subscript[u, i]) (Subscript[u, j]) (Subscript[w, i]^2) (Subscript[f, 21])[(Subscript[t, i]),(Subscript[t, j])]+(Subscript[u, i]) (Subscript[u, j]) (Subscript[w, i]) (Subscript[w, j]) (Subscript[f, 22])[(Subscript[t, i]),(Subscript[t, j])]+(Subscript[u, i]) (Subscript[u, j]) (Subscript[w, j]^2) (Subscript[f, 23])[(Subscript[t, i]),(Subscript[t, j])]+(Subscript[u, j]^2) (Subscript[w, i]^2) (Subscript[f, 24])[(Subscript[t, i]),(Subscript[t, j])]+(Subscript[u, j]^2) (Subscript[w, i]) (Subscript[w, j]) (Subscript[f, 25])[(Subscript[t, i]),(Subscript[t, j])]+(Subscript[u, j]^2) (Subscript[w, j]^2) (Subscript[f, 26])[(Subscript[t, i]),(Subscript[t, j])]]

Short[t1=Rp[1,2]Rp[3,4]Rp[5,6] // m[3,5,x]//m[1,6,y]//m[2,4,z]]

Short[t2=Rp[1,2]Rp[3,4]Rp[5,6]//m[1,3,x]//m[2,5,y]//m[4,6,z]]

Take[t1,3]==Take[t2,3]

Short[t3=t1[[4]]-t2[[4]]]

Short[Errors=CoefficientRules[Expand[t3], {(Subscript[c, x]),(Subscript[c, y]),(Subscript[c, z]),(Subscript[u, x]),(Subscript[u, y]),(Subscript[u, z]),(Subscript[w, x]),(Subscript[w, y]),(Subscript[w, z])}] /. {(_->c_) :> c} /. {(Subscript[t, x])->x, (Subscript[t, y])->y, (Subscript[t, z])->z}]

EC={}; E0=Errors; gn=hn=0;

RunMess

eq1=E0[[1]]

Select[eq1,FreeQ[x]]+Select[eq1,FreeQ[y]]-eq1

Solve[Select[eq1,FreeQ[y]]==(Subscript[g, ++gn])[z], (Subscript[f, 1])[x,z]]

rule1=((Subscript[f, 1])[x_,z_]:>(1/2) (-2 (Subscript[f, 6])[x,z]+2 x (Subscript[f, 6])[x,z]+2 (Subscript[f, 20])[x,z]-(Subscript[g, 5])[x]+4 (Subscript[g, 6])[x]-12 x (Subscript[g, 6])[x]+12 (x^2) (Subscript[g, 6])[x]-4 (x^3) (Subscript[g, 6])[x]+(Subscript[g, 8])[z]))

E0/.rule1

AppendTo[EC,rule1]

RunMess

eq2=E0[[1]]

eq3=Collect[(eq2/(2(((-1+y))^2))), (Subscript[f, _])[__],FullSimplify]

eq4=Collect[D[eq3, x],(Subscript[f, _])[__], Simplify]

(eq4/.(Subscript[f, _])[x,y]|Derivative[1,0][(Subscript[f, _])][x,y]->0)+(eq4/.(Subscript[f, _])[x,z]|Derivative[1,0][(Subscript[f, _])][x,z]->0) == eq4 // Simplify

FreeQ[(eq4/.(Subscript[f, _])[x,y]|Derivative[1,0][(Subscript[f, _])][x,y]->0),y]

FreeQ[(eq4/.(Subscript[f, _])[x,z]|Derivative[1,0][(Subscript[f, _])][x,z]->0),z]

eq5=Collect[(-eq4/.(Subscript[f, _])[x,y]|Derivative[1,0][(Subscript[f, _])][x,y]->0), (Subscript[f, 6])[x,z]|Derivative[1,0][(Subscript[f, 6])][x,z], Simplify]

D[(x-1)(Subscript[f, 6])[x,z]+(Subscript[f, 7])[x,z], x]-eq5

Solve[(x-1)(Subscript[f, 6])[x,z]+(Subscript[f, 7])[x,z]==(Subscript[g, ++gn])[x]+(Subscript[g, ++gn])[z], (Subscript[f, 7])[x,z]]

eq6=Simplify[eq2/.(Subscript[f, 7])[x_,z_]:>(Subscript[f, 6])[x,z]-x (Subscript[f, 6])[x,z]+(Subscript[g, 9])[x]+(Subscript[g, 10])[z]]

First@Solve[eq6==0,(Subscript[g, 10])[y]]

Simplify[E0 /. (Subscript[f, 7])[x_,z_]:>(Subscript[f, 6])[x,z]-x (Subscript[f, 6])[x,z]+(Subscript[g, 9])[x]+(Subscript[g, 10])[z] /. (Subscript[g, 10])[y_]:>(1/(2 (1-2 y+(y^2))))(1+(Subscript[g, 5])[y]-3 y (Subscript[g, 5])[y]+2 (y^2) (Subscript[g, 5])[y]-4 (Subscript[g, 6])[y]+16 y (Subscript[g, 6])[y]-24 (y^2) (Subscript[g, 6])[y]+16 (y^3) (Subscript[g, 6])[y]-4 (y^4) (Subscript[g, 6])[y]-2 y (Subscript[g, 7])[y]+2 (y^2) (Subscript[g, 7])[y]-(Subscript[g, 8])[y]+y (Subscript[g, 8])[y]-2 (Subscript[g, 9])[y]+4 y (Subscript[g, 9])[y]-2 (y^2) (Subscript[g, 9])[y]+y (Subscript[h, 1])[])]

EC=EC⋃{(Subscript[f, 7])[x_,z_]:>(Subscript[f, 6])[x,z]-x (Subscript[f, 6])[x,z]+(Subscript[g, 9])[x]+(Subscript[g, 10])[z], (Subscript[g, 10])[y_]:>(1/(2 (1-2 y+(y^2))))(1+(Subscript[g, 5])[y]-3 y (Subscript[g, 5])[y]+2 (y^2) (Subscript[g, 5])[y]-4 (Subscript[g, 6])[y]+16 y (Subscript[g, 6])[y]-24 (y^2) (Subscript[g, 6])[y]+16 (y^3) (Subscript[g, 6])[y]-4 (y^4) (Subscript[g, 6])[y]-2 y (Subscript[g, 7])[y]+2 (y^2) (Subscript[g, 7])[y]-(Subscript[g, 8])[y]+y (Subscript[g, 8])[y]-2 (Subscript[g, 9])[y]+4 y (Subscript[g, 9])[y]-2 (y^2) (Subscript[g, 9])[y]+y (Subscript[h, 1])[])}

Simplify[Errors //. EC]

Rp[1,2]//.EC