Wrong calculation of negative powers of variables (Sorry for not describing the error in detail).

[itumekanik]

from numexpr import evaluate

x = 2

print(evaluate("x(2.0)")) print(evaluate("x(-2.0)")) print(x**(-2))

OUTPUT:

4 0 # WRONG!!!!! 0.25

[bashtage]

Is integer truncation playing a role here?

[gdementen]

Note that this is only an issue when x is an integer:

>>> from numexpr import evaluate
>>> x = 2
>>> print(evaluate("x**(2.0)"))
4
>>> print(evaluate("x**(-2.0)"))
0
>>> x = 2.0
>>> print(evaluate("x**(2.0)"))
4.0
>>> print(evaluate("x**(-2.0)"))
0.25

[gdementen]

@bashtage : you beat me by a few seconds 😉

[itumekanik]

Hi, Then why not make "x" float automatically?

[gdementen]

Then why not make "x" float automatically?

My guess would be that nobody noticed it (AFAIK numexpr is rarely used to do scalar computations). But I just noticed the problem with negative powers is larger than I thought:

>>> import numexpr as ne
>>> a = np.array([1, 2, 3])
>>> ne.evaluate("a ** -2")            # <-- this is correct (mimics numpy)
ValueError: Integers to negative integer powers are not allowed.
>>> ne.evaluate("a ** -2.0")          # <-- this is not correct
array([1, 0, 0], dtype=int32)'

While with numpy, this gives:

>>> a ** -2
ValueError: Integers to negative integer powers are not allowed
>>> a ** -2.0
array([ 1.        ,  0.25      ,  0.11111111])

If that is your case, you may be able to work around the issue by changing your expression (I guess this is what numexpr should do behind the scene -- but I am unsure it would preserve numeric accuracy for high values):

>>> ne.evaluate("1 / (a ** 2.0)")
array([ 1.        ,  0.25      ,  0.11111111])
>>> ne.evaluate("1 / (a ** 2)")
array([ 1.        ,  0.25      ,  0.11111111])

Or by casting to float before calling numexpr:

>>> a = a.astype(float)
>>> ne.evaluate("a ** -2.0")
array([ 1.        ,  0.25      ,  0.11111111])
>>> ne.evaluate("a ** -2")
array([ 1.        ,  0.25      ,  0.11111111])

[robbmcleod]

I'll be back from vacation on the 22nd, if someone wants to remind me of this issue then I can take a deeper look at it to see if it's resolvable or not.

[gdementen]

I think it shouldn't be too hard to solve.

I don't have time to test that, but it seems to me that it would be enough to add a kind argument to the line at: https://github.com/pydata/numexpr/blob/62cc05148fd6eb8bb32ed1bbcaf102149da6445e/numexpr/expressions.py#L315

Or maybe replace that line with:

    r = truediv_op(ConstantNode(1), r)

While looking at the code, I think line 319 is also problematic for integer input (i.e. a ** -1.0 is wrong when a is an integer array) and should either use truediv_op or specify kind.

Here is an additional test case for that:

>>> import numexpr as ne
>>> a = np.array([1, 2, 3])
>>> a ** -1.0
array([ 1.        ,  0.5       ,  0.33333333])
>>> ne.evaluate("a ** -1.0")
array([1, 0, 0], dtype=int32)

[itumekanik]

Hi @gdementen ,

I tested your change proposal and replaced line 315 with the following statement

r = truediv_op(ConstantNode(1), r)

This seems to work for my case.

Still I guess it would require some more testing.

[itumekanik]

Following is a test script of my case for the calculation of two different matrices, one with numexpr the other with numpy. It seems to work fine after the change you mentioned in line 315 of the expressions.py.

import numpy as np

def _lambdifygenerated0(sym_en_3, sym_en_2, sym_en_0, sym_en_1):
    from numexpr import evaluate
    x0 = sym_en_0 - sym_en_1
    x1 = evaluate("x0**2", truediv=True)
    x2 = sym_en_2 - sym_en_3
    x3 = evaluate("x2**2", truediv=True)
    x4 = evaluate("x1 + x3", truediv=True)
    x5 = evaluate("x4**1.5", truediv=True)
    # x6 = evaluate("5040000.0/x5", truediv=True)
    x6 = 5040000.0/x5
    x7 = evaluate("x1*x6", truediv=True)
    # x8 = evaluate("sqrt(x4)", truediv=True)
    x8 = np.sqrt(x4)
    x9 = evaluate("sym_en_0**2 - 2*sym_en_0*sym_en_1 + sym_en_1**2 + sym_en_2**2 - 2*sym_en_2*sym_en_3 + sym_en_3**2", truediv=True)
    x10 = evaluate("x9**(-2.0)", truediv=True)
    x11 = evaluate("x4**1.0", truediv=True)
    x12 = evaluate("x9**(-2.5)", truediv=True)
    x13 = evaluate("5443200.0*x10*x8 - 10886400.0*x11*x12 + 7257600.0*x5*x9**(-3.0)", truediv=True)
    # x13 = 5443200.0*x10*x8 - 10886400.0*x11*x12 + 7257600.0*x5*x9**(-3.0)
    x14 = evaluate("x11**(-1.0)", truediv=True)
    x15 = evaluate("x13*x14*x3 + x7", truediv=True)
    x16 = evaluate("x0*x2", truediv=True)
    x17 = evaluate("x16*x6", truediv=True)
    x18 = evaluate("x14*x16", truediv=True)
    x19 = evaluate("-x13*x18 + x17", truediv=True)
    x20 = evaluate("x9**(-1.5)", truediv=True)
    x21 = evaluate("x20*x8", truediv=True)
    x22 = evaluate("x10*x11", truediv=True)
    x23 = evaluate("3628800.0*x12*x5", truediv=True)
    x24 = evaluate("3628800.0*x21 - 6350400.0*x22 + x23", truediv=True)
    x25 = evaluate("x8**(-1.0)", truediv=True)
    x26 = evaluate("x2*x25", truediv=True)
    x27 = evaluate("-x24*x26", truediv=True)
    x28 = evaluate("-x13", truediv=True)
    x29 = evaluate("x14*x28*x3 - x7", truediv=True)
    x30 = evaluate("-x17 - x18*x28", truediv=True)
    x31 = evaluate("1814400.0*x21 - 4536000.0*x22 + x23", truediv=True)
    x32 = evaluate("-x26*x31", truediv=True)
    x33 = evaluate("x3*x6", truediv=True)
    x34 = evaluate("x1*x14", truediv=True)
    x35 = evaluate("x13*x34 + x33", truediv=True)
    x36 = evaluate("x0*x25", truediv=True)
    x37 = evaluate("x24*x36", truediv=True)
    x38 = evaluate("x28*x34 - x33", truediv=True)
    x39 = evaluate("x31*x36", truediv=True)
    x40 = evaluate("x8*x9**(-1.0)", truediv=True)
    x41 = evaluate("x11*x20", truediv=True)
    x42 = evaluate("1814400.0*x10*x5", truediv=True)
    x43 = evaluate("-x24", truediv=True)
    x44 = evaluate("-x26*x43", truediv=True)
    x45 = evaluate("x36*x43", truediv=True)
    x46 = evaluate("1209600.0*x40 - 2721600.0*x41 + x42", truediv=True)
    x47 = evaluate("-x31", truediv=True)
    x48 = evaluate("-x26*x47", truediv=True)
    x49 = evaluate("x36*x47", truediv=True)
    return np.array([[x15, x19, x27, x29, x30, x32],
                     [x19, x35, x37, x30, x38, x39],
                     [x27, x37, 2419200.0*x40 - 3628800.0*x41 + x42, x44, x45, x46],
                     [x29, x30, x44, x15, x19, x48],
                     [x30, x38, x45, x19, x35, x49],
                     [x32, x39, x46, x48, x49, 604800.0*x40 - 1814400.0*x41 + x42]])

def _lambdifygenerated1(sym_en_0, sym_en_1, sym_en_2, sym_en_3):
    x0 = sym_en_0 - sym_en_1
    x1 = x0**2
    x2 = sym_en_2 - sym_en_3
    x3 = x2**2
    x4 = x1 + x3
    x5 = x4**1.5
    x6 = 5040000.0/x5
    x7 = x1*x6
    x8 = np.sqrt(x4)
    x9 = sym_en_0**2 - 2*sym_en_0*sym_en_1 + sym_en_1**2 + sym_en_2**2 - 2*sym_en_2*sym_en_3 + sym_en_3**2
    x10 = x9**(-2.0)
    x11 = x4**1.0
    x12 = x9**(-2.5)
    x13 = 5443200.0*x10*x8 - 10886400.0*x11*x12 + 7257600.0*x5*x9**(-3.0)
    x14 = x11**(-1.0)
    x15 = x13*x14*x3 + x7
    x16 = x0*x2
    x17 = x16*x6
    x18 = x14*x16
    x19 = -x13*x18 + x17
    x20 = x9**(-1.5)
    x21 = x20*x8
    x22 = x10*x11
    x23 = 3628800.0*x12*x5
    x24 = 3628800.0*x21 - 6350400.0*x22 + x23
    x25 = x8**(-1.0)
    x26 = x2*x25
    x27 = -x24*x26
    x28 = -x13
    x29 = x14*x28*x3 - x7
    x30 = -x17 - x18*x28
    x31 = 1814400.0*x21 - 4536000.0*x22 + x23
    x32 = -x26*x31
    x33 = x3*x6
    x34 = x1*x14
    x35 = x13*x34 + x33
    x36 = x0*x25
    x37 = x24*x36
    x38 = x28*x34 - x33
    x39 = x31*x36
    x40 = x8*x9**(-1.0)
    x41 = x11*x20
    x42 = 1814400.0*x10*x5
    x43 = -x24
    x44 = -x26*x43
    x45 = x36*x43
    x46 = 1209600.0*x40 - 2721600.0*x41 + x42
    x47 = -x31
    x48 = -x26*x47
    x49 = x36*x47
    return np.array([[x15, x19, x27, x29, x30, x32], [x19, x35, x37, x30, x38, x39], [x27, x37, 2419200.0*x40 - 3628800.0*x41 + x42, x44, x45, x46], [x29, x30, x44, x15, x19, x48], [x30, x38, x45, x19, x35, x49], [x32, x39, x46, x48, x49, 604800.0*x40 - 1814400.0*x41 + x42]])

s0, s1, s2, s3 = 1, 2, 3, 4

a = _lambdifygenerated0(sym_en_0=s0, sym_en_1=s1, sym_en_2=s2, sym_en_3=s3)
b = _lambdifygenerated1(sym_en_0=s0, sym_en_1=s1, sym_en_2=s2, sym_en_3=s3)

print(a-b)

from numexpr import evaluate
x = 2
print(evaluate("x**(2.0)"))
print(evaluate("x**(-2.0)"))
print(x**(-2))

[robbmcleod]

Should be fixed, thanks to the contributors.