Why is numpy.prod() incorrectly returning negative results, or 0, for my long lists of natural numbers?
Note that Python uses "unlimited" integers, but in numpy everything is typed, and so it is a "C"-style (probably 64-bit) integer here. You're probably experiencing an overflow.
If you look at the documentation for numpy.prod, you can see the dtype parameter:
The type of the returned array, as well as of the accumulator in which the elements are multiplied.
There are a few things you can do:
Drop back to Python, and multiply using its "unlimited integers" (see this question for how to do so).
Consider whether you actually need to find the product of such huge numbers. Often, when you're working with the product of very small or very large numbers, you switch to sums of logarithms. As @WarrenWeckesser notes, this is obviously imprecise (it's not like taking the exponent at the end will give you the exact solution) - rather, it's used to gauge whether one product is growing faster than another.
Those numbers get very big, fast.
>>> np.prod(array[:25])
7034535277573963776
>>> np.prod(array[:26])
-1569523520172457984
>>> type(_)
numpy.int64
You're actually overflowing numpy's data type here, hence the wack results. If you stick to python ints, you won't have overflow.
>>> import operator
>>> reduce(operator.mul, array, 1)
933262154439441526816992388562667004907159682643816214685929638952175999932299156089414639761565182862536979208272237582511852109168640000000000000000000000L