The point of the "u = 1.0 - self.random()" line was to prevent the case where *u* was exactly equal to zero; however, as you noted, raising a very small number to a small can round down to zero.

We could wrap the calculation in a try/except to catch the ZeroDivisionError but that is probably insufficient.  The codomain can easily spew beyond the range of a C double:

   >>> alpha = 0.01
   >>> u = 0.005
   >>> 1.0 / u ** (1.0 / alpha)
   1.2676506002282268e+230             <== Huge!

We could clamp the values as is done in gammavariate().  Or we can raise a clean Overflow or Underflow exception but that would be unpleasant for users if it were to arise unexpectedly.