76 while (i !=
seq.end()) {
83 GINAC_ASSERT(ex_to<numeric>(i->coeff) < ex_to<numeric>(ip1->coeff));
93 : seq(
std::move(ops_))
97 while (i !=
seq.end()) {
104 GINAC_ASSERT(ex_to<numeric>(i->coeff) < ex_to<numeric>(ip1->coeff));
121 inherited::read_archive(
n, sym_lst);
122 auto range =
n.find_property_range(
"coeff",
"power");
123 seq.reserve((range.end-range.begin)/2);
125 for (
auto loc = range.begin; loc < range.end;) {
128 n.find_ex_by_loc(loc++, rest, sym_lst);
129 n.find_ex_by_loc(loc++,
coeff, sym_lst);
133 n.find_ex(
"var",
var, sym_lst);
134 n.find_ex(
"point",
point, sym_lst);
139 inherited::archive(
n);
140 for (
auto & it :
seq) {
141 n.add_ex(
"coeff", it.rest);
142 n.add_ex(
"power", it.coeff);
144 n.add_ex(
"var",
var);
163 auto i =
seq.begin(), end =
seq.end();
167 if (i !=
seq.begin())
177 c.s << openbrace <<
'(';
179 c.s <<
')' << closebrace;
183 if (!i->coeff.is_zero()) {
186 c.s << openbrace <<
'(';
188 c.s <<
')' << closebrace;
191 if (i->coeff.compare(
_ex1)) {
229 c.s << std::string(level,
' ') << class_name() <<
" @" <<
this
230 << std::hex <<
", hash=0x" <<
hashvalue <<
", flags=0x" <<
flags << std::dec
232 size_t num =
seq.size();
233 for (
size_t i=0; i<num; ++i) {
234 seq[i].rest.print(
c, level +
c.delta_indent);
235 seq[i].coeff.print(
c, level +
c.delta_indent);
236 c.s << std::string(level +
c.delta_indent,
' ') <<
"-----" << std::endl;
244 c.s << class_name() <<
"(relational(";
249 size_t num =
seq.size();
250 for (
size_t i=0; i<num; ++i) {
254 seq[i].rest.print(
c);
256 seq[i].coeff.print(
c);
268 if (
seq.size()>o.
seq.size())
270 if (
seq.size()<o.
seq.size())
282 auto it =
seq.begin(), o_it = o.
seq.begin();
283 while (it!=
seq.end() && o_it!=o.
seq.end()) {
284 cmpval = it->compare(*o_it);
305 throw (std::out_of_range(
"op() out of range"));
322 return ex_to<numeric>((
seq.end()-1)->coeff).to_int();
324 int max_pow = std::numeric_limits<int>::min();
325 for (
auto & it :
seq)
326 max_pow = std::max(max_pow, it.rest.degree(s));
342 return ex_to<numeric>((
seq.begin())->coeff).to_int();
344 int min_pow = std::numeric_limits<int>::max();
345 for (
auto & it :
seq)
346 min_pow = std::min(min_pow, it.rest.degree(s));
365 int lo = 0, hi =
seq.size() - 1;
367 int mid = (lo + hi) / 2;
369 int cmp = ex_to<numeric>(
seq[mid].
coeff).compare(looking_for);
375 return seq[mid].rest;
380 throw(std::logic_error(
"pseries::coeff: compare() didn't return -1, 0 or 1"));
405 new_seq.reserve(
seq.size());
406 for (
auto & it :
seq)
407 new_seq.emplace_back(
expair(it.rest.evalf(), it.coeff));
415 return conjugate_function(*this).hold();
424 return dynallocate<pseries>(
var==newpoint, newseq ? std::move(*newseq) :
seq);
430 return real_part_function(*this).hold();
432 if(newpoint !=
point)
433 return real_part_function(*this).hold();
436 v.reserve(
seq.size());
437 for (
auto & it :
seq)
438 v.emplace_back(
expair(it.rest.real_part(), it.coeff));
439 return dynallocate<pseries>(
var==
point, std::move(v));
445 return imag_part_function(*this).hold();
447 if(newpoint !=
point)
448 return imag_part_function(*this).hold();
451 v.reserve(
seq.size());
452 for (
auto & it :
seq)
453 v.emplace_back(
expair(it.rest.imag_part(), it.coeff));
454 return dynallocate<pseries>(
var==
point, std::move(v));
459 std::unique_ptr<epvector> newseq(
nullptr);
460 for (
auto i=
seq.begin(); i!=
seq.end(); ++i) {
462 newseq->emplace_back(
expair(i->rest.eval_integ(), i->coeff));
468 newseq->reserve(
seq.size());
469 for (
auto j=
seq.begin(); j!=i; ++j)
470 newseq->push_back(*j);
477 return dynallocate<pseries>(
var==newpoint, std::move(*newseq));
485 bool something_changed =
false;
486 for (
auto i=
seq.begin(); i!=
seq.end(); ++i) {
487 if (something_changed) {
494 something_changed =
true;
495 newseq.reserve(
seq.size());
496 std::copy(
seq.begin(), i, std::back_inserter<epvector>(newseq));
502 if (something_changed)
503 return dynallocate<pseries>(
var==
point, std::move(newseq));
513 if (
m.find(
var) !=
m.end())
519 newseq.reserve(
seq.size());
520 for (
auto & it :
seq)
530 for (
auto & it :
seq) {
547 for (
auto & it :
seq) {
549 new_seq.emplace_back(
expair(it.rest, it.coeff - 1));
551 ex c = it.rest * it.coeff;
553 new_seq.emplace_back(
expair(
c, it.coeff - 1));
559 for (
auto & it :
seq) {
561 new_seq.push_back(it);
563 ex c = it.rest.diff(s);
565 new_seq.emplace_back(
expair(
c, it.coeff));
576 for (
auto & it :
seq) {
594 throw (std::out_of_range(
"coeffop() out of range"));
601 throw (std::out_of_range(
"exponop() out of range"));
615 const symbol &s = ex_to<symbol>(
r.lhs());
648 deriv = deriv.
diff(s);
660 const ex point =
r.rhs();
692 auto a =
seq.begin(), a_end =
seq.end();
693 auto b = other.
seq.begin(), b_end = other.
seq.end();
694 int pow_a = std::numeric_limits<int>::max(), pow_b = std::numeric_limits<int>::max();
699 new_seq.push_back(*b);
704 pow_a = ex_to<numeric>((*a).coeff).to_int();
709 new_seq.push_back(*a);
714 pow_b = ex_to<numeric>((*b).coeff).to_int();
719 new_seq.push_back(*a);
723 }
else if (pow_b < pow_a) {
725 new_seq.push_back(*b);
732 new_seq.emplace_back(
expair(Order(
_ex1), (*a).coeff));
735 ex sum = (*a).rest + (*b).rest;
758 for (
auto & it :
seq) {
760 if (is_exactly_a<pseries>(it.rest))
764 if (!it.coeff.is_equal(
_ex1))
765 op = ex_to<pseries>(
op).mul_const(ex_to<numeric>(it.coeff));
768 acc = ex_to<pseries>(acc).add_series(ex_to<pseries>(
op));
782 new_seq.reserve(
seq.size());
784 for (
auto & it :
seq) {
786 new_seq.emplace_back(
expair(it.rest * other, it.
coeff));
788 new_seq.push_back(it);
808 if (
seq.empty() || other.
seq.empty()) {
818 const int cdeg_min = a_min + b_min;
819 int cdeg_max = a_max + b_max;
821 int higher_order_a = std::numeric_limits<int>::max();
822 int higher_order_b = std::numeric_limits<int>::max();
824 higher_order_a = a_max + b_min;
826 higher_order_b = b_max + a_min;
827 const int higher_order_c = std::min(higher_order_a, higher_order_b);
828 if (cdeg_max >= higher_order_c)
829 cdeg_max = higher_order_c - 1;
831 std::map<int, ex> rest_map_a, rest_map_b;
832 for (
const auto& it :
seq)
833 rest_map_a[ex_to<numeric>(it.coeff).to_int()] = it.rest;
836 for (
const auto& it : other.
seq)
837 rest_map_b[ex_to<numeric>(it.coeff).to_int()] = it.rest;
839 for (
int cdeg=cdeg_min; cdeg<=cdeg_max; ++cdeg) {
842 for (
int i=a_min; cdeg-i>=b_min; ++i) {
843 const auto& ita = rest_map_a.
find(i);
844 if (ita == rest_map_a.end())
846 const auto& itb = rest_map_b.find(cdeg-i);
847 if (itb == rest_map_b.end())
850 co += ita->second * itb->second;
855 if (higher_order_c < std::numeric_limits<int>::max())
869 const ex& sym =
r.lhs();
872 std::vector<int> ldegrees;
873 std::vector<bool> ldegree_redo;
878 for (
auto & it :
seq) {
885 factor = ex_to<numeric>(expon).to_int();
890 int real_ldegree = 0;
891 bool flag_redo =
false;
894 }
catch (std::runtime_error &) {}
896 if (real_ldegree == 0) {
904 }
while (real_ldegree == orderloop);
910 if (real_ldegree == 0)
915 ldegrees.push_back(
factor * real_ldegree);
916 ldegree_redo.push_back(flag_redo);
919 int degbound =
order-std::accumulate(ldegrees.begin(), ldegrees.end(), 0);
924 for (
auto & it :
seq) {
925 if ( ldegree_redo[j] ) {
931 factor = ex_to<numeric>(expon).to_int();
935 int real_ldegree = 0;
940 }
while ((real_ldegree == orderloop)
941 && (
factor*real_ldegree < degbound));
942 ldegrees[j] =
factor * real_ldegree;
943 degbound -=
factor * real_ldegree;
948 int degsum = std::accumulate(ldegrees.begin(), ldegrees.end(), 0);
950 if (degsum >
order) {
955 auto itd = ldegrees.begin();
956 for (
auto it=
seq.begin(), itend=
seq.end(); it!=itend; ++it, ++itd) {
962 if (it ==
seq.begin())
963 acc = ex_to<pseries>(
op);
965 acc = ex_to<pseries>(acc.
mul_series(ex_to<pseries>(
op)));
1003 throw std::domain_error(
"pseries::power_const(): pow(0,I) is undefined");
1005 throw pole_error(
"pseries::power_const(): division by zero",1);
1012 throw std::runtime_error(
"pseries::power_const(): trying to assemble a Puiseux series");
1015 int numcoeff = deg - (p*ldeg).
to_int();
1016 if (numcoeff <= 0) {
1023 throw pole_error(
"pseries::power_const(): division by zero",1);
1027 co.reserve(numcoeff);
1029 for (
int i=1; i<numcoeff; ++i) {
1031 for (
int j=1; j<=i; ++j) {
1034 co.push_back(Order(
_ex1));
1037 sum += (p * j - (i - j)) * co[i - j] *
c;
1039 co.push_back(sum /
coeff(
var, ldeg) / i);
1044 bool higher_order =
false;
1045 for (
int i=0; i<numcoeff; ++i) {
1047 new_seq.emplace_back(
expair(co[i], p * ldeg + i));
1049 higher_order =
true;
1054 new_seq.emplace_back(
expair(Order(
_ex1), p * ldeg + numcoeff));
1064 for (
auto & it : newseq)
1076 if (is_exactly_a<pseries>(
basis))
1080 bool must_expand_basis =
false;
1084 must_expand_basis =
true;
1087 bool exponent_is_regular =
true;
1091 exponent_is_regular =
false;
1094 if (!exponent_is_regular) {
1127 return pseries(
r, std::move(new_seq));
1138 const ex& sym =
r.lhs();
1141 int real_ldegree = 0;
1143 real_ldegree = eb.
ldegree(sym-
r.rhs());
1144 if (real_ldegree == 0) {
1149 }
while (real_ldegree == orderloop);
1153 throw std::runtime_error(
"pseries::power_const(): trying to assemble a Puiseux series");
1158 result = ex_to<pseries>(e).power_const(numexp,
order);
1161 result =
pseries(
r, std::move(ser));
1171 const ex p =
r.rhs();
1173 const symbol &s = ex_to<symbol>(
r.lhs());
1180 for (
auto & it :
seq) {
1181 int o = ex_to<numeric>(it.coeff).to_int();
1183 new_seq.emplace_back(
expair(Order(
_ex1), o));
1186 new_seq.push_back(it);
1188 return pseries(
r, std::move(new_seq));
1197 throw std::logic_error(
"Cannot series expand wrt dummy variable");
1202 fexpansion.reserve(fseries.
nops());
1203 for (
size_t i=0; i<fseries.
nops(); ++i) {
1204 ex currcoeff = ex_to<pseries>(fseries).coeffop(i);
1205 currcoeff = (currcoeff == Order(
_ex1))
1209 fexpansion.emplace_back(
1210 expair(currcoeff, ex_to<pseries>(fseries).exponop(i)));
1214 ex result = dynallocate<pseries>(
r, std::move(fexpansion));
1217 for (
size_t i=0; i<fseries.
nops(); ++i) {
1218 ex currcoeff = ex_to<pseries>(fseries).coeffop(i);
1221 ex currexpon = ex_to<pseries>(fseries).exponop(i);
1222 int orderforf =
order-ex_to<numeric>(currexpon).to_int()-1;
1223 currcoeff = currcoeff.
series(
r, orderforf);
1224 ex term = ex_to<pseries>(aseries).power_const(ex_to<numeric>(currexpon+1),
order);
1225 term = ex_to<pseries>(term).mul_const(ex_to<numeric>(-1/(currexpon+1)));
1226 term = ex_to<pseries>(term).mul_series(ex_to<pseries>(currcoeff));
1227 result = ex_to<pseries>(result).add_series(ex_to<pseries>(term));
1233 for (
size_t i=0; i<fseries.
nops(); ++i) {
1234 ex currcoeff = ex_to<pseries>(fseries).coeffop(i);
1237 ex currexpon = ex_to<pseries>(fseries).exponop(i);
1238 int orderforf =
order-ex_to<numeric>(currexpon).to_int()-1;
1239 currcoeff = currcoeff.
series(
r, orderforf);
1240 ex term = ex_to<pseries>(bseries).power_const(ex_to<numeric>(currexpon+1),
order);
1241 term = ex_to<pseries>(term).mul_const(ex_to<numeric>(1/(currexpon+1)));
1242 term = ex_to<pseries>(term).mul_series(ex_to<pseries>(currcoeff));
1243 result = ex_to<pseries>(result).add_series(ex_to<pseries>(term));
1264 if (is_a<relational>(
r))
1265 rel_ = ex_to<relational>(
r);
1266 else if (is_a<symbol>(
r))
1269 throw (std::logic_error(
"ex::series(): expansion point has unknown type"));
Interface to GiNaC's sums of expressions.
Archiving of GiNaC expressions.
#define GINAC_ASSERT(X)
Assertion macro for checking invariances.
ex series(const relational &r, int order, unsigned options=0) const override
Implementation of ex::series() for sums.
This class stores all properties needed to record/retrieve the state of one object of class basic (or...
This class is the ABC (abstract base class) of GiNaC's class hierarchy.
unsigned hashvalue
hash value
virtual bool has(const ex &other, unsigned options=0) const
Test for occurrence of a pattern.
unsigned flags
of type status_flags
virtual void print(const print_context &c, unsigned level=0) const
Output to stream.
const basic & hold() const
Stop further evaluation.
virtual ex series(const relational &r, int order, unsigned options=0) const
Default implementation of ex::series().
virtual int compare_same_type(const basic &other) const
Returns order relation between two objects of same type.
virtual ex coeff(const ex &s, int n=1) const
Return coefficient of degree n in object s.
Wrapper template for making GiNaC classes out of STL containers.
Lightweight wrapper for GiNaC's symbolic objects.
bool find(const ex &pattern, exset &found) const
Find all occurrences of a pattern.
ex diff(const symbol &s, unsigned nth=1) const
Compute partial derivative of an expression.
ex expand(unsigned options=0) const
bool is_equal(const ex &other) const
ptr< basic > bp
pointer to basic object managed by this
ex series(const ex &r, int order, unsigned options=0) const
Compute the truncated series expansion of an expression.
ex subs(const exmap &m, unsigned options=0) const
bool info(unsigned inf) const
int compare(const ex &other) const
ex lhs() const
Left hand side of relational expression.
void print(const print_context &c, unsigned level=0) const
Print expression to stream.
int ldegree(const ex &s) const
ex rhs() const
Right hand side of relational expression.
ex coeff(const ex &s, int n=1) const
ex op(size_t i) const override
Return operand/member at position i.
integral(const ex &x_, const ex &a_, const ex &b_, const ex &f_)
ex series(const relational &r, int order, unsigned options=0) const override
Default implementation of ex::series().
ex series(const relational &s, int order, unsigned options=0) const override
Implementation of ex::series() for product.
ex recombine_pair_to_ex(const expair &p) const override
Form an ex out of an expair, using the corresponding semantics.
This class is a wrapper around CLN-numbers within the GiNaC class hierarchy.
const numeric real() const
Real part of a number.
ex coeff(const ex &s, int n=1) const override
Return coefficient of degree n in object s.
bool is_negative() const
True if object is not complex and less than zero.
bool is_zero() const
True if object is zero.
const numeric div(const numeric &other) const
Numerical division method.
Exception class thrown when a singularity is encountered.
ex series(const relational &s, int order, unsigned options=0) const override
Implementation of ex::series() for powers.
Base class for print_contexts.
Context for latex-parsable output.
Context for python-parsable output.
Context for python pretty-print output.
Context for tree-like output for debugging.
This class holds a extended truncated power series (positive and negative integer powers).
int ldegree(const ex &s) const override
Return degree of lowest power of the series.
ex evalf() const override
Evaluate coefficients numerically.
ex var
Series variable (holds a symbol)
bool is_zero() const
Check whether series has the value zero.
pseries shift_exponents(int deg) const
Return a new pseries object with the powers shifted by deg.
void do_print(const print_context &c, unsigned level) const
ex op(size_t i) const override
Return the ith term in the series when represented as a sum.
unsigned precedence() const override
Return relative operator precedence (for parenthezing output).
void do_print_latex(const print_latex &c, unsigned level) const
ex mul_const(const numeric &other) const
Multiply a pseries object with a numeric constant, producing a pseries object that represents the pro...
ex coeff(const ex &s, int n=1) const override
Return coefficient of degree n in power series if s is the expansion variable.
void do_print_tree(const print_tree &c, unsigned level) const
ex convert_to_poly(bool no_order=false) const
Convert the pseries object to an ordinary polynomial.
size_t nops() const override
Return the number of operands including a possible order term.
void read_archive(const archive_node &n, lst &syms) override
Read (a.k.a.
ex subs(const exmap &m, unsigned options=0) const override
Substitute a set of objects by arbitrary expressions.
pseries(const ex &rel_, const epvector &ops_)
Construct pseries from a vector of coefficients and powers.
ex eval_integ() const override
Evaluate integrals, if result is known.
ex expand(unsigned options=0) const override
Implementation of ex::expand() for a power series.
ex derivative(const symbol &s) const override
Implementation of ex::diff() for a power series.
ex real_part() const override
ex evalm() const override
Evaluate sums, products and integer powers of matrices.
epvector seq
Vector of {coefficient, power} pairs.
bool is_compatible_to(const pseries &other) const
Check whether series is compatible to another series (expansion variable and point are the same.
ex exponop(size_t i) const
ex coeffop(size_t i) const
Get coefficients and exponents.
void archive(archive_node &n) const override
Save (a.k.a.
void print_series(const print_context &c, const char *openbrace, const char *closebrace, const char *mul_sym, const char *pow_sym, unsigned level) const
ex series(const relational &r, int order, unsigned options=0) const override
Re-expansion of a pseries object.
void do_print_python(const print_python &c, unsigned level) const
ex imag_part() const override
void do_print_python_repr(const print_python_repr &c, unsigned level) const
bool is_terminating() const
Returns true if there is no order term, i.e.
ex conjugate() const override
ex mul_series(const pseries &other) const
Multiply one pseries object to another, producing a pseries object that represents the product.
ex power_const(const numeric &p, int deg) const
Compute the p-th power of a series.
int degree(const ex &s) const override
Return degree of highest power of the series.
ex add_series(const pseries &other) const
Add one series object to another, producing a pseries object that represents the sum.
ex collect(const ex &s, bool distributed=false) const override
Does nothing.
ex eval() const override
Perform coefficient-wise automatic term rewriting rules in this class.
This class holds a relation consisting of two expressions and a logical relation between them.
@ expanded
.expand(0) has already done its job (other expand() options ignore this flag)
@ evaluated
.eval() has already done its job
@ no_pattern
disable pattern matching
bool is_equal_same_type(const basic &other) const override
Returns true if two objects of same type are equal.
ex series(const relational &s, int order, unsigned options=0) const override
Implementation of ex::series() for symbols.
Interface to GiNaC's initially known functions.
Interface to GiNaC's symbolic integral.
Definition of GiNaC's lst.
Interface to GiNaC's products of expressions.
const numeric pow(const numeric &x, const numeric &y)
std::map< ex, ex, ex_is_less > exmap
std::vector< expair > epvector
expair-vector
epvector * conjugateepvector(const epvector &epv)
Complex conjugate every element of an epvector.
bool are_ex_trivially_equal(const ex &e1, const ex &e2)
Compare two objects of class quickly without doing a deep tree traversal.
const numeric exp(const numeric &x)
Exponential function.
print_func< print_context >(&varidx::do_print). print_func< print_latex >(&varidx
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
const numeric log(const numeric &x)
Natural logarithm.
ex factor(const ex &poly, unsigned options)
Interface function to the outside world.
int to_int(const numeric &x)
bool is_integer(const numeric &x)
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.
std::vector< ex > exvector
bool is_order_function(const ex &e)
Check whether a function is the Order (O(n)) function.
Interface to GiNaC's overloaded operators.
Interface to GiNaC's symbolic exponentiation (basis^exponent).
Interface to class for extended truncated power series.
Interface to relations between expressions.
Interface to GiNaC's symbolic objects.
Interface to several small and furry utilities needed within GiNaC but not of any interest to the use...