74 #ifdef DO_GINAC_ASSERT 76 while (i !=
seq.end()) {
83 GINAC_ASSERT(ex_to<numeric>(i->coeff) < ex_to<numeric>(ip1->coeff));
86 #endif // def DO_GINAC_ASSERT 93 : seq(
std::move(ops_))
95 #ifdef DO_GINAC_ASSERT 97 while (i !=
seq.end()) {
104 GINAC_ASSERT(ex_to<numeric>(i->coeff) < ex_to<numeric>(ip1->coeff));
107 #endif // def DO_GINAC_ASSERT 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) {
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] ) {
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"));
const numeric exp(const numeric &x)
Exponential function.
bool is_negative() const
True if object is not complex and less than zero.
ex power_const(const numeric &p, int deg) const
Compute the p-th power of a series.
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(const print_context &c, unsigned level) const
Interface to GiNaC's symbolic exponentiation (basis^exponent).
virtual ex coeff(const ex &s, int n=1) const
Return coefficient of degree n in object s.
ex expand(unsigned options=0) const override
Implementation of ex::expand() for a power series.
ex coeff(const ex &s, int n=1) const
Interface to GiNaC's symbolic objects.
unsigned hashvalue
hash value
ex mul_const(const numeric &other) const
Multiply a pseries object with a numeric constant, producing a pseries object that represents the pro...
bool is_equal_same_type(const basic &other) const override
Returns true if two objects of same type are equal.
int to_int(const numeric &x)
int compare(const ex &other) const
int degree(const ex &s) const override
Return degree of highest power of the series.
ex subs(const exmap &m, unsigned options=0) const
void print_series(const print_context &c, const char *openbrace, const char *closebrace, const char *mul_sym, const char *pow_sym, unsigned level) const
virtual bool has(const ex &other, unsigned options=0) const
Test for occurrence of a pattern.
const basic & hold() const
Stop further evaluation.
ex conjugate() const override
ex derivative(const symbol &s) const override
Implementation of ex::diff() for a power series.
void do_print_latex(const print_latex &c, unsigned level) const
bool is_equal(const ex &other) const
virtual ex series(const relational &r, int order, unsigned options=0) const
Default implementation of ex::series().
pseries(const ex &rel_, const epvector &ops_)
Construct pseries from a vector of coefficients and powers.
const basic & setflag(unsigned f) const
Set some status_flags.
ex expand(unsigned options=0) const
ex exponop(size_t i) const
integral(const ex &x_, const ex &a_, const ex &b_, const ex &f_)
ex real_part() const override
ex evalf() const override
Evaluate coefficients numerically.
This class holds a relation consisting of two expressions and a logical relation between them...
Archiving of GiNaC expressions.
Interface to GiNaC's sums of expressions.
This class is the ABC (abstract base class) of GiNaC's class hierarchy.
ex imag_part() const override
Context for python-parsable output.
bool are_ex_trivially_equal(const ex &e1, const ex &e2)
Compare two objects of class quickly without doing a deep tree traversal.
void do_print_python_repr(const print_python_repr &c, unsigned level) const
Interface to GiNaC's products of expressions.
bool is_terminating() const
Returns true if there is no order term, i.e.
ex convert_to_poly(bool no_order=false) const
Convert the pseries object to an ordinary polynomial.
pseries shift_exponents(int deg) const
Return a new pseries object with the powers shifted by deg.
ex collect(const ex &s, bool distributed=false) const override
Does nothing.
bool is_integer(const numeric &x)
bool info(unsigned inf) const
Interface to several small and furry utilities needed within GiNaC but not of any interest to the use...
Context for python pretty-print output.
ex series(const relational &r, int order, unsigned options=0) const override
Re-expansion of a pseries object.
ex recombine_pair_to_ex(const expair &p) const override
int ldegree(const ex &s) const
int compare(const basic &other) const
Compare objects syntactically to establish canonical ordering.
void do_print_python(const print_python &c, unsigned level) const
ex series(const relational &s, int order, unsigned options=0) const override
Implementation of ex::series() for product.
const_iterator begin() const noexcept
This class is a wrapper around CLN-numbers within the GiNaC class hierarchy.
Context for latex-parsable output.
bool is_order_function(const ex &e)
Check whether a function is the Order (O(n)) function.
print_func< print_context >(&varidx::do_print). print_func< print_latex >(&varidx
ex eval() const override
Perform coefficient-wise automatic term rewriting rules in this class.
void do_print_tree(const print_tree &c, unsigned level) const
bool find(const ex &pattern, exset &found) const
Find all occurrences of a pattern.
ex coeffop(size_t i) const
Get coefficients and exponents.
ex rhs() const
Right hand side of relational expression.
ex eval_integ() const override
Evaluate integrals, if result is known.
const numeric div(const numeric &other) const
Numerical division method.
#define GINAC_ASSERT(X)
Assertion macro for checking invariances.
ex var
Series variable (holds a symbol)
This class holds a extended truncated power series (positive and negative integer powers)...
epvector * conjugateepvector(const epvector &)
Complex conjugate every element of an epvector.
bool is_compatible_to(const pseries &other) const
Check whether series is compatible to another series (expansion variable and point are the same...
ex coeff(const ex &s, int n=1) const override
Return coefficient of degree n in object s.
ex evalm() const override
Evaluate sums, products and integer powers of matrices.
const numeric log(const numeric &x)
Natural logarithm.
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 real() const
Real part of a number.
ex diff(const symbol &s, unsigned nth=1) const
Compute partial derivative of an expression.
std::vector< ex > exvector
std::map< ex, ex, ex_is_less > exmap
ex add_series(const pseries &other) const
Add one series object to another, producing a pseries object that represents the sum.
Interface to GiNaC's overloaded operators.
void read_archive(const archive_node &n, lst &syms) override
Read (a.k.a.
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).
Interface to GiNaC's symbolic integral.
ex factor(const ex &poly, unsigned options)
Interface function to the outside world.
Base class for print_contexts.
ex series(const ex &r, int order, unsigned options=0) const
Compute the truncated series expansion of an expression.
Definition of GiNaC's lst.
ex op(size_t i) const override
Return operand/member at position i.
const numeric pow(const numeric &x, const numeric &y)
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.
size_t nops() const override
Return the number of operands including a possible order term.
ex series(const relational &r, int order, unsigned options=0) const override
Default implementation of ex::series().
bool is_zero() const
True if object is zero.
void print(const print_context &c, unsigned level=0) const
Print expression to stream.
This class stores all properties needed to record/retrieve the state of one object of class basic (or...
Lightweight wrapper for GiNaC's symbolic objects.
ex series(const relational &r, int order, unsigned options=0) const override
Implementation of ex::series() for sums.
Interface to GiNaC's initially known functions.
virtual void print(const print_context &c, unsigned level=0) const
Output to stream.
int ldegree(const ex &s) const override
Return degree of lowest power of the series.
Interface to class for extended truncated power series.
virtual int compare_same_type(const basic &other) const
Returns order relation between two objects of same type.
Wrapper template for making GiNaC classes out of STL containers.
.expand(0) has already done its job (other expand() options ignore this flag)
Interface to relations between expressions.
ex mul_series(const pseries &other) const
Multiply one pseries object to another, producing a pseries object that represents the product...
void archive(archive_node &n) const override
Save (a.k.a.
std::vector< expair > epvector
expair-vector
ptr< basic > bp
pointer to basic object managed by this
bool is_zero() const
Check whether series has the value zero.
epvector seq
Vector of {coefficient, power} pairs.
ex lhs() const
Left hand side of relational expression.
ex subs(const exmap &m, unsigned options=0) const override
Substitute a set of objects by arbitrary expressions.
.eval() has already done its job
ex series(const relational &s, int order, unsigned options=0) const override
Implementation of ex::series() for symbols.
Context for tree-like output for debugging.
unsigned flags
of type status_flags
GINAC_BIND_UNARCHIVER(add)