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