Why can't I raise to a negative power in numpy?

Question

I'm modelling the Riemann theta function:

import numpy as np
def theta(s, n=100):
    a_range = np.arange(2, n + 1)
    return 1 + sum(1/(a_range ** s))

It does not work for negative s; e.g. theta(-2) leads to this error:

      1 def theta(s, n=100):
      2     a_range = np.arange(1)
----> 3     return 1 + sum(1/(a_range ** s))
      4 
      5 theta(-2)

      ValueError: Integers to negative integer powers are not allowed.

Why's that? x^-1 should just be 1/x if i recall my math correctly.

Answer 1

In NumPy, the logic used to choose the output dtype of an operation like a_range ** s is based on dtypes, not values. That means that a_range ** -2 has to have the same output dtype as a_range ** 2.

It's important that something like numpy.array([2]) ** 2 give integer output, so that means that numpy.array([2]) ** -2 has to give integers or nothing. They picked nothing; raising integers to negative integer powers is an error in NumPy.

If you want floating-point output, make floating-point input:

a_range = np.arange(2, n + 1, dtype=float)

or

a_range = np.arange(2, n + 1).astype(float)

---

There are a few weird aspects of NumPy's type rules you might not expect from the above description. One is that for operations involving both scalars and arrays, the scalar's dtype may actually be "demoted" based on its value before the input dtypes are used to choose the result dtype:

>>> (numpy.array([1], dtype='int8') + numpy.int32(1)).dtype
dtype('int8')
>>> (numpy.array([1], dtype='int8') + numpy.array([1], dtype='int32')).dtype
dtype('int32')

Here, the scalar numpy.int32(1) gets "demoted" to int8, but the arrays don't get demoted. (It's actually a bit more complex than just demoting to int8, particularly for signed/unsigned handling; see the implementation for the full details.)

Second, when uint64s are involved, NumPy may suddenly appear to be okay with negative exponents:

>>> numpy.arange(5, dtype='uint64') ** -2
__main__:1: RuntimeWarning: divide by zero encountered in power
array([       inf, 1.        , 0.25      , 0.11111111, 0.0625    ])

This is because NumPy cannot find an integer dtype big enough for both uint64 values and negative values, so it gives up and coerces the inputs to floats. The same can be seen with a positive exponent of signed dtype, as long as you avoid scalar type "demotion":

>>> numpy.arange(5, dtype='uint64') ** numpy.array([2], dtype='int32')
array([ 0.,  1.,  4.,  9., 16.])

Answer 2

What I did to solve the problem is -

float_value = float(input_value)
power_value = float_value**-power
final_value = int(power_value)

Here, input_value is the input number, it is converted into float, then I calculate the negative power and then take the int value of it.