Better simplification of complex numbers and radicals.

expreduce/builtin_arithmetic.go

@@ -39,7 +39,32 @@ const (
 	FoldFnMul
 )
 
-func typedRealPart(fn FoldFn, i *Integer, r *Rational, f *Flt) Ex {
+func typedRealPart(fn FoldFn, i *Integer, r *Rational, f *Flt, c *Complex) Ex {
+	if c != nil {
+		toReturn := c
+		if f != nil {
+			if fn == FoldFnAdd {
+				toReturn.AddF(f)
+			} else if fn == FoldFnMul {
+				toReturn.MulF(f)
+			}
+		}
+		if r != nil {
+			if fn == FoldFnAdd {
+				toReturn.AddR(r)
+			} else if fn == FoldFnMul {
+				toReturn.MulR(r)
+			}
+		}
+		if i != nil {
+			if fn == FoldFnAdd {
+				toReturn.AddI(i)
+			} else if fn == FoldFnMul {
+				toReturn.MulI(i)
+			}
+		}
+		return toReturn
+	}
 	if f != nil {
 		toReturn := f
 		if r != nil {
@@ -75,10 +100,11 @@ func typedRealPart(fn FoldFn, i *Integer, r *Rational, f *Flt) Ex {
 	return nil
 }
 
-func computeRealPart(fn FoldFn, e *Expression) (Ex, int) {
+func computeNumericPart(fn FoldFn, e *Expression) (Ex, int) {
 	var foldedInt *Integer
 	var foldedRat *Rational
 	var foldedFlt *Flt
+	var foldedComp *Complex
 	for i := 1; i < len(e.Parts); i++ {
 		// TODO: implement short circuiting if we encounter a zero while
 		// multiplying.
@@ -124,18 +150,31 @@ func computeRealPart(fn FoldFn, e *Expression) (Ex, int) {
 			}
 			continue
 		}
-		return typedRealPart(fn, foldedInt, foldedRat, foldedFlt), i
+		asComp, isComp := e.Parts[i].(*Complex)
+		if isComp {
+			if foldedComp == nil {
+				foldedComp = asComp.DeepCopy().(*Complex)
+				continue
+			}
+			if fn == FoldFnAdd {
+				foldedComp.AddC(asComp)
+			} else if fn == FoldFnMul {
+				foldedComp.MulC(asComp)
+			}
+			continue
+		}
+		return typedRealPart(fn, foldedInt, foldedRat, foldedFlt, foldedComp), i
 	}
-	return typedRealPart(fn, foldedInt, foldedRat, foldedFlt), -1
+	return typedRealPart(fn, foldedInt, foldedRat, foldedFlt, foldedComp), -1
 }
 
 // Define a special NumberQ for our purposes since this logic does not support
 // complex numbers yet. TODO(corywalker): fix this.
 func numberQForTermCollection(e Ex) bool {
-	_, ok := e.(*Complex)
-	if ok {
-		return false
-	}
+	// _, ok := e.(*Complex)
+	// if ok {
+	// 	return false
+	// }
 	return numberQ(e)
 }
 
@@ -177,7 +216,7 @@ func collectedToTerm(coeffs []Ex, vars Ex, fullPart Ex) Ex {
 		return fullPart
 	}
 
-	finalC, _ := computeRealPart(FoldFnAdd, NewExpression(append([]Ex{
+	finalC, _ := computeNumericPart(FoldFnAdd, NewExpression(append([]Ex{
 		NewSymbol("System`Plus")}, coeffs...)))
 
 	toAdd := NewExpression([]Ex{NewSymbol("System`Times")})
@@ -243,7 +282,7 @@ func getArithmeticDefinitions() (defs []Definition) {
 			}
 
 			res := this
-			realPart, symStart := computeRealPart(FoldFnAdd, this)
+			realPart, symStart := computeNumericPart(FoldFnAdd, this)
 			if realPart != nil {
 				if symStart == -1 {
 					return realPart
@@ -296,7 +335,7 @@ func getArithmeticDefinitions() (defs []Definition) {
 			}
 
 			res := this
-			realPart, symStart := computeRealPart(FoldFnMul, this)
+			realPart, symStart := computeNumericPart(FoldFnMul, this)
 			if realPart != nil {
 				if symStart == -1 {
 					return realPart

expreduce/ex_complex.go

@@ -11,6 +11,8 @@ type Complex struct {
 }
 
 func (this *Complex) Eval(es *EvalState) Ex {
+	this.Re = this.Re.Eval(es)
+	this.Im = this.Im.Eval(es)
 	if IsSameQ(this.Im, NewInt(0), &es.CASLogger) {
 		return this.Re
 	}
@@ -68,3 +70,65 @@ func (this *Complex) Hash() uint64 {
 	h.Write(b)
 	return h.Sum64()
 }
+
+func (this *Complex) addReal(e Ex) {
+	a, _ := computeNumericPart(FoldFnAdd, E(S("Dummy"), this.Re, e))
+	this.Re = a
+	this.needsEval = true
+}
+
+func (this *Complex) AddI(i *Integer) {
+	this.addReal(i)
+}
+
+func (this *Complex) AddF(f *Flt) {
+	this.addReal(f)
+}
+
+func (this *Complex) AddR(r *Rational) {
+	this.addReal(r)
+}
+
+func (this *Complex) AddC(c *Complex) {
+	a, _ := computeNumericPart(FoldFnAdd, E(S("Dummy"), this.Re, c.Re))
+	b, _ := computeNumericPart(FoldFnAdd, E(S("Dummy"), this.Im, c.Im))
+	this.Re = a
+	this.Im = b
+	this.needsEval = true
+}
+
+func (this *Complex) mulReal(e Ex) {
+	a, _ := computeNumericPart(FoldFnMul, E(S("Dummy"), this.Re, e))
+	b, _ := computeNumericPart(FoldFnMul, E(S("Dummy"), this.Im, e))
+	this.Re = a
+	this.Im = b
+	this.needsEval = true
+}
+
+func (this *Complex) MulI(i *Integer) {
+	this.mulReal(i)
+}
+
+func (this *Complex) MulF(f *Flt) {
+	this.mulReal(f)
+}
+
+func (this *Complex) MulR(r *Rational) {
+	this.mulReal(r)
+}
+
+func (this *Complex) MulC(c *Complex) {
+	// HoldPattern[Complex[x_, y_]*Complex[u_, v_]*rest___] -> Complex[x*u + (y*v)*(-1), x*v + y*u]*rest)
+	// This is ugly. Need to refactor.
+	// Perhaps create "Calculator" utility??
+	// TODO(corywalker) Remove the definition that this implements in code.
+	a, _ := computeNumericPart(FoldFnMul, E(S("Dummy"), this.Re, c.Re))
+	b, _ := computeNumericPart(FoldFnMul, E(S("Dummy"), NewInt(-1), this.Im, c.Im))
+	cc, _ := computeNumericPart(FoldFnMul, E(S("Dummy"), this.Re, c.Im))
+	d, _ := computeNumericPart(FoldFnMul, E(S("Dummy"), this.Im, c.Re))
+	e, _ := computeNumericPart(FoldFnAdd, E(S("Dummy"), a, b))
+	f, _ := computeNumericPart(FoldFnAdd, E(S("Dummy"), cc, d))
+	this.Re = e
+	this.Im = f
+	this.needsEval = true
+}

expreduce/resources/arithmetic.m

@@ -122,7 +122,11 @@
 Times[den_Integer^-1, num_Integer, rest___] := Rational[num,den] * rest;
 Times[ComplexInfinity, rest___] := ComplexInfinity;
 a_Integer?Negative^b_Rational*c_Integer^d_Rational*rest___ := (-1)^b*rest /; (a == -c && b == -d);
-Verbatim[Times][beg___, a_Integer^m_Rational, a_Integer^n_Rational, end___] := beg*end /; (m == -n);
+Verbatim[Times][beg___, a_Integer^m_Rational, a_Integer^n_Rational, end___] := beg*a^(m+n)*end;
+Times[c : (Rational[_Integer, d_Integer] | 
+     Complex[_Integer, Rational[_Integer, d_Integer]]), 
+  Power[a_Integer, Rational[1, r_Integer]]] :=
+ Times[c*a, a^(1/r - 1)] /; (Mod[d, a] === 0 && a > 1)
 Sin[x_]*Cos[x_]^(-1)*rest___ := Tan[x]*rest;
 Cos[x_]*Sin[x_]^(-1)*rest___ := Cot[x]*rest;
 Attributes[Times] = {Flat, Listable, NumericFunction, OneIdentity, Orderless, Protected};
@@ -192,6 +196,17 @@
 
         ESameTest[(-1)^(1/3) a b c, (-2)^(1/3)*2^(-1/3)*a*b*c],
         ESameTest[1/12, (1/12)*2^(-2/3)*2^(2/3)],
+        ESameTest[I/(4 Sqrt[3]), (0+I/12) Sqrt[3]],
+        (* 1/12*Sqrt[3]===1/12*3^(1/2)===12^-1*3^(1/2)===(2*2*3)^-1*3^(1/2)===(2*2)^-1*3^-1*3^(1/2)===(2*2)^-1*3^(1/2-1)===1/4*3^(-1/2)===1/(4 Sqrt[3]) *)
+        ESameTest[1/(4 Sqrt[3]), 1/12*Sqrt[3]],
+        ESameTest[-(I/(4 Sqrt[3])), (0-I/4) 1/Sqrt[3]],
+        ESameTest[-(1/(4 Sqrt[3])), Times[-1/12,Sqrt[3]]],
+        ESameTest[-(5/(4 Sqrt[3])), Times[-5/12,Sqrt[3]]],
+        ESameTest[(3+I/4)/Sqrt[3], Times[1+1/12 I,Sqrt[3]]],
+        ESameTest[(I Sqrt[3])/4, Times[3/12 I,Sqrt[3]]],
+        ESameTest[-(I/(2 Sqrt[3])), Times[-2/12 I,Sqrt[3]]],
+        ESameTest[-(I/(2 Sqrt[3])), Times[2/-12 I,Sqrt[3]]],
+        ESameTest[(3+(5 I)/4)/Sqrt[3], Times[1+5/12 I,Sqrt[3]]],
     ], EKnownFailures[
         ESameTest[-2^(1/3), (-2)*2^(-2/3)],
         ESameTest[-2^(1+a), (-2)*2^(a)],

expreduce/resources/power.m

@@ -841,14 +841,14 @@
 
 ComplexExpand::usage = "`ComplexExpand[e]` returns a complex expansion of `e`.";
 Attributes[ComplexExpand] = {Protected};
-ComplexExpandInner[e_] := e;
-ComplexExpandInner[(a_Integer?Negative)^b_Rational] := 
+complexExpandInner[e_] := e;
+complexExpandInner[(a_Integer?Negative)^b_Rational] := 
   Module[{coeff, inner},
    coeff = ((a^2)^(b/2));
    inner = b*Arg[a];
    coeff*Cos[inner] + I*coeff*Sin[inner]];
 ComplexExpand[exp_] := 
-  Map[ComplexExpandInner, exp, {0, Infinity}] // Expand;
+  Map[complexExpandInner, exp, {0, Infinity}] // Expand;
 Tests`ComplexExpand = {
     ESimpleExamples[
         ESameTest[a, ComplexExpand[a]],

expreduce/resources/simplify.m

@@ -45,6 +45,8 @@
       tryVal = FixedPoint[Replace[#, expandRules]&, e];
       If[LeafCount[tryVal] < LeafCount[e], e = tryVal];
     ];
+    tryVal = ComplexExpand[e];
+    If[LeafCount[tryVal] < LeafCount[e], e = tryVal];
     (* also need to try complexexpand to simplify cases like (-1)^(1/3) (1 + I Sqrt[3]) *)
     e = Replace[e, Sqrt[inner_] :> Sqrt[FactorTerms[inner]]];
     e

expreduce/resources/trig.m

@@ -9,6 +9,7 @@
 Sin[(-5/2)*Pi] := -1;
 Sin[(-3/2)*Pi] := 1;
 Sin[(-1/2)*Pi] := -1;
+Sin[(-1/3)*Pi] := -(Sqrt[3]/2);
 Sin[(1/3)*Pi] := Sqrt[3]/2;
 Sin[(1/2)*Pi] := 1;
 Sin[(3/2)*Pi] := -1;
@@ -25,6 +26,7 @@
 Cos[(-5/2)*Pi] := 0;
 Cos[(-3/2)*Pi] := 0;
 Cos[(-1/2)*Pi] := 0;
+Cos[(-1/3)*Pi] := 1/2;
 Cos[(1/3)*Pi] := 1/2;
 Cos[(1/2)*Pi] := 0;
 Cos[(3/2)*Pi] := 0;