I checked his comment, I think he's incorrect on the approaches being different. Using my function `foo` that does the iterative approach, we can compare the distributions and they're fairly identical.

  >>> np.mean([foo(100) for _ in range(100000)])
  3.258425093913613e-33

  >>> np.mean([np.product(np.random.random(100)) for _ in range(100000)])
  8.814732867008917e-33

(There's quite a bit of variance on a log scale, so 3 vs 8 is not a huge difference. I re-ran this and got varying numbers with both approaches. But the iterated code is very slow...)

Note that the mean is actually quite a bit closer to 2^-100 even though the vast majority of numbers from either distribution fall below 2^-100. Even so, the mean for both is approximately a factor of 100 less than 2^-100. Suspicious! Although I think we've both burned enough time on this.