(* ::Package:: *)

(* ::Text:: *)
(*Initial assumptions*)


Needs["QuaternionAnalysis`"]


ScalarQ[x_]:=And@@((Head[x]=!=#)&/@{Complex,Quaternion,Paravector,Polynomial});


Quaternion[a_?ScalarQ,0,0,0]=a;


Polynomial[a_?ScalarQ]:=a;
Polynomial[a_Quaternion]:=a;
Polynomial[0..]:=0;
Polynomial[Longest[PatternSequence[0]..],x__]:=Polynomial[x];


(* ::Text:: *)
(*Basic Functions*)


Polynomial/:Plus[sc_?ScalarQ,p_Polynomial]:=Polynomial@@Plus@@@Transpose@PadLeft[{{sc},List@@p}];
Polynomial/:Plus[sc_Quaternion,p_Polynomial]:=Polynomial@@Plus@@@Transpose@PadLeft[{{sc},List@@p}];
Polynomial/:Plus[p1_Polynomial,p2_Polynomial]:=Polynomial@@Plus@@@Transpose@PadLeft[{List@@p1,List@@p2}];


Polynomial/:Times[sc_?ScalarQ,p_Polynomial]:=Map[sc #&,p];


Polynomial/:NonCommutativeMultiply[sc_Quaternion,p_Polynomial]:=Map[sc **#&,p];
Polynomial/:NonCommutativeMultiply[p_Polynomial,sc_Quaternion]:=Map[ #**sc&,p];
Polynomial/:NonCommutativeMultiply[p1_Polynomial,p2_Polynomial]:=Module[{dim1=Length@p1,dim2=Length@p2},
Polynomial@@Plus@@@Transpose[PadRight[MapThread[PadLeft,{Transpose[Outer[NonCommutativeMultiply,List@@p1,List@@p2]],Range[dim1,dim1+dim2-1]}]]]];


Polynomial/:Power[p_Polynomial,n_Integer?Positive]:=Nest[NonCommutativeMultiply[p,#]&,p,n-1];


Polynomial/:Conjugate[p_Polynomial]:=Conjugate/@p;


Unprotect[CharacteristicPolynomial];
CharacteristicPolynomial[q_?ScalarQ]:=Polynomial[1,-q]^2;
Quaternion/:CharacteristicPolynomial[q_Quaternion]:=Polynomial[1,-q]**Polynomial[1,-Conjugate[q]]
SyntaxInformation[CharacteristicPolynomial]={"ArgumentsPattern"->{_,_.}};
Protect[CharacteristicPolynomial];


(* ::Input::GrayLevel[0]:: *)
CompanionPolynomial[p_Polynomial]:=p**Conjugate[p];


Eval[p_Polynomial]:=Plus@@MapThread[NonCommutativeMultiply,{List@@p,Reverse@Function[u,NestList[(#**u)&,1,Length@p-1]][#]}]&


Eval[p_Polynomial,x_]:=Eval[p][x];


HornerEval[p_Polynomial,x_]:= Fold[NonCommutativeMultiply[#1,x]+#2&,0,p];


(* ::Text:: *)
(*Euclidean Algorithm*)


PolynomialDivisionR[p1_Polynomial,sc_?ScalarQ]:={p1/sc,0};
PolynomialDivisionR[p1_Polynomial,sc_Quaternion]:={p1**(1/sc),0};
PolynomialDivisionR[p1_Polynomial,p2_Polynomial]:=Module[{t,qq=0,rr=p1,degree=Length@p1-Length@p2},While[Head@rr===Polynomial&&degree>=0,t=Polynomial@@(PadRight[{First@rr**(1/First@p2)},degree+1]);
qq=qq+t;
rr=rr-t**p2;
degree=Length[rr]-Length[p2];];
{qq,rr}];


PNormalizeR[0]=0;
PNormalizeR[sc_?ScalarQ]:=(1/sc)**sc;
PNormalizeR[sc_Quaternion]:=1;
PNormalizeR[p_]:=p**(1/First@p);
PNormalizeL[0]=0;
PNormalizeL[sc_?ScalarQ]:=(1/sc)**sc;
PNormalizeL[sc_Quaternion]:=1;
PNormalizeL[p_]:=(1/First@p)**p;


GCDR[a_,0]= a;
GCDR[a_,b_]:=GCDR[b, Last[PolynomialDivisionR[a,b]]];
GCDR[a_,b_,c__]:=GCDR[GCDR[a,b],c];
PolynomialGCDR[a_,b_,c___]:=PNormalizeL[GCDR[a,b,c]];


PolynomialDivisionL[p1_Polynomial,sc_?ScalarQ]:={p1/sc,0};
PolynomialDivisionL[p1_Polynomial,sc_Quaternion]:={(1/sc)**p1,0};
PolynomialDivisionL[p1_Polynomial,p2_Polynomial]:=Module[{t,qq=0,rr=p1,degree=Length@p1-Length@p2},
While[Head@rr===Polynomial&&degree>=0,
t=Polynomial@@(PadRight[{(1/First@p2)**First@rr},degree+1]);
qq=qq+t;
rr=rr-t**p2;
degree=Length[rr]- Length[p2];];
{qq,rr}]


GCDL[a_,0]= a;
GCDL[a_,b_]:=GCDL[b, Last[PolynomialDivisionL[a,b]]]
GCDL[a_,b_,c__]:=GCDL[GCDL[a,b],c]
PolynomialGCDL[a_,b_,c___]:=PNormalizeR[GCDL[a,b,c]]


(* ::Text:: *)
(*Zeros and factors*)


ZeroQ[a_]:=PossibleZeroQ[a,Method->"ExactAlgebraics"];
ZeroQ[q_Quaternion]:=And@@ZeroQ/@q;
ZeroQ[p_Polynomial]:=And@@(ZeroQ/@p);


SimilarQ[q1_Quaternion,q2_Quaternion]:=ZeroQ[Re@q1-Re@q2] &&ZeroQ[Norm@q1-Norm@q2];


SphericalQ[p_Polynomial,q_Quaternion]:=ZeroQ[Eval[p][q]]&&ZeroQ[Eval[p][Conjugate@q]];


PolynomialFromChain[list_List]:=NonCommutativeMultiply@@Map[Polynomial[1,-#]&,Reverse@list];


SimilarChainQ[c1_List,c2_List]:=ZeroQ[PolynomialFromChain[c1]-PolynomialFromChain[c2]];


FactorShift[{q1_,q2_}]:=Module[{h=(Conjugate@q2-q1)},
If[ZeroQ@h,{q2,q1},{(1/h)**q2**h,(1/h)**q1**h}]];
FactorShift[l_List,m_Integer?Positive,n_Integer?Positive]:=Module[{lista=l,dim=Length@l,seq},
If[m>dim||n>dim,Message[FactorShift::args],seq=Sort/@Partition[Range[m,n,Sign[n-m]],2,1];
Map[(Part[lista,#]=FactorShift[Part[lista,#]])&,seq]];
lista];
FactorShift::args="Chain to short.";


RPol[fact_,i_]:=Module[{n=Length@fact},
Which[
i==1,Polynomial[1],
i==2,Polynomial[1,-Conjugate@#]&@@(Drop[fact,i-n-1]),
i<=n,(NonCommutativeMultiply@@(Polynomial[1,-Conjugate@#]&/@(Drop[fact,i-n-1])))
]
]


ZerosFromChain[fact_List]:=Module[{n=Length@fact,factors=fact,roots={},RPoli,RPoliz},If[n>Length@Union@Map[{Re@#,Norm@#}&,fact],Message[ZerosFromChain::args],
For[i=1,i<=n,i++,
RPoli=RPol[factors,i];
RPoliz=Eval[RPoli][factors[[i]]];
AppendTo[roots,RPoliz**factors[[i]]**(1/RPoliz)];
];
roots
]];
ZerosFromChain::args="Arguments in the same similarity class.";


ChainFromZeros[root_List]:=Module[{n=Length@root,factor,factors,QPoli=1,QPoliz},
If[n>Length@Union@Map[{Re@#,Norm@#}&,root],Message[ChainFromZeros::args1],
factor=First@root;factors={factor};
For[i=2,i<=n,i++,
QPoli=Polynomial[1,-factor]**QPoli;
QPoliz=Eval[QPoli][root[[i]]];
factor=QPoliz**root[[i]]**(1/QPoliz);
AppendTo[factors,factor];
];
factors
]];
ChainFromZeros::args1="Arguments in the same similarity class. Use an alternative syntax.";


ChainFromZeros[{},{}]:={};
ChainFromZeros[isoroots:{{_,_}..},{}]:=Module[{nisoroots=Length@isoroots,factor,multiplicity,factors,QPoli=1,QPoliz},
If[nisoroots>Length@Union@Map[{Re@#,Norm@#}&,First@Transpose@isoroots],Message[ChainFromZeros::args2],
factor=isoroots[[1,1]];multiplicity=isoroots[[1,2]];factors=Table[factor,{multiplicity}];
For[i=2,i<=nisoroots,i++,
QPoli=Power[Polynomial[1,-factor],multiplicity]**QPoli;
QPoliz=Eval[QPoli][isoroots[[i,1]]];
factor=QPoliz**isoroots[[i,1]]**(1/QPoliz);multiplicity=isoroots[[i,2]];
factors=Join[factors,Table[factor,{multiplicity}]];
];
factors
]];
ChainFromZeros[{},sphroots:{{_,_}..}]:=Module[{nsphroots=Length@sphroots},
If[nsphroots>Count[sphroots,{_Quaternion,_}],Message[ChainFromZeros::args4],
If[nsphroots>Length@Union@Map[{Re@#,Norm@#}&,First@Transpose@sphroots],Message[ChainFromZeros::args3],Null];
Flatten@(Table[{#1,Conjugate@#1},{#2}]&@@@sphroots)
]];
ChainFromZeros[isoroots:{{_,_}...},sphroots:{{_,_}...}]:=Module[{ch1=ChainFromZeros[isoroots,{}],ch2=ChainFromZeros[{},sphroots]},
If[ch1=!=Null && ch2=!= Null,Join[ch1,ch2],Null]];
ChainFromZeros::args2="Two or more isolated zeros in the same similarity class.";
ChainFromZeros::args3="Two or more spherical zeros in the same similarity class.";
ChainFromZeros::args4="Spherical zeros must be non-reals.";
