sp2 = x /. Solve[x^(37) == x + 1, x][[1]];
G = Array[1, 36];
Do[G[[i]] = 1/2 + i95/sp2^i, {i, 1, 36}];
rB = 0;
Do[hw = TimeUsed[];
P = InverseCDF[NormalDistribution[0, 1], N[FractionalPart[G]]];
Y1 = (Orthogonalize[ArrayReshape[Take[P, {1, 16}], {4, 4}]] +
     IdentityMatrix[4]).ArrayReshape[Take[P, {17, 36}], {4, 5}];
z = Partition[Y1.Transpose[Y1], {2, 2}];
If[PositiveDefiniteMatrixQ[
    ArrayFlatten@{{z[[1, 1]], z[[2, 1]]}, {z[[1, 2]], z[[2, 2]]}}] ==
   True, rB = rB + 1];
If[Mod[i95, 500000] == 0,
  Print[{TimeUsed[] - hw, i95, rB, rB/i95, N[rB/i95, 20]}];
{i95,rB} >> LowDiscrepancyBuresTwoRebitsSave], {i95, 1, 1000000000}]