[BUGFIX] Fix crash in parser.

[ginac.git] / ginac / inifcns.cpp

diff --git a/ginac/inifcns.cpp b/ginac/inifcns.cpp
index 6ac1b471abaa4b3d3c39b811054a214d0c312bec..b72d81bdfe8364c37f14df1e7a45558edfe5ace5 100644 (file)
--- a/ginac/inifcns.cpp
+++ b/ginac/inifcns.cpp
@@ -3,7 +3,7 @@
  *  Implementation of GiNaC's initially known functions. */
 
 /*
- *  GiNaC Copyright (C) 1999-2015 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2024 Johannes Gutenberg University Mainz, Germany
  *
  *  This program is free software; you can redistribute it and/or modify
  *  it under the terms of the GNU General Public License as published by
@@ -1021,6 +1021,20 @@ static ex Order_imag_part(const ex & x)
        return Order(x).hold();
 }
 
+static ex Order_power(const ex & x, const ex & e)
+{
+       // Order(x)^e -> Order(x^e) for positive integer e
+       if (is_exactly_a<numeric>(e) && e.info(info_flags::posint))
+               return Order(pow(x, e));
+       // NB: For negative exponents, the above could be wrong.
+       // This is because series() produces Order(x^n) to denote the order where
+       // it gave up. So, Order(x^n) can also be an x^(n+1) term if the x^n term
+       // vanishes. In this situation, 1/Order(x^n) can also be a x^(-n-1) term.
+       // Transforming it to Order(x^-n) would miss that.
+
+       return power(Order(x), e).hold();
+}
+
 static ex Order_expl_derivative(const ex & arg, const symbol & s)
 {
        return Order(arg.diff(s));
@@ -1032,12 +1046,36 @@ REGISTER_FUNCTION(Order, eval_func(Order_eval).
                          expl_derivative_func(Order_expl_derivative).
                          conjugate_func(Order_conjugate).
                          real_part_func(Order_real_part).
-                         imag_part_func(Order_imag_part));
+                         imag_part_func(Order_imag_part).
+                         power_func(Order_power));
 
 //////////
 // Solve linear system
 //////////
 
+class symbolset {
+       exset s;
+       void insert_symbols(const ex &e)
+       {
+               if (is_a<symbol>(e)) {
+                       s.insert(e);
+               } else {
+                       for (const ex &sube : e) {
+                               insert_symbols(sube);
+                       }
+               }
+       }
+public:
+       explicit symbolset(const ex &e)
+       {
+               insert_symbols(e);
+       }
+       bool has(const ex &e) const
+       {
+               return s.find(e) != s.end();
+       }
+};
+
 ex lsolve(const ex &eqns, const ex &symbols, unsigned options)
 {
        // solve a system of linear equations
@@ -1053,20 +1091,20 @@ ex lsolve(const ex &eqns, const ex &symbols, unsigned options)
        }
        
        // syntax checks
-       if (!eqns.info(info_flags::list)) {
-               throw(std::invalid_argument("lsolve(): 1st argument must be a list or an equation"));
+       if (!(eqns.info(info_flags::list) || eqns.info(info_flags::exprseq))) {
+               throw(std::invalid_argument("lsolve(): 1st argument must be a list, a sequence, or an equation"));
        }
        for (size_t i=0; i<eqns.nops(); i++) {
                if (!eqns.op(i).info(info_flags::relation_equal)) {
                        throw(std::invalid_argument("lsolve(): 1st argument must be a list of equations"));
                }
        }
-       if (!symbols.info(info_flags::list)) {
-               throw(std::invalid_argument("lsolve(): 2nd argument must be a list or a symbol"));
+       if (!(symbols.info(info_flags::list) || symbols.info(info_flags::exprseq))) {
+               throw(std::invalid_argument("lsolve(): 2nd argument must be a list, a sequence, or a symbol"));
        }
        for (size_t i=0; i<symbols.nops(); i++) {
                if (!symbols.op(i).info(info_flags::symbol)) {
-                       throw(std::invalid_argument("lsolve(): 2nd argument must be a list of symbols"));
+                       throw(std::invalid_argument("lsolve(): 2nd argument must be a list or a sequence of symbols"));
                }
        }
        
@@ -1077,8 +1115,11 @@ ex lsolve(const ex &eqns, const ex &symbols, unsigned options)
        
        for (size_t r=0; r<eqns.nops(); r++) {
                const ex eq = eqns.op(r).op(0)-eqns.op(r).op(1); // lhs-rhs==0
+               const symbolset syms(eq);
                ex linpart = eq;
                for (size_t c=0; c<symbols.nops(); c++) {
+                       if (!syms.has(symbols.op(c)))
+                               continue;
                        const ex co = eq.coeff(ex_to<symbol>(symbols.op(c)),1);
                        linpart -= co*symbols.op(c);
                        sys(r,c) = co;
@@ -1088,11 +1129,13 @@ ex lsolve(const ex &eqns, const ex &symbols, unsigned options)
        }
        
        // test if system is linear and fill vars matrix
+       const symbolset sys_syms(sys);
+       const symbolset rhs_syms(rhs);
        for (size_t i=0; i<symbols.nops(); i++) {
                vars(i,0) = symbols.op(i);
-               if (sys.has(symbols.op(i)))
+               if (sys_syms.has(symbols.op(i)))
                        throw(std::logic_error("lsolve: system is not linear"));
-               if (rhs.has(symbols.op(i)))
+               if (rhs_syms.has(symbols.op(i)))
                        throw(std::logic_error("lsolve: system is not linear"));
        }