On random numbers extracted from NormalDistribution by two different methods under the same seed

I think this is due to the Box–Muller transform that is employed to generate pairs of normally distributed values. This method is much faster than using the inverse commulative probability function.

The algorithm is essentially as follows:

n = 1000000;
a = Flatten[ReIm[
   Times[
    Sqrt[-2. Log[RandomReal[{0, 1}, n]]],
    Exp[(2. Pi I) RandomReal[{0, 1}, n]]
    ]
   ]]; // AbsoluteTiming //First

0.003597

As comparison:

b = RandomVariate[NormalDistribution[], 2 n]; // AbsoluteTiming // First
Histogram[{a, b}]

0.003579

Although I cannot reproduce the precise order used in Mathematica, the timings are very similar. I guess, this supports my hypothesis.

Whenever only a single sampling of the normal distribution is required, the second sampled value is just thrown away. This explains why the results of the bulk sampling equal the results one-by-one sampling, interleaved with other pseudorandom numbers.

When I change the distribution from NormalDistribution to UniformDistribution there is no difference in the two results.

SeedRandom[1234];
RandomVariate[UniformDistribution[], 10]

{0.876608, 0.521964, 0.0862234, 0.377913, 0.0116446, 0.927266, 0.543757, 0.479332, 0.245349, 0.759896}

SeedRandom[1234];
Table[First @ RandomVariate[UniformDistribution[], 1], 10]

{0.876608, 0.521964, 0.0862234, 0.377913, 0.0116446, 0.927266, 0.543757, 0.479332, 0.245349, 0.759896}

From this I conclude that the anomalous behavior is caused by the distribution function and not by RandomVariate. This is additional support for Henrik Schumacher's hypothesis.