(* Independent Wolfram Language handle for 031 — The Shape Belongs to the Rule.
   The constants come from Wolfram's Pi and E, not from the page's spigot and
   integer-series implementations.  The finite claim is about 900 displayed
   digits only; this is not a normality test. *)

ClearAll[n, seed, piDigits, eDigits, u32, mulberryDigits, numericDirections, yy, theta,
  witness, results, expectedCounts, expectedRatios, openingsHold,
  unitDirections, recordedValuesHold, finiteClaim];

n = 900;
seed = 20260828;

piDigits = First@RealDigits[Pi, 10, n];
eDigits = First@RealDigits[E, 10, n];

u32[x_] := Mod[x, 2^32];
mulberryDigits[count_Integer, initialSeed_Integer] := Module[
  {state = u32[initialSeed], value, mixed},
  Table[
    state = u32[state + 16^^6D2B79F5];
    value = state;
    value = u32[BitXor[value, BitShiftRight[value, 15]]*
      BitOr[value, 1]];
    mixed = u32[BitXor[value, BitShiftRight[value, 7]]*
      BitOr[value, 61]];
    value = u32[BitXor[value, u32[value + mixed]]];
    Floor[10 u32[BitXor[value, BitShiftRight[value, 14]]]/2^32],
    {count}
  ]
];

numericDirections = Table[
  With[
    {y = 1.0 - 2.0 i/9.0,
     theta = N[i Pi (3 - Sqrt[5]), MachinePrecision]},
    {Cos[theta] Sqrt[Max[0.0, 1.0 - y^2]], y,
     Sin[theta] Sqrt[Max[0.0, 1.0 - y^2]]}
  ],
  {i, 0, 9}
];

witness[digits_List] := Module[
  {counts, position, displacement, ratio, chiSquared},
  counts = BinCounts[digits, {-1/2, 19/2, 1}];
  position = Total[numericDirections[[# + 1]] & /@ digits];
  displacement = Norm[position];
  ratio = displacement/Sqrt[n];
  chiSquared = Total[(counts - n/10)^2/(n/10)];
  <|
    "counts" -> counts,
    "displacement over sqrt(n)" -> N[ratio, 16],
    "chi squared" -> N[chiSquared, 16]
  |>
];

openingsHold =
  FromDigits[Take[piDigits, 32]] === 31415926535897932384626433832795 &&
  FromDigits[Take[eDigits, 32]] === 27182818284590452353602874713526;
unitDirections = FullSimplify[
  (Cos[theta] Sqrt[1 - yy^2])^2 + yy^2 +
      (Sin[theta] Sqrt[1 - yy^2])^2 == 1,
  Assumptions -> Element[{yy, theta}, Reals] && -1 <= yy <= 1
];

results = <|
  "pi" -> witness[piDigits],
  "e" -> witness[eDigits],
  "control" -> witness[mulberryDigits[n, seed]]
|>;

expectedCounts = <|
  "pi" -> {86, 103, 93, 94, 89, 84, 85, 83, 88, 95},
  "e" -> {88, 83, 93, 99, 91, 72, 88, 87, 94, 105},
  "control" -> {75, 87, 86, 103, 85, 102, 94, 85, 101, 82}
|>;
expectedRatios = <|
  "pi" -> 0.577712621887,
  "e" -> 0.733539871626,
  "control" -> 1.326808687770
|>;

recordedValuesHold = And @@ Table[
  results[name]["counts"] === expectedCounts[name] &&
    Abs[results[name]["displacement over sqrt(n)"] - expectedRatios[name]] < 5 10^-12,
  {name, Keys[results]}
];

finiteClaim = And @@ Table[
  With[
    {ratio = results[name]["displacement over sqrt(n)"],
     chi = results[name]["chi squared"]},
    1/4 < ratio < 5/2 && chi < 169/10
  ],
  {name, Keys[results]}
];

Print[<|
  "N" -> n,
  "seed" -> seed,
  "known openings hold" -> openingsHold,
  "ten directions are unit vectors" -> unitDirections,
  "witnesses" -> results
|>];
Print["claim holds: ",
  TrueQ[openingsHold && unitDirections && recordedValuesHold && finiteClaim]
];
