complex: Repair auto-simplificiation for i*i=-1

There was a regression on the i*i=-1 test from the work on constants.
Add it back and extend it to any case that produces real-only results
when auto-simplify is on.

Fixes: #497

Signed-off-by: Christophe de Dinechin <christophe@dinechin.org>

src/arithmetic.cc

@@ -32,6 +32,7 @@
 #include "array.h"
 #include "bignum.h"
 #include "compare.h"
+#include "constants.h"
 #include "decimal.h"
 #include "expression.h"
 #include "fraction.h"
@@ -419,7 +420,12 @@ algebraic_p arithmetic::non_numeric<mul>(algebraic_r x, algebraic_r y)
         if (y->is_one(false))                   // X * 1 = X
             return x;
         if (x->is_symbolic() && x->is_same_as(y))
+        {
+            if (constant_p cst = x->as<constant>())
+                if (cst->is_imaginary_unit())
+                    return integer::make(-1);
             return sq::run(x);                  // X * X = X²
+        }
     }
 
     // Check multiplication of unit objects
@@ -1269,7 +1275,12 @@ algebraic_p arithmetic::evaluate(id          op,
         complex_g xc = complex_p(algebraic_p(x));
         complex_g yc = complex_p(algebraic_p(y));
         if (ops.complex_ok(xc, yc))
+        {
+            if (Settings.AutoSimplify())
+                if (algebraic_p re = xc->is_real())
+                    return re;
             return xc;
+        }
     }
 
     if (!x || !y)

src/complex.cc

@@ -125,6 +125,17 @@ algebraic_g complex::pifrac() const
 }
 
 
+algebraic_p complex::is_real() const
+// ----------------------------------------------------------------------------
+//   Check if the complex is a purely real value
+// ----------------------------------------------------------------------------
+{
+    if (type() == ID_polar)
+        return polar_p(this)->is_real();
+    return rectangular_p(this)->is_real();
+}
+
+
 complex_g complex::conjugate() const
 // ----------------------------------------------------------------------------
 //   Return complex conjugate in a format-independent way
@@ -402,6 +413,7 @@ PARSE_BODY(complex)
 //   f. 1+3ⅈ            ⅈ as a postfix
 //   g. 1-3ⅈ
 //   h. 1∡30            ∡ as a separator
+//   i. ⅈ               Imaginary unit by itself
 //   u. 1_km            _ as a separator for unit objects
 //
 //   Cases a-g generate a rectangular form, case i generates a polar form
@@ -594,6 +606,17 @@ PARSE_BODY(complex)
         last += utf8_size(cp);
     }
 
+    // If we just have the imaginary unit
+    if (type == ID_rectangular && !xlen)
+    {
+        // Build the imaginary unit
+        rectangular_g result = rectangular::make(integer::make(0),
+                                                 integer::make(1));
+        p.out = +result;
+        p.end = last - first + 2*paren;
+        return OK;
+    }
+
     // If we did not find the necessary structure, just skip
     if (type == ID_object || !xlen || !ybeg)
         return SKIP;
@@ -730,6 +753,15 @@ bool rectangular::is_one() const
 }
 
 
+algebraic_p rectangular::is_real() const
+// ----------------------------------------------------------------------------
+//   Check if the complex is a purely real value
+// ----------------------------------------------------------------------------
+{
+    return y()->is_zero(false) ? x() : nullptr;
+}
+
+
 RENDER_BODY(rectangular)
 // ----------------------------------------------------------------------------
 //   Render a complex number in rectangular form
@@ -814,6 +846,21 @@ bool polar::is_one() const
 }
 
 
+algebraic_p polar::is_real() const
+// ----------------------------------------------------------------------------
+//   Check if the complex is a purely real value
+// ----------------------------------------------------------------------------
+{
+    polar_g o = this;
+    algebraic_g pifrac = o->pifrac();
+    if (pifrac->is_zero(false))
+        return o->mod();
+    if (pifrac->is_one(false))
+        return -o->mod();
+    return nullptr;
+}
+
+
 polar_p polar::make(algebraic_r mr, algebraic_r ar, angle_unit unit)
 // ----------------------------------------------------------------------------
 //   Build a normalized polar from given modulus and argument

src/complex.h

@@ -83,6 +83,7 @@ struct complex : algebraic
     algebraic_g         arg(angle_unit unit) const;
     algebraic_g         pifrac() const;
     complex_g           conjugate() const;
+    algebraic_p         is_real() const;
 
     polar_g             as_polar() const;
     rectangular_g       as_rectangular() const;
@@ -166,6 +167,7 @@ struct rectangular : complex
     algebraic_g pifrac() const;
     bool        is_zero() const;
     bool        is_one()  const;
+    algebraic_p is_real() const;
 
     static rectangular_p make(algebraic_r r, algebraic_r i)
     {
@@ -196,6 +198,7 @@ struct polar : complex
     algebraic_g pifrac() const  { return y(); }
     bool        is_zero() const;
     bool        is_one()  const;
+    algebraic_p is_real() const;
 
     static polar_p make(algebraic_r mod, algebraic_r arg, angle_unit unit);
 

src/constants.h

@@ -59,6 +59,15 @@ struct constant : symbol
         return text::value(size);
     }
     algebraic_p value() const;
+    bool is_imaginary_unit() const              { return matches("i"); }
+    bool is_pi() const                          { return matches("π"); }
+    bool matches(cstring ref) const
+    {
+        size_t nlen = strlen(ref);
+        size_t len = 0;
+        utf8 txt = name(&len);
+        return len == nlen && memcmp(ref, txt, len) == 0;
+    }
 
 public:
     OBJECT_DECL(constant);

src/tests.cc

@@ -2804,7 +2804,8 @@ void tests::fraction_decimal_conversions()
     test("→Q", ENTER).expect("[ 1/4-1/2ⅈ 3/4 ]");
 
     step("Expressions");
-    test(CLEAR, "355 113 / pi -", ENTER) .expect("'355/113-π'");
+    test(CLEAR, "355 113 /",
+         LSHIFT, I, F1, F1, "-", ENTER) .expect("'355/113-π'");
     test("→Num", ENTER).expect("2.66764 18906 2⁳⁻⁷");
 
     step("Restoring small fraction mode")
@@ -4029,6 +4030,17 @@ void tests::auto_simplification()
     step("i*i == -1");
     test(CLEAR, "ⅈ", ENTER, ENTER, MUL).expect("-1");
 
+    step("i*i == -1 (symbolic constant)");
+    test(CLEAR, LSHIFT, I, F1, F3, ENTER, MUL).expect("-1");
+
+    step("Simplification of rectangular real-only results");
+    test(CLEAR, "0ⅈ3 0ⅈ5", ENTER, MUL).expect("-15");
+    test(CLEAR, "0ⅈ3 0-ⅈ5", ENTER, MUL).expect("15");
+
+    step("Simplification of polar real-only results");
+    test(CLEAR, "2∡90 3∡90", ENTER, MUL).expect("-6");
+    test(CLEAR, "2∡90 3∡-90", ENTER, MUL).expect("6");
+
     step("Applies when building a matrix");
     test(CLEAR, "[[3 0 2][2 0 -2][ 0 1 1 ]] [x y z] *", ENTER)
         .expect("[ '3·x+2·z' '2·x+-2·z' 'y+z' ]");