(* ::Package:: *)

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


Needs["Coquaternions`"]

(*<<Coquaternions.m*)


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


(* Coquaternion[a_?ScalarQ,0,0,0]=a; *)



Polynomial[a___, Coquaternion[b_, 0, 0, 0], c___] := Polynomial[a, b, c];
Polynomial[a_?ScalarQ]:=a;
Polynomial[a_Coquaternion]:=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_Coquaternion,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_Coquaternion,p_Polynomial]:=Map[sc **#&,p];
Polynomial/:NonCommutativeMultiply[p_Polynomial,sc_Coquaternion]:=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;
Coquaternion/:CharacteristicPolynomial[q_Coquaternion]:=Polynomial[1,-q]**Polynomial[1,-Conjugate[q]]
SyntaxInformation[CharacteristicPolynomial]={"ArgumentsPattern"->{_,_.}};
Protect[CharacteristicPolynomial];


CompanionPolynomial[p_Polynomial]:=Simplify[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]//Simplify//QSimplify;


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

Polynomial /: PolynomialZeroQ[P_Polynomial, q_] := QSimplify[Eval[P, q]] === 0;


(*Euclidean Algorithm*)

PolynomialDivisionR[p1_Polynomial,sc_?ScalarQ]:={p1/sc,0};
PolynomialDivisionR[p1_Polynomial,sc_Coquaternion]:={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_Coquaternion]:=1;
PNormalizeR[p_]:=p**(1/First@p);
PNormalizeL[0]=0;
PNormalizeL[sc_?ScalarQ]:=(1/sc)**sc;
PNormalizeL[sc_Coquaternion]:=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_Coquaternion]:={(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]]


(*Zeros and factors*)


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


SimilarQ[q1_Coquaternion,q2_Coquaternion]:=ZeroQ[Re@q1-Re@q2] &&ZeroQ[Norm@q1-Norm@q2];


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


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])))
]
]


(* ::Input:: *)
(**)


\[Phi][q_Coquaternion]:=ArcSinh[Re@q/Norm@q];
\[Omega][q_Coquaternion]:=Vec@q/Norm[Vec@q];
\[Psi][q_Coquaternion]:=ArcSinh[Norm[Vec@q]/Norm@q];
\[Theta][q_Coquaternion]:=ArcCos[Re@q/Norm@q];


Clear@PolarForm


(* ::Input:: *)
(**)


PolarForm[q_Coquaternion?LightLikeQ]:=Message[PolarForm::LightLikeCoq]
PolarForm[q_Coquaternion?SpaceLikeQ,{nat_,r_,\[Phi]_,\[Omega]_}]:=Module[{normq=Norm[q],q0=Re[q],vecq=Vec[q]},If[LightLikeQ[vecq],Message[PolarForm::LightLikeVecCoq],{nat->Nature[q],r-> normq,\[Phi]->ArcSinh[q0/normq],\[Omega]->vecq/Norm[vecq]}]]
PolarForm[q_Coquaternion?TimeLikeQ,{nat_,r_,\[Psi]_,\[Omega]_,sgn_}]:=Module[{normq=Norm[q],q0=Re[q],vecq=Vec[q]},Switch[Nature[vecq],"Spacelike",{nat->"Timelike-Spacelike",r-> normq,\[Psi]->ArcSinh[Norm[vecq]/normq],\[Omega]->vecq/Norm[vecq],sgn->Sign[Re[q]]},"Timelike",{nat->"Timelike-Timelike",r-> normq,\[Psi]->ArcCos[Re[q]/normq],\[Omega]->vecq/Norm[vecq]},_,Message[PolarForm::LightLikeVecCoq]]]



(* ::Input:: *)
(**)


PolarForm[q_Coquaternion?LightLikeQ]:=Message[PolarForm::LightLikeCoq]
PolarForm[q_Coquaternion?SpaceLikeQ]:=Module[{normq=Norm[q],q0=Re[q],vecq=Vec[q]},If[LightLikeQ[vecq],Message[PolarForm::LightLikeVecCoq],{nat->Nature[q],r-> normq,\[Phi]->ArcSinh[q0/normq],\[Omega]->vecq/Norm[vecq]}]]
PolarForm[q_Coquaternion?TimeLikeQ]:=Module[{normq=Norm[q],q0=Re[q],vecq=Vec[q]},Switch[Nature[vecq],"Spacelike",{nat->"Timelike-Spacelike",r-> normq,\[Psi]->ArcSinh[Norm[vecq]/normq],\[Omega]->vecq/Norm[vecq],sgn->Sign[Re[q]]},"Timelike",{nat->"Timelike-Timelike",r-> normq,\[Psi]->ArcCos[Re[q]/normq],\[Omega]->vecq/Norm[vecq]},_,Message[PolarForm::LightLikeVecCoq]]]


(* ::Input:: *)
(**)


FromPolarForm[{"Spacelike",r_,\[Phi]_,\[Omega]_,sgn_:1}]:=r(Sinh[\[Phi]]+Cosh[\[Phi]]\[Omega])


(* ::Input:: *)
(**)


PolyForm[p_, x_] := Plus @@ MapThread[Times, {List @@ p, Append[Table[Power[x, n], {n, Length@p - 1, 1, -1}], 1]}] 


ToCoquaternion[a_?(Element[#, Complexes] &)] := Coquaternion[Re[a], Im[a], 0, 0] 
ToCoquaternion[a_?ScalarQ]:=Coquaternion[a,0,0,0]
ToCoquaternion[a_Coquaternion]:=a



(* ::Funcoes definidas pela Joana:: *)
(**)


Clear[ClassRealRoots]
ClassRealRoots[a_,b_]:=Coquaternion[(a+b)/2,0,(b-a)/2,0];

Clear[AdmissibleClasses];
AdmissibleClasses[pol_]:=Module[{CP,rootsCP,realRoots,imagRoots,imagClasses,realClasses,x},
CP=CompanionPolynomial[pol];
(*rootsCP=Simplify[ComplexExpand[x/.Solve[PolyForm[CP]\[Equal]0]]];*)
rootsCP=x/.Solve[PolyForm[CP,x]==0];
realRoots=Select[rootsCP,(Im[#]==0)&];
imagRoots=Select[rootsCP,(Im[#]>0)&];
If [imagRoots!={},
imagClasses=Union[Map[Coquaternion[Re[#],Im[#],0,0]&,imagRoots]],imagClasses={}];
If[realRoots!={},
realClasses=Union[Flatten[Table[ClassRealRoots[Part[realRoots,i],Part[realRoots,j]],{i,1,Length[realRoots]-1},{j,i+1,Length[realRoots]}]]],
    realClasses={}];
classes=Join[imagClasses,realClasses];
If[classes=={},
    Print["Unable to compute classes"];,{}];
 Return[classes]]


Clear[Zeros];

Zeros[pol_]:=Module[{coefs,degree,isolatedZeros,linearZeros,hyperboloidalZeros,classes,nClasses,qC,qC0,qCV,dqC,tqC,alfas,betas,ax,bx,a,b,detb,z,z0,z1,z2,z3,matb,realSol,complexSol,sol,condition,gamma0,gamma1,x0,x1,bv,b0,b1,b2,b3,k1,k2,denominator,beta},
coefs=Reverse[Map[ToCoquaternion,List@@pol]];
degree=Length[coefs]-1;
isolatedZeros={};
linearZeros={};
hyperboloidalZeros={};
classes=AdmissibleClasses[pol];
nClasses=Length[classes];
Do[{
qC=classes[[j]];
dqC=Det[qC];
qC0=Re[qC];
qCV=Vec[qC];
tqC=2*qC0;
alfas[0]=1;
betas[0]=0;
alfas[k_]:=-dqC*betas[k-1];
betas[k_]:=alfas[k-1]+tqC*betas[k-1];
ax=Table[alfas[i],{i,0,degree}];
bx=Table[betas[i],{i,0,degree}];
b=Simplify[bx . coefs];
a=Simplify[ax . coefs];
detb=Simplify[Det[b]];
If[QSimplify[b] === 0,
     If[QSimplify[a] === 0,
         hyperboloidalZeros=Append[hyperboloidalZeros,qC], {}], 
     If[detb!=0,
             {z=Simplify[-(Conjugate[b]**a)/detb];
            isolatedZeros=Append[isolatedZeros,z];    
            },
            {matb=CoquaternionTo4DMatrixL[b];
           realSol=Solve[matb . {x0,0,0,0}==-List@@a];
         If[realSol!={},
           {gamma0=x0/.realSol[[1]];
          condition=(gamma0-qC0)^2==-Det[qCV];
          If[condition,
                    {{b0,b1,b2,b3}=List@@b;
		  denominator=b0^2+b1^2;
		  k1=-(b0*b2+b1*b3)/denominator;
		  k2=(b1*b2-b0*b3)/denominator;
		  z=Simplify[Coquaternion[qC0,C[1],k2*C[1]+k1(qC0-gamma0),-k1*C[1]+k2(qC0-gamma0)]];
		 linearZeros=Append[linearZeros,z];               
     },{}]},
          {
           matb=CoquaternionTo4DMatrixL[b];
           complexSol=Solve[matb . {x0,x1,0,0}==-List@@a];
           If[complexSol!={},
                            {
                            gamma0=x0/.complexSol[[1]];
                           gamma1=x1/.complexSol[[1]];
                            {b0,b1,b2,b3}=List@@b;
		         denominator=b0^2+b1^2;
		         k1=-(b0*b2+b1*b3)/denominator;
		         k2=(b1*b2-b0*b3)/denominator;
		         beta=(Det[qCV]+(qC0-gamma0)^2-gamma1^2)/(2*gamma1);
                           z1=beta+gamma1;
                            z2=k2*beta+k1*(qC0-gamma0);
                            z3=  -k1*beta+k2*(qC0-gamma0);
                            z=Simplify[Coquaternion[qC0,z1,z2,z3]];
                            isolatedZeros=Append[isolatedZeros,z];
                            },{}]}]}]]},{j,nClasses}];
zeros={isolatedZeros,linearZeros,hyperboloidalZeros};
Return[zeros]]


Clear[NumZeros];

NumZeros[pol_]:=Module[{coefs,degree,nisolatedZeros,nlinearZeros,nhyperboloidalZeros,classes,nClasses,qC,qC0,qCV,dqC,tqC,alfas,betas,ax,bx,a,b,detb,z,z0,z1,z2,z3,matb,realSol,complexSol,sol,condition,gamma0,gamma1,x0,x1,bv,b0,b1,b2,b3,k1,k2,denominator,beta,nzeros},
coefs=Reverse[Map[ToCoquaternion,List@@pol]];
degree=Length[coefs]-1;
nisolatedZeros=0;
nlinearZeros=0;
nhyperboloidalZeros=0;
classes=AdmissibleClasses[pol];
nClasses=Length[classes];
Do[{
qC=classes[[j]];
dqC=Det[qC];
qC0=Re[qC];
qCV=Vec[qC];
tqC=2*qC0; 
alfas[0]=1;
betas[0]=0;
alfas[k_]:=-dqC*betas[k-1];
betas[k_]:=alfas[k-1]+tqC*betas[k-1];
ax=Table[alfas[i],{i,0,degree}];
bx=Table[betas[i],{i,0,degree}];
b=Simplify[bx . coefs];
a=Simplify[ax . coefs];
detb=Simplify[Det[b]];
If[QSimplify[b] === 0,
     If[QSimplify[a] === 0,
         nhyperboloidalZeros=nhyperboloidalZeros+1, {}], 
     If[detb!=0,
             {
            nisolatedZeros=nisolatedZeros+1;  
            (*Print["Isolated-type1"];  *)
            },
            {matb=CoquaternionTo4DMatrixL[b];
           realSol=Solve[matb . {x0,0,0,0}==-List@@a];
         If[realSol!={},
           {gamma0=x0/.realSol[[1]];
          condition=(gamma0-qC0)^2==-Det[qCV];
          If[condition,
                 nlinearZeros=nlinearZeros+1;               
     ,{}]},
          {
           matb=CoquaternionTo4DMatrixL[b];
           complexSol=Solve[matb . {x0,x1,0,0}==-List@@a];
           If[complexSol!={},
                            {
                           nisolatedZeros=nisolatedZeros+1;
(*Print["Isolated-type2"];*)
                            },{}]}]}]]},{j,nClasses}];
nzeros={nisolatedZeros,nlinearZeros,nhyperboloidalZeros};
Return[nzeros]]



(* ::Funcoes definidas pela Irene:: *)
(**)


RootsN[Coquaternion[x0_?Positive, 0, 0, 0], 
  n_Integer?(OddQ[#] && # > 1 &)] := {{Surd[x0, n] Coquaternion[1, 0, 
     0, 0]}, Map[ToCoquaternion, 
   Table[Surd[x0, n] ComplexExpand[
      Cos[(2 k \[Pi])/n] + I Sin[(2 k \[Pi])/n]] , {k, 1, (n - 1)/
     2} ]]}
RootsN[Coquaternion[x0_?Positive, 0, 0, 0], 
   n_Integer?(EvenQ[#] && # > 2 &)] := {{-Surd[x0, n], 
     Surd[x0, n]} Coquaternion[1, 0, 0, 0], 
   Union[{Surd[x0, n] Coquaternion[0, 0, 1, 0]}, 
    Map[ToCoquaternion, 
     Table[Surd[x0, n] ComplexExpand[
        Cos[(2 k \[Pi])/n] + I Sin[(2 k \[Pi])/n]] , {k, 1, 
       n/2 - 1} ]]]};
RootsN[Coquaternion[x0_?Positive, 0, 0, 0], 
   2] := {{-Surd[x0, 2], Surd[x0, 2]} Coquaternion[1, 0, 0, 
     0], {Surd[x0, 2] Coquaternion[0, 0, 1, 0]}};
RootsN[Coquaternion[x0_?Negative, 0, 0, 0], 
  n_Integer?(OddQ[#] && # > 1 &)] := {{Surd[x0, n]} Coquaternion[1, 0,
     0, 0], Map[ToCoquaternion, 
   Table[Surd[-x0, n] ComplexExpand[
      Cos[((2 k + 1) \[Pi])/n] + I Sin[((2 k + 1) \[Pi])/n]] , {k, 
     0, (n - 3)/2} ]]}
RootsN[Coquaternion[x0_?Negative, 0, 0, 0], 
   n_Integer?(EvenQ[#] && # > 1 &)] := {{}, 
   Map[ToCoquaternion, 
    Table[Surd[-x0, n] ComplexExpand[
       Cos[((2 k + 1) \[Pi])/n] + I Sin[((2 k + 1) \[Pi])/n]] , {k, 0,
       n/2 - 1} ]]};
RootsN[Coquaternion[0, 0, 0, 0], n_Integer?(# > 1 &)] := {{}, {0}};
RootsN[q_Coquaternion, 
   n_Integer?(# > 1 &)] := {RootsN[q, n, Nature[q], 
    Nature[Vec@q]], {}};
RootsN[q_Coquaternion, n_?OddQ, "Spacelike", 
  "Spacelike"] := {Norm[q]^(1/
   n) (Sinh[\[Phi][q]/n] + \[Omega][q] Cosh[\[Phi][q]/n])}
RootsN[q_Coquaternion, n_?EvenQ, "Spacelike", "Spacelike"] := {}
RootsN[q_Coquaternion, n_?OddQ, "Timelike", 
   "Spacelike"] := {Norm[q]^(1/
    n) (Sign[Re@q] Cosh[\[Psi][q]/n] + \[Omega][q] Sinh[\[Psi][q]/n])};
RootsN[q_Coquaternion, n_?EvenQ, "Timelike", "Spacelike"] := 
  If[Re[q] < 
    0, {}, {Norm[q]^(1/
      n) (Cosh[\[Psi][q]/n] + \[Omega][q] Sinh[\[Psi][q]/n]) // 
     FullSimplify, 
    Norm[q]^(1/
      n) (-Cosh[\[Psi][q]/n] - \[Omega][q] Sinh[\[Psi][q]/n]) // 
     FullSimplify, 
    Norm[q]^(1/
      n) (Sinh[\[Psi][q]/n] + \[Omega][q] Cosh[\[Psi][q]/n]) // 
     FullSimplify, 
    Norm[q]^(1/
      n) (-Sinh[\[Psi][q]/n] - \[Omega][q] Cosh[\[Psi][q]/n]) // 
     FullSimplify}];
RootsN[q_Coquaternion, n_, "Timelike", "Timelike"] := 
  Table[Norm[q]^(1/
    n) (Cos[(\[Theta][q] + 2 k \[Pi])/
       n] + \[Omega][q] Sin[(\[Theta][q] + 2 k \[Pi])/n]), {k, 0, 
    n - 1}];
RootsN[q_Coquaternion, n_, "Lightlike", "Lightlike"] := {}
RootsN[q_Coquaternion, n_?OddQ, "Timelike", 
   "Lightlike"] := {Surd[Re@q, n] + Vec[q]/(
    n (Surd[Re@q, n])^(n - 1))};
RootsN[q_Coquaternion, n_?EvenQ, "Timelike", "Lightlike"] := 
  If[Re[q] < 0, 
   Null, {-Surd[Re@q, n] - Vec[q]/(n (Surd[Re@q, n])^(n - 1)) , 
    Surd[Re@q, n] + Vec[q]/(n (Surd[Re@q, n])^(n - 1)) }];
RootsN[q_Coquaternion, n_?OddQ, "Lightlike", 
   "Spacelike"] := {Surd[Re@q/2^(n - 1), n] + 
    Vec[q]/(Surd[2 Re@q, n])^(n - 1)};
RootsN[q_Coquaternion, n_?EvenQ, "Lightlike", "Spacelike"] := 
  If[Re[q] < 
    0, {}, {-Surd[Re@q/2^(n - 1), n] - Vec[q]/(Surd[2 Re@q, n])^(
     n - 1), Surd[Re@q/2^(n - 1), n] + Vec[q]/(Surd[2 Re@q, n])^(
     n - 1)}];
