[GiNaC-list] matrix::solve()-related problems

[Vitaly Magerya]

On 01/25/2018 10:10 PM, Richard B. Kreckel wrote:
>> On another related note: there's a flag that tells you if an expression
>> is expanded (status_flag::expanded), which presumably exists so that
>> expand() repeated twice would be a quick no-op. There's no such flag for
>> normal() though, and I'd argue it's much more needed there than in
>> expand(). With such a flag you could pretty much put normal() all over
>> the place without fearing for wasted computation: performance-minded
>> code would have already normal()ized whatever it's passing into
>> is_zero() and the like, so no performance penalty there, and more
>> careless code would get the benefit of always obtaining correct results
>> even when it forgot to normal()ize it's expressions.
> 
> This sounds very reasonable. Would you like to try your hand at a patch?

I've been testing the attached patch for a bit, but I do not see as big
of a performance gain on my workload as I hoped for. This is because
normalizing an already normal expression is already much, much faster
than normalizing it in the first place, so this optimization is nice,
but not as big of a deal.

I haven't covered all the places status_flag::normal could be exploited
though, so don't commit it just yet.
-------------- next part --------------
diff --git a/ginac/flags.h b/ginac/flags.h
index c97229af..c657a4f8 100644
--- a/ginac/flags.h
+++ b/ginac/flags.h
@@ -200,7 +200,8 @@ public:
 		has_no_indices	= 0x0040, // ! (has_indices || has_no_indices) means "don't know"
 		is_positive	= 0x0080,
 		is_negative	= 0x0100,
-		purely_indefinite = 0x0200  // If set in a mul, then it does not contains any terms with determined signs, used in power::expand()
+		purely_indefinite = 0x0200, // If set in a mul, then it does not contains any terms with determined signs, used in power::expand()
+		normal			= 0x0400 ///< .normal() was used to make this expression
 	};
 };
 
diff --git a/ginac/mul.cpp b/ginac/mul.cpp
index 68ea6851..aa6ad700 100644
--- a/ginac/mul.cpp
+++ b/ginac/mul.cpp
@@ -1080,7 +1080,7 @@ ex mul::expand(unsigned options) const
 		}
 	}
 	if (monomial_case) {
-		setflag(status_flags::expanded);
+		setflag(status_flags::expanded | status_flags::normal);
 		return *this;
 	}
 
diff --git a/ginac/normal.cpp b/ginac/normal.cpp
index 231ef7e8..ce0fc1f6 100644
--- a/ginac/normal.cpp
+++ b/ginac/normal.cpp
@@ -2306,6 +2306,10 @@ ex pseries::normal(exmap & repl, exmap & rev_lookup) const
  *  @return normalized expression */
 ex ex::normal() const
 {
+	if (bp->flags & status_flags::normal) {
+		return *this;
+	}
+
 	exmap repl, rev_lookup;
 
 	ex e = bp->normal(repl, rev_lookup);
@@ -2316,7 +2320,9 @@ ex ex::normal() const
 		e = e.subs(repl, subs_options::no_pattern);
 
 	// Convert {numerator, denominator} form back to fraction
-	return e.op(0) / e.op(1);
+	ex r = e.op(0) / e.op(1);
+	r.bp->setflag(status_flags::normal);
+	return r;
 }
 
 /** Get numerator of an expression. If the expression is not of the normal
@@ -2333,10 +2339,10 @@ ex ex::numer() const
 	GINAC_ASSERT(is_a<lst>(e));
 
 	// Re-insert replaced symbols
-	if (repl.empty())
-		return e.op(0);
-	else
-		return e.op(0).subs(repl, subs_options::no_pattern);
+    ex r =
+        repl.empty() ? e.op(0) : e.op(0).subs(repl, subs_options::no_pattern);
+	r.bp->setflag(status_flags::normal);
+	return r;
 }
 
 /** Get denominator of an expression. If the expression is not of the normal
@@ -2353,10 +2359,10 @@ ex ex::denom() const
 	GINAC_ASSERT(is_a<lst>(e));
 
 	// Re-insert replaced symbols
-	if (repl.empty())
-		return e.op(1);
-	else
-		return e.op(1).subs(repl, subs_options::no_pattern);
+    ex r =
+        repl.empty() ? e.op(1) : e.op(1).subs(repl, subs_options::no_pattern);
+	r.bp->setflag(status_flags::normal);
+	return r;
 }
 
 /** Get numerator and denominator of an expression. If the expression is not
@@ -2373,10 +2379,10 @@ ex ex::numer_denom() const
 	GINAC_ASSERT(is_a<lst>(e));
 
 	// Re-insert replaced symbols
-	if (repl.empty())
-		return e;
-	else
-		return e.subs(repl, subs_options::no_pattern);
+    ex r = repl.empty() ? e : e.subs(repl, subs_options::no_pattern);
+	r.let_op(0).bp->setflag(status_flags::normal);
+	r.let_op(1).bp->setflag(status_flags::normal);
+    return r;
 }
 
 
diff --git a/ginac/power.cpp b/ginac/power.cpp
index 62fc3a8d..334fecd3 100644
--- a/ginac/power.cpp
+++ b/ginac/power.cpp
@@ -779,7 +779,7 @@ ex power::expand(unsigned options) const
 {
 	if (is_a<symbol>(basis) && exponent.info(info_flags::integer)) {
 		// A special case worth optimizing.
-		setflag(status_flags::expanded);
+		setflag(status_flags::expanded | status_flags::normal);
 		return *this;
 	}