int ldegree(const symbol & s) const;
ex coeff(const symbol & s, int n=1) const;
ex eval(int level=0) const;
- ex series(const relational & r, int order, bool branchcut = true) const;
+ ex series(const relational & r, int order, unsigned options = 0) const;
ex normal(lst &sym_lst, lst &repl_lst, int level=0) const;
numeric integer_content(void) const;
ex smod(const numeric &xi) const;
virtual ex collect(const symbol & s) const;
virtual ex eval(int level = 0) const;
virtual ex evalf(int level = 0) const;
- virtual ex series(const relational & r, int order, bool branchcut = true) const;
+ virtual ex series(const relational & r, int order, unsigned options = 0) const;
virtual ex subs(const lst & ls, const lst & lr) const;
virtual ex normal(lst &sym_lst, lst &repl_lst, int level = 0) const;
virtual ex to_rational(lst &repl_lst) const;
ex eval(int level = 0) const;
ex evalf(int level = 0) const;
ex diff(const symbol & s, unsigned nth = 1) const;
- ex series(const ex & r, int order, bool branchcut = true) const;
+ ex series(const ex & r, int order, unsigned options = 0) const;
ex subs(const lst & ls, const lst & lr) const;
ex subs(const ex & e) const;
exvector get_indices(void) const;
inline ex diff(const ex & thisex, const symbol & s, unsigned nth = 1)
{ return thisex.diff(s, nth); }
-inline ex series(const ex & thisex, const ex & r, int order, bool branchcut = true)
-{ return thisex.series(r, order, branchcut); }
+inline ex series(const ex & thisex, const ex & r, int order, unsigned options = 0)
+{ return thisex.series(r, order, options); }
inline ex subs(const ex & thisex, const ex & e)
{ return thisex.subs(e); }
};
};
+class series_options {
+public:
+ enum { suppress_branchcut = 0x0001
+ };
+};
+
+class determinant_options {
+public:
+ enum { laplace,
+ bareiss
+ };
+};
+
class status_flags {
public:
enum { dynallocated = 0x0001,
enum { delete_never, // let table grow undefinitely, not recommmended, but currently default
delete_lru, // least recently used
delete_lfu, // least frequently used
- delete_cyclic // first one in list (oldest)
+ delete_cyclic // first (oldest) one in list
};
};
'const ex &','');
$typedef_series_funcp=generate(
-'typedef ex (* series_funcp_${N})(${SEQ1}, const relational &, int, bool);'."\n",
+'typedef ex (* series_funcp_${N})(${SEQ1}, const relational &, int, unsigned);'."\n",
'const ex &','');
$eval_func_interface=generate(' function_options & eval_func(eval_funcp_${N} e);'."\n",'','');
<<'END_OF_SERIES_SWITCH_STATEMENT','seq[${N}-1]','');
case ${N}:
try {
- res = ((series_funcp_${N})(registered_functions()[serial].series_f))(${SEQ1},r,order,branchcut);
+ res = ((series_funcp_${N})(registered_functions()[serial].series_f))(${SEQ1},r,order,options);
} catch (do_taylor) {
- res = basic::series(r, order, branchcut);
+ res = basic::series(r, order, options);
}
return res;
break;
ex expand(unsigned options=0) const;
ex eval(int level=0) const;
ex evalf(int level=0) const;
- ex series(const relational & r, int order, bool branchcut = true) const;
+ ex series(const relational & r, int order, unsigned options = 0) const;
ex thisexprseq(const exvector & v) const;
ex thisexprseq(exvector * vp) const;
protected:
/** Implementation of ex::series for functions.
* \@see ex::series */
-ex function::series(const relational & r, int order, bool branchcut = true) const
+ex function::series(const relational & r, int order, unsigned options = 0) const
{
GINAC_ASSERT(serial<registered_functions().size());
static ex csgn_series(const ex & arg,
const relational & rel,
int order,
- bool branchcut)
+ unsigned options)
{
const ex arg_pt = arg.subs(rel);
if (arg_pt.info(info_flags::numeric)) {
return -log(1-x)/x;
}
-static ex Li2_series(const ex &x, const relational &rel, int order, bool branchcut)
+static ex Li2_series(const ex &x, const relational &rel, int order, unsigned options)
{
const ex x_pt = x.subs(rel);
if (x_pt.info(info_flags::numeric)) {
return ser.series(rel, order);
}
// third special case: x real, >=1 (branch cut)
- if (ex_to_numeric(x_pt).is_real() && ex_to_numeric(x_pt)>1) {
+ if (!(options & series_options::suppress_branchcut) &&
+ ex_to_numeric(x_pt).is_real() && ex_to_numeric(x_pt)>1) {
// method:
// This is the branch cut: assemble the primitive series manually
// and then add the corresponding complex step function.
// compute the intermediate terms:
ex replarg = series(Li2(x), *s==foo, order);
for (unsigned i=1; i<replarg.nops()-1; ++i)
- seq.push_back(expair((replarg.op(i)/power(*s-foo,i)).series(foo==point,1,branchcut).op(0).subs(foo==*s),i));
+ seq.push_back(expair((replarg.op(i)/power(*s-foo,i)).series(foo==point,1,options).op(0).subs(foo==*s),i));
// append an order term:
seq.push_back(expair(Order(_ex1()), replarg.nops()-1));
return pseries(rel, seq);
return Order(x).hold();
}
-static ex Order_series(const ex & x, const relational & r, int order, bool branchcut)
+static ex Order_series(const ex & x, const relational & r, int order, unsigned options)
{
// Just wrap the function into a pseries object
epvector new_seq;
static ex lgamma_series(const ex & arg,
const relational & rel,
int order,
- bool branchcut)
+ unsigned options)
{
// method:
// Taylor series where there is no pole falls back to psi function
static ex tgamma_series(const ex & arg,
const relational & rel,
int order,
- bool branchcut)
+ unsigned options)
{
// method:
// Taylor series where there is no pole falls back to psi function
const ex & arg2,
const relational & rel,
int order,
- bool branchcut)
+ unsigned options)
{
// method:
// Taylor series where there is no pole of one of the tgamma functions
static ex psi1_series(const ex & arg,
const relational & rel,
int order,
- bool branchcut)
+ unsigned options)
{
// method:
// Taylor series where there is no pole falls back to polygamma function
const ex & arg,
const relational & rel,
int order,
- bool branchcut)
+ unsigned options)
{
// method:
// Taylor series where there is no pole falls back to polygamma function
static ex log_series(const ex &arg,
const relational &rel,
int order,
- bool branchcut)
+ unsigned options)
{
GINAC_ASSERT(is_ex_exactly_of_type(rel.lhs(),symbol));
ex arg_pt;
// Return a plain n*log(x) for the x^n part and series expand the
// other part. Add them together and reexpand again in order to have
// one unnested pseries object. All this also works for negative n.
- const pseries argser = ex_to_pseries(arg.series(rel, order, branchcut));
+ const pseries argser = ex_to_pseries(arg.series(rel, order, options));
const symbol *s = static_cast<symbol *>(rel.lhs().bp);
const ex point = rel.rhs();
const int n = argser.ldegree(*s);
seq.push_back(expair(n*log(*s-point), _ex0()));
if (!argser.is_terminating() || argser.nops()!=1) {
// in this case n more terms are needed
- ex newarg = ex_to_pseries(arg.series(rel, order+n, branchcut)).shift_exponents(-n).convert_to_poly(true);
- return pseries(rel, seq).add_series(ex_to_pseries(log(newarg).series(rel, order, branchcut)));
+ ex newarg = ex_to_pseries(arg.series(rel, order+n, options)).shift_exponents(-n).convert_to_poly(true);
+ return pseries(rel, seq).add_series(ex_to_pseries(log(newarg).series(rel, order, options)));
} else // it was a monomial
return pseries(rel, seq);
}
- if (branchcut && arg_pt.info(info_flags::negative)) {
+ if (!(options & series_options::suppress_branchcut) &&
+ arg_pt.info(info_flags::negative)) {
// method:
// This is the branch cut: assemble the primitive series manually and
// then add the corresponding complex step function.
static ex tan_series(const ex &x,
const relational &rel,
int order,
- bool branchcut)
+ unsigned options)
{
// method:
// Taylor series where there is no pole falls back to tan_deriv.
static ex tanh_series(const ex &x,
const relational &rel,
int order,
- bool branchcut)
+ unsigned options)
{
// method:
// Taylor series where there is no pole falls back to tanh_deriv.
exvector prod(row*other.col);
- for (unsigned r1=0; r1<row; ++r1) {
- for (unsigned c=0; c<col; ++c) {
+ for (unsigned r1=0; r1<rows(); ++r1) {
+ for (unsigned c=0; c<cols(); ++c) {
if (m[r1*col+c].is_zero())
continue;
for (unsigned r2=0; r2<other.col; ++r2)
/** operator() to access elements.
*
* @param ro row of element
- * @param co column of element
+ * @param co column of element
* @exception range_error (index out of range) */
const ex & matrix::operator() (unsigned ro, unsigned co) const
{
if (ro<0 || ro>=row || co<0 || co>=col)
throw (std::range_error("matrix::operator(): index out of range"));
-
+
return m[ro*col+co];
}
throw (std::runtime_error("matrix::inverse(): singular matrix"));
}
if (indx != 0) { // swap rows r and indx of matrix tmp
- for (unsigned i=0; i<col; ++i) {
+ for (unsigned i=0; i<col; ++i)
tmp.m[r1*col+i].swap(tmp.m[indx*col+i]);
- }
}
ex a1 = cpy.m[r1*col+r1];
for (unsigned c=0; c<col; ++c) {
return tmp;
}
+// superfluous helper function, to be removed:
+void matrix::swap(unsigned r1, unsigned c1, unsigned r2 ,unsigned c2)
+{
+ ensure_if_modifiable();
+
+ ex tmp = (*this)(r1,c1);
+ set(r1,c1,(*this)(r2,c2));
+ set(r2,c2,tmp);
+}
+
/** Solve a set of equations for an m x n matrix by fraction-free Gaussian
* elimination. Based on algorithm 9.1 from 'Algorithms for Computer Algebra'
for (unsigned k=0; (k<n)&&(r<m); ++k) {
// find a nonzero pivot
unsigned p;
- for (p=r; (p<m)&&(a.m[p*a.cols()+k].is_zero()); ++p) {}
+ for (p=r; (p<m)&&(a(p,k).is_zero()); ++p) {}
// pivot is in row p
if (p<m) {
if (p!=r) {
// swap rows p and r
for (unsigned j=k; j<n; ++j)
- a.m[p*a.cols()+j].swap(a.m[r*a.cols()+j]);
- a.m[p*a.cols()].swap(a.m[r*a.cols()]);
+ a.swap(p,j,r,j);
+ b.swap(p,0,r,0);
// keep track of sign changes due to row exchange
- sign = -sign;
+ sign *= -1;
}
for (unsigned i=r+1; i<m; ++i) {
for (unsigned j=k+1; j<n; ++j) {
- a.set(i,j,(a.m[r*a.cols()+k]*a.m[i*a.cols()+j]
- -a.m[r*a.cols()+j]*a.m[i*a.cols()+k])/divisor);
- a.set(i,j,a.m[i*a.cols()+j].normal());
+ a.set(i,j,(a(r,k)*a(i,j)
+ -a(r,j)*a(i,k))/divisor);
+ a.set(i,j,a(i,j).normal());
}
- b.set(i,0,(a.m[r*a.cols()+k]*b.m[i*b.cols()]
- -b.m[r*b.cols()]*a.m[i*a.cols()+k])/divisor);
- b.set(i,0,b.m[i*b.cols()].normal());
+ b.set(i,0,(a(r,k)*b(i,0)
+ -b(r,0)*a(i,k))/divisor);
+ b.set(i,0,b(i,0).normal());
a.set(i,k,0);
}
- divisor = a.m[r*a.cols()+k];
- r++;
+ divisor = a(r,k);
+ ++r;
}
}
for (unsigned r=0; r<m; ++r) {
int zero_in_this_row=0;
for (unsigned c=0; c<n; ++c) {
- if (a.m[r*a.cols()+c].is_zero())
+ if (a(r,c).is_zero())
zero_in_this_row++;
else
break;
unsigned last_assigned_sol = n+1;
for (int r=m-1; r>=0; --r) {
unsigned first_non_zero = 1;
- while ((first_non_zero<=n)&&(a.m[r*a.cols()+(first_non_zero-1)].is_zero()))
+ while ((first_non_zero<=n)&&(a(r,first_non_zero-1).is_zero()))
first_non_zero++;
if (first_non_zero>n) {
// row consists only of zeroes, corresponding rhs must be 0 as well
- if (!b.m[r*b.cols()].is_zero()) {
+ if (!b(r,0).is_zero()) {
throw (std::runtime_error("matrix::fraction_free_elim(): singular matrix"));
}
} else {
// assign solutions for vars between first_non_zero+1 and
// last_assigned_sol-1: free parameters
for (unsigned c=first_non_zero; c<last_assigned_sol-1; ++c)
- sol.set(c,0,vars.m[c*vars.cols()]);
- ex e = b.m[r*b.cols()];
+ sol.set(c,0,vars(c,0));
+ ex e = b(r,0);
for (unsigned c=first_non_zero; c<n; ++c)
- e -= a.m[r*a.cols()+c]*sol.m[c*sol.cols()];
+ e -= a(r,c)*sol(c,0);
sol.set(first_non_zero-1,0,
- (e/a.m[r*a.cols()+(first_non_zero-1)]).normal());
+ (e/a(r,first_non_zero-1)).normal());
last_assigned_sol = first_non_zero;
}
}
// assign solutions for vars between 1 and
// last_assigned_sol-1: free parameters
for (unsigned c=0; c<last_assigned_sol-1; ++c)
- sol.set(c,0,vars.m[c*vars.cols()]);
+ sol.set(c,0,vars(c,0));
#ifdef DO_GINAC_ASSERT
// test solution with echelon matrix
for (unsigned r=0; r<m; ++r) {
ex e = 0;
for (unsigned c=0; c<n; ++c)
- e += a.m[r*a.cols()+c]*sol.m[c*sol.cols()];
- if (!(e-b.m[r*b.cols()]).normal().is_zero()) {
+ e += a(r,c)*sol(c,0);
+ if (!(e-b(r,0)).normal().is_zero()) {
cout << "e=" << e;
- cout << "b(" << r <<",0)=" << b.m[r*b.cols()] << endl;
- cout << "diff=" << (e-b.m[r*b.cols()]).normal() << endl;
+ cout << "b(" << r <<",0)=" << b(r,0) << endl;
+ cout << "diff=" << (e-b(r,0)).normal() << endl;
}
- GINAC_ASSERT((e-b.m[r*b.cols()]).normal().is_zero());
+ GINAC_ASSERT((e-b(r,0)).normal().is_zero());
}
// test solution with original matrix
for (unsigned r=0; r<m; ++r) {
ex e = 0;
for (unsigned c=0; c<n; ++c)
- e += this->m[r*cols()+c]*sol.m[c*sol.cols()];
+ e += this->m[r*cols()+c]*sol(c,0);
try {
- if (!(e-rhs.m[r*rhs.cols()]).normal().is_zero()) {
+ if (!(e-rhs(r,0)).normal().is_zero()) {
cout << "e==" << e << endl;
e.printtree(cout);
ex en = e.normal();
cout << "e.normal()=" << en << endl;
en.printtree(cout);
- cout << "rhs(" << r <<",0)=" << rhs.m[r*rhs.cols()] << endl;
- cout << "diff=" << (e-rhs.m[r*rhs.cols()]).normal() << endl;
+ cout << "rhs(" << r <<",0)=" << rhs(r,0) << endl;
+ cout << "diff=" << (e-rhs(r,0)).normal() << endl;
}
} catch (...) {
- ex xxx = e - rhs.m[r*rhs.cols()];
+ ex xxx = e - rhs(r,0);
cerr << "xxx=" << xxx << endl << endl;
}
- GINAC_ASSERT((e-rhs.m[r*rhs.cols()]).normal().is_zero());
+ GINAC_ASSERT((e-rhs(r,0)).normal().is_zero());
}
#endif // def DO_GINAC_ASSERT
#include "basic.h"
#include "ex.h"
+namespace std {
+ // forward declaration, so <stdexcept> need not be included:
+ class range_error;
+}
+
#ifndef NO_NAMESPACE_GINAC
namespace GiNaC {
#endif // ndef NO_NAMESPACE_GINAC
int division_free_elimination(void);
int fraction_free_elimination(bool det = false);
int pivot(unsigned ro, bool symbolic=true);
+ void swap(unsigned r1, unsigned c1, unsigned r2 ,unsigned c2);
// member variables
protected:
static unsigned precedence;
};
+
// global constants
extern const matrix some_matrix;
extern const type_info & typeid_matrix;
+
// wrapper functions around member functions
inline unsigned nops(const matrix & m)
ex coeff(const symbol & s, int n = 1) const;
ex eval(int level=0) const;
ex evalf(int level=0) const;
- ex series(const relational & s, int order, bool branchcut = true) const;
+ ex series(const relational & s, int order, unsigned options = 0) const;
ex normal(lst &sym_lst, lst &repl_lst, int level = 0) const;
numeric integer_content(void) const;
ex smod(const numeric &xi) const;
ex coeff(const symbol & s, int n = 1) const;
ex eval(int level=0) const;
ex evalf(int level=0) const;
- ex series(const relational & s, int order, bool branchcut = true) const;
+ ex series(const relational & s, int order, unsigned options = 0) const;
ex subs(const lst & ls, const lst & lr) const;
ex normal(lst &sym_lst, lst &repl_lst, int level = 0) const;
ex to_rational(lst &repl_lst) const;
/** Default implementation of ex::series(). This performs Taylor expansion.
* @see ex::series */
-ex basic::series(const relational & r, int order, bool branchcut) const
+ex basic::series(const relational & r, int order, unsigned options) const
{
epvector seq;
numeric fac(1);
/** Implementation of ex::series() for symbols.
* @see ex::series */
-ex symbol::series(const relational & r, int order, bool branchcut) const
+ex symbol::series(const relational & r, int order, unsigned options) const
{
epvector seq;
const ex point = r.rhs();
/** Implementation of ex::series() for sums. This performs series addition when
* adding pseries objects.
* @see ex::series */
-ex add::series(const relational & r, int order, bool branchcut) const
+ex add::series(const relational & r, int order, unsigned options) const
{
ex acc; // Series accumulator
// Get first term from overall_coeff
- acc = overall_coeff.series(r, order, branchcut);
+ acc = overall_coeff.series(r, order, options);
// Add remaining terms
epvector::const_iterator it = seq.begin();
if (is_ex_exactly_of_type(it->rest, pseries))
op = it->rest;
else
- op = it->rest.series(r, order, branchcut);
+ op = it->rest.series(r, order, options);
if (!it->coeff.is_equal(_ex1()))
op = ex_to_pseries(op).mul_const(ex_to_numeric(it->coeff));
/** Implementation of ex::series() for product. This performs series
* multiplication when multiplying series.
* @see ex::series */
-ex mul::series(const relational & r, int order, bool branchcut) const
+ex mul::series(const relational & r, int order, unsigned options) const
{
ex acc; // Series accumulator
// Get first term from overall_coeff
- acc = overall_coeff.series(r, order, branchcut);
+ acc = overall_coeff.series(r, order, options);
// Multiply with remaining terms
epvector::const_iterator it = seq.begin();
acc = ex_to_pseries(acc).mul_const(ex_to_numeric(f));
continue;
} else if (!is_ex_exactly_of_type(op, pseries))
- op = op.series(r, order, branchcut);
+ op = op.series(r, order, options);
if (!it->coeff.is_equal(_ex1()))
op = ex_to_pseries(op).power_const(ex_to_numeric(it->coeff), order);
/** Implementation of ex::series() for powers. This performs Laurent expansion
* of reciprocals of series at singularities.
* @see ex::series */
-ex power::series(const relational & r, int order, bool branchcut) const
+ex power::series(const relational & r, int order, unsigned options) const
{
ex e;
if (!is_ex_exactly_of_type(basis, pseries)) {
// Basis is not a series, may there be a singulary?
if (!exponent.info(info_flags::negint))
- return basic::series(r, order, branchcut);
+ return basic::series(r, order, options);
// Expression is of type something^(-int), check for singularity
if (!basis.subs(r).is_zero())
- return basic::series(r, order, branchcut);
+ return basic::series(r, order, options);
// Singularity encountered, expand basis into series
- e = basis.series(r, order, branchcut);
+ e = basis.series(r, order, options);
} else {
// Basis is a series
e = basis;
/** Re-expansion of a pseries object. */
-ex pseries::series(const relational & r, int order, bool branchcut) const
+ex pseries::series(const relational & r, int order, unsigned options) const
{
const ex p = r.rhs();
GINAC_ASSERT(is_ex_exactly_of_type(r.lhs(),symbol));
return pseries(r, new_seq);
}
} else
- return convert_to_poly().series(r, order, branchcut);
+ return convert_to_poly().series(r, order, options);
}
*
* @param r expansion relation, lhs holds variable and rhs holds point
* @param order truncation order of series calculations
- * @param branchcut when set to false, branch cuts are not honored
+ * @param options of class series_options
* @return an expression holding a pseries object */
-ex ex::series(const ex & r, int order, bool branchcut) const
+ex ex::series(const ex & r, int order, unsigned options) const
{
GINAC_ASSERT(bp!=0);
ex e;
throw (std::logic_error("ex::series(): expansion point has unknown type"));
try {
- e = bp->series(rel_, order, branchcut);
+ e = bp->series(rel_, order, options);
} catch (std::exception &x) {
throw (std::logic_error(std::string("unable to compute series (") + x.what() + ")"));
}
ex collect(const symbol &s) const;
ex eval(int level=0) const;
ex evalf(int level=0) const;
- ex series(const relational & r, int order, bool branchcut = true) const;
+ ex series(const relational & r, int order, unsigned options = 0) const;
ex subs(const lst & ls, const lst & lr) const;
ex normal(lst &sym_lst, lst &repl_lst, int level = 0) const;
ex expand(unsigned options = 0) const;
int ldegree(const symbol & s) const;
ex coeff(const symbol & s, int n = 1) const;
ex eval(int level = 0) const;
- ex series(const relational & s, int order, bool branchcut = true) const;
+ ex series(const relational & s, int order, unsigned options = 0) const;
ex normal(lst &sym_lst, lst &repl_lst, int level = 0) const;
ex to_rational(lst &repl_lst) const;
ex subs(const lst & ls, const lst & lr) const;
git://www.ginac.de / ginac.git / commitdiff