X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fpower.cpp;h=491f2d332a86a181184d9abad1b2a97bea796de8;hp=72806768dea4564b754b912f3734651331559e22;hb=345f0c93051c01cdb5ab879b57f70118978fa9d0;hpb=acae7ab5a4dc94d1f54ba794f32f5764cdb4d704 diff --git a/ginac/power.cpp b/ginac/power.cpp index 72806768..491f2d33 100644 --- a/ginac/power.cpp +++ b/ginac/power.cpp @@ -139,7 +139,7 @@ void power::do_print_latex(const print_latex & c, unsigned level) const static void print_sym_pow(const print_context & c, const symbol &x, int exp) { // Optimal output of integer powers of symbols to aid compiler CSE. - // C.f. ISO/IEC 14882:1998, section 1.9 [intro execution], paragraph 15 + // C.f. ISO/IEC 14882:2011, section 1.9 [intro execution], paragraph 15 // to learn why such a parenthesation is really necessary. if (exp == 1) { x.print(c); @@ -524,8 +524,8 @@ ex power::eval(int level) const addp->setflag(status_flags::dynallocated); addp->clearflag(status_flags::hash_calculated); addp->overall_coeff = ex_to(addp->overall_coeff).div_dyn(icont); - for (epvector::iterator i = addp->seq.begin(); i != addp->seq.end(); ++i) - i->coeff = ex_to(i->coeff).div_dyn(icont); + for (auto & i : addp->seq) + i.coeff = ex_to(i.coeff).div_dyn(icont); const numeric c = icont.power(*num_exponent); if (likely(c != *_num1_p)) @@ -653,12 +653,12 @@ ex power::subs(const exmap & m, unsigned options) const if (!(options & subs_options::algebraic)) return subs_one_level(m, options); - for (exmap::const_iterator it = m.begin(); it != m.end(); ++it) { + for (auto & it : m) { int nummatches = std::numeric_limits::max(); exmap repls; - if (tryfactsubs(*this, it->first, nummatches, repls)) { - ex anum = it->second.subs(repls, subs_options::no_pattern); - ex aden = it->first.subs(repls, subs_options::no_pattern); + if (tryfactsubs(*this, it.first, nummatches, repls)) { + ex anum = it.second.subs(repls, subs_options::no_pattern); + ex aden = it.first.subs(repls, subs_options::no_pattern); ex result = (*this)*power(anum/aden, nummatches); return (ex_to(result)).subs_one_level(m, options); } @@ -695,51 +695,71 @@ ex power::conjugate() const ex power::real_part() const { + // basis == a+I*b, exponent == c+I*d + const ex a = basis.real_part(); + const ex c = exponent.real_part(); + if (basis.is_equal(a) && exponent.is_equal(c)) { + // Re(a^c) + return *this; + } + + const ex b = basis.imag_part(); if (exponent.info(info_flags::integer)) { - ex basis_real = basis.real_part(); - if (basis_real == basis) - return *this; - realsymbol a("a"),b("b"); - ex result; - if (exponent.info(info_flags::posint)) - result = power(a+I*b,exponent); - else - result = power(a/(a*a+b*b)-I*b/(a*a+b*b),-exponent); - result = result.expand(); - result = result.real_part(); - result = result.subs(lst( a==basis_real, b==basis.imag_part() )); + // Re((a+I*b)^c) w/ c ∈ ℤ + long N = ex_to(c).to_long(); + // Use real terms in Binomial expansion to construct + // Re(expand(power(a+I*b, N))). + long NN = N > 0 ? N : -N; + ex numer = N > 0 ? _ex1 : power(power(a,2) + power(b,2), NN); + ex result = 0; + for (long n = 0; n <= NN; n += 2) { + ex term = binomial(NN, n) * power(a, NN-n) * power(b, n) / numer; + if (n % 4 == 0) { + result += term; // sign: I^n w/ n == 4*m + } else { + result -= term; // sign: I^n w/ n == 4*m+2 + } + } return result; } - - ex a = basis.real_part(); - ex b = basis.imag_part(); - ex c = exponent.real_part(); - ex d = exponent.imag_part(); + + // Re((a+I*b)^(c+I*d)) + const ex d = exponent.imag_part(); return power(abs(basis),c)*exp(-d*atan2(b,a))*cos(c*atan2(b,a)+d*log(abs(basis))); } ex power::imag_part() const { + const ex a = basis.real_part(); + const ex c = exponent.real_part(); + if (basis.is_equal(a) && exponent.is_equal(c)) { + // Im(a^c) + return 0; + } + + const ex b = basis.imag_part(); if (exponent.info(info_flags::integer)) { - ex basis_real = basis.real_part(); - if (basis_real == basis) - return 0; - realsymbol a("a"),b("b"); - ex result; - if (exponent.info(info_flags::posint)) - result = power(a+I*b,exponent); - else - result = power(a/(a*a+b*b)-I*b/(a*a+b*b),-exponent); - result = result.expand(); - result = result.imag_part(); - result = result.subs(lst( a==basis_real, b==basis.imag_part() )); + // Im((a+I*b)^c) w/ c ∈ ℤ + long N = ex_to(c).to_long(); + // Use imaginary terms in Binomial expansion to construct + // Im(expand(power(a+I*b, N))). + long p = N > 0 ? 1 : 3; // modulus for positive sign + long NN = N > 0 ? N : -N; + ex numer = N > 0 ? _ex1 : power(power(a,2) + power(b,2), NN); + ex result = 0; + for (long n = 1; n <= NN; n += 2) { + ex term = binomial(NN, n) * power(a, NN-n) * power(b, n) / numer; + if (n % 4 == p) { + result += term; // sign: I^n w/ n == 4*m+p + } else { + result -= term; // sign: I^n w/ n == 4*m+2+p + } + } return result; } - - ex a=basis.real_part(); - ex b=basis.imag_part(); - ex c=exponent.real_part(); - ex d=exponent.imag_part(); + + // Im((a+I*b)^(c+I*d)) + const ex d = exponent.imag_part(); return power(abs(basis),c)*exp(-d*atan2(b,a))*sin(c*atan2(b,a)+d*log(abs(basis))); } @@ -802,21 +822,18 @@ ex power::expand(unsigned options) const epvector powseq; prodseq.reserve(m.seq.size() + 1); powseq.reserve(m.seq.size() + 1); - epvector::const_iterator last = m.seq.end(); - epvector::const_iterator cit = m.seq.begin(); bool possign = true; // search for positive/negative factors - while (cit!=last) { - ex e=m.recombine_pair_to_ex(*cit); + for (auto & cit : m.seq) { + ex e=m.recombine_pair_to_ex(cit); if (e.info(info_flags::positive)) prodseq.push_back(pow(e, exponent).expand(options)); else if (e.info(info_flags::negative)) { prodseq.push_back(pow(-e, exponent).expand(options)); possign = !possign; } else - powseq.push_back(*cit); - ++cit; + powseq.push_back(cit); } // take care on the numeric coefficient @@ -847,11 +864,8 @@ ex power::expand(unsigned options) const const add &a = ex_to(expanded_exponent); exvector distrseq; distrseq.reserve(a.seq.size() + 1); - epvector::const_iterator last = a.seq.end(); - epvector::const_iterator cit = a.seq.begin(); - while (cit!=last) { - distrseq.push_back(power(expanded_basis, a.recombine_pair_to_ex(*cit))); - ++cit; + for (auto & cit : a.seq) { + distrseq.push_back(power(expanded_basis, a.recombine_pair_to_ex(cit))); } // Make sure that e.g. (x+y)^(2+a) expands the (x+y)^2 factor @@ -1006,8 +1020,8 @@ private: element *head, *i, *after_i; // NB: Partition must be sorted in non-decreasing order. explicit coolmulti(const std::vector& partition) + : head(nullptr), i(nullptr), after_i(nullptr) { - head = nullptr; for (unsigned n = 0; n < partition.size(); ++n) { head = new element(partition[n], head); if (n <= 1) @@ -1077,14 +1091,12 @@ public: * where n = p1+p2+...+pk, i.e. p is a partition of n. */ const numeric -multinomial_coefficient(const std::vector p) +multinomial_coefficient(const std::vector & p) { numeric n = 0, d = 1; - std::vector::const_iterator it = p.begin(), itend = p.end(); - while (it != itend) { - n += numeric(*it); - d *= factorial(numeric(*it)); - ++it; + for (auto & it : p) { + n += numeric(it); + d *= factorial(numeric(it)); } return factorial(numeric(n)) / d; } @@ -1195,7 +1207,6 @@ ex power::expand_add(const add & a, int n, unsigned options) const numeric factor = coeff; for (unsigned i = 0; i < exponent.size(); ++i) { const ex & r = a.seq[i].rest; - const ex & c = a.seq[i].coeff; GINAC_ASSERT(!is_exactly_a(r)); GINAC_ASSERT(!is_exactly_a(r) || !is_exactly_a(ex_to(r).exponent) || @@ -1203,15 +1214,19 @@ ex power::expand_add(const add & a, int n, unsigned options) const !is_exactly_a(ex_to(r).basis) || !is_exactly_a(ex_to(r).basis) || !is_exactly_a(ex_to(r).basis)); + GINAC_ASSERT(is_exactly_a(a.seq[i].coeff)); + const numeric & c = ex_to(a.seq[i].coeff); if (exponent[i] == 0) { // optimize away } else if (exponent[i] == 1) { // optimized term.push_back(r); - factor = factor.mul(ex_to(c)); + if (c != *_num1_p) + factor = factor.mul(c); } else { // general case exponent[i] > 1 term.push_back((new power(r, exponent[i]))->setflag(status_flags::dynallocated)); - factor = factor.mul(ex_to(c).power(exponent[i])); + if (c != *_num1_p) + factor = factor.mul(c.power(exponent[i])); } } result.push_back(a.combine_ex_with_coeff_to_pair(mul(term).expand(options), factor)); @@ -1327,17 +1342,14 @@ ex power::expand_mul(const mul & m, const numeric & n, unsigned options, bool fr distrseq.reserve(m.seq.size()); bool need_reexpand = false; - epvector::const_iterator last = m.seq.end(); - epvector::const_iterator cit = m.seq.begin(); - while (cit!=last) { - expair p = m.combine_pair_with_coeff_to_pair(*cit, n); - if (from_expand && is_exactly_a(cit->rest) && ex_to(p.coeff).is_pos_integer()) { + for (auto & cit : m.seq) { + expair p = m.combine_pair_with_coeff_to_pair(cit, n); + if (from_expand && is_exactly_a(cit.rest) && ex_to(p.coeff).is_pos_integer()) { // this happens when e.g. (a+b)^(1/2) gets squared and // the resulting product needs to be reexpanded need_reexpand = true; } distrseq.push_back(p); - ++cit; } const mul & result = static_cast((new mul(distrseq, ex_to(m.overall_coeff).power_dyn(n)))->setflag(status_flags::dynallocated));