Work at SourceForge, help us to make it a better place! We have an immediate need for a Support Technician in our San Francisco or Denver office.

Close

#68 real/imag/conjugate fails on some Complex forms

1.3.0
open-accepted
algebra (9)
5
2009-07-01
2009-06-17
Arnold Doray
No

I get the wrong answer for real/imag for some forms of complex numbers.
For example:
a := (-1*%i)^(1/2)
b := (%i)^(3/2)
are both equivalent, but imag/real gives complex (and incorrect) answers
for "a" (ie, real(a) and imag(a) are both complex), but I get the right
answers for "b". Also "conjugate" fails on both with an error message.

Discussion

    • milestone: --> 1.3.0
    • assigned_to: nobody --> dos-reis
     
  • This is indeed a bug in all versions of OpenAxiom.

    Furthermore, OpenAxiom reports that a and b are not equal as
    Expression Complex Integer

    (3) -> (a=b)@Boolean
    (3) false
    Type: Boolean

     
  • Trying another equivalent form:

    c := %e^(%i*%pi*3/4)

    gives real(c) = real(b) and imag(c) = imag(b), but yet (b=c)@Boolean is FALSE

     
  • Bill Page
    Bill Page
    2009-06-18

    fix for real and imag, improve trigs simplification

     
    Attachments
  • Bill Page
    Bill Page
    2009-06-18

    I have attached a patch that seems to correct this problem. Simplification of complex-valued expressions is at least in part handled by the same internal routines involved in complex trigonometric functions. The perhaps rather poorly name function 'trigs' performs this simplification but is was not able to properly handle roots of complex values. This patch uses De Moivre's method to compute the roots.

     
  • Bill Page
    Bill Page
    2009-06-18

    Note that contrary to the original claim 'a' and 'b' are not equal.

    (10) -> (-1.0*%i)^(1/2)

    (10) 0.7071067811 865475244 - 0.7071067811 865475244 %i
    Type: Complex Float
    (11) -> (0.0+%i)^(3/2)

    (11) - 0.7071067811 865475244 + 0.7071067811 865475244 %i
    Type: Complex Float

    Also I would like to comment that I do not think it is too strange that

    (b=c)@Boolean

    does not always return true when two complex numbers are in fact equal. A more reliable test would make use of the "normalization" done by the 'trigs' operator:

    (trigs(b-c)=0)@Boolean

    Here is a transcript showing the results following the patch.

    OpenAxiom: The Open Scientific Computation Platform
    Version: OpenAxiom 1.3.0-2009-06-07
    Built on Monday June 8, 2009 at 13:21:11
    -----------------------------------------------------------------------------
    Issue )copyright to view copyright notices.
    Issue )summary for a summary of useful system commands.
    Issue )quit to leave OpenAxiom and return to shell.
    -----------------------------------------------------------------------------

    (1) -> )co itrigmnp.spad
    Compiling OpenAxiom source code from file /home/wspage/itrigmnp.spad
    using Spad compiler.
    ITRIGMNP abbreviates package InnerTrigonometricManipulations
    ------------------------------------------------------------------------
    InnerTrigonometricManipulations will be automatically loaded when
    needed from /home/wspage/ITRIGMNP.NRLIB/code.o

    (1) -> a := (-1*%i)^(1/2)

    +----+
    (1) \|- %i
    Type: Expression Complex Integer
    (2) -> b := (%i)^(3/2)

    +--+
    (2) %i\|%i
    Type: Expression Complex Integer
    (3) -> trigs a

    +-+
    \|2
    (3) ------
    1 + %i
    Type: Expression Complex Integer
    (4) -> trigs b

    +-+
    \|2
    (4) - ------
    1 + %i
    Type: Expression Complex Integer
    (5) -> [real a,imag a]

    +-+ +-+
    \|2 \|2
    (5) [----,- ----]
    2 2
    Type: List Expression Integer
    (6) -> [real b,imag b]

    +-+ +-+
    \|2 \|2
    (6) [- ----,----]
    2 2
    Type: List Expression Integer
    (7) -> trigs(a-b)

    +-+
    (7) (1 - %i)\|2
    Type: Expression Complex Integer
    (8) -> c := %e^(%i*%pi*3/4)

    +-+
    \|2
    (8) - ------
    1 + %i
    Type: Expression Complex Integer
    (9) -> trigs(b-c)

    (9) 0
    Type: Expression Complex Integer
    (10) ->

     
    • status: open --> open-accepted