Possible bug when computing powers of complex numbers

[Florian Hanisch]

Hi,

this might be a duplicate (sorry) but I could not find this issue described anywhere else. I have been running the code

k=1.0
(I*k/(2*pi*5))^(6/2)

and got the result

(1.00000046096651e-4294967297*I)/pi^3

This is clearly way too small, adding CC(...) to get a numerical value does not resolve 
the problem. However, simply replacing 6/2 by 3 yields the correct result: 

(I*k/(2*pi*5))^(3)
result: -0.00100000000000000*I/pi^3

Did I make some stupid and obvious error or is this just a strange bug ? I am using 
sage 8.6 (jupyter notebook) on a linux computer (Centos 7, 64 bit). I tried these code 
snippets on cocalc and ran into the same problems. Friends of mine have also reported the 
problem on sage 8.8. 

Thanks a lot,

Florian

 

[Dima Pasechnik]

yes, it appears to be a bug. here is a slightly easier example:

sage: t
(0.100000000000000*I/pi)^(3/2)
sage: t.n()^2
-6.77626357803440e-21 - 0.0000322515344331995*I
sage: t.n()
# a correct numerical value
-0.00401569012955429 + 0.00401569012955429*I
sage: t^2
# pynac bug ?
(1.00000046096651e-4294967297*I)/pi^3

maxima is doing it right:

sage: maxima(t)
((0.1*%i)/%pi)^(3/2)
sage: maxima(t)^2
-(0.001*%i)/%pi^3




>
> Thanks a lot,
>
> Florian
>
>
>
> --
> You received this message because you are subscribed to the Google Groups "sage-devel" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sage-devel+...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/sage-devel/ba2015a4-6249-4d4b-98ce-5d8e115fd95a%40googlegroups.com.

[Emmanuel Charpentier]

It's worse than that...

sage: t=((1/10)*I/pi)^(3/2)

sage: t

(1/10*I/pi)^(3/2)

sage: t^2

This doesn't return. An I had to insist to interrupt (had to strike C-C 6 times in emacs' sage-shell-mode...)

^C^C^C^C^C^C---------------------------------------------------------------------------

KeyboardInterrupt                         Traceback (most recent call last)

<ipython-input-31-f00da4fd424a> in <module>()

----> 1 t**Integer(2)

/usr/local/sage-P3-2/local/lib/python3.7/site-packages/sage/structure/element.pyx in sage.structure.element.Element.__pow__ (build/cythonized/sage/structure/element.c:14088)()

   2038             return (<Element>left)._pow_(right)

   2039         if BOTH_ARE_ELEMENT(cl):

-> 2040             return coercion_model.bin_op(left, right, pow)

   2041 

   2042         cdef long value

/usr/local/sage-P3-2/local/lib/python3.7/site-packages/sage/structure/coerce.pyx in sage.structure.coerce.CoercionModel.bin_op (build/cythonized/sage/structure/coerce.c:9800)()

   1153                 return (<Action>action)._act_(x, y)

   1154             else:

-> 1155                 return (<Action>action)._act_(y, x)

   1156 

   1157         # Now coerce to a common parent and do the operation there

/usr/local/sage-P3-2/local/lib/python3.7/site-packages/sage/structure/coerce_actions.pyx in sage.structure.coerce_actions.IntegerPowAction._act_ (build/cythonized/sage/structure/coerce_actions.c:10173)()

    861         integer_check_long(n, &value, &err)

    862         if not err:

--> 863             return e._pow_long(value)

    864         return e._pow_int(n)

    865 

/usr/local/sage-P3-2/local/lib/python3.7/site-packages/sage/structure/element.pyx in sage.structure.element.Element._pow_long (build/cythonized/sage/structure/element.c:14635)()

   2115         Generic path for powering with a C long.

   2116         """

-> 2117         return self._pow_int(n)

   2118 

   2119 

/usr/local/sage-P3-2/local/lib/python3.7/site-packages/sage/symbolic/expression.pyx in sage.symbolic.expression.Expression._pow_int (build/cythonized/sage/symbolic/expression.cpp:25386)()

   4049             pi^4

   4050         """

-> 4051         return self._pow_(self._parent(other))

   4052 

   4053     def derivative(self, *args):

/usr/local/sage-P3-2/local/lib/python3.7/site-packages/sage/symbolic/expression.pyx in sage.symbolic.expression.Expression._pow_ (build/cythonized/sage/symbolic/expression.cpp:25211)()

   4033                            relational_operator(self._gobj))

   4034         else:

-> 4035             x = g_pow(self._gobj, nexp._gobj)

   4036         return new_Expression_from_GEx(self._parent, x)

   4037 

/usr/local/sage-P3-2/local/lib/python3.7/site-packages/sage/rings/number_field/number_field_element.pyx in sage.rings.number_field.number_field_element.NumberFieldElement.__pow__ (build/cythonized/sage/rings/number_field/number_field_element.cpp:22076)()

   2316         if (isinstance(base, NumberFieldElement) and

   2317             (isinstance(exp, Integer) or type(exp) is int or exp in ZZ)):

-> 2318             return generic_power(base, exp)

   2319 

   2320         if (isinstance(base, NumberFieldElement) and exp in QQ):

/usr/local/sage-P3-2/local/lib/python3.7/site-packages/sage/arith/power.pyx in sage.arith.power.generic_power (build/cythonized/sage/arith/power.c:2376)()

     81         raise NotImplementedError("non-integral exponents not supported")

     82     if not err:

---> 83         return generic_power_long(a, value)

     84 

     85     if n < 0:

/usr/local/sage-P3-2/local/lib/python3.7/site-packages/sage/arith/power.pyx in sage.arith.power.generic_power_long (build/cythonized/sage/arith/power.c:2663)()

    100         u = -u

    101         a = invert(a)

--> 102     return generic_power_pos(a, u)

    103 

    104 

/usr/local/sage-P3-2/local/lib/python3.7/site-packages/sage/arith/power.pyx in sage.arith.power.generic_power_pos (build/cythonized/sage/arith/power.c:2796)()

    118     n >>= 1

    119     while n:

--> 120         sig_check()

    121         apow *= apow

    122         if n & 1:

KeyboardInterrupt: 

sage: 

This, IMNSHO, is a (quite serious) bug, affecting a very basic functionality of Sage.

That's why I have opened Trac#28620 against "symbolics", and flagged it as "critical"

HTH,

> To unsubscribe from this group and stop receiving emails from it, send an email to sage-...@googlegroups.com.

[Emmanuel Charpentier]

And, BTW:

sage: maxima(t)

(%i/(10*%pi))^(3/2)

sage: maxima(t)^2

-%i/(1000*%pi^3)

HTH,...