-- SocDeg3-8-25-24.package
-- {\small \begin{verbatim}
needsPackage "PieriMaps"

-- The Package SocDeg3Package

-- The purpose of this package is to calculate the 
-- resolutions described in the paper 
--     "Artinian Gorenstein algebras of socle degree 
--      three which have the weak Lefschetz property"

-- The main script in SocDeg3Package is called SocDeg3.

-- Some of the other scripts of interest in this 
-- package are: 
--     BasisXY(r,d) (which gives the bases used in 
--       SocDeg3), 
--     tau(T) and theta(T) (which turn Tableaux into
--        algebraic objects),
--     EonE (which gives the action of an exterior 
--        algebra on its dual), 
--     SonD (which gives the action of the symmteric 
--        algebra on its dual), 
--     MakeResAnnXPhi and MakeResAnnPhi (which give 
--        the whole resolution), and
--     TestResAnnXPhi and TestResAnnPhi (which test 
--        the resolution).
-- All of these scripts are explained below.

-- **************************** 
-- *** The Script SocDeg3.  ***
-- ****************************

--- INPUT (Phi,r): where Phi is a homogeneous 
--  polynomial of degree three in a polynomial
--- ring over a field and r is an integer 
--  between 1 and d, whered is the number of variables. 

--- OUTPUT: If r = 1, then the output is (f11,f12).
---         If 2<= r <= d-1, then the output is 
---            Four matrices (A_r, B_r, C_r, D_r).
---         If r=d, then the output is (fd1,fd2).

---         The output is explained below.

-- Let Phi be a homogeneous polynomial of degree three in 
-- a polynomial ring P with d variables P_0,...,P_(d-1). 
-- View Phi as an element of the 
-- divided power algebra in the variables 
-- P_0^*,...,P_(d-1)^*.
--
-- The script requires that ann Phi does not contain any 
-- linear elements and that P_0 is a weak Lefschetz 
-- element on P/ann(Phi). An error message will be given 
-- if either of these requirments fail to hold. 

-- If r is an integer with 2<= r<= d-1, then  
-- the script SocDeg3 makes the matrices A_r, B_r, C_r, and 
-- D_r. The minimal resolution of P/ ann P_0(Phi)  has the 
-- matrix
--        A_r      P_0*B_r
--        P_0*C_r  D_r
-- from P\otimes K_1^{r-1}U_0 \oplus P\t L_2^{r-1} U_0 
-- to   P\otimes K_1^{r-2}U_0 \oplus P\t L_2^{r-2} U_0,
-- in position r for 2<= r <= d-1. 
--
-- The resolution (not necessarily minimal) of P/ann(Phi) 
-- has the matrix 
--        A_r         B_r                      (*)
--        P_0^2*C_r   D_r
-- in position r for 2<= r <= d-1.
-- The entries of A_r and D_r are homogeneous of degree one. 
-- The entries of B_r and C_r are homogeneous of degree 
-- zero. In particular, the minimal resolution of P/ann(Phi)
-- is obtained from (*) by splitting off a rank B_r summand 
-- of P\otimes L_2^{r-1} U_0  and  P\otimes K_1^{r-2}U_0, 
-- for each r.

-- If r is 1, then SocDeg3(Phi,1) produces (f11,f12).
-- The first map in the resolution of P/ ann P_0(Phi) is
--      P_0*f11    f12. 
-- The first map in the resolution of P/ ann(Phi) is 
--      P_0^2*f11  f12.

-- If r is d, then SocDeg3(Phi,d) produces (fd1,fd2).
-- The last map in the resolution of P/ ann P_0(Phi) is
--       fd1
--       P_0* fd2
-- The last map in the resolution of P/ ann (Phi) is
--       fd1
--       P_0^2*fd2

-- To produce the entire resolution all at once, use the 
-- commands MakeResAnnXPhi and/or MakeResAnnPhi. These 
-- commands are described below.



-- *******************************************
-- ** The function modifiedstandardTableaux **
-- *******************************************

-- INPUT (di,mu)
-- di is an Integer. It represents the number of labels 
-- to be used. mu is a partition, which is a list of 
-- positive integers written in descending order.
-- OUTPUT: 
-- All standard tableaux of shape mu and with labels from
-- 2..(di+1)
--
-- The function standardTableaux(di,mu) of the package 
-- PieriMaps produces all standard Tableaux of shape mu 
-- with labels from 0..(di-1). The function 
-- modifiedstandardTableaux applies standardTableaux and 
-- then adds 2 to each label.

-- A Tableaux with weakly increasing rows and strictly 
-- increasing columns iscalled a standard Tableaux.

modifiedstandardTableaux=(di,mu)->(
Plus2:={};
for i from 0 to #mu-1 do(
Plus2=append(Plus2,toList(mu_i:2)));
T:=standardTableaux(di,mu);
MST:={};
for i from 0 to #T-1 do (MST=append(MST,T_i+Plus2));
return(MST)
);

-- ****************************************
-- ** The bases of the source and target **
-- ** of A_r, B_r, C_r, and D_r          ** 
-- ****************************************

-- There are eight commands. They all have the form 
-- BasisXY(r,d), 
-- where X can be "A", "B", "C", or "D" and Y can be 
-- "Source" or "Target". In particular, we describe 
-- BasisASource(r,d).

-- INPUT: (r,d) where d is the number of variables 
-- and r is any integer between 2 and d-1.
--
-- OUTPUT: is the list of tableaux T for which  
--         {eta(k_T)| T in OUTPUT} form the basis for 
--         source A_r = K_1^{r-1} U_0.  

BasisASource=(r,d)->(
sourceShapeAr:=flatten(append({2}, toList ((d-r-1):1))); 
return (modifiedstandardTableaux(d-1,sourceShapeAr))
);

BasisATarget=(r,d)->(
targetShapeAr:=flatten(append({2}, toList ((d-r):1)));
return (modifiedstandardTableaux(d-1,targetShapeAr))
);

BasisBSource=(r,d)->(
sourceShapeBr:=flatten(append({2}, toList ((r-1):1))); 
return (modifiedstandardTableaux(d-1,sourceShapeBr))
);

BasisBTarget=(r,d)->(
targetShapeBr:=flatten(append({2}, toList ((d-r):1)));
return (modifiedstandardTableaux(d-1,targetShapeBr))
);

BasisCSource=(r,d)->(
sourceShapeCr:=flatten(append({2}, toList ((d-r-1):1))); 
return (modifiedstandardTableaux(d-1,sourceShapeCr))
);

BasisCTarget=(r,d)->(
targetShapeCr:=flatten(append({2}, toList ((r-2):1)));
return modifiedstandardTableaux(d-1,targetShapeCr))

BasisDSource=(r,d)->(
sourceShapeDr:=flatten(append({2}, toList ((r-1):1))); 
return(modifiedstandardTableaux(d-1,sourceShapeDr))
);

BasisDTarget=(r,d)->(
targetShapeDr:=flatten(append({2}, toList ((r-2):1)));
return(modifiedstandardTableaux(d-1,targetShapeDr))
);

 Basisf11source=(d)->BasisATarget(2,d);
 Basisf12source=(d)->BasisCTarget(2,d);

-- ***************************************
-- ** The functions tau(T) and theta(T) **
-- ***************************************


-- If T is a Hook Tableau with a rows and b columns, 
-- then k_T = tauT(omega)\otimes nuT in 
-- \bw^{d-a}U_0\otimes D_bU_0^*. The function tau(T) 
-- creates tauT. The answer tau(T) is a list of integers:
-- i_1,...,i_(a-1). This list represents 
-- X_{i_1}\wedge ... \wedge X_{i_(a-1)}.
-- The i's are in descending order.

-- In a similar manner, ell_T= thetaT\otimes uT in 
-- \bw^aU_0\otimes Sym_{b-1}U_0. The function theta(T) 
-- creates thetaT. The OUTPUT of theta(T) is a list 
-- of integers: i_1,...,i_a. This list represents
-- X_{i_1}\wedge ... \wedge X_{i_a}.
-- The i's are in ascending order.

tau=(T)->(
ans:={};
for i from 1 to #T-1 do ans=append(ans,T#i#0);
return(rsort ans)
);

theta=(T)->(
ans:={};
for i from 0 to #T-1 do ans=append(ans,T#i#0);
return(ans)
);


-- *************************
-- *** The function EonE ***
-- *************************

-- The function EonE computes the action of 
-- \bigwedge F on \bigwedge F^* for any free module F.
-- Fix a pair of dual bases x_1,...,x_N and x_1^*,...,x_N^*
-- for F and F^*, respectively. 
-- The input is two lists L1 and L2. If L1 is {a_1,...,a_s}
-- and L2 is {b_1,...,b_t}, then the output represents
--
-- (x_{a_1}^*\wedge ... \wedge x_{a_s}^*)(x_{b_1}\wedge 
-- ... \wedge x_{b_t}). (*)
--
-- If (*) is zero, then the output is ({},0).
-- If (*) is not zero, then the output is (L,\epsilon) where 
-- L ={c_1,...,c_u} is a list and \epsilon is \pm 1 and
-- (*) equals \epsilon \cdot x_{c_1}\wedge...\wedge x_{c_u}.

EonE=(L1,L2)->(
if not isSubset(L1,L2) then return({},0);
if not #L1 == #(unique L1) then return({},0);
if not #L2 == #(unique L2) then return({},0);
c:=1;
A:=L1;
B:=L2;
for i from 0 to (#L1-1) do(
pos:=position(B,i-> i==A_(#A-1));
A=delete(A_(#A-1),A);
B=delete(B_pos,B);
c=(-1)^pos*c);
return(B,c)
);
 
-- **************************
-- *** The function SonD  ***
-- **************************

-- The function SonD mimics the action of 
-- Sym_1U on D_2U^*.
-- INPUT: (a_1, a_2, a_3). Three integers
-- OUTPUT: Either 
--      ( {},0) if a_1 is not equal to a_2 or a_3, or
--      (a_3,1) if a_1=a_2, or
--      (a_2,1) if a_1=a_3.
-- In other words, SonD(a_1,a_2,a_3) gives information about
--      X_{a_1}(X_{a_2}X_{a_3})^*.
-- The Macaulay2 function contract(X_(a_1), X_(a_2)*X_(a_3))
-- gives the same information in a different form   

SonD=(ell1,ell2,ell3)->(
if ell1==ell2 then return ({ell3},1);
if ell1==ell3 then return ({ell2},1);
return({},0)
);

-- *************************************
-- *************************************
-- ** THE MAIN SCRIPT: SocDeg3(Phi,r) **
-- *************************************
-- *************************************

SocDeg3=(Phi,r)->(

-- ********************
-- ** The Basic Data **
-- ********************
  
PP:=ring(Phi);
d:=dim PP;   
PE:=coefficientRing(PP)[X_1..X_d];
map2PE:=map(PE,PP,{X_1..X_d});
varPP:={};
for i from 0 to d-1 do varPP=append(varPP,PP_i);
map2PP:=map(PP,PE, varPP); 
PHI:=map2PE(Phi);

if not isHomogeneous(PHI) then
error "Phi must be homogeneous.";

if not (degree(PHI))_0==3 then
error "Phi must have degree 3.";

v:=matrix{{X_1..X_d}};
v':=matrix{{0,X_2..X_d}};
p:=contract (X_1*(transpose v)*v,PHI);


if det p==0 then
error "One must choose Phi so that the only linear
       element in ann(Phi) is zero and the variable 
       P_0 is a weak Lefschetz element on P/ann Phi.";

pinverse:=inverse (p);

-- ********************
-- *****  If r=1 ******
-- ********************

 if (r==1 or r==d) then (
 sourcef11=Basisf11source(d);
 sourcef12=Basisf12source(d);

 f11:=(i)-> (
 Ts:=sourcef11#i;
 lTs1:=Ts#0#1;
 return(v*submatrix(pinverse,,{lTs1-1}))
 );

 f11mat:=f11(0);
 for i from 1 to d-2 do f11mat=f11mat|f11(i);
 f11mat=(-1)^d*f11mat;

 f12:=(i)-> (
 Ts:=sourcef12#i;
 lTs0:=Ts#0#0;
 lTs1:=Ts#0#1;
 muTs:=X_lTs0*X_lTs1;
 re:=muTs-X_1*v*pinverse*( matrix{{contract(muTs,PHI)}}//v);
 return(re)
 );
 f12mat:=f12(0);
 for i from 1 to (#sourcef12-1) do f12mat=f12mat|f12(i);
 if r==1 then return(map2PP(f11mat),map2PP(f12mat));
 fd1mat:=transpose f12mat;
 fd2mat:=(-1)^(d+1)* transpose f11mat;   
 if r==d then return(map2PP(fd1mat),map2PP(fd2mat))
 );

-- ***********************************
-- ** The functions blop, ip, and G **
-- *********************************** 

-- If u is an element of Sym_2U, (written as a degree two 
-- homogeneous polynomial in the variables X_1,...,X_d),  
-- then blop(u) is  u(Phi_3), written as a column vector. 
-- The basis for this vector is X_1,X_2, ... ,X_d.

-- If u and u' are homogeneous forms of degree two in 
-- X_1,..., X_d, then ip(u,u') 
-- is u(Phi_3)[p^{-1}(u'(Phi_3))]. The answer ip(u,u') 
-- is in the field k. (Of course, ip stands for inner 
-- product.)

-- In the paper G:\Sym(\bw^2 U)--> k is the homomorphism 
-- which sends (l1\wedge l2)(l3\w l4) to 
-- \langle l1*l3,l2*l4\rangle - \langle l1*l4,l2*l3\rangle.  
-- This function is accomplished in the program by 
-- G(l1,l2,l3,l4). The l's are linear elements in U.  

blop:=(u)-> (contract(matrix{{u}},matrix{{PHI}})//v);
ip:=(u,u')-> (
(transpose( blop(u))*pinverse*blop(u'))_(0,0));
G:=(l1,l2,l3,l4)->(ip(l1*l3,l2*l4)-ip(l1*l4,l2*l3));
 
zmat:=(r,c)->(map(PE^r,PE^c,(i,j)->0));

-- ****************************************
-- ** The functions coefD1, lala, and GG **
-- ****************************************


-- If nu1 is an element of U^* (written as a linear form 
-- in X_1,...,X_d), then coefD1(nu1) is the row vector of 
-- coefficients of nu1.

-- If nu1 and nu1' are elements of U^* (written as linear 
-- forms in X_1,...,X_d), then 
-- lala(nu1,nu1')=nu1(p^{-1}(nu1')). In the paper, the 
-- inner product "ip" is more important than the inner 
-- product "lala". In the paper, "ip" is denoted 
-- \langle -,-\rangle; and "lala" is denoted
-- \langle\langle -,- \rangle\rangle. This explains the 
-- name "lala".

-- If u1 is a linear form in X_1,...,X_d and nu2 is an 
-- element of D_2U^*, (with nu2 written as a homogeneous 
-- form of degree 2 in X_1,...,X_d), then 
-- contract(u1,nu2)=u1(nu2). (Recall that DU^* is a 
-- module over Sym(U).) The answer is an element of 
-- D_1U^*. The answer given as a linear form in  
-- X_1,...,X_d.

-- If a1 and a2 are integers between 2 and d and nu2 and 
-- nu2' are elements of D_2U. (The nu's should be written 
-- as homogeneous forms of degree 2 in X_1,..., X_d), then
-- GG(a1,a2,nu2,nu2') is equal to 
-- lala( X_{a1}(nu2),X_{a2}(nu2')) 
-- - lala( X_{a1}(nu2'),X_{a2}(nu2)) 



coefD1=(nu)-> (matrix{{nu}}//v); 
lala=(nu1,nu1')-> (
(transpose(coefD1(nu1)))*pinverse*coefD1(nu1'));

GG=(a1,a2,nu2,nu2')->(
lala(contract(X_a1,nu2),contract(X_a2,nu2'))-
lala(contract(X_a1,nu2'),contract(X_a2,nu2)));

-- ***************************
-- *** Make the matrix A_r *** 
-- ***************************

-- The basis for the source of A_r consists of all k_Ts 
-- as Ts roams over the set of standard Hooks with d-r 
-- rows and 2 columns.

-- The basis for the target of A_r consists of all k_Tt 
-- as Tt roams over the set of standard Hooks with d-r+1 
-- rows and 2 columns.

sourceAr:=BasisASource(r,d);
targetAr:=BasisATarget(r,d);

-- The function funcA(i,j) is equal to (-1)^{d+r} times 
-- (  [uTt(nuTs)\w tauTs](thetaTt)+
-- X_1* GG(tauTs(thetaTt)\t (uTt Phi\w nuTs))  ),
-- where Ts is the j-th element of the basis of the 
-- source of Ar, Tt is the i-th element of the basis of 
-- target of Ar, k_Ts is 
-- tauTs(omega_{U_0})\otimes nuTs \in \bw^rU_0\t  D_2U_0,
-- ell_Tt is thetaTt\otimes uTs \in \bw^{d+1-r}\t Sym_1 U_0,
-- omega is the basis X_{2}\w ...\w X_{d} of \bw^{d-1}U_0, 
-- and GG : (\bw^2 U_0)\t (D_2 U^*) -> k, with 
-- GG((u1\w u1')\t (nu2\w nu2'))  equal to 
-- [u1(nu2)][p^{-1}(u2(nu2'))]
-- -[u1(nu2')][p^{-1}(u1'(nu2))]]. 

funcA=(i,j)->(
Ts:=(sourceAr)#j;
Tt:=(targetAr)#i; 
lTs0:=Ts#0#0;
lTs1:=Ts#0#1; 
lTt1:=Tt#0#1; 
uTt:=X_lTt1;
nuTs:=X_lTs0*X_lTs1;

fst:=SonD(lTt1,lTs0,lTs1);
snd:=EonE(fst#0|tau(Ts),theta(Tt));
func1=(A,B)->(
  if (A#1==0 or B#1==0) then return(0_PE);
  if (not A#1==0 and not B#1==0) then 
     return(B#1*X_(snd#0#0))
);

ANS:=func1(fst,snd);

 temp:=EonE(tau(Ts),theta(Tt));
 if temp#1==0 then  return(ANS);
 if not temp#1==0 then (
 guts:=GG(temp#0#0,temp#0#1,contract(uTt,PHI),nuTs);
 ans':=temp#1*guts;
 return(ANS+X_1*ans')) 
);

row= (a)-> (TEMPA:=funcA(a,0),
for i from 1 to (#sourceAr-1) do 
TEMPA=TEMPA|funcA(a,i); return TEMPA);

preAr:=row(0);
for i from 1 to (#targetAr-1) do preAr=preAr||row(i);

Ar:=(-1)^(d+r)*preAr;
Armat:=map(PE^#targetAr,PE^#sourceAr,Ar);

-- ******************
-- ** Make B_r     **
-- ******************

-- The basis for the source of B_r consists of all 
-- ell_Ts as Ts roams over the set of standard Hooks 
-- with r rows and 2 columns.

-- The basis for the target of B_r consists of all 
-- k_Tt as Tt roams over the set of standard Hooks 
-- with d+1-r rows and 2 columns.

sourceBr:=BasisBSource(r,d);
targetBr:=BasisBTarget(r,d);

-- The function funcB(i,j) is equal to 
-- G((uTt\w uTs) [thetaTt(omegastar)](thetaTs) ),
-- where Ts is the j-th element of the basis of 
-- the source of Br, Tt is the i-th element of 
-- the basis of target of Br,
-- ell_Ts is thetaTs\t uTs in \bw^r U_0\t U_0, 
-- ell_{Tt} is thetaTt\t uTt in 
-- \bw^{d-1-r} U_0 \t U_0, omegastar is the 
-- basis X_{d}^*\w ...\w X_2^* of \bw^{d-1}U_0^*,
-- and G : Sym_2 (\bw^2 U_0) -> k is given by
-- G( (u_1\w u_1')\cdot (u_1''\t u_1''') ) equals 
-- <u_1u_1'',u_1'u_1'''>-<u_1u_1''',u_1''u_1'> .

omega:=toList(2..d);
omegastar:=rsort toList(2..d);
funcB=(i,j)->(
Ts:=(sourceBr)#j;
Tt:=(targetBr)#i;
ls:=Ts#0#1;
lt:=Tt#0#1;
uTs:=X_ls;
uTt:=X_lt;
thetaTs:= toList (apply(0..(r-1), c-> Ts#c#0 ));
thetaTt:= toList (apply(0..(d-r), c-> Tt#c#0 ));
imp:=(theta,theta')-> (
fst:=EonE(theta,omegastar);
scd:=EonE(fst#0,theta');
return(scd#0,fst#1*scd#1)
);
timp:=imp(thetaTt,thetaTs);
if timp#1==0 then return (0);
ans:=timp#1*G(uTt,uTs,X_(timp#0#0),X_(timp#0#1) );
return(ans)
);


row= (a) -> toList (apply(0..(#sourceBr-1), 
     c-> funcB(a,c) ));
MB:=matrix (toList apply (0..(#targetBr-1), 
     a-> row(a)));
c:=(-1)^((r-1)*(d-1));
Brmat:=c*map(PE^#targetBr,PE^#sourceBr,MB);

--- ***************************
--- *** Make the matrix C_r ***
--- ***************************

-- The basis for the source of C_r consists of all 
-- k_Ts as Ts roams over the set of standard Hooks 
-- with d-r rows and 2 columns.

-- The basis for the target of C_r consists of all 
-- k_Tt as Tt roams over the set of standard Hooks 
-- with r-1 rows and 2 columns.

sourceCr:=BasisCSource(r,d);
targetCr:=BasisCTarget(r,d);

-- The function funcC(i,j) is equal to 
-- GG( ((tauTt\w tauTs)(omega_{U_0}) \t (nuTt\w nuTs)),               
-- where Ts is the j-th element of the basis of the 
-- source of Cr, Tt is the i-th element of the basis 
-- of target of Cr, k_Ts is 
-- tauTs(omega_{U_0})\otimes nuTs in 
-- \bw^rU_0\t D_2 U_0^*, k_Tt is 
-- tauTt(omega_{U_0})\otimes nuTs in 
-- \bw^{d+1-r}\t D_2 U_0^*, omega is the basis 
-- X_{2}\w ...\w X_{d} of \bw^{d-1}U_0, and
-- GG : (\bw^2 U_0)\t (D_2 U^*) -> k, with 
-- GG((u1\w u1')\t (nu2\w nu2'))  equal to 
-- [u1(nu2)][p^{-1}(u2(nu2'))]
-- -  [u1(nu2')][p^{-1}(u1'(nu2))]]. 

funcC=(i,j)->(
Ts:=(sourceCr)#j;
Tt:=(targetCr)#i;
lTs0:=Ts#0#0;
lTs1:=Ts#0#1; 
lTt0:=Tt#0#0;
lTt1:=Tt#0#1; 
nuTs:=X_lTs0*X_lTs1;
nuTt:=X_lTt0*X_lTt1;
impC:=EonE((tau(Tt))|(tau(Ts)),omega);
if impC#1==0 then return (zmat(1,1));
ans:=impC#1*GG(impC#0#0,impC#0#1,nuTt,nuTs);
return(ans)
);

row=(a)-> (TEMPC:=funcC(a,0),
for i from 1 to (#sourceCr-1) do TEMPC=TEMPC|funcC(a,i); 
     return TEMPC);

MC:=row(0);
for i from 1 to (#targetCr-1) do MC=MC||row(i);

Crmat:=map(PE^#targetCr,PE^#sourceCr,MC);

--- ***************************
--- *** Make the matrix D_r ***
--- ***************************

-- The basis for the source of D_r consists of all 
-- ell_Ts as Ts roams over the set of standard Hooks
-- with r rows and 2 columns.

-- The basis for the target of D_r consists of all
-- k_Tt as Tt roams over the set of standard Hooks 
-- with r-1 rows and 2 columns.

sourceDr:=BasisDSource(r,d);
targetDr:=BasisDTarget(r,d);

-- The function funcD(i,j) is equal to 
-- [uTs(nuTt)\w tauTt](thetaTs)
-- +X_1* GG(tauTt(thetaTs)\t (uTs Phi\w nuTt)),
              
-- where Ts is the j-th element of the basis of the
-- source of Dr, Tt is the i-th element of the basis 
-- of target of Dr, ell_Ts is thetaTs\otimes uTs 
-- \in \bw^rU_0\t  U_0, k_Tt is 
-- tauTt(omega_{U_0})\otimes nuTs 
-- \in \bw^{d+1-r}\t D_2 U_0^*,
-- omega is the basis X_{2}\w ...\w X_{d} of 
-- \bw^{d-1}U_0, and
-- GG : (\bw^2 U_0)\t (D_2 U^*) -> k, with 
-- GG((u1\w u1')\t (nu2\w nu2'))  equal to 
-- [u1(nu2)][p^{-1}(u2(nu2'))] 
-- -  [u1(nu2')][p^{-1}(u1'(nu2))]]. 

funcD=(i,j)->(
Ts:=(sourceDr)#j;
Tt:=(targetDr)#i;
lTs1:=Ts#0#1; 
lTt0:=Tt#0#0;
lTt1:=Tt#0#1; 
uTs:=X_lTs1;
nuTt:=X_lTt0*X_lTt1;

fst:=SonD(lTs1,lTt0,lTt1);
snd:=EonE(fst#0|tau(Tt),theta(Ts));
func1=(A,B)->(
  if (A#1==0 or B#1==0) then return(0_PE);
  if (not A#1==0 and not B#1==0) then 
     return(B#1*X_(snd#0#0))
);
ANS:=func1(fst,snd);

 temp:=EonE(tau(Tt),theta(Ts));
 if temp#1==0 then  return(ANS); 
 if not temp#1==0 then (
 ans':=temp#1*GG(temp#0#0,temp#0#1,contract(uTs,PHI),nuTt);
 return(ANS+X_1*ans'))
);


row= (a) -> (TEMPD:=funcD(a,0),
   for i from 1 to (#sourceDr-1) do 
                    TEMPD=TEMPD|funcD(a,i); return TEMPD);

   MD:=row(0);
   for i from 1 to (#targetDr-1) do MD=MD||row(i);


Drmat:=map(PE^#targetDr,PE^#sourceDr,MD);

TA:=map2PP(Armat);
TB:=map2PP(Brmat);
TC:=map2PP(Crmat);
TD:=map2PP(Drmat);

return(TA,TB,TC,TD)
);

-- ********************************
-- ** The script MakeResAnnXPhi  **
-- ********************************

-- INPUT: Phi

-- OUTPUT: The list  {f_1,f_2,...,f_d} of matrices in 
--         the resolution (F,f) of P/ann (xPhi). 

MakeResAnnXPhi=(Phi)->(
 d:=dim(ring(Phi));
 XX:=(ring(Phi))_0;
 (f11,f12):=SocDeg3(Phi,1);
Flast:=map((ring(Phi))^1,,(XX*f11)|f12); 
MANS:={Flast};
for r from 2 to d-1 do (
(A,B,C,D):=SocDeg3(Phi,r);
Flast=map(source Flast,,(A|(XX*B))||((XX*C)|D));
MANS=append(MANS,Flast);
                        ); 
Flast=map(source Flast,,
      (transpose f12)||((-1)^(d+1)*XX*transpose f11));
MANS=append(MANS,Flast);
 return (MANS)
);


-- *******************************
-- ** The script MakeResAnnPhi  **
-- *******************************
-- INPUT: Phi

-- OUTPUT: The list of {f_1,f_2,...,f_d} of matrices 
--         in the resolution (F,f) of P/ann (Phi). 

MakeResAnnPhi=(Phi)->(
 d:=dim(ring(Phi));
 XX:=(ring(Phi))_0;
 (f11,f12):=SocDeg3(Phi,1);
Flast:=map((ring(Phi))^1,,(XX^2*f11)|f12); 
MANS:={Flast};
for r from 2 to d-1 do (
(A,B,C,D):=SocDeg3(Phi,r);
Flast=map(source Flast,,(A|B)||((XX^2*C)|D));
MANS=append(MANS,Flast);
                        ); 
Flast=map(source Flast,,
(transpose f12)||((-1)^(d+1)*XX^2*transpose f11));
MANS=append(MANS,Flast);
 return (MANS)
);

-- **********************************
-- **  The script TestResAnnXPhi  *** 
-- **********************************

-- INPUT: (Phi, Res), where Res is the resolution 
-- of P/ann xPhi which has previously been constructed.

-- The script checks that the ideal generated by the 
-- first differential really is equal to the the 
-- annihilator of x Phi. The script also checks that 
-- the proposed resolution does indeed resolve P/ann xPhi. 

-- Output: "pass" if everything is in order, 
--         "fail" otherwise. 
 

TestResAnnXPhi=(Phi,ResAnnXPhi)->(
XX:=(ring(Phi))_0;
d:=#ResAnnXPhi;
RlAn:=ideal(fromDual(contract(XX,Phi)));
if not (ideal(ResAnnXPhi#0))==(RlAn) then return("Fail");
for i from 1 to d-1 do (
if not image(ResAnnXPhi#i)==kernel(ResAnnXPhi#(i-1)) 
then return("Fail"));
return ("Pass")
);

-- ******************************
-- ** The script TestResAnnPhi **
-- ******************************


-- INPUT: (Phi, Res), where Res is the resolution of 
-- P/ann Phi which has previously been constructed.

-- The script checks that the ideal generated by 
-- the first differential really is equal to the 
-- annihilator of  Phi. The script also checks 
-- that the proposed resolution does indeed resolve 
-- P/ann Phi. 

-- Output: "pass" if everything is in order, 
--         "fail" otherwise. 
 

TestResAnnPhi=(Phi,ResAnnPhi)->(
XX:=(ring(Phi))_0;
d:=#ResAnnPhi;
if not (ideal (ResAnnPhi#0))==(ideal(fromDual(Phi))) 
            then return("Fail");
for i from 1 to d-1 do (
if not image(ResAnnPhi#i)==kernel(ResAnnPhi#(i-1)) 
            then return("Fail"));
return ("Pass")
);

-- \end{verbatim}
--  }