Updated ConstantMatrix.Multiply to accept more types + added tests for more Polynomial functions

Author: Kwesi Rutledge

symbolic/constant_matrix.go

@@ -274,6 +274,32 @@ func (km KMatrix) Multiply(e interface{}) Expression {
 
 	case K:
 		return km.Multiply(float64(right)) // Reuse float64 case
+	case Polynomial:
+		// Choose the correct output type based on the size of km
+		nR, nC := km.Dims()[0], km.Dims()[1]
+		switch {
+		case (nR == 1) && (nC == 1):
+			// If the output is a scalar, return a scalar
+			return km[0][0].Multiply(right)
+		case nC == 1:
+			// If the output is a vector, return a vector
+			var outputVec PolynomialVector = make([]Polynomial, nR)
+			for rIndex := 0; rIndex < nR; rIndex++ {
+				outputVec[rIndex] = km[rIndex][0].Multiply(right.Copy()).(Polynomial)
+			}
+			return outputVec
+		default:
+			// If the output is a matrix, return a matrix
+			var outputMat PolynomialMatrix = make([][]Polynomial, nR)
+			for rIndex := 0; rIndex < nR; rIndex++ {
+				outputMat[rIndex] = make([]Polynomial, nC)
+				for cIndex := 0; cIndex < nC; cIndex++ {
+					outputMat[rIndex][cIndex] = km[rIndex][cIndex].Multiply(right.Copy()).(Polynomial)
+				}
+			}
+			return outputMat
+		}
+
 	case *mat.VecDense:
 		// Use gonum's built-in multiplication function
 		var product mat.VecDense

symbolic/polynomial.go

@@ -239,12 +239,6 @@ func (p Polynomial) Minus(e interface{}) Expression {
 		return Minus(p, eAsE)
 	}
 
-	// Constants
-	switch right := e.(type) {
-	case float64:
-		return p.Copy().Minus(K(right))
-	}
-
 	// If the function has reached this point, then
 	// the input is not recognized
 	panic(
@@ -429,7 +423,13 @@ Description:
 	The transpose operator when applied to a scalar is just the same scalar object.
 */
 func (p Polynomial) Transpose() Expression {
-	return p
+	// Input Processing
+	err := p.Check()
+	if err != nil {
+		panic(err)
+	}
+
+	return p.Copy()
 }
 
 /*
@@ -624,8 +624,19 @@ func (p Polynomial) DerivativeWrt(vIn Variable) Expression {
 		}
 
 		// Append
-		components := monomial.DerivativeWrt(vIn)
-		derivative.Monomials = append(derivative.Monomials, components.(Monomial))
+		dMonomial := monomial.DerivativeWrt(vIn)
+		switch component := dMonomial.(type) {
+		case Monomial:
+			derivative.Monomials = append(derivative.Monomials, component)
+		case K:
+			// Skip zero monomials
+			if float64(component) == 0.0 {
+				continue
+			}
+			derivative.Monomials = append(derivative.Monomials, component.ToMonomial())
+		default:
+			panic(fmt.Errorf("Unexpected type in Polynomial.Derivative: %T", component))
+		}
 	}
 
 	// If the derivative is empty, then return 0.0
@@ -780,7 +791,7 @@ func (p Polynomial) Substitute(vIn Variable, eIn ScalarExpression) Expression {
 	var out Expression = K(0.0)
 	for _, monomial := range p.Monomials {
 		newMonomial := monomial.Substitute(vIn, eIn)
-		out = out.Plus(newMonomial)
+		out = out.Plus(newMonomial).(Polynomial).Simplify()
 	}
 
 	return out