Matt Cao
Matt Cao

Reputation: 63

numpy large integer failed

I recently work on some project Euler problems

Smallest multiple

Problem 5

2520 is the smallest number that can be divided by each of the numbers from 1 to 10 without any remainder.

What is the smallest positive number that is evenly divisible by all of the numbers from 1 to 20?

I wrote my code it works great

def factor_finder(n, j=2):

    factor_list = []

    if n == 2:
        return [2]
    elif n == 3:
        return [3]
    else:
        while n >= j * 2:
            while n % j == 0:
                n = int(n / j)
                factor_list.append(j)
            j += 1

    if n > 1:
        factor_list.append(n)

    return factor_list



def smallest_multiples(n):

    from functools import reduce

    factor_list = []
    final_list = []

    for i in range(2, n + 1):
        factor_list += factor_finder(i)
    # print(factor_list)

    for i in set(factor_list):
        l1 = []
        l2 = []
        for j in factor_list:
            if j == i:
                l1.append(j)
            else:
                if len(l1) > len(l2):
                    l2 = l1
                    l1 = []
                else:
                    l1 = []
        # print(l2)
        final_list += l2
    # print(final_list)

    return (
        np.array(final_list).cumprod()[-1],
        reduce((lambda x, y: x * y), final_list),
    )

The result is:

%time

smallest_multiples(1000)

CPU times: user 5 µs, sys: 0 ns, total: 5 µs Wall time: 32.4 µs

(-4008056434385126912, 7128865274665093053166384155714272920668358861885893040452001991154324087581111499476444151913871586911717817019575256512980264067621009251465871004305131072686268143200196609974862745937188343705015434452523739745298963145674982128236956232823794011068809262317708861979540791247754558049326475737829923352751796735248042463638051137034331214781746850878453485678021888075373249921995672056932029099390891687487672697950931603520000)

My question is why numpy.cumprod() failed to get the right number. I thought numpy is the very number tool. Can Somebody give me some idea?

Upvotes: 0

Views: 335

Answers (2)

rLevv
rLevv

Reputation: 508

The problem is that the number reached a size that meant that it was no longer representable by ints in Python. If you look here, you'll see that ints max out in size around 19 digits (i.e. 2^63 from 63 bits + sign bit) and then go into overflow. Numpy is based in C which uses fixed-precision for much faster computations with the trade off that it is limited by the 64bit integer and will overflow. Some functions in numpy even guard against this by converting to floats to do calculations which can hold even more digits.

If you tell numpy to use "object" as your datatype, there is a significant time penalty but it'll let you use the arbitrary-precision that you're used to in Python. For your code, it would look like:

return (
    np.cumprod(final_list, dtype="object")[-1],
    reduce((lambda x, y: x * y), final_list),
)

More about overflow in numpy.

Upvotes: 1

Paul Panzer
Paul Panzer

Reputation: 53029

Numerical analysis is not number theory. Correctness is not the only goal but has to be weighed against efficiency. Arbitrary precision numbers (like large integers) are slow, so numpy defaults to using fixed length integers. These just overflow when they become too large. You can instruct numpy to use arbitrary precision integers, but you will lose much of its speed:

np.arange(1, 100).prod() # fast but wrong
# 0
np.arange(1, 100, dtype=object).prod() # slow but correct
# 933262154439441526816992388562667004907159682643816214685929638952175999932299156089414639761565182862536979208272237582511852109168640000000000000000000000

Upvotes: 1

Related Questions