GiNaC  1.8.3
power.cpp
Go to the documentation of this file.
1 
5 /*
6  * GiNaC Copyright (C) 1999-2022 Johannes Gutenberg University Mainz, Germany
7  *
8  * This program is free software; you can redistribute it and/or modify
9  * it under the terms of the GNU General Public License as published by
10  * the Free Software Foundation; either version 2 of the License, or
11  * (at your option) any later version.
12  *
13  * This program is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  * GNU General Public License for more details.
17  *
18  * You should have received a copy of the GNU General Public License
19  * along with this program; if not, write to the Free Software
20  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
21  */
22 
23 #include "power.h"
24 #include "expairseq.h"
25 #include "add.h"
26 #include "mul.h"
27 #include "ncmul.h"
28 #include "numeric.h"
29 #include "constant.h"
30 #include "operators.h"
31 #include "inifcns.h" // for log() in power::derivative()
32 #include "matrix.h"
33 #include "indexed.h"
34 #include "symbol.h"
35 #include "lst.h"
36 #include "archive.h"
37 #include "utils.h"
38 #include "relational.h"
39 #include "compiler.h"
40 
41 #include <iostream>
42 #include <limits>
43 #include <stdexcept>
44 #include <vector>
45 #include <algorithm>
46 
47 namespace GiNaC {
48 
51  print_func<print_latex>(&power::do_print_latex).
52  print_func<print_csrc>(&power::do_print_csrc).
53  print_func<print_python>(&power::do_print_python).
54  print_func<print_python_repr>(&power::do_print_python_repr).
55  print_func<print_csrc_cl_N>(&power::do_print_csrc_cl_N))
56 
57 // default constructor
60 
61 power::power() { }
62 
64 // other constructors
66 
67 // all inlined
68 
70 // archiving
72 
73 void power::read_archive(const archive_node &n, lst &sym_lst)
74 {
75  inherited::read_archive(n, sym_lst);
76  n.find_ex("basis", basis, sym_lst);
77  n.find_ex("exponent", exponent, sym_lst);
78 }
79 
81 {
82  inherited::archive(n);
83  n.add_ex("basis", basis);
84  n.add_ex("exponent", exponent);
85 }
86 
88 // functions overriding virtual functions from base classes
90 
91 // public
92 
93 void power::print_power(const print_context & c, const char *powersymbol, const char *openbrace, const char *closebrace, unsigned level) const
94 {
95  // Ordinary output of powers using '^' or '**'
96  if (precedence() <= level)
97  c.s << openbrace << '(';
98  basis.print(c, precedence());
99  c.s << powersymbol;
100  c.s << openbrace;
102  c.s << closebrace;
103  if (precedence() <= level)
104  c.s << ')' << closebrace;
105 }
106 
107 void power::do_print_dflt(const print_dflt & c, unsigned level) const
108 {
109  if (exponent.is_equal(_ex1_2)) {
110 
111  // Square roots are printed in a special way
112  c.s << "sqrt(";
113  basis.print(c);
114  c.s << ')';
115 
116  } else
117  print_power(c, "^", "", "", level);
118 }
119 
120 void power::do_print_latex(const print_latex & c, unsigned level) const
121 {
122  if (is_exactly_a<numeric>(exponent) && ex_to<numeric>(exponent).is_negative()) {
123 
124  // Powers with negative numeric exponents are printed as fractions
125  c.s << "\\frac{1}{";
126  power(basis, -exponent).eval().print(c);
127  c.s << '}';
128 
129  } else if (exponent.is_equal(_ex1_2)) {
130 
131  // Square roots are printed in a special way
132  c.s << "\\sqrt{";
133  basis.print(c);
134  c.s << '}';
135 
136  } else
137  print_power(c, "^", "{", "}", level);
138 }
139 
140 static void print_sym_pow(const print_context & c, const symbol &x, int exp)
141 {
142  // Optimal output of integer powers of symbols to aid compiler CSE.
143  // C.f. ISO/IEC 14882:2011, section 1.9 [intro.execution], paragraph 15
144  // to learn why such a parenthesation is really necessary.
145  if (exp == 1) {
146  x.print(c);
147  } else if (exp == 2) {
148  x.print(c);
149  c.s << "*";
150  x.print(c);
151  } else if (exp & 1) {
152  x.print(c);
153  c.s << "*";
154  print_sym_pow(c, x, exp-1);
155  } else {
156  c.s << "(";
157  print_sym_pow(c, x, exp >> 1);
158  c.s << ")*(";
159  print_sym_pow(c, x, exp >> 1);
160  c.s << ")";
161  }
162 }
163 
164 void power::do_print_csrc_cl_N(const print_csrc_cl_N& c, unsigned level) const
165 {
166  if (exponent.is_equal(_ex_1)) {
167  c.s << "recip(";
168  basis.print(c);
169  c.s << ')';
170  return;
171  }
172  c.s << "expt(";
173  basis.print(c);
174  c.s << ", ";
175  exponent.print(c);
176  c.s << ')';
177 }
178 
179 void power::do_print_csrc(const print_csrc & c, unsigned level) const
180 {
181  // Integer powers of symbols are printed in a special, optimized way
183  (is_a<symbol>(basis) || is_a<constant>(basis))) {
184  int exp = ex_to<numeric>(exponent).to_int();
185  if (exp > 0)
186  c.s << '(';
187  else {
188  exp = -exp;
189  c.s << "1.0/(";
190  }
191  print_sym_pow(c, ex_to<symbol>(basis), exp);
192  c.s << ')';
193 
194  // <expr>^-1 is printed as "1.0/<expr>" or with the recip() function of CLN
195  } else if (exponent.is_equal(_ex_1)) {
196  c.s << "1.0/(";
197  basis.print(c);
198  c.s << ')';
199 
200  // Otherwise, use the pow() function
201  } else {
202  c.s << "pow(";
203  basis.print(c);
204  c.s << ',';
205  exponent.print(c);
206  c.s << ')';
207  }
208 }
209 
210 void power::do_print_python(const print_python & c, unsigned level) const
211 {
212  print_power(c, "**", "", "", level);
213 }
214 
215 void power::do_print_python_repr(const print_python_repr & c, unsigned level) const
216 {
217  c.s << class_name() << '(';
218  basis.print(c);
219  c.s << ',';
220  exponent.print(c);
221  c.s << ')';
222 }
223 
224 bool power::info(unsigned inf) const
225 {
226  switch (inf) {
232  return basis.info(inf) && exponent.info(info_flags::nonnegint);
234  return basis.info(inf) && exponent.info(info_flags::integer);
235  case info_flags::real:
236  return basis.info(inf) && exponent.info(info_flags::integer);
238  return (flags & status_flags::expanded);
246  return true;
248  return false;
249  else if (basis.info(info_flags::has_indices)) {
252  return true;
253  } else {
256  return false;
257  }
258  }
259  }
260  return inherited::info(inf);
261 }
262 
263 size_t power::nops() const
264 {
265  return 2;
266 }
267 
268 ex power::op(size_t i) const
269 {
270  GINAC_ASSERT(i<2);
271 
272  return i==0 ? basis : exponent;
273 }
274 
276 {
277  const ex &mapped_basis = f(basis);
278  const ex &mapped_exponent = f(exponent);
279 
280  if (!are_ex_trivially_equal(basis, mapped_basis)
281  || !are_ex_trivially_equal(exponent, mapped_exponent))
282  return dynallocate<power>(mapped_basis, mapped_exponent);
283  else
284  return *this;
285 }
286 
287 bool power::is_polynomial(const ex & var) const
288 {
289  if (basis.is_polynomial(var)) {
290  if (basis.has(var))
291  // basis is non-constant polynomial in var
293  else
294  // basis is constant in var
295  return !exponent.has(var);
296  }
297  // basis is a non-polynomial function of var
298  return false;
299 }
300 
301 int power::degree(const ex & s) const
302 {
303  if (is_equal(ex_to<basic>(s)))
304  return 1;
305  else if (is_exactly_a<numeric>(exponent) && ex_to<numeric>(exponent).is_integer()) {
306  if (basis.is_equal(s))
307  return ex_to<numeric>(exponent).to_int();
308  else
309  return basis.degree(s) * ex_to<numeric>(exponent).to_int();
310  } else if (basis.has(s))
311  throw(std::runtime_error("power::degree(): undefined degree because of non-integer exponent"));
312  else
313  return 0;
314 }
315 
316 int power::ldegree(const ex & s) const
317 {
318  if (is_equal(ex_to<basic>(s)))
319  return 1;
320  else if (is_exactly_a<numeric>(exponent) && ex_to<numeric>(exponent).is_integer()) {
321  if (basis.is_equal(s))
322  return ex_to<numeric>(exponent).to_int();
323  else
324  return basis.ldegree(s) * ex_to<numeric>(exponent).to_int();
325  } else if (basis.has(s))
326  throw(std::runtime_error("power::ldegree(): undefined degree because of non-integer exponent"));
327  else
328  return 0;
329 }
330 
331 ex power::coeff(const ex & s, int n) const
332 {
333  if (is_equal(ex_to<basic>(s)))
334  return n==1 ? _ex1 : _ex0;
335  else if (!basis.is_equal(s)) {
336  // basis not equal to s
337  if (n == 0)
338  return *this;
339  else
340  return _ex0;
341  } else {
342  // basis equal to s
343  if (is_exactly_a<numeric>(exponent) && ex_to<numeric>(exponent).is_integer()) {
344  // integer exponent
345  int int_exp = ex_to<numeric>(exponent).to_int();
346  if (n == int_exp)
347  return _ex1;
348  else
349  return _ex0;
350  } else {
351  // non-integer exponents are treated as zero
352  if (n == 0)
353  return *this;
354  else
355  return _ex0;
356  }
357  }
358 }
359 
375 {
377  return *this;
378 
379  const numeric *num_basis = nullptr;
380  const numeric *num_exponent = nullptr;
381 
382  if (is_exactly_a<numeric>(basis)) {
383  num_basis = &ex_to<numeric>(basis);
384  }
385  if (is_exactly_a<numeric>(exponent)) {
386  num_exponent = &ex_to<numeric>(exponent);
387  }
388 
389  // ^(x,0) -> 1 (0^0 also handled here)
390  if (exponent.is_zero()) {
391  if (basis.is_zero())
392  throw (std::domain_error("power::eval(): pow(0,0) is undefined"));
393  else
394  return _ex1;
395  }
396 
397  // ^(x,1) -> x
398  if (exponent.is_equal(_ex1))
399  return basis;
400 
401  // ^(0,c1) -> 0 or exception (depending on real value of c1)
402  if (basis.is_zero() && num_exponent) {
403  if ((num_exponent->real()).is_zero())
404  throw (std::domain_error("power::eval(): pow(0,I) is undefined"));
405  else if ((num_exponent->real()).is_negative())
406  throw (pole_error("power::eval(): division by zero",1));
407  else
408  return _ex0;
409  }
410 
411  // ^(1,x) -> 1
412  if (basis.is_equal(_ex1))
413  return _ex1;
414 
415  // power of a function calculated by separate rules defined for this function
416  if (is_exactly_a<function>(basis))
417  return ex_to<function>(basis).power(exponent);
418 
419  // Turn (x^c)^d into x^(c*d) in the case that x is positive and c is real.
420  if (is_exactly_a<power>(basis) && basis.op(0).info(info_flags::positive) && basis.op(1).info(info_flags::real))
421  return dynallocate<power>(basis.op(0), basis.op(1) * exponent);
422 
423  if ( num_exponent ) {
424 
425  // ^(c1,c2) -> c1^c2 (c1, c2 numeric(),
426  // except if c1,c2 are rational, but c1^c2 is not)
427  if ( num_basis ) {
428  const bool basis_is_crational = num_basis->is_crational();
429  const bool exponent_is_crational = num_exponent->is_crational();
430  if (!basis_is_crational || !exponent_is_crational) {
431  // return a plain float
432  return dynallocate<numeric>(num_basis->power(*num_exponent));
433  }
434 
435  const numeric res = num_basis->power(*num_exponent);
436  if (res.is_crational()) {
437  return res;
438  }
439  GINAC_ASSERT(!num_exponent->is_integer()); // has been handled by now
440 
441  // ^(c1,n/m) -> *(c1^q,c1^(n/m-q)), 0<(n/m-q)<1, q integer
442  if (basis_is_crational && exponent_is_crational
443  && num_exponent->is_real()
444  && !num_exponent->is_integer()) {
445  const numeric n = num_exponent->numer();
446  const numeric m = num_exponent->denom();
447  numeric r;
448  numeric q = iquo(n, m, r);
449  if (r.is_negative()) {
450  r += m;
451  --q;
452  }
453  if (q.is_zero()) { // the exponent was in the allowed range 0<(n/m)<1
454  if (num_basis->is_rational() && !num_basis->is_integer()) {
455  // try it for numerator and denominator separately, in order to
456  // partially simplify things like (5/8)^(1/3) -> 1/2*5^(1/3)
457  const numeric bnum = num_basis->numer();
458  const numeric bden = num_basis->denom();
459  const numeric res_bnum = bnum.power(*num_exponent);
460  const numeric res_bden = bden.power(*num_exponent);
461  if (res_bnum.is_integer())
462  return dynallocate<mul>(dynallocate<power>(bden,-*num_exponent),res_bnum).setflag(status_flags::evaluated);
463  if (res_bden.is_integer())
464  return dynallocate<mul>(dynallocate<power>(bnum,*num_exponent),res_bden.inverse()).setflag(status_flags::evaluated);
465  }
466  return this->hold();
467  } else {
468  // assemble resulting product, but allowing for a re-evaluation,
469  // because otherwise we'll end up with something like
470  // (7/8)^(4/3) -> 7/8*(1/2*7^(1/3))
471  // instead of 7/16*7^(1/3).
472  return pow(basis, r.div(m)) * pow(basis, q);
473  }
474  }
475  }
476 
477  // ^(^(x,c1),c2) -> ^(x,c1*c2)
478  // (c1, c2 numeric(), c2 integer or -1 < c1 <= 1 or (c1=-1 and c2>0),
479  // case c1==1 should not happen, see below!)
480  if (is_exactly_a<power>(basis)) {
481  const power & sub_power = ex_to<power>(basis);
482  const ex & sub_basis = sub_power.basis;
483  const ex & sub_exponent = sub_power.exponent;
484  if (is_exactly_a<numeric>(sub_exponent)) {
485  const numeric & num_sub_exponent = ex_to<numeric>(sub_exponent);
486  GINAC_ASSERT(num_sub_exponent!=numeric(1));
487  if (num_exponent->is_integer() || (abs(num_sub_exponent) - (*_num1_p)).is_negative() ||
488  (num_sub_exponent == *_num_1_p && num_exponent->is_positive())) {
489  return dynallocate<power>(sub_basis, num_sub_exponent.mul(*num_exponent));
490  }
491  }
492  }
493 
494  // ^(*(x,y,z),c1) -> *(x^c1,y^c1,z^c1) (c1 integer)
495  if (num_exponent->is_integer() && is_exactly_a<mul>(basis)) {
496  return expand_mul(ex_to<mul>(basis), *num_exponent, false);
497  }
498 
499  // (2*x + 6*y)^(-4) -> 1/16*(x + 3*y)^(-4)
500  if (num_exponent->is_integer() && is_exactly_a<add>(basis)) {
501  numeric icont = basis.integer_content();
502  const numeric lead_coeff =
503  ex_to<numeric>(ex_to<add>(basis).seq.begin()->coeff).div(icont);
504 
505  const bool canonicalizable = lead_coeff.is_integer();
506  const bool unit_normal = lead_coeff.is_pos_integer();
507  if (canonicalizable && (! unit_normal))
508  icont = icont.mul(*_num_1_p);
509 
510  if (canonicalizable && (icont != *_num1_p)) {
511  const add& addref = ex_to<add>(basis);
512  add & addp = dynallocate<add>(addref);
514  addp.overall_coeff = ex_to<numeric>(addp.overall_coeff).div_dyn(icont);
515  for (auto & i : addp.seq)
516  i.coeff = ex_to<numeric>(i.coeff).div_dyn(icont);
517 
518  const numeric c = icont.power(*num_exponent);
519  if (likely(c != *_num1_p))
520  return dynallocate<mul>(dynallocate<power>(addp, *num_exponent), c);
521  else
522  return dynallocate<power>(addp, *num_exponent);
523  }
524  }
525 
526  // ^(*(...,x;c1),c2) -> *(^(*(...,x;1),c2),c1^c2) (c1, c2 numeric(), c1>0)
527  // ^(*(...,x;c1),c2) -> *(^(*(...,x;-1),c2),(-c1)^c2) (c1, c2 numeric(), c1<0)
528  if (is_exactly_a<mul>(basis)) {
529  GINAC_ASSERT(!num_exponent->is_integer()); // should have been handled above
530  const mul & mulref = ex_to<mul>(basis);
531  if (!mulref.overall_coeff.is_equal(_ex1)) {
532  const numeric & num_coeff = ex_to<numeric>(mulref.overall_coeff);
533  if (num_coeff.is_real()) {
534  if (num_coeff.is_positive()) {
535  mul & mulp = dynallocate<mul>(mulref);
536  mulp.overall_coeff = _ex1;
538  return dynallocate<mul>(dynallocate<power>(mulp, exponent),
539  dynallocate<power>(num_coeff, *num_exponent));
540  } else {
541  GINAC_ASSERT(num_coeff.compare(*_num0_p)<0);
542  if (!num_coeff.is_equal(*_num_1_p)) {
543  mul & mulp = dynallocate<mul>(mulref);
544  mulp.overall_coeff = _ex_1;
546  return dynallocate<mul>(dynallocate<power>(mulp, exponent),
547  dynallocate<power>(abs(num_coeff), *num_exponent));
548  }
549  }
550  }
551  }
552  }
553 
554  // ^(nc,c1) -> ncmul(nc,nc,...) (c1 positive integer, unless nc is a matrix)
555  if (num_exponent->is_pos_integer() &&
557  !is_a<matrix>(basis)) {
558  return ncmul(exvector(num_exponent->to_int(), basis));
559  }
560  }
561 
562  return this->hold();
563 }
564 
566 {
567  ex ebasis = basis.evalf();
568  ex eexponent;
569 
570  if (!is_exactly_a<numeric>(exponent))
571  eexponent = exponent.evalf();
572  else
573  eexponent = exponent;
574 
575  return dynallocate<power>(ebasis, eexponent);
576 }
577 
579 {
580  const ex ebasis = basis.evalm();
581  const ex eexponent = exponent.evalm();
582  if (is_a<matrix>(ebasis)) {
583  if (is_exactly_a<numeric>(eexponent)) {
584  return dynallocate<matrix>(ex_to<matrix>(ebasis).pow(eexponent));
585  }
586  }
587  return dynallocate<power>(ebasis, eexponent);
588 }
589 
590 bool power::has(const ex & other, unsigned options) const
591 {
593  return basic::has(other, options);
594  if (!is_a<power>(other))
595  return basic::has(other, options);
597  !other.op(1).info(info_flags::integer))
598  return basic::has(other, options);
600  other.op(1).info(info_flags::posint) &&
601  ex_to<numeric>(exponent) > ex_to<numeric>(other.op(1)) &&
602  basis.match(other.op(0)))
603  return true;
605  other.op(1).info(info_flags::negint) &&
606  ex_to<numeric>(exponent) < ex_to<numeric>(other.op(1)) &&
607  basis.match(other.op(0)))
608  return true;
609  return basic::has(other, options);
610 }
611 
612 // from mul.cpp
613 extern bool tryfactsubs(const ex &, const ex &, int &, exmap&);
614 
615 ex power::subs(const exmap & m, unsigned options) const
616 {
617  const ex &subsed_basis = basis.subs(m, options);
618  const ex &subsed_exponent = exponent.subs(m, options);
619 
620  if (!are_ex_trivially_equal(basis, subsed_basis)
621  || !are_ex_trivially_equal(exponent, subsed_exponent))
622  return dynallocate<power>(subsed_basis, subsed_exponent);
623 
625  return subs_one_level(m, options);
626 
627  for (auto & it : m) {
628  int nummatches = std::numeric_limits<int>::max();
629  exmap repls;
630  if (tryfactsubs(*this, it.first, nummatches, repls)) {
631  ex anum = it.second.subs(repls, subs_options::no_pattern);
632  ex aden = it.first.subs(repls, subs_options::no_pattern);
633  ex result = (*this) * pow(anum/aden, nummatches);
634  return (ex_to<basic>(result)).subs_one_level(m, options);
635  }
636  }
637 
638  return subs_one_level(m, options);
639 }
640 
641 ex power::eval_ncmul(const exvector & v) const
642 {
643  return inherited::eval_ncmul(v);
644 }
645 
647 {
648  // conjugate(pow(x,y))==pow(conjugate(x),conjugate(y)) unless on the
649  // branch cut which runs along the negative real axis.
651  ex newexponent = exponent.conjugate();
652  if (are_ex_trivially_equal(exponent, newexponent)) {
653  return *this;
654  }
655  return dynallocate<power>(basis, newexponent);
656  }
658  ex newbasis = basis.conjugate();
659  if (are_ex_trivially_equal(basis, newbasis)) {
660  return *this;
661  }
662  return dynallocate<power>(newbasis, exponent);
663  }
664  return conjugate_function(*this).hold();
665 }
666 
668 {
669  // basis == a+I*b, exponent == c+I*d
670  const ex a = basis.real_part();
671  const ex c = exponent.real_part();
672  if (basis.is_equal(a) && exponent.is_equal(c) &&
674  // Re(a^c)
675  return *this;
676  }
677 
678  const ex b = basis.imag_part();
680  // Re((a+I*b)^c) w/ c ∈ ℤ
681  long N = ex_to<numeric>(c).to_long();
682  // Use real terms in Binomial expansion to construct
683  // Re(expand(pow(a+I*b, N))).
684  long NN = N > 0 ? N : -N;
685  ex numer = N > 0 ? _ex1 : pow(pow(a,2) + pow(b,2), NN);
686  ex result = 0;
687  for (long n = 0; n <= NN; n += 2) {
688  ex term = binomial(NN, n) * pow(a, NN-n) * pow(b, n) / numer;
689  if (n % 4 == 0) {
690  result += term; // sign: I^n w/ n == 4*m
691  } else {
692  result -= term; // sign: I^n w/ n == 4*m+2
693  }
694  }
695  return result;
696  }
697 
698  // Re((a+I*b)^(c+I*d))
699  const ex d = exponent.imag_part();
700  return pow(abs(basis),c) * exp(-d*atan2(b,a)) * cos(c*atan2(b,a)+d*log(abs(basis)));
701 }
702 
704 {
705  // basis == a+I*b, exponent == c+I*d
706  const ex a = basis.real_part();
707  const ex c = exponent.real_part();
708  if (basis.is_equal(a) && exponent.is_equal(c) &&
710  // Im(a^c)
711  return 0;
712  }
713 
714  const ex b = basis.imag_part();
716  // Im((a+I*b)^c) w/ c ∈ ℤ
717  long N = ex_to<numeric>(c).to_long();
718  // Use imaginary terms in Binomial expansion to construct
719  // Im(expand(pow(a+I*b, N))).
720  long p = N > 0 ? 1 : 3; // modulus for positive sign
721  long NN = N > 0 ? N : -N;
722  ex numer = N > 0 ? _ex1 : pow(pow(a,2) + pow(b,2), NN);
723  ex result = 0;
724  for (long n = 1; n <= NN; n += 2) {
725  ex term = binomial(NN, n) * pow(a, NN-n) * pow(b, n) / numer;
726  if (n % 4 == p) {
727  result += term; // sign: I^n w/ n == 4*m+p
728  } else {
729  result -= term; // sign: I^n w/ n == 4*m+2+p
730  }
731  }
732  return result;
733  }
734 
735  // Im((a+I*b)^(c+I*d))
736  const ex d = exponent.imag_part();
737  return pow(abs(basis),c) * exp(-d*atan2(b,a)) * sin(c*atan2(b,a)+d*log(abs(basis)));
738 }
739 
740 // protected
741 
744 ex power::derivative(const symbol & s) const
745 {
746  if (is_a<numeric>(exponent)) {
747  // D(b^r) = r * b^(r-1) * D(b) (faster than the formula below)
748  const epvector newseq = {expair(basis, exponent - _ex1), expair(basis.diff(s), _ex1)};
749  return dynallocate<mul>(std::move(newseq), exponent);
750  } else {
751  // D(b^e) = b^e * (D(e)*ln(b) + e*D(b)/b)
752  return *this * (exponent.diff(s)*log(basis) + exponent*basis.diff(s)*pow(basis, _ex_1));
753  }
754 }
755 
756 int power::compare_same_type(const basic & other) const
757 {
758  GINAC_ASSERT(is_exactly_a<power>(other));
759  const power &o = static_cast<const power &>(other);
760 
761  int cmpval = basis.compare(o.basis);
762  if (cmpval)
763  return cmpval;
764  else
765  return exponent.compare(o.exponent);
766 }
767 
768 unsigned power::return_type() const
769 {
770  return basis.return_type();
771 }
772 
774 {
775  return basis.return_type_tinfo();
776 }
777 
778 ex power::expand(unsigned options) const
779 {
780  if (is_a<symbol>(basis) && exponent.info(info_flags::integer)) {
781  // A special case worth optimizing.
783  return *this;
784  }
785 
786  // (x*p)^c -> x^c * p^c, if p>0
787  // makes sense before expanding the basis
788  if (is_exactly_a<mul>(basis) && !basis.info(info_flags::indefinite)) {
789  const mul &m = ex_to<mul>(basis);
790  exvector prodseq;
791  epvector powseq;
792  prodseq.reserve(m.seq.size() + 1);
793  powseq.reserve(m.seq.size() + 1);
794  bool possign = true;
795 
796  // search for positive/negative factors
797  for (auto & cit : m.seq) {
798  ex e=m.recombine_pair_to_ex(cit);
799  if (e.info(info_flags::positive))
800  prodseq.push_back(pow(e, exponent).expand(options));
801  else if (e.info(info_flags::negative)) {
802  prodseq.push_back(pow(-e, exponent).expand(options));
803  possign = !possign;
804  } else
805  powseq.push_back(cit);
806  }
807 
808  // take care on the numeric coefficient
809  ex coeff=(possign? _ex1 : _ex_1);
810  if (m.overall_coeff.info(info_flags::positive) && m.overall_coeff != _ex1)
811  prodseq.push_back(pow(m.overall_coeff, exponent));
812  else if (m.overall_coeff.info(info_flags::negative) && m.overall_coeff != _ex_1) {
813  prodseq.push_back(pow(-m.overall_coeff, exponent));
814  coeff = -coeff;
815  } else
816  coeff *= m.overall_coeff;
817 
818  // If positive/negative factors are found, then extract them.
819  // In either case we set a flag to avoid the second run on a part
820  // which does not have positive/negative terms.
821  if (prodseq.size() > 0) {
822  ex newbasis = dynallocate<mul>(std::move(powseq), coeff);
823  ex_to<basic>(newbasis).setflag(status_flags::purely_indefinite);
824  return dynallocate<mul>(std::move(prodseq)) * pow(newbasis, exponent);
825  } else
826  ex_to<basic>(basis).setflag(status_flags::purely_indefinite);
827  }
828 
829  const ex expanded_basis = basis.expand(options);
830  const ex expanded_exponent = exponent.expand(options);
831 
832  // x^(a+b) -> x^a * x^b
833  if (is_exactly_a<add>(expanded_exponent)) {
834  const add &a = ex_to<add>(expanded_exponent);
835  exvector distrseq;
836  distrseq.reserve(a.seq.size() + 1);
837  for (auto & cit : a.seq) {
838  distrseq.push_back(pow(expanded_basis, a.recombine_pair_to_ex(cit)));
839  }
840 
841  // Make sure that e.g. (x+y)^(2+a) expands the (x+y)^2 factor
842  if (ex_to<numeric>(a.overall_coeff).is_integer()) {
843  const numeric &num_exponent = ex_to<numeric>(a.overall_coeff);
844  long int_exponent = num_exponent.to_int();
845  if (int_exponent > 0 && is_exactly_a<add>(expanded_basis))
846  distrseq.push_back(expand_add(ex_to<add>(expanded_basis), int_exponent, options));
847  else
848  distrseq.push_back(pow(expanded_basis, a.overall_coeff));
849  } else
850  distrseq.push_back(pow(expanded_basis, a.overall_coeff));
851 
852  // Make sure that e.g. (x+y)^(1+a) -> x*(x+y)^a + y*(x+y)^a
853  ex r = dynallocate<mul>(distrseq);
854  return r.expand(options);
855  }
856 
857  if (!is_exactly_a<numeric>(expanded_exponent) ||
858  !ex_to<numeric>(expanded_exponent).is_integer()) {
859  if (are_ex_trivially_equal(basis,expanded_basis) && are_ex_trivially_equal(exponent,expanded_exponent)) {
860  return this->hold();
861  } else {
862  return dynallocate<power>(expanded_basis, expanded_exponent).setflag(options == 0 ? status_flags::expanded : 0);
863  }
864  }
865 
866  // integer numeric exponent
867  const numeric & num_exponent = ex_to<numeric>(expanded_exponent);
868  long int_exponent = num_exponent.to_long();
869 
870  // (x+y)^n, n>0
871  if (int_exponent > 0 && is_exactly_a<add>(expanded_basis))
872  return expand_add(ex_to<add>(expanded_basis), int_exponent, options);
873 
874  // (x*y)^n -> x^n * y^n
875  if (is_exactly_a<mul>(expanded_basis))
876  return expand_mul(ex_to<mul>(expanded_basis), num_exponent, options, true);
877 
878  // cannot expand further
879  if (are_ex_trivially_equal(basis,expanded_basis) && are_ex_trivially_equal(exponent,expanded_exponent))
880  return this->hold();
881  else
882  return dynallocate<power>(expanded_basis, expanded_exponent).setflag(options == 0 ? status_flags::expanded : 0);
883 }
884 
886 // new virtual functions which can be overridden by derived classes
888 
889 // none
890 
892 // non-virtual functions in this class
894 
897 ex power::expand_add(const add & a, long n, unsigned options)
898 {
899  // The special case power(+(x,...y;x),2) can be optimized better.
900  if (n==2)
901  return expand_add_2(a, options);
902 
903  // method:
904  //
905  // Consider base as the sum of all symbolic terms and the overall numeric
906  // coefficient and apply the binomial theorem:
907  // S = power(+(x,...,z;c),n)
908  // = power(+(+(x,...,z;0);c),n)
909  // = sum(binomial(n,k)*power(+(x,...,z;0),k)*c^(n-k), k=1..n) + c^n
910  // Then, apply the multinomial theorem to expand all power(+(x,...,z;0),k):
911  // The multinomial theorem is computed by an outer loop over all
912  // partitions of the exponent and an inner loop over all compositions of
913  // that partition. This method makes the expansion a combinatorial
914  // problem and allows us to directly construct the expanded sum and also
915  // to re-use the multinomial coefficients (since they depend only on the
916  // partition, not on the composition).
917  //
918  // multinomial power(+(x,y,z;0),3) example:
919  // partition : compositions : multinomial coefficient
920  // [0,0,3] : [3,0,0],[0,3,0],[0,0,3] : 3!/(3!*0!*0!) = 1
921  // [0,1,2] : [2,1,0],[1,2,0],[2,0,1],... : 3!/(2!*1!*0!) = 3
922  // [1,1,1] : [1,1,1] : 3!/(1!*1!*1!) = 6
923  // => (x + y + z)^3 =
924  // x^3 + y^3 + z^3
925  // + 3*x^2*y + 3*x*y^2 + 3*y^2*z + 3*y*z^2 + 3*x*z^2 + 3*x^2*z
926  // + 6*x*y*z
927  //
928  // multinomial power(+(x,y,z;0),4) example:
929  // partition : compositions : multinomial coefficient
930  // [0,0,4] : [4,0,0],[0,4,0],[0,0,4] : 4!/(4!*0!*0!) = 1
931  // [0,1,3] : [3,1,0],[1,3,0],[3,0,1],... : 4!/(3!*1!*0!) = 4
932  // [0,2,2] : [2,2,0],[2,0,2],[0,2,2] : 4!/(2!*2!*0!) = 6
933  // [1,1,2] : [2,1,1],[1,2,1],[1,1,2] : 4!/(2!*1!*1!) = 12
934  // (no [1,1,1,1] partition since it has too many parts)
935  // => (x + y + z)^4 =
936  // x^4 + y^4 + z^4
937  // + 4*x^3*y + 4*x*y^3 + 4*y^3*z + 4*y*z^3 + 4*x*z^3 + 4*x^3*z
938  // + 6*x^2*y^2 + 6*y^2*z^2 + 6*x^2*z^2
939  // + 12*x^2*y*z + 12*x*y^2*z + 12*x*y*z^2
940  //
941  // Summary:
942  // r = 0
943  // for k from 0 to n:
944  // f = c^(n-k)*binomial(n,k)
945  // for p in all partitions of n with m parts (including zero parts):
946  // h = f * multinomial coefficient of p
947  // for c in all compositions of p:
948  // t = 1
949  // for e in all elements of c:
950  // t = t * a[e]^e
951  // r = r + h*t
952  // return r
953 
954  epvector result;
955  // The number of terms will be the number of combinatorial compositions,
956  // i.e. the number of unordered arrangements of m nonnegative integers
957  // which sum up to n. It is frequently written as C_n(m) and directly
958  // related with binomial coefficients: binomial(n+m-1,m-1).
959  size_t result_size = binomial(numeric(n+a.nops()-1), numeric(a.nops()-1)).to_long();
960  if (!a.overall_coeff.is_zero()) {
961  // the result's overall_coeff is one of the terms
962  --result_size;
963  }
964  result.reserve(result_size);
965 
966  // Iterate over all terms in binomial expansion of
967  // S = power(+(x,...,z;c),n)
968  // = sum(binomial(n,k)*power(+(x,...,z;0),k)*c^(n-k), k=1..n) + c^n
969  for (int k = 1; k <= n; ++k) {
970  numeric binomial_coefficient; // binomial(n,k)*c^(n-k)
971  if (a.overall_coeff.is_zero()) {
972  // degenerate case with zero overall_coeff:
973  // apply multinomial theorem directly to power(+(x,...z;0),n)
974  binomial_coefficient = 1;
975  if (k < n) {
976  continue;
977  }
978  } else {
979  binomial_coefficient = binomial(numeric(n), numeric(k)) * pow(ex_to<numeric>(a.overall_coeff), numeric(n-k));
980  }
981 
982  // Multinomial expansion of power(+(x,...,z;0),k)*c^(n-k):
983  // Iterate over all partitions of k with exactly as many parts as
984  // there are symbolic terms in the basis (including zero parts).
985  partition_with_zero_parts_generator partitions(k, a.seq.size());
986  do {
987  const std::vector<unsigned>& partition = partitions.get();
988  // All monomials of this partition have the same number of terms and the same coefficient.
989  const unsigned msize = std::count_if(partition.begin(), partition.end(), [](int i) { return i > 0; });
990  const numeric coeff = multinomial_coefficient(partition) * binomial_coefficient;
991 
992  // Iterate over all compositions of the current partition.
993  composition_generator compositions(partition);
994  do {
995  const std::vector<unsigned>& exponent = compositions.get();
996  epvector monomial;
997  monomial.reserve(msize);
998  numeric factor = coeff;
999  for (unsigned i = 0; i < exponent.size(); ++i) {
1000  const ex & r = a.seq[i].rest;
1001  GINAC_ASSERT(!is_exactly_a<add>(r));
1002  GINAC_ASSERT(!is_exactly_a<power>(r) ||
1003  !is_exactly_a<numeric>(ex_to<power>(r).exponent) ||
1004  !ex_to<numeric>(ex_to<power>(r).exponent).is_pos_integer() ||
1005  !is_exactly_a<add>(ex_to<power>(r).basis) ||
1006  !is_exactly_a<mul>(ex_to<power>(r).basis) ||
1007  !is_exactly_a<power>(ex_to<power>(r).basis));
1008  GINAC_ASSERT(is_exactly_a<numeric>(a.seq[i].coeff));
1009  const numeric & c = ex_to<numeric>(a.seq[i].coeff);
1010  if (exponent[i] == 0) {
1011  // optimize away
1012  } else if (exponent[i] == 1) {
1013  // optimized
1014  monomial.emplace_back(expair(r, _ex1));
1015  if (c != *_num1_p)
1016  factor = factor.mul(c);
1017  } else { // general case exponent[i] > 1
1018  monomial.emplace_back(expair(r, exponent[i]));
1019  if (c != *_num1_p)
1020  factor = factor.mul(c.power(exponent[i]));
1021  }
1022  }
1023  result.emplace_back(expair(mul(std::move(monomial)).expand(options), factor));
1024  } while (compositions.next());
1025  } while (partitions.next());
1026  }
1027 
1028  GINAC_ASSERT(result.size() == result_size);
1029  if (a.overall_coeff.is_zero()) {
1030  return dynallocate<add>(std::move(result)).setflag(status_flags::expanded);
1031  } else {
1032  return dynallocate<add>(std::move(result), ex_to<numeric>(a.overall_coeff).power(n)).setflag(status_flags::expanded);
1033  }
1034 }
1035 
1036 
1039 ex power::expand_add_2(const add & a, unsigned options)
1040 {
1041  epvector result;
1042  size_t result_size = (a.nops() * (a.nops()+1)) / 2;
1043  if (!a.overall_coeff.is_zero()) {
1044  // the result's overall_coeff is one of the terms
1045  --result_size;
1046  }
1047  result.reserve(result_size);
1048 
1049  auto last = a.seq.end();
1050 
1051  // power(+(x,...,z;c),2)=power(+(x,...,z;0),2)+2*c*+(x,...,z;0)+c*c
1052  // first part: ignore overall_coeff and expand other terms
1053  for (auto cit0=a.seq.begin(); cit0!=last; ++cit0) {
1054  const ex & r = cit0->rest;
1055  const ex & c = cit0->coeff;
1056 
1057  GINAC_ASSERT(!is_exactly_a<add>(r));
1058  GINAC_ASSERT(!is_exactly_a<power>(r) ||
1059  !is_exactly_a<numeric>(ex_to<power>(r).exponent) ||
1060  !ex_to<numeric>(ex_to<power>(r).exponent).is_pos_integer() ||
1061  !is_exactly_a<add>(ex_to<power>(r).basis) ||
1062  !is_exactly_a<mul>(ex_to<power>(r).basis) ||
1063  !is_exactly_a<power>(ex_to<power>(r).basis));
1064 
1065  if (c.is_equal(_ex1)) {
1066  if (is_exactly_a<mul>(r)) {
1067  result.emplace_back(expair(expand_mul(ex_to<mul>(r), *_num2_p, options, true),
1068  _ex1));
1069  } else {
1070  result.emplace_back(expair(dynallocate<power>(r, _ex2),
1071  _ex1));
1072  }
1073  } else {
1074  if (is_exactly_a<mul>(r)) {
1075  result.emplace_back(expair(expand_mul(ex_to<mul>(r), *_num2_p, options, true),
1076  ex_to<numeric>(c).power_dyn(*_num2_p)));
1077  } else {
1078  result.emplace_back(expair(dynallocate<power>(r, _ex2),
1079  ex_to<numeric>(c).power_dyn(*_num2_p)));
1080  }
1081  }
1082 
1083  for (auto cit1=cit0+1; cit1!=last; ++cit1) {
1084  const ex & r1 = cit1->rest;
1085  const ex & c1 = cit1->coeff;
1086  result.emplace_back(expair(mul(r,r1).expand(options),
1087  _num2_p->mul(ex_to<numeric>(c)).mul_dyn(ex_to<numeric>(c1))));
1088  }
1089  }
1090 
1091  // second part: add terms coming from overall_coeff (if != 0)
1092  if (!a.overall_coeff.is_zero()) {
1093  for (auto & i : a.seq)
1094  result.push_back(a.combine_pair_with_coeff_to_pair(i, ex_to<numeric>(a.overall_coeff).mul_dyn(*_num2_p)));
1095  }
1096 
1097  GINAC_ASSERT(result.size() == result_size);
1098 
1099  if (a.overall_coeff.is_zero()) {
1100  return dynallocate<add>(std::move(result)).setflag(status_flags::expanded);
1101  } else {
1102  return dynallocate<add>(std::move(result), ex_to<numeric>(a.overall_coeff).power(2)).setflag(status_flags::expanded);
1103  }
1104 }
1105 
1108 ex power::expand_mul(const mul & m, const numeric & n, unsigned options, bool from_expand)
1109 {
1110  GINAC_ASSERT(n.is_integer());
1111 
1112  if (n.is_zero()) {
1113  return _ex1;
1114  }
1115 
1116  // do not bother to rename indices if there are no any.
1118  m.info(info_flags::has_indices))
1120  // Leave it to multiplication since dummy indices have to be renamed
1122  (get_all_dummy_indices(m).size() > 0) && n.is_positive()) {
1123  ex result = m;
1125  sort(va.begin(), va.end(), ex_is_less());
1126 
1127  for (int i=1; i < n.to_int(); i++)
1128  result *= rename_dummy_indices_uniquely(va, m);
1129  return result;
1130  }
1131 
1132  epvector distrseq;
1133  distrseq.reserve(m.seq.size());
1134  bool need_reexpand = false;
1135 
1136  for (auto & cit : m.seq) {
1137  expair p = m.combine_pair_with_coeff_to_pair(cit, n);
1138  if (from_expand && is_exactly_a<add>(cit.rest) && ex_to<numeric>(p.coeff).is_pos_integer()) {
1139  // this happens when e.g. (a+b)^(1/2) gets squared and
1140  // the resulting product needs to be reexpanded
1141  need_reexpand = true;
1142  }
1143  distrseq.push_back(p);
1144  }
1145 
1146  const mul & result = dynallocate<mul>(std::move(distrseq), ex_to<numeric>(m.overall_coeff).power_dyn(n));
1147  if (need_reexpand)
1148  return ex(result).expand(options);
1149  if (from_expand)
1150  return result.setflag(status_flags::expanded);
1151  return result;
1152 }
1153 
1155 
1156 } // namespace GiNaC
Interface to GiNaC's sums of expressions.
Archiving of GiNaC expressions.
#define GINAC_ASSERT(X)
Assertion macro for checking invariances.
Definition: assertion.h:33
Sum of expressions.
Definition: add.h:32
ex recombine_pair_to_ex(const expair &p) const override
Definition: add.cpp:564
expair combine_pair_with_coeff_to_pair(const expair &p, const ex &c) const override
Definition: add.cpp:550
This class stores all properties needed to record/retrieve the state of one object of class basic (or...
Definition: archive.h:49
This class is the ABC (abstract base class) of GiNaC's class hierarchy.
Definition: basic.h:105
virtual bool has(const ex &other, unsigned options=0) const
Test for occurrence of a pattern.
Definition: basic.cpp:280
friend class ex
Definition: basic.h:108
unsigned flags
of type status_flags
Definition: basic.h:302
ex subs_one_level(const exmap &m, unsigned options) const
Helper function for subs().
Definition: basic.cpp:585
const basic & hold() const
Stop further evaluation.
Definition: basic.cpp:887
bool is_equal(const basic &other) const
Test for syntactic equality.
Definition: basic.cpp:863
const basic & setflag(unsigned f) const
Set some status_flags.
Definition: basic.h:288
virtual int compare_same_type(const basic &other) const
Returns order relation between two objects of same type.
Definition: basic.cpp:719
const basic & clearflag(unsigned f) const
Clear some status_flags.
Definition: basic.h:291
Generate all compositions of a partition of an integer n, starting with the compositions which has no...
Definition: utils.h:395
const std::vector< unsigned > & get() const
Definition: utils.h:460
Wrapper template for making GiNaC classes out of STL containers.
Definition: container.h:73
Lightweight wrapper for GiNaC's symbolic objects.
Definition: ex.h:72
numeric integer_content() const
Compute the integer content (= GCD of all numeric coefficients) of an expanded polynomial.
Definition: normal.cpp:318
bool match(const ex &pattern) const
Check whether expression matches a specified pattern.
Definition: ex.cpp:95
bool is_polynomial(const ex &vars) const
Check whether expression is a polynomial.
Definition: ex.cpp:241
ex diff(const symbol &s, unsigned nth=1) const
Compute partial derivative of an expression.
Definition: ex.cpp:86
ex expand(unsigned options=0) const
Definition: ex.cpp:73
bool is_equal(const ex &other) const
Definition: ex.h:345
int degree(const ex &s) const
Definition: ex.h:173
bool has(const ex &pattern, unsigned options=0) const
Definition: ex.h:151
ex evalf() const
Definition: ex.h:121
ex conjugate() const
Definition: ex.h:146
unsigned return_type() const
Definition: ex.h:230
return_type_t return_type_tinfo() const
Definition: ex.h:231
ex imag_part() const
Definition: ex.h:148
ex subs(const exmap &m, unsigned options=0) const
Definition: ex.h:826
bool info(unsigned inf) const
Definition: ex.h:132
int compare(const ex &other) const
Definition: ex.h:322
bool is_zero() const
Definition: ex.h:213
void print(const print_context &c, unsigned level=0) const
Print expression to stream.
Definition: ex.cpp:56
ex op(size_t i) const
Definition: ex.h:136
int ldegree(const ex &s) const
Definition: ex.h:174
ex real_part() const
Definition: ex.h:147
ex evalm() const
Definition: ex.h:122
ex coeff(const ex &s, int n=1) const
Definition: ex.h:175
A pair of expressions.
Definition: expair.h:38
ex coeff
second member of pair, must be numeric
Definition: expair.h:91
size_t nops() const override
Number of operands/members.
epvector seq
Definition: expairseq.h:127
@ expand_rename_idx
used internally by mul::expand()
Definition: flags.h:34
@ algebraic
enable algebraic matching
Definition: flags.h:43
@ integer_polynomial
Definition: flags.h:256
@ cinteger_polynomial
Definition: flags.h:257
@ crational_polynomial
Definition: flags.h:259
@ rational_polynomial
Definition: flags.h:258
Product of expressions.
Definition: mul.h:32
Non-commutative product of expressions.
Definition: ncmul.h:33
This class is a wrapper around CLN-numbers within the GiNaC class hierarchy.
Definition: numeric.h:82
bool is_pos_integer() const
True if object is an exact integer greater than zero.
Definition: numeric.cpp:1161
const numeric & mul_dyn(const numeric &other) const
Numerical multiplication method.
Definition: numeric.cpp:957
bool is_crational() const
True if object is an exact rational number, may even be complex (denominator may be unity).
Definition: numeric.cpp:1243
bool is_rational() const
True if object is an exact rational number, may even be complex (denominator may be unity).
Definition: numeric.cpp:1201
const numeric real() const
Real part of a number.
Definition: numeric.cpp:1339
long to_long() const
Converts numeric types to machine's long.
Definition: numeric.cpp:1313
int compare(const numeric &other) const
This method establishes a canonical order on all numbers.
Definition: numeric.cpp:1104
bool is_positive() const
True if object is not complex and greater than zero.
Definition: numeric.cpp:1136
bool is_real() const
True if object is a real integer, rational or float (but not complex).
Definition: numeric.cpp:1208
const numeric numer() const
Numerator.
Definition: numeric.cpp:1356
bool is_integer() const
True if object is a non-complex integer.
Definition: numeric.cpp:1154
const numeric power(const numeric &other) const
Numerical exponentiation.
Definition: numeric.cpp:900
const numeric denom() const
Denominator.
Definition: numeric.cpp:1387
const numeric mul(const numeric &other) const
Numerical multiplication method.
Definition: numeric.cpp:880
bool is_equal(const numeric &other) const
Definition: numeric.cpp:1122
int to_int() const
Converts numeric types to machine's int.
Definition: numeric.cpp:1303
const numeric inverse() const
Inverse of a number.
Definition: numeric.cpp:1053
bool is_zero() const
True if object is zero.
Definition: numeric.cpp:1129
const numeric div(const numeric &other) const
Numerical division method.
Definition: numeric.cpp:890
Generate all bounded combinatorial partitions of an integer n with exactly m parts (including zero pa...
Definition: utils.h:327
const std::vector< unsigned > & get() const
Definition: utils.h:337
Exception class thrown when a singularity is encountered.
Definition: numeric.h:70
This class holds a two-component object, a basis and and exponent representing exponentiation.
Definition: power.h:39
static ex expand_mul(const mul &m, const numeric &n, unsigned options, bool from_expand=false)
Expand factors of m in m^n where m is a mul and n is an integer.
Definition: power.cpp:1108
void do_print_dflt(const print_dflt &c, unsigned level) const
Definition: power.cpp:107
void do_print_csrc(const print_csrc &c, unsigned level) const
Definition: power.cpp:179
static ex expand_add(const add &a, long n, unsigned options)
expand a^n where a is an add and n is a positive integer.
Definition: power.cpp:897
int degree(const ex &s) const override
Return degree of highest power in object s.
Definition: power.cpp:301
void read_archive(const archive_node &n, lst &syms) override
Read (a.k.a.
Definition: power.cpp:73
ex subs(const exmap &m, unsigned options=0) const override
Substitute a set of objects by arbitrary expressions.
Definition: power.cpp:615
int ldegree(const ex &s) const override
Return degree of lowest power in object s.
Definition: power.cpp:316
static ex expand_add_2(const add &a, unsigned options)
Special case of power::expand_add.
Definition: power.cpp:1039
ex real_part() const override
Definition: power.cpp:667
void archive(archive_node &n) const override
Save (a.k.a.
Definition: power.cpp:80
ex derivative(const symbol &s) const override
Implementation of ex::diff() for a power.
Definition: power.cpp:744
friend class mul
Definition: power.h:42
ex map(map_function &f) const override
Construct new expression by applying the specified function to all sub-expressions (one level only,...
Definition: power.cpp:275
ex eval() const override
Perform automatic term rewriting rules in this class.
Definition: power.cpp:374
ex basis
Definition: power.h:105
void do_print_python(const print_python &c, unsigned level) const
Definition: power.cpp:210
void print_power(const print_context &c, const char *powersymbol, const char *openbrace, const char *closebrace, unsigned level) const
Definition: power.cpp:93
void do_print_python_repr(const print_python_repr &c, unsigned level) const
Definition: power.cpp:215
ex exponent
Definition: power.h:106
unsigned precedence() const override
Return relative operator precedence (for parenthezing output).
Definition: power.h:53
ex op(size_t i) const override
Return operand/member at position i.
Definition: power.cpp:268
power(const ex &lh, const ex &rh)
Definition: power.h:48
void do_print_latex(const print_latex &c, unsigned level) const
Definition: power.cpp:120
ex conjugate() const override
Definition: power.cpp:646
ex coeff(const ex &s, int n=1) const override
Return coefficient of degree n in object s.
Definition: power.cpp:331
bool has(const ex &other, unsigned options=0) const override
Test for occurrence of a pattern.
Definition: power.cpp:590
ex evalm() const override
Evaluate sums, products and integer powers of matrices.
Definition: power.cpp:578
return_type_t return_type_tinfo() const override
Definition: power.cpp:773
size_t nops() const override
Number of operands/members.
Definition: power.cpp:263
ex expand(unsigned options=0) const override
Expand expression, i.e.
Definition: power.cpp:778
bool info(unsigned inf) const override
Information about the object.
Definition: power.cpp:224
ex eval_ncmul(const exvector &v) const override
Definition: power.cpp:641
void do_print_csrc_cl_N(const print_csrc_cl_N &c, unsigned level) const
Definition: power.cpp:164
unsigned return_type() const override
Definition: power.cpp:768
ex imag_part() const override
Definition: power.cpp:703
bool is_polynomial(const ex &var) const override
Check whether this is a polynomial in the given variables.
Definition: power.cpp:287
ex evalf() const override
Evaluate object numerically.
Definition: power.cpp:565
Base class for print_contexts.
Definition: print.h:103
Context for C source output using CLN numbers.
Definition: print.h:182
Base context for C source output.
Definition: print.h:158
Context for default (ginsh-parsable) output.
Definition: print.h:115
Context for latex-parsable output.
Definition: print.h:123
Context for python-parsable output.
Definition: print.h:139
Context for python pretty-print output.
Definition: print.h:131
@ expanded
.expand(0) has already done its job (other expand() options ignore this flag)
Definition: flags.h:204
@ evaluated
.eval() has already done its job
Definition: flags.h:203
@ hash_calculated
.calchash() has already done its job
Definition: flags.h:205
@ no_pattern
disable pattern matching
Definition: flags.h:51
@ algebraic
enable algebraic substitutions
Definition: flags.h:53
Basic CAS symbol.
Definition: symbol.h:39
Definition of optimizing macros.
#define likely(cond)
Definition: compiler.h:32
Interface to GiNaC's constant types and some special constants.
Interface to sequences of expression pairs.
unsigned options
Definition: factor.cpp:2480
vector< int > k
Definition: factor.cpp:1466
ex x
Definition: factor.cpp:1641
size_t n
Definition: factor.cpp:1463
size_t c
Definition: factor.cpp:770
size_t r
Definition: factor.cpp:770
mvec m
Definition: factor.cpp:771
size_t last
Definition: factor.cpp:1465
Interface to GiNaC's indexed expressions.
Interface to GiNaC's initially known functions.
Definition of GiNaC's lst.
Interface to symbolic matrices.
Interface to GiNaC's products of expressions.
Definition: add.cpp:38
const numeric * _num_1_p
Definition: utils.cpp:351
const ex _ex2
Definition: utils.cpp:389
const numeric pow(const numeric &x, const numeric &y)
Definition: numeric.h:251
const ex _ex1_2
Definition: utils.cpp:381
std::map< ex, ex, ex_is_less > exmap
Definition: basic.h:50
std::vector< expair > epvector
expair-vector
Definition: expairseq.h:33
const numeric abs(const numeric &x)
Absolute value.
Definition: numeric.cpp:2315
bool is_negative(const numeric &x)
Definition: numeric.h:269
const ex _ex1
Definition: utils.cpp:385
bool are_ex_trivially_equal(const ex &e1, const ex &e2)
Compare two objects of class quickly without doing a deep tree traversal.
Definition: ex.h:684
const numeric binomial(const numeric &n, const numeric &k)
The Binomial coefficients.
Definition: numeric.cpp:2143
print_func< print_dflt >(&diracone::do_print). print_func< print_latex >(&diracone
Definition: clifford.cpp:51
const numeric * _num2_p
Definition: utils.cpp:388
const numeric exp(const numeric &x)
Exponential function.
Definition: numeric.cpp:1439
const numeric cos(const numeric &x)
Numeric cosine (trigonometric function).
Definition: numeric.cpp:1470
const ex _ex_1
Definition: utils.cpp:352
GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(add, expairseq, print_func< print_context >(&add::do_print). print_func< print_latex >(&add::do_print_latex). print_func< print_csrc >(&add::do_print_csrc). print_func< print_tree >(&add::do_print_tree). print_func< print_python_repr >(&add::do_print_python_repr)) add
Definition: add.cpp:40
const numeric iquo(const numeric &a, const numeric &b)
Numeric integer quotient.
Definition: numeric.cpp:2404
GINAC_IMPLEMENT_REGISTERED_CLASS_OPT_T(lst, basic, print_func< print_context >(&lst::do_print). print_func< print_tree >(&lst::do_print_tree)) template<> bool lst GINAC_BIND_UNARCHIVER(lst)
Specialization of container::info() for lst.
Definition: lst.cpp:42
bool is_pos_integer(const numeric &x)
Definition: numeric.h:275
const numeric log(const numeric &x)
Natural logarithm.
Definition: numeric.cpp:1450
const numeric sin(const numeric &x)
Numeric sine (trigonometric function).
Definition: numeric.cpp:1461
const numeric * _num1_p
Definition: utils.cpp:384
ex numer(const ex &thisex)
Definition: ex.h:745
ex factor(const ex &poly, unsigned options)
Interface function to the outside world.
Definition: factor.cpp:2581
bool is_integer(const numeric &x)
Definition: numeric.h:272
static void print_sym_pow(const print_context &c, const symbol &x, int exp)
Definition: power.cpp:140
lst rename_dummy_indices_uniquely(const exvector &va, const exvector &vb)
Similar to above, where va and vb are the same and the return value is a list of two lists for substi...
Definition: indexed.cpp:1460
const ex _ex0
Definition: utils.cpp:369
bool tryfactsubs(const ex &origfactor, const ex &patternfactor, int &nummatches, exmap &repls)
Definition: mul.cpp:672
std::vector< ex > exvector
Definition: basic.h:46
exvector get_all_dummy_indices(const ex &e)
Returns all dummy indices from the exvector.
Definition: indexed.cpp:1435
const numeric multinomial_coefficient(const std::vector< unsigned > &p)
Compute the multinomial coefficient n!/(p1!*p2!*...*pk!) where n = p1+p2+...+pk, i....
Definition: utils.cpp:60
const numeric * _num0_p
Definition: utils.cpp:367
Interface to GiNaC's non-commutative products of expressions.
Makes the interface to the underlying bignum package available.
Interface to GiNaC's overloaded operators.
Interface to GiNaC's symbolic exponentiation (basis^exponent).
Interface to relations between expressions.
Function object for map().
Definition: basic.h:85
To distinguish between different kinds of non-commutative objects.
Definition: registrar.h:44
Interface to GiNaC's symbolic objects.
Interface to several small and furry utilities needed within GiNaC but not of any interest to the use...

This page is part of the GiNaC developer's reference. It was generated automatically by doxygen. For an introduction, see the tutorial.