Mutually Unbiased Bases for prime-power dimensions d=p^n // v1.0 for mathematica // 04.10.2013

(c)	Christoph Spengler 2013, Theoretical Physics, University of Innsbruck
	Contact : Christoph.Spengler@uibk.ac.at

Reference : 	-- PLEASE CITE THIS PAPER WHEN USING THIS FILE --

[1]	C. Spengler, B. Kraus
	'A graph state formalism for mutually unbiased bases'
	arXiv:1309.6557
	
Don’t hesitate to contact me if you have any questions, comments or problems!

Note: Irreducible polynomials may be taken from http://theory.cs.uvic.ca/gen/poly.html or the tables and references given in Ref.[1].

Code: The code is hidden in this cell. Run it!

QMatrix::usage = "For a given monic irreducible polynomial f(x)=c(0)+c(1)*x+c(2)*x^2+...+c(n-1)*x^(n-1)+x^n of degree n with coefficients c={c(0),...,c(n-1)} from the prime field Z_p={0,...,p-1}, the function Q=Qmatrix[c,p] symmetrizes the associated companion matrix. Specifically, it holds that char(Q)=det(x*Id-Q)=f(x) where Q=Transpose(Q).";
QMatrix[c0_,p0_]:=If[p0==2,SYM2[Mod[c0,2]],SYMP[Mod[c0,p0],p0]];

SYM2::Error1="Used coefficients c={c(0),...,c(n-1)} do not correspond to an irreducible polynomial.";
SYM2[c0_]:=Module[{c=c0,CREATECf,CREATEB,n,Cf,s,b0,b,B,Π,Λ,Ω,Ωinv,M,P,Pinv,TEMP1,STEP1,STEP2,STEP3,SYMMETRIZE,G},CREATECf:={n=Dimensions[c][[1]],Cf=0*IdentityMatrix[n],Do[Cf[[k-1,k]]=1,{k,2,n}],Do[Cf[[n,k]]=c[[k]],{k,1,n}],Cf}[[5]];
CREATEB:={n=Dimensions[c][[1]],b=Array[b0,n-1],b[[1]]=c[[1]],Do[b[[i]]=Mod[Sum[c[[n-i+k+1]]*b[[k]],{k,1,i-1}],2],{i,2,n-1}],B=0*IdentityMatrix[n],B[[1,1]]=1,Do[Do[B[[i,n-i+k+1]]=b[[k]],{i,k+1,n}],{k,1,n-1}],B}[[8]];
Π[i_,j_]:={M=IdentityMatrix[n],M[[i,i]]=0,M[[j,j]]=0,M[[i,j]]=1,M[[j,i]]=1,M}[[6]];
Λ[i_,j_]:={M=IdentityMatrix[n],M[[j,i]]=1,M}[[3]];
Ω[i_,j_,k_]:={M=IdentityMatrix[n],M[[j,j]]=0,M[[i,j]]=1,M[[j,i]]=1,M[[k,i]]=1,M[[j,k]]=1,M[[k,j]]=1,M}[[8]];
Ωinv[i_,j_,k_]:={M=IdentityMatrix[n],M[[k,i]]=1,M[[i,j]]=1,M[[i,k]]=1,M[[j,k]]=1,M[[k,j]]=0,M}[[7]];
STEP1[s_]:={K=Drop[Diagonal[B],s-1],If[K[[Ordering[K,-1]]][[1]]==1&&Ordering[K,-1][[1]]≠ 1,{P=Mod[Π[s,Ordering[K,-1]+s-1].P,2],Pinv=Mod[Pinv.Π[s,Ordering[K,-1]+s-1],2],B=Π[s,Ordering[K,-1]+s-1].B.Π[s,Ordering[K,-1]+s-1]}]};
STEP2[s_]:={TEMP1=Drop[Take[B,{s}][[1]],s],If[TEMP1[[Ordering[TEMP1,-1]]][[1]]==1,{P=Mod[Λ[s,Ordering[TEMP1,-1]+s].P,2],Pinv=Mod[Pinv.Λ[s,Ordering[TEMP1,-1]+s],2],B=Mod[Λ[s,Ordering[TEMP1,-1]+s].B.Transpose[Λ[s,Ordering[TEMP1,-1]+s]],2]}],B}[[3]];
STEP3[s_]:={TEMP1=Drop[Take[B,{s}][[1]],s],If[TEMP1[[Ordering[TEMP1,-1]]][[1]]==1,{P=Mod[Ω[s-1,s,Ordering[TEMP1,-1]+s].P,2],Pinv=Mod[Pinv.Ωinv[s-1,s,Ordering[TEMP1,-1]+s],2],B=Mod[Ω[s-1,s,Ordering[TEMP1,-1]+s].B.Transpose[Ω[s-1,s,Ordering[TEMP1,-1]+s]],2]}],B}[[2]];
SYMMETRIZE:={s=2,CREATECf,CREATEB,P=IdentityMatrix[n],Pinv=IdentityMatrix[n],While[s≠n,If[Max[Drop[Diagonal[B],s-1]]≠0,{STEP1[s],While[Max[Drop[Take[B,{s}][[1]],s]]≠0,STEP2[s]],s=s+1},{STEP3[s],s=s-1}]],G=Mod[P.Cf.Pinv,2]}[[7]];If[c[[1]]==0,G=Message[SYM2::Error1],SYMMETRIZE];G]

SYMP::Error1="Please use coefficients c={c(0),...,c(n-1)} of an irreducible polynomial";
SYMP::Error2="Please use coefficients c={c(0),...,c(n-1)} of a primitive polynomial for (n mod 4)=2.";
SYMP::Error3="Second argument 'p' has to be a prime number.";
SYMP[c0_,p0_]:=Module[{c=c0,p=p0,IntNums,Int,QR,QN,NumInv,a,n,CREATECf,Cf,CREATEB,b,B,MQ,findsquare,temps,Π,Λ,Λinv,Ω,Ωinv,M,Φ,Φinv,ctemp,b0,STEP1,STEP2,STEP3,K,P,Pinv,s,ss,TEMP1,SYMMETRIZE,trafo,Q,k,nonres},
IntNums:={Do[Int[i]=i-1,{i,1,p}],Array[Int,p]}[[2]];
QR:=Complement[IntNums,Complement[IntNums,Mod[IntNums^2,p]],{0}];
QN:=Complement[IntNums,Mod[IntNums^2,p]];
NumInv[x_]:={a=1,While[Mod[x*a,p]≠1,a=a+1],a}[[3]];
CREATECf:={n=Dimensions[c][[1]],Cf=0*IdentityMatrix[n],Do[Cf[[k-1,k]]=1,{k,2,n}],Do[Cf[[n,k]]=-c[[k]],{k,1,n}],Cf=Mod[Cf,p]}[[5]];
CREATEB:={n=Dimensions[c][[1]],b=Array[b0,n],b[[1]]=1,Do[b[[i]]=Mod[Sum[-c[[n-i+k+1]]*b[[k]],{k,1,i-1}],p],{i,2,n}],B=0*IdentityMatrix[n],Do[Do[B[[i,n-i+k]]=b[[k]],{i,k,n}],{k,1,n}],MQ=MemberQ[QR,Mod[Det[B],p]],If[MQ,B=B,Switch[Mod[n,4],3,B=QN[[1]]*B,2,B=Cf.B]],B=Mod[B,p]}[[9]];
findsquare:={temps=1,While[Mod[1+temps,p]≠ QN[[1]],temps=temps+1],temps}[[3]];
Π[i_,j_]:={M=IdentityMatrix[n],M[[i,i]]=0,M[[j,j]]=0,M[[i,j]]=1,M[[j,i]]=1,M=Mod[M,p]}[[6]];
Λ[i_,j_,a_]:={M=IdentityMatrix[n],M[[j,i]]=a,M=Mod[M,p]}[[3]];
Λinv[i_,j_,a_]:={M=IdentityMatrix[n],M[[j,i]]=-a,M=Mod[M,p]}[[3]];
Ω[i_,j_]:={M=IdentityMatrix[n],M[[j,j]]=-1,M[[i,j]]=1,M[[j,i]]=1,M=Mod[M,p]}[[5]];
Ωinv[i_,j_]:={M=IdentityMatrix[n],M[[i,i]]=1*(p+1)/2,M[[j,j]]=-1*(p+1)/2,M[[i,j]]=1*(p+1)/2,M[[j,i]]=1*(p+1)/2,M=Mod[M,p]}[[6]];
Φ[i_,j_,b_]:={M=IdentityMatrix[n],M[[j,i]]=-b,M[[i,j]]=b,M=Mod[M,p]}[[4]];
Φinv[i_,j_,b_]:={M=IdentityMatrix[n],ctemp=NumInv[Mod[1+b^2,p]],M[[i,i]]=ctemp,M[[j,j]]=ctemp,M[[j,i]]=ctemp*b,M[[i,j]]=-b*ctemp,M=Mod[M,p]}[[7]];
STEP1[s_]:={K=Drop[Diagonal[B],s-1],If[K[[Ordering[K,-1]]][[1]]≠ 0&&B[[s,s]]==0,{P=Mod[Π[s,Ordering[K,-1]+s-1].P,p],Pinv=Mod[Pinv.Π[s,Ordering[K,-1]+s-1],p],B=Π[s,Ordering[K,-1]+s-1].B.Π[s,Ordering[K,-1]+s-1]}],B}[[3]];
STEP2[s_]:={TEMP1=Drop[Take[B,{s}][[1]],s],If[TEMP1[[Ordering[TEMP1,-1]]][[1]]≠ 0,{a=1,While[Mod[a*B[[s,s]]+TEMP1[[Ordering[TEMP1,-1]]][[1]],p]≠0,a=a+1],P=Mod[Λ[s,Ordering[TEMP1,-1][[1]]+s,a].P,p],Pinv=Mod[Pinv.Λinv[s,Ordering[TEMP1,-1][[1]]+s,a],p],B=Mod[Λ[s,Ordering[TEMP1,-1][[1]]+s,a].B.Transpose[Λ[s,Ordering[TEMP1,-1][[1]]+s,a]],p]}],B}[[3]];
STEP3[s_]:={TEMP1=Drop[Take[B,{s}][[1]],s],If[TEMP1[[Ordering[TEMP1,-1]]][[1]]≠0,{P=Mod[Ω[s,Ordering[TEMP1,-1]+s].P,p],Pinv=Mod[Pinv.Ωinv[s,Ordering[TEMP1,-1]+s],p],B=Mod[Ω[s,Ordering[TEMP1,-1]+s].B.Transpose[Ω[s,Ordering[TEMP1,-1]+s]],p]}],B}[[3]];
SYMMETRIZE:={{s=1,CREATECf,CREATEB,IntNums,QR,QN,P=IdentityMatrix[n],Pinv=IdentityMatrix[n],If[(MemberQ[Union[QR,{0}],Mod[Det[Cf],p]]&&Mod[n,4]==2)||c[[1]]==0,If[MemberQ[Union[QR,{0}],Mod[Det[Cf],p]]&&Mod[n,4]==2,Q=Message[SYMP::Error2],Q=Message[SYMP::Error1]],{While[s≠n,If[Max[Drop[Diagonal[B],s-1]]≠0,{STEP1[s],While[Max[Drop[Take[B,{s}][[1]],s]]≠0,STEP2[s]],s=s+1},{STEP3[s],s=s-1}]],Clear[nonres],k=1,ss=findsquare,Do[If[MemberQ[QN,B[[i,i]]],{a=1,While[Mod[a^2*B[[i,i]],p]≠QN[[1]],a=a+1],trafo=IdentityMatrix[n],trafo[[i,i]]=a,P=trafo.P,B=Mod[trafo.B.trafo,p],trafo=IdentityMatrix[n],trafo[[i,i]]=NumInv[a],Pinv=Mod[Pinv.trafo,p],nonres[k]=i,k=k+1}],{i,1,n}],nonres=Array[nonres,k-1],Do[{P=Mod[Φ[nonres[[i]],nonres[[i+1]],ss].P,p],Pinv=Mod[Pinv.Φinv[nonres[[i]],nonres[[i+1]],ss],p],B=Mod[Φ[nonres[[i]],nonres[[i+1]],ss].B.Transpose[Φ[nonres[[i]],nonres[[i+1]],ss]],p]},{i,1,(k-1),2}],Do[{a=1,While[Mod[a^2*B[[i,i]],p]≠1,a=a+1],trafo=IdentityMatrix[n],trafo[[i,i]]=a,P=Mod[trafo.P,p],B=Mod[trafo.B.trafo,p],trafo=IdentityMatrix[n],trafo[[i,i]]=NumInv[a],Pinv=Mod[Pinv.trafo,p]},{i,1,n}],B,Pinv=Mod[Pinv,p],Q=Mod[P.Cf.Pinv,p]}]},{Q}}[[2,1]];If[Element[p,Primes],SYMMETRIZE,Q=Message[SYMP::Error3]];Q]

AdjMatrix::Error="Coefficient vector and matrix have incompatible shapes.";
AdjMatrix::usage = "AdjMatrix[a,Q,p] yields a linear combination of the powers i={0,...,n-1} of the nxn matrix Q using the coefficients a={a(0),...,a(n-1)} over Z_p, i.e. AdjMatrix[a,Q,p]=a(0)*Id+a(1)*Q+a(2)*Q^2+...+a(n-1)*Q^(n-1) mod p.";
AdjMatrix[c_,Q_,p_]:=Module[{c0=c,Q0=Q,p0=p},
If[Dimensions[c0][[1]]==Dimensions[Q0][[1]],Mod[Sum[c0[[i]]*MatrixPower[Q0,i-1],{i,2,Dimensions[Q0][[1]]}]+c0[[1]]*IdentityMatrix[Dimensions[Q0][[1]]],p0],Message[AdjMatrix::Error]]]

Graphstatebasis::Error1="Dimension p must be a prime number.";
Graphstatebasis::Error2="Adjacency matrix A must be symmetric.";
Graphstatebasis::usage = "Graphstatebasis[A,p] yields the graph state basis corresponding to the (n x n) adjacency matrix A for a multi-qupit system H=(C^p)^n where p is a prime number.";
Graphstatebasis[A0_,p0_]:=Module[{P,p1,NLPos,TEMP,Z,n,LPos,w,L2Pos,U,FB,TEMP2,generatebasis,TEMPBASIS,V},
P[i_]:=Transpose[{UnitVector[p1,i+1]}].{UnitVector[p1,i+1]};
NLPos[i_,k_,l_,r_]:={TEMP=P[i],TEMP=KroneckerProduct[IdentityMatrix[p1^Max[0,k-1]],TEMP],TEMP=KroneckerProduct[TEMP,IdentityMatrix[p1^Max[0,l-k-1]]],TEMP=KroneckerProduct[TEMP,MatrixPower[Z,(i*r)]],TEMP=KroneckerProduct[TEMP,IdentityMatrix[p1^Max[0,n-l]]]}[[5]];
LPos[i_,k_,r_]:={TEMP=w^(i*(i-1)/2*r)*P[i],TEMP=KroneckerProduct[IdentityMatrix[p1^Max[0,k-1]],TEMP],TEMP=KroneckerProduct[TEMP,IdentityMatrix[p1^Max[0,n-k]]]}[[3]];
L2Pos[i_,k_,r_]:={TEMP=(I^(i*r))*P[i],TEMP=KroneckerProduct[IdentityMatrix[p1^Max[0,k-1]],TEMP],TEMP=KroneckerProduct[TEMP,IdentityMatrix[p1^Max[0,n-k]]]}[[3]];
U[k_,l_,r_]:=If[k==l,If[p1==2,Sum[L2Pos[i,k,r],{i,0,p1-1}],Sum[LPos[i,k,r],{i,0,p1-1}]],Sum[NLPos[i,k,l,r],{i,0,p1-1}]];
FB[p_]:={TEMP=IdentityMatrix[p],Do[TEMP[[k+1,l+1]]=1/Sqrt[p]*w^(k*l),{k,0,p-1},{l,0,p-1}],TEMP2=IdentityMatrix[1],Do[TEMP2=KroneckerProduct[TEMP2,TEMP],{z,1,n}],TEMP2}[[5]];
generatebasis[A_,p_]:={n=Dimensions[A][[1]],p1=p,w=Exp[2*Pi*I/p],Z=IdentityMatrix[p],Do[Z[[i,i]]=w^(i-1),{i,1,p}],TEMPBASIS=FB[p],Do[If[A[[k,l]]≠0,TEMPBASIS=U[k,l,A[[k,l]]].TEMPBASIS],{k,1,n},{l,k,n}],TEMPBASIS}[[8]];If[Element[p0,Primes],If[A0==Transpose[A0],V=generatebasis[A0,p0],V=Message[Graphstatebasis::Error2]],V=Message[Graphstatebasis::Error1]];V]

MUBCheck::Error1="Bases have incompatible shapes";
MUBCheck::usage ="Given two bases in the form of two (d x d) matrices Bx and By, the function MUBCheck[Bx,By] checks if they are mutually unbiased.";
MUBCheck[Bx_,By_]:=If[Dimensions[Bx]==Dimensions[By],Module[{T=Chop[Abs[ConjugateTranspose[N[Bx]].N[By]]^2],Maxi,Mini},Maxi=Max[T];Mini=Min[T];Abs[Maxi-Mini]<10^(-6)&&Abs[Maxi-Chop[N[1/(Dimensions[Bx][[1]])]]]<10^(-6)],Message[MUBCheck::Error1]]

AdjCheck::usage ="Given two adjacency matrices Ax and Ay over Z_p, the function MUBCheck[Ax,Ay] checks if the corresponding graph state bases are mutually unbiased.";
AdjCheck[Ax_,Ay_,p_]:=Mod[Det[Ax-Ay],p]≠0

Usage:

1) Run the code in the cell above.
2) For a given monic irreducible polynomial f(x)=c(0)+c(1)*x+c(2)*x^2+...+c(n-1)*x^(n-1)+x^n of degree n with coefficients c={c(0),...,c(n-1)} in Z_p={0,...,p-1} where p is a prime number, use the function Q=QMatrix[c,p] to obtain the associated symmetrized companion matrix Q. (see definition and algorithm introduced in Ref.[1])
3) Use the function A=AdjMatrix[a,Q,p] with coefficients a={a(0),...,a(n-1)} from Z_p to obtain the linear combination a(0)*Id+a(1)*Q+a(2)*Q^2+...+a(n-1)*Q^(n-1), which for each setting of a={a(0),...,a(n-1)} is a valid adjacency matrix from the set S as described in Ref.[1].
4) Use the function B=Graphstatebasis[A,p] to generate the graph state basis corresponding to the adjacency matrix A over Z_p.
5) The set of bases {B_0,...,B_(p^n-1)} for the p^n settings of a={a(0),...,a(n-1)} together with the computational Basis B_C=IdentityMatrix[p^n] form a complete set of MUBs!
6) Optional: Use MUBCheck[Bx,By] to check if two basis, say Bx and By, are mutually unbiased. Use AdjCheck[Ax,Ay,p] to check if two matrices, Ax and Ay, satisfy det(Ax-Ay) unequal 0.

Example 1 (p=2,n=3) ==> please open:

Over Z_2 the polynomial f(x)=x^3+x^1+1 is irreducible (see Ref.[1] and link above). This polynomial has the coefficients (in ascending order) c={1,1,0}. First, we use Q=QMatrix to obtain a symmetric matrix Q which has f(x) as its characteristic polynomial, i.e. it holds char(Q)=det(x*Id-Q) mod 2=f(x) where Q=Transpose[Q].

Q=QMatrix[{1,1,0},2]

The possible 2^3=8 linear combinations of the powers of Q^i with i=0,..,2 are obtained via AdjMatrix

A0=AdjMatrix[{0,0,0},Q,2]

A1=AdjMatrix[{1,0,0},Q,2]

A2=AdjMatrix[{0,1,0},Q,2]

A3=AdjMatrix[{1,1,0},Q,2]

A4=AdjMatrix[{0,0,1},Q,2]

A5=AdjMatrix[{1,0,1},Q,2]

A6=AdjMatrix[{0,1,1},Q,2]

A7=AdjMatrix[{1,1,1},Q,2]

The function Graphstatebasis yields the graph state basis corresponding to each adjacency matrix A.

B0=Graphstatebasis[A0,2]

B1=Graphstatebasis[A1,2]

B2=Graphstatebasis[A2,2]

B3=Graphstatebasis[A3,2]

B4=Graphstatebasis[A4,2]

B5=Graphstatebasis[A5,2]

B6=Graphstatebasis[A6,2]

B7=Graphstatebasis[A7,2]

Together with the computational basis BC, this is a complete set of 9 MUBs for d=2^3=8

BC=IdentityMatrix[2^3]

Let’s see if they are really mutually unbiased:

{MUBCheck[BC,B0],MUBCheck[BC,B1],MUBCheck[BC,B2],MUBCheck[BC,B3],MUBCheck[BC,B4],MUBCheck[BC,B5],MUBCheck[BC,B6],MUBCheck[BC,B7],MUBCheck[B1,B2],MUBCheck[B1,B3],MUBCheck[B1,B4],MUBCheck[B1,B5],MUBCheck[B1,B6],MUBCheck[B1,B7],MUBCheck[B2,B3],MUBCheck[B2,B4],MUBCheck[B2,B5],MUBCheck[B2,B6],MUBCheck[B2,B7],MUBCheck[B3,B4],MUBCheck[B3,B5],MUBCheck[B3,B6],MUBCheck[B3,B7],MUBCheck[B4,B5],MUBCheck[B4,B6],MUBCheck[B4,B7],MUBCheck[B5,B6],MUBCheck[B5,B7],MUBCheck[B6,B7]}

...or check if MUB condition on adjacency matrices is satisfied:

{AdjCheck[A0,A1,2],AdjCheck[A0,A2,2],AdjCheck[A0,A3,2],AdjCheck[A0,A4,2],AdjCheck[A0,A5,2],AdjCheck[A0,A6,2],AdjCheck[A0,A7,2],AdjCheck[A1,A2,2],AdjCheck[A1,A3,2],AdjCheck[A1,A4,2],AdjCheck[A1,A5,2],AdjCheck[A1,A6,2],AdjCheck[A1,A7,2],AdjCheck[A2,A3,2],AdjCheck[A2,A4,2],AdjCheck[A2,A5,2],AdjCheck[A2,A6,2],AdjCheck[A2,A7,2],AdjCheck[A3,A4,2],AdjCheck[A3,A5,2],AdjCheck[A3,A6,2],AdjCheck[A3,A7,2],AdjCheck[A4,A5,2],AdjCheck[A4,A6,2],AdjCheck[A4,A7,2],AdjCheck[A5,A6,2],AdjCheck[A5,A7,2],AdjCheck[A6,A7,2]}

Example 2 (p=3,n=2):

Over Z_3 the polynomial f(x)=x^2+x+2 is irreducible (and primitive, see Ref.[1] and link above). This polynomial has the coefficients (in ascending order) c={2,1}. First, we use Q=QMatrix to obtain a symmetric matrix Q which has f(x) as its characteristic polynomial, i.e. it holds char(Q)=det(x*Id-Q) mod 3=f(x) where Q=Transpose[Q].

Q=QMatrix[{2,1},3]

The possible 3^2=9 linear combinations of the powers of Q^i with i=0,1 are obtained via AdjMatrix

A0=AdjMatrix[{0,0},Q,3]

A1=AdjMatrix[{1,0},Q,3]

A2=AdjMatrix[{2,0},Q,3]

A3=AdjMatrix[{0,1},Q,3]

A4=AdjMatrix[{1,1},Q,3]

A5=AdjMatrix[{2,1},Q,3]

A6=AdjMatrix[{0,2},Q,3]

A7=AdjMatrix[{1,2},Q,3]

A8=AdjMatrix[{2,2},Q,3]

The function Graphstatebasis yields the graph state basis corresponding to each adjacency matrix A.

B0=Graphstatebasis[A0,3]

B1=Graphstatebasis[A1,3]

B2=Graphstatebasis[A2,3]

B3=Graphstatebasis[A3,3]

B4=Graphstatebasis[A4,3]

B5=Graphstatebasis[A5,3]

B6=Graphstatebasis[A6,3]

B7=Graphstatebasis[A7,3]

B8=Graphstatebasis[A8,3]

Together with the computational basis BC, this is a complete set of 10 MUBs for d=3^2=9

BC=IdentityMatrix[3^2]

Let’s see if they are really mutually unbiased:

{MUBCheck[BC,B0],MUBCheck[BC,B1],MUBCheck[BC,B2],MUBCheck[BC,B3],MUBCheck[BC,B4],MUBCheck[BC,B5],MUBCheck[BC,B6],MUBCheck[BC,B7],MUBCheck[BC,B8],MUBCheck[B1,B2],MUBCheck[B1,B3],MUBCheck[B1,B4],MUBCheck[B1,B5],MUBCheck[B1,B6],MUBCheck[B1,B7],MUBCheck[B1,B8],MUBCheck[B2,B3],MUBCheck[B2,B4],MUBCheck[B2,B5],MUBCheck[B2,B6],MUBCheck[B2,B7],MUBCheck[B2,B8],MUBCheck[B3,B4],MUBCheck[B3,B5],MUBCheck[B3,B6],MUBCheck[B3,B7],MUBCheck[B3,B8],MUBCheck[B4,B5],MUBCheck[B4,B6],MUBCheck[B4,B7],MUBCheck[B4,B8],MUBCheck[B5,B6],MUBCheck[B5,B7],MUBCheck[B5,B8],MUBCheck[B6,B7],MUBCheck[B6,B8],MUBCheck[B7,B8]}

...or check if MUB condition on adjacency matrices is satisfied:

{AdjCheck[A0,A1,3],AdjCheck[A0,A2,3],AdjCheck[A0,A3,3],AdjCheck[A0,A4,3],AdjCheck[A0,A5,3],AdjCheck[A0,A6,3],AdjCheck[A0,A7,3],AdjCheck[A0,A8,3],AdjCheck[A1,A2,3],AdjCheck[A1,A3,3],AdjCheck[A1,A4,3],AdjCheck[A1,A5,3],AdjCheck[A1,A6,3],AdjCheck[A1,A7,3],AdjCheck[A1,A8,3],AdjCheck[A2,A3,3],AdjCheck[A2,A4,3],AdjCheck[A2,A5,3],AdjCheck[A2,A6,3],AdjCheck[A2,A7,3],AdjCheck[A2,A8,3],AdjCheck[A3,A4,3],AdjCheck[A3,A5,3],AdjCheck[A3,A6,3],AdjCheck[A3,A7,3],AdjCheck[A3,A8,3],AdjCheck[A4,A5,3],AdjCheck[A4,A6,3],AdjCheck[A4,A7,3],AdjCheck[A4,A8,3],AdjCheck[A5,A6,3],AdjCheck[A5,A7,3],AdjCheck[A5,A8,3],AdjCheck[A6,A7,3],AdjCheck[A6,A8,3],AdjCheck[A7,A8,3]}

Example 3 (p=5,n=1) :

This also works for a single qupit of prime dimension p. Just use the p different 1x1 (adjacency) matrices {{r}} with r=0,...,p-1 as the argument of Graphstatebasis

B0=Graphstatebasis[{{0}},5]

B1=Graphstatebasis[{{1}},5]

B2=Graphstatebasis[{{2}},5]

B3=Graphstatebasis[{{3}},5]

B4=Graphstatebasis[{{4}},5]

BC=IdentityMatrix[5]

Check mutually unbiasedness:

{MUBCheck[BC,B0],MUBCheck[BC,B1],MUBCheck[BC,B2],MUBCheck[BC,B3],MUBCheck[BC,B4],MUBCheck[B0,B1],MUBCheck[B0,B2],MUBCheck[B0,B3],MUBCheck[B0,B4],MUBCheck[B1,B2],MUBCheck[B1,B3],MUBCheck[B1,B4],MUBCheck[B2,B3],MUBCheck[B2,B4],MUBCheck[B3,B4]}