|
Needs["Quaternions`"];
(* Initialize variables *)
R = 20;
NN = 1010;
(* Quaternion operations *)
test[q_Quaternion] :=
Module[{unit, res, a, b, c, u, v, w, p},
If[Round[Norm[q]] > R, Return[]];
If[q == Quaternion[0, 0, 0, 0], Return[]];
unit = Quaternion[0, 1, 0, 0];
res = q ** unit ** Conjugate[q];
a = Abs[res[[2]]] + Abs[res[[3]]] + Abs[res[[4]]];
unit = Quaternion[0, 0, 1, 0];
res = q ** unit ** Conjugate[q];
b = Abs[res[[2]]] + Abs[res[[3]]] + Abs[res[[4]]];
unit = Quaternion[0, 0, 0, 1];
res = q ** unit ** Conjugate[q];
c = Abs[res[[2]]] + Abs[res[[3]]] + Abs[res[[4]]];
For[i = 1, i <= (R - 1)/Max[a, b, c], i++,
If[SquareFreeQ[i], {u = a*i;
v = b*i;
w = c*i;
p = Max[u, v, w] + 1;
coe[[p + 1, 4]] += (1);
coe[[p + 1, 3]] -= (u + v + w);
coe[[p + 1, 2]] += (u*v + v*w + w*u);
coe[[p + 1, 1]] -= (u*v*w)}]]];
(* Set up coefficient matrix *)
coe = ConstantArray[0, {NN, 4}];
(* Loop through quaternions *)
rt = Ceiling[Sqrt[R]] + 1;
For[s = -rt, s <= rt, s++,
For[x = -rt, x <= rt, x++,
For[y = -rt, y <= rt, y++,
For[z = -rt, z <= rt, z++, test[Quaternion[s, x, y, z]];
test[Quaternion[s + 0.5, x + 0.5, y + 0.5, z + 0.5]]; ]]]];
newCoe = coe;
newCoe[[2 ;; ;; 2]] = coe[[2 ;; ;; 2]]/2;
(* Calculate and output results *)
For[i = 2, i <= R + 1, i++, ans = 0;
For[j = 4, j >= 1, j--, newCoe[[i, j]] += newCoe[[i - 1, j]];
ans = ans*(i - 1) + newCoe[[i, j]];
];
Print[i - 1, " ", ans/24]; ];
| 