(GiNaC is Not a CAS (Computer Algebra System)) is a C++ library for symbolic calculations.
| normal.cpp (ginac-1.7.11.tar.bz2) | : | normal.cpp (ginac-1.8.0.tar.bz2) |
|---|
| | |
| skipping to change at line 1776 | | skipping to change at line 1776 |
|---|
| | |
| /** Compute square-free factorization of multivariate polynomial a(x) using | | /** Compute square-free factorization of multivariate polynomial a(x) using |
| * Yun's algorithm. Used internally by sqrfree(). | | * Yun's algorithm. Used internally by sqrfree(). |
| * | | * |
| * @param a multivariate polynomial over Z[X], treated here as univariate | | * @param a multivariate polynomial over Z[X], treated here as univariate |
| * polynomial in x (needs not be expanded). | | * polynomial in x (needs not be expanded). |
| * @param x variable to factor in | | * @param x variable to factor in |
| * @return vector of expairs (factor, exponent), sorted by exponents */ | | * @return vector of expairs (factor, exponent), sorted by exponents */ |
| static epvector sqrfree_yun(const ex &a, const symbol &x) | | static epvector sqrfree_yun(const ex &a, const symbol &x) |
| { | | { |
|
| // make polynomial primitive w.r.t. x | | ex w = a; |
| ex w = a.primpart(x); | | |
| ex z = w.diff(x); | | ex z = w.diff(x); |
| ex g = gcd(w, z); | | ex g = gcd(w, z); |
| if (g.is_zero()) { | | if (g.is_zero()) { |
| // manifest zero or hidden zero | | // manifest zero or hidden zero |
| return {}; | | return {}; |
| } | | } |
| if (g.is_equal(_ex1)) { | | if (g.is_equal(_ex1)) { |
| // w(x) and w'(x) share no factors: w(x) is square-free | | // w(x) and w'(x) share no factors: w(x) is square-free |
| return {expair(a, _ex1)}; | | return {expair(a, _ex1)}; |
| } | | } |
| | |
| epvector factors; | | epvector factors; |
| ex i = 0; // exponent | | ex i = 0; // exponent |
| do { | | do { |
| w = quo(w, g, x); | | w = quo(w, g, x); |
|
| | if (w.is_zero()) { |
| | // hidden zero |
| | break; |
| | } |
| z = quo(z, g, x) - w.diff(x); | | z = quo(z, g, x) - w.diff(x); |
| i += 1; | | i += 1; |
| if (w.is_equal(x)) { | | if (w.is_equal(x)) { |
| // shortcut for x^n with n ∈ ℕ | | // shortcut for x^n with n ∈ ℕ |
| i += quo(z, w.diff(x), x); | | i += quo(z, w.diff(x), x); |
| factors.push_back(expair(w, i)); | | factors.push_back(expair(w, i)); |
| break; | | break; |
| } | | } |
| g = gcd(w, z); | | g = gcd(w, z); |
| if (!g.is_equal(_ex1)) { | | if (!g.is_equal(_ex1)) { |
| factors.push_back(expair(g, i)); | | factors.push_back(expair(g, i)); |
| } | | } |
| } while (!z.is_zero()); | | } while (!z.is_zero()); |
| | |
| // correct for lost factor | | // correct for lost factor |
|
| // (NB: This factor is not necessarily the unit and content lost when | | // (being based on GCDs, Yun's algorithm only finds factors up to a unit) |
| // we made the polynomial primitive since Yun's algorithm only finds | | |
| // factors up to a unit.) | | |
| const ex lost_factor = quo(a, mul{factors}, x); | | const ex lost_factor = quo(a, mul{factors}, x); |
| if (lost_factor.is_equal(_ex1)) { | | if (lost_factor.is_equal(_ex1)) { |
| // trivial lost factor | | // trivial lost factor |
| return factors; | | return factors; |
| } | | } |
|
| if (factors[0].coeff.is_equal(1)) { | | if (!factors.empty() && factors[0].coeff.is_equal(1)) { |
| // multiply factor^1 with lost_factor | | // multiply factor^1 with lost_factor |
| factors[0].rest *= lost_factor; | | factors[0].rest *= lost_factor; |
| return factors; | | return factors; |
| } | | } |
|
| // only factors with higher exponent: prepend lost_factor^1 to the result
s | | // no factor^1: prepend lost_factor^1 to the results |
| epvector results = {expair(lost_factor, 1)}; | | epvector results = {expair(lost_factor, 1)}; |
| std::move(factors.begin(), factors.end(), std::back_inserter(results)); | | std::move(factors.begin(), factors.end(), std::back_inserter(results)); |
| return results; | | return results; |
| } | | } |
| | |
| /** Compute a square-free factorization of a multivariate polynomial in Q[X]. | | /** Compute a square-free factorization of a multivariate polynomial in Q[X]. |
| * | | * |
| * @param a multivariate polynomial over Q[X] (needs not be expanded) | | * @param a multivariate polynomial over Q[X] (needs not be expanded) |
| * @param l lst of variables to factor in, may be left empty for autodetection | | * @param l lst of variables to factor in, may be left empty for autodetection |
| * @return a square-free factorization of \p a. | | * @return a square-free factorization of \p a. |
| | |
| skipping to change at line 1994 | | skipping to change at line 1995 |
|---|
| * Note: The internal normal() functions (= basic::normal() and overloaded | | * Note: The internal normal() functions (= basic::normal() and overloaded |
| * functions) all return lists of the form {numerator, denominator}. This | | * functions) all return lists of the form {numerator, denominator}. This |
| * is to get around mul::eval()'s automatic expansion of numeric coefficients. | | * is to get around mul::eval()'s automatic expansion of numeric coefficients. |
| * E.g. (a+b)/3 is automatically converted to a/3+b/3 but we want to keep | | * E.g. (a+b)/3 is automatically converted to a/3+b/3 but we want to keep |
| * the information that (a+b) is the numerator and 3 is the denominator. | | * the information that (a+b) is the numerator and 3 is the denominator. |
| */ | | */ |
| | |
| /** Create a symbol for replacing the expression "e" (or return a previously | | /** Create a symbol for replacing the expression "e" (or return a previously |
| * assigned symbol). The symbol and expression are appended to repl, for | | * assigned symbol). The symbol and expression are appended to repl, for |
| * a later application of subs(). | | * a later application of subs(). |
|
| | * An entry in the replacement table repl can be changed in some cases. |
| | * If it was altered, we need to provide the modifier for the previously build |
| | expressions. |
| | * The modifier is an (ordered) list, because those substitutions need to be do |
| | ne in the |
| | * incremental order. |
| | * As an example let us consider a rationalisation of the expression |
| | * e = exp(2*x)*cos(exp(2*x)+1)*exp(x) |
| | * The first factor GiNaC denotes by something like symbol1 and will record: |
| | * e =symbol1*cos(symbol1 + 1)*exp(x) |
| | * repl = {symbol1 : exp(2*x)} |
| | * Similarly, the second factor would be denoted as symbol2 and we will have |
| | * e =symbol1*symbol2*exp(x) |
| | * repl = {symbol1 : exp(2*x), symbol2 : cos(symbol1 + 1)} |
| | * Denoting the third term as symbol3 GiNaC is willing to re-think exp(2*x) as |
| | * symbol3^2 rather than just symbol1. Here are two issues: |
| | * 1) The replacement "symbol1 -> symbol3^2" in the previous part of the expres |
| | sion |
| | * needs to be done outside of the present routine; |
| | * 2) The pair "symbol1 : exp(2*x)" shall be deleted from the replacement table |
| | repl. |
| | * However, this will create illegal substitution "symbol2 : cos(symbol1 + |
| | 1)" with |
| | * undefined symbol1. |
| | * These both problems are mitigated through the additions of the record |
| | * "symbol1==symbol3^2" to the list modifier. Changed length of the modifier si |
| | gnals |
| | * to the calling code that the previous portion of the expression needs to be |
| | * altered (it solves 1). Thus GiNaC can record now |
| | * e =symbol3^2*symbol2*symbol3 |
| | * repl = {symbol2 : cos(symbol1 + 1), symbol3 : exp(x)} |
| | * modifier = {symbol1==symbol3^2} |
| | * Then, doing the backward substitutions the list modifier will be used to res |
| | tore |
| | * such iterative substitutions in the right way (this solves 2). |
| * @see ex::normal */ | | * @see ex::normal */ |
|
| static ex replace_with_symbol(const ex & e, exmap & repl, exmap & rev_lookup) | | static ex replace_with_symbol(const ex & e, exmap & repl, exmap & rev_lookup, ls
t & modifier) |
| { | | { |
| // Since the repl contains replaced expressions we should search for them | | // Since the repl contains replaced expressions we should search for them |
| ex e_replaced = e.subs(repl, subs_options::no_pattern); | | ex e_replaced = e.subs(repl, subs_options::no_pattern); |
| | |
| // Expression already replaced? Then return the assigned symbol | | // Expression already replaced? Then return the assigned symbol |
| auto it = rev_lookup.find(e_replaced); | | auto it = rev_lookup.find(e_replaced); |
| if (it != rev_lookup.end()) | | if (it != rev_lookup.end()) |
| return it->second; | | return it->second; |
| | |
|
| | // The expression can be the base of substituted power, which requires a |
| | more careful search |
| | if (! is_a<numeric>(e_replaced)) |
| | for (auto & it : repl) |
| | if (is_a<power>(it.second) && e_replaced.is_equal(it.seco |
| | nd.op(0))) { |
| | ex degree = pow(it.second.op(1), _ex_1); |
| | if (is_a<numeric>(degree) && ex_to<numeric>(degre |
| | e).is_integer()) |
| | return pow(it.first, degree); |
| | } |
| | |
| | // We treat powers and the exponent functions differently because |
| | // they can be rationalised more efficiently |
| | if (is_a<function>(e_replaced) && is_ex_the_function(e_replaced, exp)) { |
| | for (auto & it : repl) { |
| | if (is_a<function>(it.second) && is_ex_the_function(e_rep |
| | laced, exp)) { |
| | ex ratio = normal(e_replaced.op(0) / it.second.op |
| | (0)); |
| | if (is_a<numeric>(ratio) && ex_to<numeric>(ratio) |
| | .is_rational()) { |
| | // Different exponents can be treated as |
| | powers of the same basic equation |
| | if (ex_to<numeric>(ratio).is_integer()) { |
| | // If ratio is an integer then th |
| | is is simply the power of the existing symbol. |
| | // std::clog << e_replaced << " i |
| | s a " << ratio << " power of " << it.first << std::endl; |
| | return dynallocate<power>(it.firs |
| | t, ratio); |
| | } else { |
| | // otherwise we need to give the |
| | replacement pattern to change |
| | // the previous expression... |
| | ex es = dynallocate<symbol>(); |
| | ex Num = numer(ratio); |
| | modifier.append(it.first == power |
| | (es, denom(ratio))); |
| | // std::clog << e_replaced << " i |
| | s power " << Num << " and " |
| | // << it.first << |
| | " is power " << denom(ratio) << " of the common base " |
| | // << exp(e_replac |
| | ed.op(0)/Num) << std::endl; |
| | // ... and modify the replacemen |
| | t tables |
| | rev_lookup.erase(it.second); |
| | rev_lookup.insert({exp(e_replaced |
| | .op(0)/Num), es}); |
| | repl.erase(it.first); |
| | repl.insert({es, exp(e_replaced.o |
| | p(0)/Num)}); |
| | return dynallocate<power>(es, Num |
| | ); |
| | } |
| | } |
| | } |
| | } |
| | } else if (is_a<power>(e_replaced) && !is_a<numeric>(e_replaced.op(0)) // |
| | We do not replace simple monomials like x^3 or sqrt(2) |
| | && ! (is_a<symbol>(e_replaced.op(0)) |
| | && is_a<numeric>(e_replaced.op(1)) && ex_to<numeric>(e_rep |
| | laced.op(1)).is_integer())) { |
| | for (auto & it : repl) { |
| | if (e_replaced.op(0).is_equal(it.second) // The base is a |
| | n allocated symbol or base of power |
| | || (is_a<power>(it.second) && e_replaced.op(0).is_equ |
| | al(it.second.op(0)))) { |
| | ex ratio; // We bind together two above cases |
| | if (is_a<power>(it.second)) |
| | ratio = normal(e_replaced.op(1) / it.seco |
| | nd.op(1)); |
| | else |
| | ratio = e_replaced.op(1); |
| | if (is_a<numeric>(ratio) && ex_to<numeric>(ratio) |
| | .is_rational()) { |
| | // Different powers can be treated as pow |
| | ers of the same basic equation |
| | if (ex_to<numeric>(ratio).is_integer()) { |
| | // If ratio is an integer then th |
| | is is simply the power of the existing symbol. |
| | //std::clog << e_replaced << " is |
| | a " << ratio << " power of " << it.first << std::endl; |
| | return dynallocate<power>(it.firs |
| | t, ratio); |
| | } else { |
| | // otherwise we need to give the |
| | replacement pattern to change |
| | // the previous expression... |
| | ex es = dynallocate<symbol>(); |
| | ex Num = numer(ratio); |
| | modifier.append(it.first == power |
| | (es, denom(ratio))); |
| | //std::clog << e_replaced << " is |
| | power " << Num << " and " |
| | // << it.first << |
| | " is power " << denom(ratio) << " of the common base " |
| | // << pow(e_replac |
| | ed.op(0), e_replaced.op(1)/Num) << std::endl; |
| | // ... and modify the replacemen |
| | t tables |
| | rev_lookup.erase(it.second); |
| | rev_lookup.insert({pow(e_replaced |
| | .op(0), e_replaced.op(1)/Num), es}); |
| | repl.erase(it.first); |
| | repl.insert({es, pow(e_replaced.o |
| | p(0), e_replaced.op(1)/Num)}); |
| | return dynallocate<power>(es, Num |
| | ); |
| | } |
| | } |
| | } |
| | } |
| | // There is no existing substitution, thus we are creating a new |
| | one. |
| | // This needs to be done separately to treat possible occurrences |
| | of |
| | // b = e_replaced.op(0) elsewhere in the expression as pow(b, 1). |
| | ex degree = pow(e_replaced.op(1), _ex_1); |
| | if (is_a<numeric>(degree) && ex_to<numeric>(degree).is_integer()) |
| | { |
| | ex es = dynallocate<symbol>(); |
| | modifier.append(e_replaced.op(0) == power(es, degree)); |
| | repl.insert({es, e_replaced}); |
| | rev_lookup.insert({e_replaced, es}); |
| | return es; |
| | } |
| | } |
| | |
| // Otherwise create new symbol and add to list, taking care that the | | // Otherwise create new symbol and add to list, taking care that the |
| // replacement expression doesn't itself contain symbols from repl, | | // replacement expression doesn't itself contain symbols from repl, |
| // because subs() is not recursive | | // because subs() is not recursive |
| ex es = dynallocate<symbol>(); | | ex es = dynallocate<symbol>(); |
| repl.insert(std::make_pair(es, e_replaced)); | | repl.insert(std::make_pair(es, e_replaced)); |
| rev_lookup.insert(std::make_pair(e_replaced, es)); | | rev_lookup.insert(std::make_pair(e_replaced, es)); |
| return es; | | return es; |
| } | | } |
| | |
| /** Create a symbol for replacing the expression "e" (or return a previously | | /** Create a symbol for replacing the expression "e" (or return a previously |
| | |
| skipping to change at line 2045 | | skipping to change at line 2163 |
|---|
| } | | } |
| | |
| /** Function object to be applied by basic::normal(). */ | | /** Function object to be applied by basic::normal(). */ |
| struct normal_map_function : public map_function { | | struct normal_map_function : public map_function { |
| ex operator()(const ex & e) override { return normal(e); } | | ex operator()(const ex & e) override { return normal(e); } |
| }; | | }; |
| | |
| /** Default implementation of ex::normal(). It normalizes the children and | | /** Default implementation of ex::normal(). It normalizes the children and |
| * replaces the object with a temporary symbol. | | * replaces the object with a temporary symbol. |
| * @see ex::normal */ | | * @see ex::normal */ |
|
| ex basic::normal(exmap & repl, exmap & rev_lookup) const | | ex basic::normal(exmap & repl, exmap & rev_lookup, lst & modifier) const |
| { | | { |
| if (nops() == 0) | | if (nops() == 0) |
|
| return dynallocate<lst>({replace_with_symbol(*this, repl, rev_loo
kup), _ex1}); | | return dynallocate<lst>({replace_with_symbol(*this, repl, rev_loo
kup, modifier), _ex1}); |
| | |
| normal_map_function map_normal; | | normal_map_function map_normal; |
|
| return dynallocate<lst>({replace_with_symbol(map(map_normal), repl, rev_l | | int nmod = modifier.nops(); // To watch new modifiers to the replacement |
| ookup), _ex1}); | | list |
| | lst result = dynallocate<lst>({replace_with_symbol(map(map_normal), repl, |
| | rev_lookup, modifier), _ex1}); |
| | for (int imod = nmod; imod < modifier.nops(); ++imod) { |
| | exmap this_repl; |
| | this_repl.insert(std::make_pair(modifier.op(imod).op(0), modifier |
| | .op(imod).op(1))); |
| | result = ex_to<lst>(result.subs(this_repl, subs_options::no_patte |
| | rn)); |
| | } |
| | |
| | return result; |
| } | | } |
| | |
| /** Implementation of ex::normal() for symbols. This returns the unmodified symb
ol. | | /** Implementation of ex::normal() for symbols. This returns the unmodified symb
ol. |
| * @see ex::normal */ | | * @see ex::normal */ |
|
| ex symbol::normal(exmap & repl, exmap & rev_lookup) const | | ex symbol::normal(exmap & repl, exmap & rev_lookup, lst & modifier) const |
| { | | { |
| return dynallocate<lst>({*this, _ex1}); | | return dynallocate<lst>({*this, _ex1}); |
| } | | } |
| | |
| /** Implementation of ex::normal() for a numeric. It splits complex numbers | | /** Implementation of ex::normal() for a numeric. It splits complex numbers |
| * into re+I*im and replaces I and non-rational real numbers with a temporary | | * into re+I*im and replaces I and non-rational real numbers with a temporary |
| * symbol. | | * symbol. |
| * @see ex::normal */ | | * @see ex::normal */ |
|
| ex numeric::normal(exmap & repl, exmap & rev_lookup) const | | ex numeric::normal(exmap & repl, exmap & rev_lookup, lst & modifier) const |
| { | | { |
| numeric num = numer(); | | numeric num = numer(); |
| ex numex = num; | | ex numex = num; |
| | |
| if (num.is_real()) { | | if (num.is_real()) { |
| if (!num.is_integer()) | | if (!num.is_integer()) |
|
| numex = replace_with_symbol(numex, repl, rev_lookup); | | numex = replace_with_symbol(numex, repl, rev_lookup, modi
fier); |
| } else { // complex | | } else { // complex |
| numeric re = num.real(), im = num.imag(); | | numeric re = num.real(), im = num.imag(); |
|
| ex re_ex = re.is_rational() ? re : replace_with_symbol(re, repl, | | ex re_ex = re.is_rational() ? re : replace_with_symbol(re, repl, |
| rev_lookup); | | rev_lookup, modifier); |
| ex im_ex = im.is_rational() ? im : replace_with_symbol(im, repl, | | ex im_ex = im.is_rational() ? im : replace_with_symbol(im, repl, |
| rev_lookup); | | rev_lookup, modifier); |
| numex = re_ex + im_ex * replace_with_symbol(I, repl, rev_lookup); | | numex = re_ex + im_ex * replace_with_symbol(I, repl, rev_lookup, |
| | modifier); |
| } | | } |
| | |
| // Denominator is always a real integer (see numeric::denom()) | | // Denominator is always a real integer (see numeric::denom()) |
| return dynallocate<lst>({numex, denom()}); | | return dynallocate<lst>({numex, denom()}); |
| } | | } |
| | |
| /** Fraction cancellation. | | /** Fraction cancellation. |
| * @param n numerator | | * @param n numerator |
| * @param d denominator | | * @param d denominator |
| * @return cancelled fraction {n, d} as a list */ | | * @return cancelled fraction {n, d} as a list */ |
| | |
| skipping to change at line 2147 | | skipping to change at line 2273 |
|---|
| } | | } |
| | |
| // Return result as list | | // Return result as list |
| //std::clog << " returns num = " << num << ", den = " << den << ", pre_factor =
" << pre_factor << std::endl; | | //std::clog << " returns num = " << num << ", den = " << den << ", pre_factor =
" << pre_factor << std::endl; |
| return dynallocate<lst>({num * pre_factor.numer(), den * pre_factor.denom
()}); | | return dynallocate<lst>({num * pre_factor.numer(), den * pre_factor.denom
()}); |
| } | | } |
| | |
| /** Implementation of ex::normal() for a sum. It expands terms and performs | | /** Implementation of ex::normal() for a sum. It expands terms and performs |
| * fractional addition. | | * fractional addition. |
| * @see ex::normal */ | | * @see ex::normal */ |
|
| ex add::normal(exmap & repl, exmap & rev_lookup) const | | ex add::normal(exmap & repl, exmap & rev_lookup, lst & modifier) const |
| { | | { |
| // Normalize children and split each one into numerator and denominator | | // Normalize children and split each one into numerator and denominator |
| exvector nums, dens; | | exvector nums, dens; |
| nums.reserve(seq.size()+1); | | nums.reserve(seq.size()+1); |
| dens.reserve(seq.size()+1); | | dens.reserve(seq.size()+1); |
|
| | int nmod = modifier.nops(); // To watch new modifiers to the replacement
list |
| for (auto & it : seq) { | | for (auto & it : seq) { |
|
| ex n = ex_to<basic>(recombine_pair_to_ex(it)).normal(repl, rev_lo
okup); | | ex n = ex_to<basic>(recombine_pair_to_ex(it)).normal(repl, rev_lo
okup, modifier); |
| nums.push_back(n.op(0)); | | nums.push_back(n.op(0)); |
| dens.push_back(n.op(1)); | | dens.push_back(n.op(1)); |
| } | | } |
|
| ex n = ex_to<numeric>(overall_coeff).normal(repl, rev_lookup); | | ex n = ex_to<numeric>(overall_coeff).normal(repl, rev_lookup, modifier); |
| nums.push_back(n.op(0)); | | nums.push_back(n.op(0)); |
| dens.push_back(n.op(1)); | | dens.push_back(n.op(1)); |
| GINAC_ASSERT(nums.size() == dens.size()); | | GINAC_ASSERT(nums.size() == dens.size()); |
| | |
| // Now, nums is a vector of all numerators and dens is a vector of | | // Now, nums is a vector of all numerators and dens is a vector of |
| // all denominators | | // all denominators |
| //std::clog << "add::normal uses " << nums.size() << " summands:\n"; | | //std::clog << "add::normal uses " << nums.size() << " summands:\n"; |
| | |
| // Add fractions sequentially | | // Add fractions sequentially |
| auto num_it = nums.begin(), num_itend = nums.end(); | | auto num_it = nums.begin(), num_itend = nums.end(); |
| auto den_it = dens.begin(), den_itend = dens.end(); | | auto den_it = dens.begin(), den_itend = dens.end(); |
| //std::clog << " num = " << *num_it << ", den = " << *den_it << std::endl; | | //std::clog << " num = " << *num_it << ", den = " << *den_it << std::endl; |
|
| | for (int imod = nmod; imod < modifier.nops(); ++imod) { |
| | while (num_it != num_itend) { |
| | *num_it = num_it->subs(modifier.op(imod), subs_options::n |
| | o_pattern); |
| | ++num_it; |
| | *den_it = den_it->subs(modifier.op(imod), subs_options::n |
| | o_pattern); |
| | ++den_it; |
| | } |
| | // Reset iterators for the next round |
| | num_it = nums.begin(); |
| | den_it = dens.begin(); |
| | } |
| | |
| ex num = *num_it++, den = *den_it++; | | ex num = *num_it++, den = *den_it++; |
| while (num_it != num_itend) { | | while (num_it != num_itend) { |
| //std::clog << " num = " << *num_it << ", den = " << *den_it << std::endl; | | //std::clog << " num = " << *num_it << ", den = " << *den_it << std::endl; |
| ex next_num = *num_it++, next_den = *den_it++; | | ex next_num = *num_it++, next_den = *den_it++; |
| | |
| // Trivially add sequences of fractions with identical denominato
rs | | // Trivially add sequences of fractions with identical denominato
rs |
| while ((den_it != den_itend) && next_den.is_equal(*den_it)) { | | while ((den_it != den_itend) && next_den.is_equal(*den_it)) { |
| next_num += *num_it; | | next_num += *num_it; |
| num_it++; den_it++; | | num_it++; den_it++; |
| } | | } |
| | |
| skipping to change at line 2198 | | skipping to change at line 2337 |
|---|
| } | | } |
| //std::clog << " common denominator = " << den << std::endl; | | //std::clog << " common denominator = " << den << std::endl; |
| | |
| // Cancel common factors from num/den | | // Cancel common factors from num/den |
| return frac_cancel(num, den); | | return frac_cancel(num, den); |
| } | | } |
| | |
| /** Implementation of ex::normal() for a product. It cancels common factors | | /** Implementation of ex::normal() for a product. It cancels common factors |
| * from fractions. | | * from fractions. |
| * @see ex::normal() */ | | * @see ex::normal() */ |
|
| ex mul::normal(exmap & repl, exmap & rev_lookup) const | | ex mul::normal(exmap & repl, exmap & rev_lookup, lst & modifier) const |
| { | | { |
| // Normalize children, separate into numerator and denominator | | // Normalize children, separate into numerator and denominator |
| exvector num; num.reserve(seq.size()); | | exvector num; num.reserve(seq.size()); |
| exvector den; den.reserve(seq.size()); | | exvector den; den.reserve(seq.size()); |
| ex n; | | ex n; |
|
| | int nmod = modifier.nops(); // To watch new modifiers to the replacement
list |
| for (auto & it : seq) { | | for (auto & it : seq) { |
|
| n = ex_to<basic>(recombine_pair_to_ex(it)).normal(repl, rev_looku
p); | | n = ex_to<basic>(recombine_pair_to_ex(it)).normal(repl, rev_looku
p, modifier); |
| num.push_back(n.op(0)); | | num.push_back(n.op(0)); |
| den.push_back(n.op(1)); | | den.push_back(n.op(1)); |
| } | | } |
|
| n = ex_to<numeric>(overall_coeff).normal(repl, rev_lookup); | | n = ex_to<numeric>(overall_coeff).normal(repl, rev_lookup, modifier); |
| num.push_back(n.op(0)); | | num.push_back(n.op(0)); |
| den.push_back(n.op(1)); | | den.push_back(n.op(1)); |
|
| | auto num_it = num.begin(), num_itend = num.end(); |
| | auto den_it = den.begin(), den_itend = den.end(); |
| | for (int imod = nmod; imod < modifier.nops(); ++imod) { |
| | while (num_it != num_itend) { |
| | *num_it = num_it->subs(modifier.op(imod), subs_options::n |
| | o_pattern); |
| | ++num_it; |
| | *den_it = den_it->subs(modifier.op(imod), subs_options::n |
| | o_pattern); |
| | ++den_it; |
| | } |
| | num_it = num.begin(); |
| | den_it = den.begin(); |
| | } |
| | |
| // Perform fraction cancellation | | // Perform fraction cancellation |
| return frac_cancel(dynallocate<mul>(num), dynallocate<mul>(den)); | | return frac_cancel(dynallocate<mul>(num), dynallocate<mul>(den)); |
| } | | } |
| | |
|
| /** Implementation of ex::normal([B) for powers. It normalizes the basis, | | /** Implementation of ex::normal() for powers. It normalizes the basis, |
| * distributes integer exponents to numerator and denominator, and replaces | | * distributes integer exponents to numerator and denominator, and replaces |
| * non-integer powers by temporary symbols. | | * non-integer powers by temporary symbols. |
| * @see ex::normal */ | | * @see ex::normal */ |
|
| ex power::normal(exmap & repl, exmap & rev_lookup) const | | ex power::normal(exmap & repl, exmap & rev_lookup, lst & modifier) const |
| { | | { |
| // Normalize basis and exponent (exponent gets reassembled) | | // Normalize basis and exponent (exponent gets reassembled) |
|
| ex n_basis = ex_to<basic>(basis).normal(repl, rev_lookup); | | int nmod = modifier.nops(); // To watch new modifiers to the replacement |
| ex n_exponent = ex_to<basic>(exponent).normal(repl, rev_lookup); | | list |
| | ex n_basis = ex_to<basic>(basis).normal(repl, rev_lookup, modifier); |
| | for (int imod = nmod; imod < modifier.nops(); ++imod) |
| | n_basis = n_basis.subs(modifier.op(imod), subs_options::no_patter |
| | n); |
| | |
| | nmod = modifier.nops(); |
| | ex n_exponent = ex_to<basic>(exponent).normal(repl, rev_lookup, modifier) |
| | ; |
| | for (int imod = nmod; imod < modifier.nops(); ++imod) |
| | n_exponent = n_exponent.subs(modifier.op(imod), subs_options::no_ |
| | pattern); |
| n_exponent = n_exponent.op(0) / n_exponent.op(1); | | n_exponent = n_exponent.op(0) / n_exponent.op(1); |
| | |
| if (n_exponent.info(info_flags::integer)) { | | if (n_exponent.info(info_flags::integer)) { |
| | |
| if (n_exponent.info(info_flags::positive)) { | | if (n_exponent.info(info_flags::positive)) { |
| | |
| // (a/b)^n -> {a^n, b^n} | | // (a/b)^n -> {a^n, b^n} |
| return dynallocate<lst>({pow(n_basis.op(0), n_exponent),
pow(n_basis.op(1), n_exponent)}); | | return dynallocate<lst>({pow(n_basis.op(0), n_exponent),
pow(n_basis.op(1), n_exponent)}); |
| | |
| } else if (n_exponent.info(info_flags::negative)) { | | } else if (n_exponent.info(info_flags::negative)) { |
| | |
| // (a/b)^-n -> {b^n, a^n} | | // (a/b)^-n -> {b^n, a^n} |
| return dynallocate<lst>({pow(n_basis.op(1), -n_exponent),
pow(n_basis.op(0), -n_exponent)}); | | return dynallocate<lst>({pow(n_basis.op(1), -n_exponent),
pow(n_basis.op(0), -n_exponent)}); |
| } | | } |
| | |
| } else { | | } else { |
| | |
| if (n_exponent.info(info_flags::positive)) { | | if (n_exponent.info(info_flags::positive)) { |
| | |
| // (a/b)^x -> {sym((a/b)^x), 1} | | // (a/b)^x -> {sym((a/b)^x), 1} |
|
| return dynallocate<lst>({replace_with_symbol(pow(n_basis.
op(0) / n_basis.op(1), n_exponent), repl, rev_lookup), _ex1}); | | return dynallocate<lst>({replace_with_symbol(pow(n_basis.
op(0) / n_basis.op(1), n_exponent), repl, rev_lookup, modifier), _ex1}); |
| | |
| } else if (n_exponent.info(info_flags::negative)) { | | } else if (n_exponent.info(info_flags::negative)) { |
| | |
| if (n_basis.op(1).is_equal(_ex1)) { | | if (n_basis.op(1).is_equal(_ex1)) { |
| | |
| // a^-x -> {1, sym(a^x)} | | // a^-x -> {1, sym(a^x)} |
|
| return dynallocate<lst>({_ex1, replace_with_symbo
l(pow(n_basis.op(0), -n_exponent), repl, rev_lookup)}); | | return dynallocate<lst>({_ex1, replace_with_symbo
l(pow(n_basis.op(0), -n_exponent), repl, rev_lookup, modifier)}); |
| | |
| } else { | | } else { |
| | |
| // (a/b)^-x -> {sym((b/a)^x), 1} | | // (a/b)^-x -> {sym((b/a)^x), 1} |
|
| return dynallocate<lst>({replace_with_symbol(pow(
n_basis.op(1) / n_basis.op(0), -n_exponent), repl, rev_lookup), _ex1}); | | return dynallocate<lst>({replace_with_symbol(pow(
n_basis.op(1) / n_basis.op(0), -n_exponent), repl, rev_lookup, modifier), _ex1})
; |
| } | | } |
| } | | } |
| } | | } |
| | |
| // (a/b)^x -> {sym((a/b)^x, 1} | | // (a/b)^x -> {sym((a/b)^x, 1} |
|
| return dynallocate<lst>({replace_with_symbol(pow(n_basis.op(0) / n_basis.
op(1), n_exponent), repl, rev_lookup), _ex1}); | | return dynallocate<lst>({replace_with_symbol(pow(n_basis.op(0) / n_basis.
op(1), n_exponent), repl, rev_lookup, modifier), _ex1}); |
| } | | } |
| | |
| /** Implementation of ex::normal() for pseries. It normalizes each coefficient | | /** Implementation of ex::normal() for pseries. It normalizes each coefficient |
| * and replaces the series by a temporary symbol. | | * and replaces the series by a temporary symbol. |
| * @see ex::normal */ | | * @see ex::normal */ |
|
| ex pseries::normal(exmap & repl, exmap & rev_lookup) const | | ex pseries::normal(exmap & repl, exmap & rev_lookup, lst & modifier) const |
| { | | { |
| epvector newseq; | | epvector newseq; |
| for (auto & it : seq) { | | for (auto & it : seq) { |
| ex restexp = it.rest.normal(); | | ex restexp = it.rest.normal(); |
| if (!restexp.is_zero()) | | if (!restexp.is_zero()) |
| newseq.push_back(expair(restexp, it.coeff)); | | newseq.push_back(expair(restexp, it.coeff)); |
| } | | } |
| ex n = pseries(relational(var,point), std::move(newseq)); | | ex n = pseries(relational(var,point), std::move(newseq)); |
|
| return dynallocate<lst>({replace_with_symbol(n, repl, rev_lookup), _ex1})
; | | return dynallocate<lst>({replace_with_symbol(n, repl, rev_lookup, modifie
r), _ex1}); |
| } | | } |
| | |
| /** Normalization of rational functions. | | /** Normalization of rational functions. |
| * This function converts an expression to its normal form | | * This function converts an expression to its normal form |
| * "numerator/denominator", where numerator and denominator are (relatively | | * "numerator/denominator", where numerator and denominator are (relatively |
| * prime) polynomials. Any subexpressions which are not rational functions | | * prime) polynomials. Any subexpressions which are not rational functions |
| * (like non-rational numbers, non-integer powers or functions like sin(), | | * (like non-rational numbers, non-integer powers or functions like sin(), |
| * cos() etc.) are replaced by temporary symbols which are re-substituted by | | * cos() etc.) are replaced by temporary symbols which are re-substituted by |
| * the (normalized) subexpressions before normal() returns (this way, any | | * the (normalized) subexpressions before normal() returns (this way, any |
| * expression can be treated as a rational function). normal() is applied | | * expression can be treated as a rational function). normal() is applied |
| * recursively to arguments of functions etc. | | * recursively to arguments of functions etc. |
| * | | * |
| * @return normalized expression */ | | * @return normalized expression */ |
| ex ex::normal() const | | ex ex::normal() const |
| { | | { |
| exmap repl, rev_lookup; | | exmap repl, rev_lookup; |
|
| | lst modifier; |
| | |
|
| ex e = bp->normal(repl, rev_lookup); | | ex e = bp->normal(repl, rev_lookup, modifier); |
| GINAC_ASSERT(is_a<lst>(e)); | | GINAC_ASSERT(is_a<lst>(e)); |
| | |
| // Re-insert replaced symbols | | // Re-insert replaced symbols |
|
| if (!repl.empty()) | | if (!repl.empty()) { |
| | for(int i=0; i < modifier.nops(); ++i) |
| | e = e.subs(modifier.op(i), subs_options::no_pattern); |
| e = e.subs(repl, subs_options::no_pattern); | | e = e.subs(repl, subs_options::no_pattern); |
|
| | } |
| | |
| // Convert {numerator, denominator} form back to fraction | | // Convert {numerator, denominator} form back to fraction |
| return e.op(0) / e.op(1); | | return e.op(0) / e.op(1); |
| } | | } |
| | |
| /** Get numerator of an expression. If the expression is not of the normal | | /** Get numerator of an expression. If the expression is not of the normal |
| * form "numerator/denominator", it is first converted to this form and | | * form "numerator/denominator", it is first converted to this form and |
| * then the numerator is returned. | | * then the numerator is returned. |
| * | | * |
| * @see ex::normal | | * @see ex::normal |
| * @return numerator */ | | * @return numerator */ |
| ex ex::numer() const | | ex ex::numer() const |
| { | | { |
| exmap repl, rev_lookup; | | exmap repl, rev_lookup; |
|
| | lst modifier; |
| | |
|
| ex e = bp->normal(repl, rev_lookup); | | ex e = bp->normal(repl, rev_lookup, modifier); |
| GINAC_ASSERT(is_a<lst>(e)); | | GINAC_ASSERT(is_a<lst>(e)); |
| | |
| // Re-insert replaced symbols | | // Re-insert replaced symbols |
| if (repl.empty()) | | if (repl.empty()) |
| return e.op(0); | | return e.op(0); |
|
| else | | else { |
| | for(int i=0; i < modifier.nops(); ++i) |
| | e = e.subs(modifier.op(i), subs_options::no_pattern); |
| | |
| return e.op(0).subs(repl, subs_options::no_pattern); | | return e.op(0).subs(repl, subs_options::no_pattern); |
|
| | } |
| } | | } |
| | |
| /** Get denominator of an expression. If the expression is not of the normal | | /** Get denominator of an expression. If the expression is not of the normal |
| * form "numerator/denominator", it is first converted to this form and | | * form "numerator/denominator", it is first converted to this form and |
| * then the denominator is returned. | | * then the denominator is returned. |
| * | | * |
| * @see ex::normal | | * @see ex::normal |
| * @return denominator */ | | * @return denominator */ |
| ex ex::denom() const | | ex ex::denom() const |
| { | | { |
| exmap repl, rev_lookup; | | exmap repl, rev_lookup; |
|
| | lst modifier; |
| | |
|
| ex e = bp->normal(repl, rev_lookup); | | ex e = bp->normal(repl, rev_lookup, modifier); |
| GINAC_ASSERT(is_a<lst>(e)); | | GINAC_ASSERT(is_a<lst>(e)); |
| | |
| // Re-insert replaced symbols | | // Re-insert replaced symbols |
| if (repl.empty()) | | if (repl.empty()) |
| return e.op(1); | | return e.op(1); |
|
| else | | else { |
| | for(int i=0; i < modifier.nops(); ++i) |
| | e = e.subs(modifier.op(i), subs_options::no_pattern); |
| | |
| return e.op(1).subs(repl, subs_options::no_pattern); | | return e.op(1).subs(repl, subs_options::no_pattern); |
|
| | } |
| } | | } |
| | |
| /** Get numerator and denominator of an expression. If the expression is not | | /** Get numerator and denominator of an expression. If the expression is not |
| * of the normal form "numerator/denominator", it is first converted to this | | * of the normal form "numerator/denominator", it is first converted to this |
| * form and then a list [numerator, denominator] is returned. | | * form and then a list [numerator, denominator] is returned. |
| * | | * |
| * @see ex::normal | | * @see ex::normal |
| * @return a list [numerator, denominator] */ | | * @return a list [numerator, denominator] */ |
| ex ex::numer_denom() const | | ex ex::numer_denom() const |
| { | | { |
| exmap repl, rev_lookup; | | exmap repl, rev_lookup; |
|
| | lst modifier; |
| | |
|
| ex e = bp->normal(repl, rev_lookup); | | ex e = bp->normal(repl, rev_lookup, modifier); |
| GINAC_ASSERT(is_a<lst>(e)); | | GINAC_ASSERT(is_a<lst>(e)); |
| | |
| // Re-insert replaced symbols | | // Re-insert replaced symbols |
| if (repl.empty()) | | if (repl.empty()) |
| return e; | | return e; |
|
| else | | else { |
| | for(int i=0; i < modifier.nops(); ++i) |
| | e = e.subs(modifier.op(i), subs_options::no_pattern); |
| | |
| return e.subs(repl, subs_options::no_pattern); | | return e.subs(repl, subs_options::no_pattern); |
|
| | } |
| } | | } |
| | |
| /** Rationalization of non-rational functions. | | /** Rationalization of non-rational functions. |
| * This function converts a general expression to a rational function | | * This function converts a general expression to a rational function |
| * by replacing all non-rational subexpressions (like non-rational numbers, | | * by replacing all non-rational subexpressions (like non-rational numbers, |
| * non-integer powers or functions like sin(), cos() etc.) to temporary | | * non-integer powers or functions like sin(), cos() etc.) to temporary |
| * symbols. This makes it possible to use functions like gcd() and divide() | | * symbols. This makes it possible to use functions like gcd() and divide() |
| * on non-rational functions by applying to_rational() on the arguments, | | * on non-rational functions by applying to_rational() on the arguments, |
| * calling the desired function and re-substituting the temporary symbols | | * calling the desired function and re-substituting the temporary symbols |
| * in the result. To make the last step possible, all temporary symbols and | | * in the result. To make the last step possible, all temporary symbols and |
| | | |
| End of changes. 50 change blocks. |
|---|
| 44 lines changed or deleted | | 282 lines changed or added |
|---|
"Fossies" - the Fresh Open Source Software Archive 