QuatDotMP[x_,y_,p_]:= (*  This function returns the quaternion product of  x  followed by  y  (mod p).  *)(If[Length[x]≠4||Length[y]≠4,Return["Undefined"],Return[Mod[{x[[1]]y[[1]]-x[[2]]y[[2]]-x[[3]]y[[3]]-x[[4]]y[[4]],x[[1]]y[[2]]+x[[2]]y[[1]]+x[[3]]y[[4]]-x[[4]]y[[3]],x[[1]]y[[3]]-x[[2]]y[[4]]+x[[3]]y[[1]]+x[[4]]y[[2]],x[[1]]y[[4]]+x[[2]]y[[3]]-x[[3]]y[[2]]+x[[4]]y[[1]]},p]]]);

(Subscript[TransMatL, 4])[q_]:=  (*  This is the Transformation Matrix corresponding to quaternion multiplication from the Left.   *)({{q[[1]],-q[[2]],-q[[3]],-q[[4]]},{q[[2]],q[[1]],-q[[4]],q[[3]]},{q[[3]],q[[4]],q[[1]],-q[[2]]},{q[[4]],-q[[3]],q[[2]],q[[1]]}});

(Subscript[TransMatR, 4])[q_]:= (*  The Transformation Matrix below corresponds to quaternion multiplication from the Right.   *)({{q[[1]],-q[[2]],-q[[3]],-q[[4]]},{q[[2]],q[[1]],q[[4]],-q[[3]]},{q[[3]],-q[[4]],q[[1]],q[[2]]},{q[[4]],q[[3]],-q[[2]],q[[1]]}});

BuildUp[n_,(Subscript[Q, 1])_,(Subscript[Q, 2])_,(Subscript[Q, 3])_,(Subscript[Q, 4])_]:=(*  This patches together 4 (n x n) quarters to return a (2n x 2n) matrix.  *)(Return[Join[Table[Join[Table[(Subscript[Q, 1])[[i,j]],{j,n}],Table[(Subscript[Q, 2])[[i,j]],{j,n}]],{i,n}],Table[Join[Table[(Subscript[Q, 3])[[i,j]],{j,n}],Table[(Subscript[Q, 4])[[i,j]],{j,n}]],{i,n}]]]);

PowerExtendedPartition[Mtx_]:=(Module[{d,n,tmp,NumRows,NumCols,RowsQ,ColsQ},d=Dimensions[Mtx];NumRows=d[[1]];NumCols=d[[2]];n=1;Label[GetPower];n=n+1;If[NumRows>(2^n),Goto[GetPower]];tmp=Mtx;If[NumRows<(2^n),tmp=IdentityMatrix[(2^n)];Do[tmp[[i,j]]=Mtx[[i,j]],{i,NumRows},{j,NumCols}]];RowsQ=(2^(n - 1));ColsQ=(2^(n - 1));(Subscript[Q, 1])=Table[tmp[[i,j]],{i,1,RowsQ},{j,1,ColsQ}];(Subscript[Q, 2])=Table[tmp[[i,j]],{i,1,RowsQ},{j,ColsQ+1,(2^n)}];(Subscript[Q, 3])=Table[tmp[[i,j]],{i,RowsQ+1,(2^n)},{j,1,ColsQ}];(Subscript[Q, 4])=Table[tmp[[i,j]],{i,RowsQ+1,(2^n)},{j,ColsQ+1,(2^n)}]]);

SwitNeg[M_]:=(Return[{{M[[2,2]],-M[[1,2]]},{-M[[2,1]],M[[1,1]]}}]);
Inv[x_,p_]:=(Module[{n},If[Floor[x]≠x,Return["Undefined"]];n=Mod[x,p];If[n==0,Return["Undefined"]];If[n==1||n==p-1,Return[n]];st=1;Label[continue];st=st+1;If[Mod[n*st,p]==1,Return[st],Goto[continue]]]);

(*  The Eight-Step Method with 2x2 matrices, each step being reduced modulo a prime, p.  *)
ModPEP2x2[A_,B_,C_,D_,p_]:=    (*  The first appropriate PEP Method is chosen.  *)(Module[{W,X,Y,Z},If[Mod[Det[A],p]≠0,Goto[(Subscript[PEP, 1])]];If[Mod[Det[B],p]≠0,Goto[(Subscript[PEP, 2])]];If[Mod[Det[C],p]≠0,Goto[(Subscript[PEP, 3])]];If[Mod[Det[D],p]≠0,Goto[(Subscript[PEP, 4])],Print["All 2x2 determinants = 0."];Return["Use IP."]];

Label[(Subscript[PEP, 1])];ivA=Inv[Det[A],p];R1=Mod[-ivA*SwitNeg[A],p];S1=Mod[R1.B,p];T1=Mod[C.S1+D,p];iv=Inv[Det[T1],p];If[iv=="Undefined",Return["No Inverse exists"]];Z=Mod[iv*SwitNeg[T1],p];X=Mod[S1.Z,p];V1=Mod[C.R1,p];Y=Mod[Z.V1,p];W=Mod[S1.Y-R1,p];Return[BuildUp[2,W,X,Y,Z]];

Label[(Subscript[PEP, 2])];ivB=Inv[Det[B],p];R2=Mod[-ivB*SwitNeg[B],p];S2=Mod[R2.A,p];T2=Mod[D.S2+C,p];iv=Inv[Det[T2],p];If[iv=="Undefined",Return["No Inverse exists"]];X=Mod[iv*SwitNeg[T2],p];Z=Mod[S2.X,p];V2=Mod[D.R2,p];W=Mod[X.V2,p];Y=Mod[S2.W-R2,p];Return[BuildUp[2,W,X,Y,Z]];

Label[(Subscript[PEP, 3])];ivC=Inv[Det[C],p];R3=Mod[-ivC*SwitNeg[C],p];S3=Mod[R3.D,p];T3=Mod[A.S3+B,p];iv=Inv[Det[T3],p];If[iv=="Undefined",Return["No Inverse exists"]];Y=Mod[iv*SwitNeg[T3],p];W=Mod[S3.Y,p];V3=Mod[A.R3,p];Z=Mod[Y.V3,p];X=Mod[S3.Z-R3,p];Return[BuildUp[2,W,X,Y,Z]];

Label[(Subscript[PEP, 4])];ivD=Inv[Det[D],p];R4=Mod[-ivD*SwitNeg[D],p];S4=Mod[R4.C,p];T4=Mod[B.S4+A,p];iv=Inv[Det[T4],p];If[iv=="Undefined",Return["No Inverse exists"]];W=Mod[iv*SwitNeg[T4],p];Y=Mod[S4.W,p];(Subscript[V, 4])=Mod[B.R4,p];X=Mod[W.V4,p];Z=Mod[S4.X-R4,p];Return[BuildUp[2,W,X,Y,Z]]]);

(*  The Eight-Step Method with 4x4 matrices, each step being reduced modulo a prime, p.  *)
ModPEP4x4[A_,B_,C_,D_,p_]:=(    (*  The first appropriate PEP Method is chosen.  *)If[Mod[Det[A],p]≠0,Goto[(Subscript[PEP, 1])]];If[Mod[Det[B],p]≠0,Goto[(Subscript[PEP, 2])]];If[Mod[Det[C],p]≠0,Goto[(Subscript[PEP, 3])]];If[Mod[Det[D],p]≠0,Goto[(Subscript[PEP, 4])]];Return["All 4x4 determinants = 0."];
 
Label[(Subscript[PEP, 1])];PowerExtendedPartition[A];(Subscript[R, 1])=Mod[-ModPEP2x2[(Subscript[Q, 1]),(Subscript[Q, 2]),(Subscript[Q, 3]),(Subscript[Q, 4]),p],p];(Subscript[S, 1])=Mod[(Subscript[R, 1]).B,p];(Subscript[T, 1])=Mod[C.(Subscript[S, 1])+D,p];PowerExtendedPartition[(Subscript[T, 1])];Z=ModPEP2x2[(Subscript[Q, 1]),(Subscript[Q, 2]),(Subscript[Q, 3]),(Subscript[Q, 4]),p];If[Z=="No Inverse exists",Print["There is no solution to the given system."];Return[]];X=Mod[(Subscript[S, 1]).Z,p];(Subscript[V, 1])=Mod[C.(Subscript[R, 1]),p];Y=Mod[Z.(Subscript[V, 1]),p];W=Mod[(Subscript[S, 1]).Y-(Subscript[R, 1]),p];Return[];

Label[(Subscript[PEP, 2])];PowerExtendedPartition[B];(Subscript[R, 2])=Mod[-ModPEP2x2[(Subscript[Q, 1]),(Subscript[Q, 2]),(Subscript[Q, 3]),(Subscript[Q, 4]),p],p];(Subscript[S, 2])=Mod[(Subscript[R, 2]).A,p];(Subscript[T, 2])=Mod[D.(Subscript[S, 2])+C,p];PowerExtendedPartition[(Subscript[T, 2])];X=ModPEP2x2[(Subscript[Q, 1]),(Subscript[Q, 2]),(Subscript[Q, 3]),(Subscript[Q, 4]),p];If[X=="No Inverse exists",Print["There is no solution to the given system."];Return[]];Z=Mod[(Subscript[S, 2]).X,p];(Subscript[V, 2])=Mod[D.(Subscript[R, 2]),p];W=Mod[X.(Subscript[V, 2]),p];Y=Mod[(Subscript[S, 2]).W-(Subscript[R, 2]),p];Return[];

Label[(Subscript[PEP, 3])];PowerExtendedPartition[C];(Subscript[R, 3])=Mod[-ModPEP2x2[(Subscript[Q, 1]),(Subscript[Q, 2]),(Subscript[Q, 3]),(Subscript[Q, 4]),p],p];(Subscript[S, 3])=Mod[(Subscript[R, 3]).D,p];(Subscript[T, 3])=Mod[A.(Subscript[S, 3])+B,p];PowerExtendedPartition[(Subscript[T, 3])];Y=ModPEP2x2[(Subscript[Q, 1]),(Subscript[Q, 2]),(Subscript[Q, 3]),(Subscript[Q, 4]),p];If[Y=="No Inverse exists",Print["There is no solution to the given system."];eturn[]];W=Mod[(Subscript[S, 3]).Y,p];(Subscript[V, 3])=Mod[A.(Subscript[R, 3]),p];Z=Mod[Y.(Subscript[V, 3]),p];X=Mod[(Subscript[S, 3]).Z-(Subscript[R, 3]),p];Return[];

Label[(Subscript[PEP, 4])];PowerExtendedPartition[D];(Subscript[R, 4])=Mod[-ModPEP2x2[(Subscript[Q, 1]),(Subscript[Q, 2]),(Subscript[Q, 3]),(Subscript[Q, 4]),p],p];(Subscript[S, 4])=Mod[(Subscript[R, 4]).D,p];(Subscript[T, 4])=Mod[B.(Subscript[S, 4])+A,p];PowerExtendedPartition[(Subscript[T, 4])];W=ModPEP2x2[(Subscript[Q, 1]),(Subscript[Q, 2]),(Subscript[Q, 3]),(Subscript[Q, 4]),p];If[W=="No Inverse exists",Print["There is no solution to the given system."];Return[]];Y=Mod[(Subscript[S, 4]).W,p];(Subscript[V, 4])=Mod[B.(Subscript[R, 4]),p];X=Mod[W.(Subscript[V, 4]),p];Z=Mod[(Subscript[S, 4]).X-(Subscript[R, 4]),p];Return[]);

RandQuatModP[p_]:=(If[PrimeQ[p],RQMP={};Do[RQMP=Append[RQMP,Floor[p*Random[]]],{i,4}];Return[RQMP],"p not prime, joker."])

(*   To quote from the wikipedia web site,"In mathematics, the general linear group of degree n is the set of n×n invertible matrices, together with the operation of ordinary matrix multiplication. This forms a group because the product of two invertible matrices is again invertible and the inverse of an invertible matrix is invertible. The name is because the columns of an invertible matrix are linearly independent, hence the vectors/points they define are in general linear position, and matrices in the general linear group take points in general linear position to points in general linear position."  The number of symmetries of a General Linear Group of degree n over any field, F -> GL(n,F)  is defined below.  *)
GL[n_,F_]:=(Module[{NumSyms,m,factor},If[F<2||n<1,Return["Undefined"]];NumSyms=1;m=0;Label[continue];factor=(F^n)-(F^m);NumSyms=NumSyms*factor;m=m+1;If[m==n,Return[NumSyms],Goto[continue]]])

(*  Choosing a small prime modulus makes the probability that no solution exists higher than if one chooses a 2 or 
3-dight prime. Let p = the prime of your choice. Since there are p numbers in the system mod p, it would
 seem that the probability that the determinant is zero of a randomly picked system is roughly (( 1 )/p). I believe the formula for the probability of a determinant being zero = 1 - GL[2,p] / (p^4).        *)

SeedRandom[9791];p=7;(Subscript[a, 1])=RandQuatModP[p];(Subscript[a, 2])=RandQuatModP[p];(Subscript[b, 1])=RandQuatModP[p];(Subscript[b, 2])=RandQuatModP[p];(Subscript[c, 1])=RandQuatModP[p];(Subscript[c, 2])=RandQuatModP[p];(Subscript[d, 1])=RandQuatModP[p];(Subscript[d, 2])=RandQuatModP[p];(Subscript[k, 1])=RandQuatModP[p];(Subscript[k, 2])=RandQuatModP[p]

(*    The system we solve for x and y is:

    ((Subscript[a, 1])x)(Subscript[b, 1]) + ((Subscript[c, 1])y)(Subscript[d, 1]) = (Subscript[k, 1])  &
   ((Subscript[a, 2])x)(Subscript[b, 2]) + ((Subscript[c, 2])y)(Subscript[d, 2]) = (Subscript[k, 2]). 
	
	Quaternions are noncommutative so that multiplication from the left is different than multiplication 
from the right. For the most general case we choose an example that illustrates multiplication of the 
unknkowns from both directions. The grouping given above is not important because quaternions are associative.   *)

(*   AA corresponds to ((Subscript[a, 1])x)(Subscript[b, 1])   *)

AA=Mod[TransMat4L[(Subscript[a, 1])].TransMat4R[(Subscript[b, 1])],p]

(*   BB corresponds to ((Subscript[c, 1])y)(Subscript[d, 1])   *)

BB=Mod[TransMat4L[(Subscript[c, 1])].TransMat4R[(Subscript[d, 1])],p]

(*   CC corresponds to ((Subscript[a, 2])x)(Subscript[b, 2])   *)

CC=Mod[TransMat4L[(Subscript[a, 2])].TransMat4R[(Subscript[b, 2])],p]

(*   DD corresponds to ((Subscript[c, 2])y)(Subscript[d, 2])   *)

DD=Mod[TransMat4L[(Subscript[c, 2])].TransMat4R[(Subscript[d, 2])],p]

ModPEP4x4[AA,BB,CC,DD,p]

x=Mod[W.(Subscript[k, 1])+X.(Subscript[k, 2]),p]

y=Mod[Y.(Subscript[k, 1])+Z.(Subscript[k, 2]),p]

Mod[QuatDotMP[QuatDotMP[(Subscript[a, 1]),x,p],(Subscript[b, 1]),p]+QuatDotMP[QuatDotMP[(Subscript[c, 1]),y,p],(Subscript[d, 1]),p],p]==(Subscript[k, 1])

True      (*  Checks that:      ((Subscript[a, 1])x)(Subscript[b, 1]) + ((Subscript[c, 1])y)(Subscript[d, 1]) = (Subscript[k, 1]).  *)

Mod[QuatDotMP[QuatDotMP[(Subscript[a, 2]),x,p],(Subscript[b, 2]),p]+QuatDotMP[QuatDotMP[(Subscript[c, 2]),y,p],(Subscript[d, 2]),p],p]==(Subscript[k, 2])

True      (*  Checks that:      ((Subscript[a, 2])x)(Subscript[b, 2]) + ((Subscript[c, 2])y)(Subscript[d, 1]) = (Subscript[k, 2]).  *)

Mod[QuatDotMP[(Subscript[a, 2]),QuatDotMP[x,(Subscript[b, 2]),p],p]+QuatDotMP[(Subscript[c, 2]),QuatDotMP[y,(Subscript[d, 2]),p],p],p]==(Subscript[k, 2])

True       (*  Checks that:      (Subscript[a, 2])((Subscript[xb, 2])) + (Subscript[c, 2])((Subscript[yd, 1])) = (Subscript[k, 2]), as well. It was mentioned that quaternions are associative.
        *)

(*    Now let us illustrate what happens when the coeficients are reversed. We again solve for x and y in:

    ((Subscript[b, 1])x)(Subscript[a, 1]) + ((Subscript[d, 1])y)(Subscript[c, 1]) = (Subscript[k, 1])  &
   ((Subscript[b, 2])x)(Subscript[a, 2]) + ((Subscript[d, 2])y)(Subscript[c, 2]) = (Subscript[k, 2]). 
	
	Quaternions are noncommutative so that multiplication from the left is different than multiplication 
from the right. For the most general case we choose an example that illustrates multiplication of the 
unknkowns from both directions. The grouping given above is not important because quaternions are associative.   *)

(*   AA corresponds to ((Subscript[a, 1])x)(Subscript[b, 1])   *)

AA=Mod[TransMat4L[(Subscript[b, 1])].TransMat4R[(Subscript[a, 1])],p]

(*   BB corresponds to ((Subscript[c, 1])y)(Subscript[d, 1])   *)

BB=Mod[TransMat4L[(Subscript[d, 1])].TransMat4R[(Subscript[c, 1])],p]

(*   CC corresponds to ((Subscript[a, 2])x)(Subscript[b, 2])   *)

CC=Mod[TransMat4L[(Subscript[b, 2])].TransMat4R[(Subscript[a, 2])],p]

(*   DD corresponds to ((Subscript[c, 2])y)(Subscript[d, 2])   *)

DD=Mod[TransMat4L[(Subscript[d, 2])].TransMat4R[(Subscript[c, 2])],p]

ModPEP4x4[AA,BB,CC,DD,p]

x=Mod[W.(Subscript[k, 1])+X.(Subscript[k, 2]),p]

{6,1,6,0}          (*  "old x" was = {0,6,0,2}  *)

y=Mod[Y.(Subscript[k, 1])+Z.(Subscript[k, 2]),p]

{6,4,6,1}          (*  "old y" was = {1,6,6,1}  *)

Mod[QuatDotMP[QuatDotMP[(Subscript[b, 1]),x,p],(Subscript[a, 1]),p]+QuatDotMP[QuatDotMP[(Subscript[d, 1]),y,p],(Subscript[c, 1]),p],p]==(Subscript[k, 1])

True       (*  Checks that:     ((Subscript[b, 1])x)(Subscript[a, 1]) + ((Subscript[d, 1])y)(Subscript[c, 1]) = (Subscript[k, 1]).   *)

Mod[QuatDotMP[QuatDotMP[(Subscript[b, 2]),x,p],(Subscript[a, 2]),p]+QuatDotMP[QuatDotMP[(Subscript[d, 2]),y,p],(Subscript[c, 2]),p],p]==(Subscript[k, 2])

True        (*  Checks that:     ((Subscript[b, 2])x)(Subscript[a, 2]) + ((Subscript[d, 2])y)(Subscript[c, 2]) = (Subscript[k, 2]).  *)

(*  Choosing p = 97, ...  *)

SeedRandom[9791];p=97;(Subscript[a, 1])=RandQuatModP[p];(Subscript[a, 2])=RandQuatModP[p];(Subscript[b, 1])=RandQuatModP[p];(Subscript[b, 2])=RandQuatModP[p];(Subscript[c, 1])=RandQuatModP[p];(Subscript[c, 2])=RandQuatModP[p];(Subscript[d, 1])=RandQuatModP[p];(Subscript[d, 2])=RandQuatModP[p];(Subscript[k, 1])=RandQuatModP[p];(Subscript[k, 2])=RandQuatModP[p]

(*    The system we solve for x and y is:

    ((Subscript[a, 1])x)(Subscript[b, 1]) + ((Subscript[c, 1])y)(Subscript[d, 1]) = (Subscript[k, 1])  &
   ((Subscript[a, 2])x)(Subscript[b, 2]) + ((Subscript[c, 2])y)(Subscript[d, 2]) = (Subscript[k, 2]). 
	
	Quaternions are noncommutative so that multiplication from the left is different than multiplication 
from the right. For the most general case we choose an example that illustrates multiplication of the 
unknkowns from both directions. The grouping given above is not important because quaternions are associative.   *)

(*   AA corresponds to ((Subscript[a, 1])x)(Subscript[b, 1])   *)

AA=Mod[TransMat4L[(Subscript[a, 1])].TransMat4R[(Subscript[b, 1])],p]

(*   BB corresponds to ((Subscript[c, 1])y)(Subscript[d, 1])   *)

BB=Mod[TransMat4L[(Subscript[c, 1])].TransMat4R[(Subscript[d, 1])],p]

(*   CC corresponds to ((Subscript[a, 2])x)(Subscript[b, 2])   *)

CC=Mod[TransMat4L[(Subscript[a, 2])].TransMat4R[(Subscript[b, 2])],p]

(*   DD corresponds to ((Subscript[c, 2])y)(Subscript[d, 2])   *)

DD=Mod[TransMat4L[(Subscript[c, 2])].TransMat4R[(Subscript[d, 2])],p]

ModPEP4x4[AA,BB,CC,DD,p]

x=Mod[W.(Subscript[k, 1])+X.(Subscript[k, 2]),p]

y=Mod[Y.(Subscript[k, 1])+Z.(Subscript[k, 2]),p]

Mod[QuatDotMP[QuatDotMP[(Subscript[a, 1]),x,p],(Subscript[b, 1]),p]+QuatDotMP[QuatDotMP[(Subscript[c, 1]),y,p],(Subscript[d, 1]),p],p]==(Subscript[k, 1])

True     (*  Checks that:      ((Subscript[a, 1])x)(Subscript[b, 1]) + ((Subscript[c, 1])y)(Subscript[d, 1]) = (Subscript[k, 1]).  *)

Mod[QuatDotMP[QuatDotMP[(Subscript[a, 2]),x,p],(Subscript[b, 2]),p]+QuatDotMP[QuatDotMP[(Subscript[c, 2]),y,p],(Subscript[d, 2]),p],p]==(Subscript[k, 2])

True     (*  Checks that:      ((Subscript[a, 2])x)(Subscript[b, 2]) + ((Subscript[c, 2])y)(Subscript[d, 1]) = (Subscript[k, 2]).  *)

Mod[QuatDotMP[(Subscript[a, 2]),QuatDotMP[x,(Subscript[b, 2]),p],p]+QuatDotMP[(Subscript[c, 2]),QuatDotMP[y,(Subscript[d, 2]),p],p],p]==(Subscript[k, 2])

True       (*  Checks that:      (Subscript[a, 2])((Subscript[xb, 2])) + (Subscript[c, 2])((Subscript[yd, 1])) = (Subscript[k, 2]), as well. It was mentioned that quaternions are associative.
        *)

(*    Now let us illustrate what happens when the coeficients are reversed. We again solve for x and y in:

    ((Subscript[b, 1])x)(Subscript[a, 1]) + ((Subscript[d, 1])y)(Subscript[c, 1]) = (Subscript[k, 1])  &
   ((Subscript[b, 2])x)(Subscript[a, 2]) + ((Subscript[d, 2])y)(Subscript[c, 2]) = (Subscript[k, 2]). 
	
	Quaternions are noncommutative so that multiplication from the left is different than multiplication 
from the right. For the most general case we choose an example that illustrates multiplication of the 
unknkowns from both directions. The grouping given above is not important because quaternions are associative.   *)

(*   AA corresponds to ((Subscript[a, 1])x)(Subscript[b, 1])   *)

AA=Mod[TransMat4L[(Subscript[b, 1])].TransMat4R[(Subscript[a, 1])],p]

(*   BB corresponds to ((Subscript[c, 1])y)(Subscript[d, 1])   *)

BB=Mod[TransMat4L[(Subscript[d, 1])].TransMat4R[(Subscript[c, 1])],p]

(*   CC corresponds to ((Subscript[a, 2])x)(Subscript[b, 2])   *)

CC=Mod[TransMat4L[(Subscript[b, 2])].TransMat4R[(Subscript[a, 2])],p]

(*   DD corresponds to ((Subscript[c, 2])y)(Subscript[d, 2])   *)

DD=Mod[TransMat4L[(Subscript[d, 2])].TransMat4R[(Subscript[c, 2])],p]

ModPEP4x4[AA,BB,CC,DD,p]

x=Mod[W.(Subscript[k, 1])+X.(Subscript[k, 2]),p]

{43,5,25,35}      (*  "old x" was = {36,37,93,83}  *)

y=Mod[Y.(Subscript[k, 1])+Z.(Subscript[k, 2]),p]

{15,28,88,5}              (*  "old y" was = {51,17,44,74}  *)

Mod[QuatDotMP[QuatDotMP[(Subscript[b, 1]),x,p],(Subscript[a, 1]),p]+QuatDotMP[QuatDotMP[(Subscript[d, 1]),y,p],(Subscript[c, 1]),p],p]==(Subscript[k, 1])

True      (*  Checks that:     ((Subscript[b, 1])x)(Subscript[a, 1]) + ((Subscript[d, 1])y)(Subscript[c, 1]) = (Subscript[k, 1]).   *)

Mod[QuatDotMP[QuatDotMP[(Subscript[b, 2]),x,p],(Subscript[a, 2]),p]+QuatDotMP[QuatDotMP[(Subscript[d, 2]),y,p],(Subscript[c, 2]),p],p]==(Subscript[k, 2])

True      (*  Checks that:     ((Subscript[b, 2])x)(Subscript[a, 2]) + ((Subscript[d, 2])y)(Subscript[c, 2]) = (Subscript[k, 2]).  *)

(*  Below is an example of where no solution exists.  *)

SeedRandom[9797];p=7;(Subscript[a, 1])=RandQuatModP[p];(Subscript[a, 2])=RandQuatModP[p];(Subscript[b, 1])=RandQuatModP[p];(Subscript[b, 2])=RandQuatModP[p];(Subscript[c, 1])=RandQuatModP[p];(Subscript[c, 2])=RandQuatModP[p];(Subscript[d, 1])=RandQuatModP[p];(Subscript[d, 2])=RandQuatModP[p];(Subscript[k, 1])=RandQuatModP[p];(Subscript[k, 2])=RandQuatModP[p]

Mod[Sum[(((Subscript[d, 2])[[i]])^2),{i,4}],7]

AA=Mod[TransMat4L[(Subscript[a, 1])].TransMat4R[(Subscript[b, 1])],p]

(*   BB corresponds to ((Subscript[c, 1])y)(Subscript[d, 1])   *)

BB=Mod[TransMat4L[(Subscript[c, 1])].TransMat4R[(Subscript[d, 1])],p]

(*   CC corresponds to ((Subscript[a, 2])x)(Subscript[b, 2])   *)

CC=Mod[TransMat4L[(Subscript[a, 2])].TransMat4R[(Subscript[b, 2])],p]

(*   DD corresponds to ((Subscript[c, 2])y)(Subscript[d, 2])   *)

DD=Mod[TransMat4L[(Subscript[c, 2])].TransMat4R[(Subscript[d, 2])],p]

ModPEP4x4[AA,BB,CC,DD,p]

AA=Mod[TransMat4L[(Subscript[b, 1])].TransMat4R[(Subscript[a, 1])],p]

(*   BB corresponds to ((Subscript[c, 1])y)(Subscript[d, 1])   *)

BB=Mod[TransMat4L[(Subscript[d, 1])].TransMat4R[(Subscript[c, 1])],p]

(*   CC corresponds to ((Subscript[a, 2])x)(Subscript[b, 2])   *)

CC=Mod[TransMat4L[(Subscript[b, 2])].TransMat4R[(Subscript[a, 2])],p]

(*   DD corresponds to ((Subscript[c, 2])y)(Subscript[d, 2])   *)

DD=Mod[TransMat4L[(Subscript[d, 2])].TransMat4R[(Subscript[c, 2])],p]

ModPEP4x4[AA,BB,CC,DD,p]

x=Mod[W.(Subscript[k, 1])+X.(Subscript[k, 2]),p]

y=Mod[Y.(Subscript[k, 1])+Z.(Subscript[k, 2]),p]

Mod[QuatDotMP[QuatDotMP[(Subscript[b, 1]),x,p],(Subscript[a, 1]),p]+QuatDotMP[QuatDotMP[(Subscript[d, 1]),y,p],(Subscript[c, 1]),p],p]==(Subscript[k, 1])

True      (*  Checks that:     ((Subscript[b, 1])x)(Subscript[a, 1]) + ((Subscript[d, 1])y)(Subscript[c, 1]) = (Subscript[k, 1]).   *)

Mod[QuatDotMP[QuatDotMP[(Subscript[b, 2]),x,p],(Subscript[a, 2]),p]+QuatDotMP[QuatDotMP[(Subscript[d, 2]),y,p],(Subscript[c, 2]),p],p]==(Subscript[k, 2])

True         (*  Checks that:     ((Subscript[b, 2])x)(Subscript[a, 2]) + ((Subscript[d, 2])y)(Subscript[c, 2]) = (Subscript[k, 2]).  *)

(*         Thus, we see that the order of the constants is crucial. One way there was a solution
& the other way there was not. There are 14 other possibilities to check out. Your homework is ...      *)