CODEWITHSUNDEEP
pythonnumpymath
Ben Quigley
Ben Quigley

Reputation: 727

Why is numpy.prod() incorrectly returning negative results, or 0, for my long lists of natural numbers?

I'm just working on Project Euler problem 12, so I need to do some testing against numbers that are multiples of over 500 unique factors.

I figured that the array [1, 2, 3... 500] would be a good starting point, since the product of that array is the lowest possible such number. However, numpy.prod() returns zero for this array. I'm sure I'm missing something obvious, but what the hell is it?

>>> import numpy as np
>>> array = []
>>> for i in range(1,100):
...   array.append(i)
... 
>>> np.prod(array)
0
>>> array.append(501)
>>> np.prod(array)
0
>>> array.append(5320934)
>>> np.prod(array)
0

Upvotes: 13

Views: 9078

Answers (3)

Lutz Lehmann
Lutz Lehmann

Reputation: 25992

You get the result 0 due to the large number of factors 2 in the product, there are more than 450 of those factors. Thus in a reduction modulo 2^64, the result is zero.

Why the data type forces this reduction is explained in the other answers.

250+125+62+31+15+7+3+1 = 494 is the multiplicity of 2 in 500!

added 12/2020: or, in closer reading the question and its code,

49+24+12+6+3+1 = 95 as the multiplicity of 2 in 99!

which is the product of the first part of your list. Still enough binary zeros at the end of the number to fill all the bit positions of a 64bit integer. Just to compare, you get

19+3 = 22 factors of 5 in 99!

which is also the number of trailing zeros in the decimal expression of this factorial.

Upvotes: 1

Ami Tavory
Ami Tavory

Reputation: 76297

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:

  1. Drop back to Python, and multiply using its "unlimited integers" (see this question for how to do so).

  2. 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.

Upvotes: 9

wim
wim

Reputation: 362786

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

Upvotes: 6

Related Questions