Project Euler #101 - how to work around numpy polynomial overflow?

Question

Project Euler #101

I just started learning Numpy and it so far looks pretty straightforward to me.

One thing I ran into is that when I evaluate the polynomial, the result is a int32, so an overflow would occur.

u = numpy.poly1d([1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1])
for i in xrange(1, 11):
    print(i, u(i))

The results are:

(1, 1)
(2, 683)
(3, 44287)
(4, 838861)
(5, 8138021)
(6, 51828151)
(7, 247165843)
(8, 954437177)
(9, -1156861335)
(10, 500974499)

The last two items are clearly incorrect.

The work around I can think of is factoring the coefficients by 100

u = numpy.poly1d([0.01, -0.01, 0.01, -0.01, 0.01, -0.01, 0.01, -0.01, 0.01, -0.01, 0.01])
for i in xrange(1, 11):
    print(i, int(u(i) * 100))

This time the results are correct

(1, 1)
(2, 682)
(3, 44286)
(4, 838860)
(5, 8138020)
(6, 51828151)
(7, 247165843)
(8, 954437177)
(9, 3138105961L)
(10, 9090909091L)

Is there a better way? Does Numpy allow me to change the data type? Thanks.

Answer 1

It wasn't the scaling by 100 that helped, but the fact that the numbers given were floats instead of ints, and thus had a higher range. Due to the floating-point calculations, there are some inaccuracies introduced to the calculations as you have seen.

You can specify the type manually like this:

u = numpy.poly1d(numpy.array([1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1], dtype=numpy.int64))

The calculations for this problem fit in 64-bit ints, so this will work.

The supported types are listed here.

Answer 2

Interjay posted a better answer while I was writing this up, but I figured you might want an alternative anyway.

Here's a simple implementation for the examples you showed:

class Poly(object):
    def __init__(self, coefficients):
        assert len(coefficients) > 0
        self.coefficients = coefficients
    def __call__(self, value):
        total = self.coefficients[0]
        for c in self.coefficients[1:]:
            total = total * value + c
        return total

along with some tests

assert Poly([5])(1) == 5
assert Poly([7])(1) == 7
assert Poly([2,3])(5) == 13
assert Poly([1,0,0,0,0])(-2) == 16
u = Poly([1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1])
for i in range(1, 11):
    print (i, u(i))

and the rather useless

assert Poly([2,"!"])("Hello ") == "Hello Hello !"