Checks of the analytic solution

Code generating solution

For convenience, we define a function which creates charges in a list. Charges take the form of objects where (Subscript[P, Q1])returns this charge of the object. l[P] returns a list of charges. 
This is so our formula match the those in the paper.

l[X_]:=Evaluate[Thread[(Subscript[X, {Q1,Q2,Q3,ν1,ν2,ν3,e1,e2,e3,u1,u2,u3,d1,d2,d3,L1,L2,L3,Hd,Hu}])]]

Special points (code assigns e.g. (Subscript[Y3, Q1]) to 1 etc)

{(Subscript[Y3, Q1]),(Subscript[Y3, Q2]),(Subscript[Y3, Q3]),(Subscript[Y3, ν1]),(Subscript[Y3, ν2]),(Subscript[Y3, ν3]),(Subscript[Y3, e1]),(Subscript[Y3, e2]),(Subscript[Y3, e3]),(Subscript[Y3, u1]),(Subscript[Y3, u2]),(Subscript[Y3, u3]),(Subscript[Y3, d1]),(Subscript[Y3, d2]),(Subscript[Y3, d3]),(Subscript[Y3, L1]),(Subscript[Y3, L2]),(Subscript[Y3, L3]),(Subscript[Y3, Hd]),(Subscript[Y3, Hu])}={1,1,-1,  0,0,0,  -6,-6,6,  4,4,-4,  -2,-2,2,  -3,-3,3 , -3,3};
{(Subscript[B3, Q1]),(Subscript[B3, Q2]),(Subscript[B3, Q3]),(Subscript[B3, ν1]),(Subscript[B3, ν2]),(Subscript[B3, ν3]),(Subscript[B3, e1]),(Subscript[B3, e2]),(Subscript[B3, e3]),(Subscript[B3, u1]),(Subscript[B3, u2]),(Subscript[B3, u3]),(Subscript[B3, d1]),(Subscript[B3, d2]),(Subscript[B3, d3]),(Subscript[B3, L1]),(Subscript[B3, L2]),(Subscript[B3, L3]),(Subscript[B3, Hd]),(Subscript[B3, Hu])}={1,1,-1, -3,-3,3, -3,-3,3, 1,1,-1, 1,1,-1, -3,-3,3, -3,3};

Quadratic function

q[X_,Z_]:=((Subscript[X, Q1])(Subscript[Z, Q1])+(Subscript[X, Q2])(Subscript[Z, Q2])+(Subscript[X, Q3])(Subscript[Z, Q3]))+((Subscript[X, e1])(Subscript[Z, e1])+(Subscript[X, e2])(Subscript[Z, e2])+(Subscript[X, e3])(Subscript[Z, e3]))-2((Subscript[X, u1])(Subscript[Z, u1])+(Subscript[X, u2])(Subscript[Z, u2])+(Subscript[X, u3])(Subscript[Z, u3]))+((Subscript[X, d1])(Subscript[Z, d1])+(Subscript[X, d2])(Subscript[Z, d2])+(Subscript[X, d3])(Subscript[Z, d3]))-((Subscript[X, L1])(Subscript[Z, L1])+(Subscript[X, L2])(Subscript[Z, L2])+(Subscript[X, L3])(Subscript[Z, L3]))-(Subscript[X, Hd])(Subscript[Z, Hd])+(Subscript[X, Hu])(Subscript[Z, Hu]);

Cubic function

c[W_,X_,Z_]:=6((Subscript[W, Q1])(Subscript[X, Q1])(Subscript[Z, Q1])+(Subscript[W, Q2])(Subscript[X, Q2])(Subscript[Z, Q2])+(Subscript[W, Q3])(Subscript[X, Q3])(Subscript[Z, Q3]))-((Subscript[W, ν1])(Subscript[X, ν1])(Subscript[Z, ν1])+(Subscript[W, ν2])(Subscript[X, ν2])(Subscript[Z, ν2])+(Subscript[W, ν3])(Subscript[X, ν3])(Subscript[Z, ν3]))-((Subscript[W, e1])(Subscript[X, e1])(Subscript[Z, e1])+(Subscript[W, e2])(Subscript[X, e2])(Subscript[Z, e2])+(Subscript[W, e3])(Subscript[X, e3])(Subscript[Z, e3]))-3((Subscript[W, u1])(Subscript[X, u1])(Subscript[Z, u1])+(Subscript[W, u2])(Subscript[X, u2])(Subscript[Z, u2])+(Subscript[W, u3])(Subscript[X, u3])(Subscript[Z, u3]))-3((Subscript[W, d1])(Subscript[X, d1])(Subscript[Z, d1])+(Subscript[W, d2])(Subscript[X, d2])(Subscript[Z, d2])+(Subscript[W, d3])(Subscript[X, d3])(Subscript[Z, d3]))+2((Subscript[W, L1])(Subscript[X, L1])(Subscript[Z, L1])+(Subscript[W, L2])(Subscript[X, L2])(Subscript[Z, L2])+(Subscript[W, L3])(Subscript[X, L3])(Subscript[Z, L3]))+2(Subscript[W, Hd])(Subscript[X, Hd])(Subscript[Z, Hd])+2(Subscript[W, Hu])(Subscript[X, Hu])(Subscript[Z, Hu]);

The accs

accs[X_]:={2((Subscript[X, Q1])+(Subscript[X, Q2])+(Subscript[X, Q3]))-((Subscript[X, u1])+(Subscript[X, u2])+(Subscript[X, u3]))-((Subscript[X, d1])+(Subscript[X, d2])+(Subscript[X, d3])),3((Subscript[X, Q1])+(Subscript[X, Q2])+(Subscript[X, Q3]))+((Subscript[X, L1])+(Subscript[X, L2])+(Subscript[X, L3]))+(Subscript[X, Hd])+(Subscript[X, Hu]),6((Subscript[X, Q1])+(Subscript[X, Q2])+(Subscript[X, Q3]))-((Subscript[X, ν1])+(Subscript[X, ν2])+(Subscript[X, ν3]))-((Subscript[X, e1])+(Subscript[X, e2])+(Subscript[X, e3]))-3((Subscript[X, u1])+(Subscript[X, u2])+(Subscript[X, u3]))-3((Subscript[X, d1])+(Subscript[X, d2])+(Subscript[X, d3]))+2((Subscript[X, L1])+(Subscript[X, L2])+(Subscript[X, L3]))+2(Subscript[X, Hd])+2(Subscript[X, Hu]),((Subscript[X, Q1])+(Subscript[X, Q2])+(Subscript[X, Q3]))-6((Subscript[X, e1])+(Subscript[X, e2])+(Subscript[X, e3]))-8((Subscript[X, u1])+(Subscript[X, u2])+(Subscript[X, u3]))-2((Subscript[X, d1])+(Subscript[X, d2])+(Subscript[X, d3]))+3((Subscript[X, L1])+(Subscript[X, L2])+(Subscript[X, L3]))+3(Subscript[X, Hd])+3(Subscript[X, Hu]),q[X,X],c[X,X,X]};

δ[x_,y_]:= If[x==y,1,0];

soln[parameters_]:=Module[{R,a,b,tc,P},
{(Subscript[R, Q1]),(Subscript[R, ν1]),(Subscript[R, e1]),(Subscript[R, e2]),(Subscript[R, u1]),(Subscript[R, d1]),(Subscript[R, d2]),(Subscript[R, d3]),(Subscript[R, L1]),(Subscript[R, L2]),(Subscript[R, Hd]),(Subscript[R, Hu]),(Subscript[R, 1]),(Subscript[R, 2]),(Subscript[a, 1]),(Subscript[a, 2]),(Subscript[a, 3]),(Subscript[b, 1]),(Subscript[b, 2]),(Subscript[b, 3]),(Subscript[tc, 1]),(Subscript[tc, 2]),(Subscript[tc, 3])}=parameters;

(Subscript[R, Q2])=3((Subscript[R, 1])+(Subscript[R, Hd])+(Subscript[R, L1])+(Subscript[R, L2])+(Subscript[R, Q1]))+4(Subscript[R, 2])+2(Subscript[R, d3])+(Subscript[R, e1])+(Subscript[R, e2]);
(Subscript[R, Q3])=−4((Subscript[R, 1])+(Subscript[R, 2])+(Subscript[R, Q1]))−2((Subscript[R, d3])+(Subscript[R, e1])+(Subscript[R, e2]))−3((Subscript[R, Hd])+(Subscript[R, L1])+(Subscript[R, L2]));
(Subscript[R, ν2])=4(Subscript[R, 1])+(Subscript[R, 2])+3((Subscript[R, e1])+(Subscript[R, e2]))−(Subscript[R, ν1])+(Subscript[R, Q1]);
(Subscript[R, ν3])=−(2(Subscript[R, 1])+(Subscript[R, 2])+(Subscript[R, d1])+(Subscript[R, d2])+(Subscript[R, d3])+(Subscript[R, e1])+(Subscript[R, e2])+(Subscript[R, Q1]));
(Subscript[R, e3])=((Subscript[R, d1])+(Subscript[R, d2])+(Subscript[R, d3]))+3((Subscript[R, e1])+(Subscript[R, e2]))+4(Subscript[R, 1]);
(Subscript[R, u2])=−(6(Subscript[R, 1])+(Subscript[R, 2])+(Subscript[R, d1])+(Subscript[R, d2])+2(Subscript[R, d3])+4((Subscript[R, e1])+(Subscript[R, e2]))+(Subscript[R, Q1])+(Subscript[R, u1]));
(Subscript[R, u3])=4(Subscript[R, 1])+(Subscript[R, 2])+(Subscript[R, d3])+2((Subscript[R, e1])+(Subscript[R, e2]))+(Subscript[R, Q1]);
(Subscript[R, L3])=3((Subscript[R, 1])+(Subscript[R, e1])+(Subscript[R, e2]))−((Subscript[R, L1])+(Subscript[R, L2]))−((Subscript[R, Hu])+(Subscript[R, Hd]));


Γ=(3q[R,R]c[B3,R,R]-2q[B3,R]c[R,R,R])+δ[q[Y3,R],0]δ[q[B3,R],0]δ[q[R,R],0]((Subscript[a, 2])c[R,R,R]-3(Subscript[a, 3])c[B3,R,R]+δ[c[Y3,R,R],0]δ[c[B3,R,R],0]δ[q[R,R,R],0] (Subscript[b, 1]))+δ[2q[B3,R]c[R,R,R],3q[R,R]c[B3,R,R]]δ[2q[Y3,R]c[R,R,R],3q[R,R]c[Y3,R,R]]δ[q[B3,R]c[Y3,R,R],q[Y3,R]c[B3,R,R]]((Subscript[tc, 2])q[R,R]-2(Subscript[tc, 3])q[B3,R]);
Σ=(2q[Y3,R]c[R,R,R]-3q[R,R]c[Y3,R,R])+δ[q[Y3,R],0]δ[q[B3,R],0]δ[q[R,R],0](3(Subscript[a, 3])c[Y3,R,R]-(Subscript[a, 1])c[R,R,R]+δ[c[Y3,R,R],0]δ[c[B3,R,R],0]δ[c[R,R,R],0] (Subscript[b, 2]))+δ[2q[B3,R]c[R,R,R],3q[R,R]c[B3,R,R]]δ[2q[Y3,R]c[R,R,R],3q[R,R]c[Y3,R,R]]δ[q[B3,R]c[Y3,R,R],q[Y3,R]c[B3,R,R]](2(Subscript[tc, 3])q[Y3,R]-(Subscript[tc, 1])q[R,R]);
Λ=6(q[B3,R]c[Y3,R,R]-q[Y3,R]c[B3,R,R])+δ[q[Y3,R],0]δ[q[B3,R],0]δ[q[R,R],0](3((Subscript[a, 1])c[B3,R,R]-(Subscript[a, 2])c[Y3,R,R])+δ[c[Y3,R,R],0]δ[c[B3,R,R],0]δ[c[R,R,R],0] (Subscript[b, 3]))+δ[2q[B3,R]c[R,R,R],3q[R,R]c[B3,R,R]]δ[2q[Y3,R]c[R,R,R],3q[R,R]c[Y3,R,R]]δ[q[B3,R]c[Y3,R,R],q[Y3,R]c[B3,R,R]](2(Subscript[tc, 1])q[B3,R]-2(Subscript[tc, 2])q[Y3,R]);

(Subscript[P, Q1])=Γ+Σ+Λ (Subscript[R, Q1]);
(Subscript[P, Q2])=Γ+Σ+Λ (Subscript[R, Q2]);
(Subscript[P, Q3])=-Γ-Σ+Λ (Subscript[R, Q3]);
(Subscript[P, ν1])=-3Σ+Λ (Subscript[R, ν1]);
(Subscript[P, ν2])=-3Σ+Λ (Subscript[R, ν2]);
(Subscript[P, ν3])=3 Σ+Λ (Subscript[R, ν3]);
(Subscript[P, e1])=-6Γ-3Σ+Λ (Subscript[R, e1]);
(Subscript[P, e2])=-6Γ-3Σ+Λ (Subscript[R, e2]);
(Subscript[P, e3])=6Γ+3Σ+Λ (Subscript[R, e3]);
(Subscript[P, u1])=4Γ+Σ+Λ (Subscript[R, u1]);
(Subscript[P, u2])=4Γ+Σ+Λ (Subscript[R, u2]);
(Subscript[P, u3])=-4Γ-Σ+Λ (Subscript[R, u3]);
(Subscript[P, d1])=-2Γ+Σ+Λ (Subscript[R, d1]);
(Subscript[P, d2])=-2Γ+Σ+Λ (Subscript[R, d2]);
(Subscript[P, d3])=2Γ-Σ+Λ (Subscript[R, d3]);
(Subscript[P, L1])=-3Γ-3Σ+Λ (Subscript[R, L1]);
(Subscript[P, L2])=-3Γ-3Σ+Λ (Subscript[R, L2]);
(Subscript[P, L3])=3Γ+3Σ+Λ (Subscript[R, L3]);
(Subscript[P, Hd])=-3Γ-3 Σ+Λ (Subscript[R, Hd]);
(Subscript[P, Hu])=3Γ+3 Σ+Λ (Subscript[R, Hu]);



P];

Inverse and check

For the inverse we need the dot product

dot[X_,Z_]:=((Subscript[X, Q1])(Subscript[Z, Q1])+(Subscript[X, Q2])(Subscript[Z, Q2])+(Subscript[X, Q3])(Subscript[Z, Q3]))+((Subscript[X, ν1])(Subscript[Z, ν1])+(Subscript[X, ν2])(Subscript[Z, ν2])+(Subscript[X, ν3])(Subscript[Z, ν3]))+((Subscript[X, e1])(Subscript[Z, e1])+(Subscript[X, e2])(Subscript[Z, e2])+(Subscript[X, e3])(Subscript[Z, e3]))+((Subscript[X, u1])(Subscript[Z, u1])+(Subscript[X, u2])(Subscript[Z, u2])+(Subscript[X, u3])(Subscript[Z, u3]))+((Subscript[X, d1])(Subscript[Z, d1])+(Subscript[X, d2])(Subscript[Z, d2])+(Subscript[X, d3])(Subscript[Z, d3]))+((Subscript[X, L1])(Subscript[Z, L1])+(Subscript[X, L2])(Subscript[Z, L2])+(Subscript[X, L3])(Subscript[Z, L3]))+(Subscript[X, Hd])(Subscript[Z, Hd])+(Subscript[X, Hu])(Subscript[Z, Hu]);

inverse[T_]:=Module[{parameters,G,R,a,b,tc},
(Subscript[G, x_]):=108 (Subscript[T, x])-(dot[Y3,T]-dot[B3,T])(Subscript[Y3, x])-(2dot[B3,T]-dot[Y3,T])(Subscript[B3, x]);
(Subscript[R, Q1])=(Subscript[G, Q1]);
(Subscript[R, ν1])=(Subscript[G, ν1]);
(Subscript[R, e1])=(Subscript[G, e1]);
(Subscript[R, e2])=(Subscript[G, e2]);
(Subscript[R, u1])=(Subscript[G, u1]);
(Subscript[R, d1])=(Subscript[G, d1]);
(Subscript[R, d2])=(Subscript[G, d2]);
(Subscript[R, d3])=(Subscript[G, d3]);
(Subscript[R, L1])=(Subscript[G, L1]);
(Subscript[R, L2])=(Subscript[G, L2]);
(Subscript[R, Hd])=(Subscript[G, Hd]);
(Subscript[R, Hu])=(Subscript[G, Hu]);
(Subscript[R, 1])=-((Subscript[G, d1])+(Subscript[G, d2])+(Subscript[G, d3]))+(Subscript[G, e3])-(Subscript[G, Hu])-(Subscript[G, Hd])-((Subscript[G, L1])+(Subscript[G, L2])+(Subscript[G, L3]));
(Subscript[R, 2])=((Subscript[G, d1])+(Subscript[G, d2])+(Subscript[G, d3]))-(Subscript[G, e1])-(Subscript[G, e2])-2(Subscript[G, e3])+2(Subscript[G, Hu])+2(Subscript[G, Hd])+2((Subscript[G, L1])+(Subscript[G, L2])+(Subscript[G, L3]))-(Subscript[G, Q1])-(Subscript[G, ν3]);
(Subscript[a, 1])=c[B3,T,T];
(Subscript[a, 2])=-c[Y3,T,T];
(Subscript[a, 3])=-c[B3,T,T](dot[Y3,T]-dot[B3,T])+c[Y3,T,T](2dot[B3,T]-dot[Y3,T]);
(Subscript[b, 1])=(dot[Y3,T]-dot[B3,T]);
(Subscript[b, 2])=(2dot[B3,T]-dot[Y3,T]);
(Subscript[b, 3])=1;
(Subscript[tc, 1])=q[B3,T];
(Subscript[tc, 2])=-q[Y3,T];
(Subscript[tc, 3])=-q[B3,T](dot[Y3,T]-dot[B3,T])+q[Y3,T](2dot[B3,T]-dot[Y3,T]);

{(Subscript[R, Q1]),(Subscript[R, ν1]),(Subscript[R, e1]),(Subscript[R, e2]),(Subscript[R, u1]),(Subscript[R, d1]),(Subscript[R, d2]),(Subscript[R, d3]),(Subscript[R, L1]),(Subscript[R, L2]),(Subscript[R, Hd]),(Subscript[R, Hu]),(Subscript[R, 1]),(Subscript[R, 2]),(Subscript[a, 1]),(Subscript[a, 2]),(Subscript[a, 3]),(Subscript[b, 1]),(Subscript[b, 2]),(Subscript[b, 3]),(Subscript[tc, 1]),(Subscript[tc, 2]),(Subscript[tc, 3])}
]

Multiplicative factor:

m[T_]:=6*(108^4)(q[B3,T]c[Y3,T,T]-q[Y3,T]c[B3,T,T])+δ[q[Y3,T],0]δ[q[B3,T],0](3*(108^3)(((c[B3,T,T])^2)+((c[Y3,T,T])^2))+108 δ[c[Y3,T,T],0]δ[c[B3,T,T],0])+2*(108^2)δ[q[B3,T]c[Y3,T,T],q[Y3,T]c[B3,T,T]](((q[B3,T])^2)+((q[Y3,T])^2))

Checking solution on large range of inputs.

We randomly generate solutions, and check whether they satisfy the ACCs

accscheck[X_]:=(accs[soln[X]]=={0,0,0,0,0,0})

inversecheck[X_]:=Module[{s}, s=soln[X];
l[soln[inverse[s]]]/m[s]-l[s]=={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}]

Monitor[Do[
p1=RandomInteger[{-10,10},23];
If[!accscheck[p1],Print["Fail"]];

p2=RandomInteger[{-10,10},23];
If[!accscheck[p2],Print["Fail"]];
,{i,1,100000}],i]