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.

Since, the error message explicitly says: "Integers to negative integer powers are not allowed.", why not try with: `a_range = np.arange(2, n+1, dtype=np.float)`? — , Apr 07 '17 at 21:30
Here is a relevant link: https://mail.scipy.org/pipermail/numpy-discussion/2016-October/076087.html — juanpa.arrivillaga, Apr 07 '17 at 21:33
@Evert -- `numpy.arange` has some issues with floating point values though that can arise from floating point precision errors etc (IIRC). — mgilson, Apr 07 '17 at 21:34
Obvious thing, thanks! I'll now continue trying to prove it :-) — shredding, Apr 07 '17 at 21:34
@Evert, Why not add your comment as an answer? I think that it will be accepted. — GurstTavo, Apr 07 '17 at 21:44

user2357112 · Accepted Answer · 2018-05-15T22:11:50.960

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.])

"_the logic used to choose the output dtype of an operation_ […] _is based on dtypes, not values_". This makes much sense to me, but why then is an integer divided by another integer converted to a float in Python 3? — pela, Nov 06 '19 at 09:29
@pela: Ordinary int/int division doesn't involve any NumPy logic or dtypes at all, so none of this is relevant. That said, dividing NumPy integers by NumPy integers on Python 3 produces NumPy floats, but that's no problem. The output dtype chosen for division with integer input dtypes is a floating point dtype, still regardless of input values. There's nothing forcing the output dtype to match any of the input dtypes. — user2357112, Nov 06 '19 at 09:55
Ah, sorry, you were referring just to numpy, not to Python. That said, I tend to disagree on "that's no problem". Ensuring that 3/2 = 1 and –3/2 = –1 (as in Fortran, IDL, and most other languages I think) is quite cumbersome in Python. Anyway, thanks for the reply :) — pela, Nov 06 '19 at 10:19
@pela: It sounds like you may be unaware of `//`. `//` doesn't *quite* behave like other languages' integer division operators, since it floors instead of rounding toward zero, but this tends to be more useful behavior, especially combined with `%` keeping the sign of the right operand instead of the left. (Maybe you were aware of `//`, but I would have expected your complaint to be worded differently if you were aware.) — user2357112, Nov 06 '19 at 10:49
Thanks, I _am_ aware of `//`, but as you say it floors, which may not be desired for negative numbers. But you're probably right that in general it's more useful behavior. — pela, Nov 06 '19 at 12:17

score 2 · Answer 2 · answered Nov 25 '18 at 19:57

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.

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

2 Answers2

Linked

Related