author Richard Kreckel Wed, 15 Aug 2001 21:16:19 +0000 (21:16 +0000) committer Richard Kreckel Wed, 15 Aug 2001 21:16:19 +0000 (21:16 +0000)
- log_series(): check if the argument expands simply to Order(x^foo) and
if so expand somewhat further in order not to produce rubbish like
log(Order(x)^(-1)... and such.

 ginac/inifcns.cpp patch | blob | history ginac/inifcns_gamma.cpp patch | blob | history ginac/inifcns_trans.cpp patch | blob | history ginac/pseries.cpp patch | blob | history

index 50ec3ec8a97dd7f92cfd241f03a1d41021abd0df..ebf04514350a06646d4676a6c1253864c9b27d32 100644 (file)
@@ -315,16 +315,16 @@ static ex Li2_series(const ex &x, const relational &rel, int order, unsigned opt
// method:
// This is the branch cut: assemble the primitive series manually
// and then add the corresponding complex step function.
-                       const symbol *s = static_cast<symbol *>(rel.lhs().bp);
+                       const symbol &s = ex_to<symbol>(rel.lhs());
const ex point = rel.rhs();
const symbol foo;
epvector seq;
// zeroth order term:
seq.push_back(expair(Li2(x_pt), _ex0()));
// compute the intermediate terms:
-                       ex replarg = series(Li2(x), *s==foo, order);
+                       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,options).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);
@@ -420,8 +420,8 @@ static ex Order_series(const ex & x, const relational & r, int order, unsigned o
// Just wrap the function into a pseries object
epvector new_seq;
GINAC_ASSERT(is_ex_exactly_of_type(r.lhs(),symbol));
-       const symbol *s = static_cast<symbol *>(r.lhs().bp);
-       new_seq.push_back(expair(Order(_ex1()), numeric(std::min(x.ldegree(*s), order))));
+       const symbol &s = ex_to<symbol>(r.lhs());
+       new_seq.push_back(expair(Order(_ex1()), numeric(std::min(x.ldegree(s), order))));
return pseries(r, new_seq);
}

index 9c1da65157426b2cb960fd29f7f491c2012b57f6..ef0169b5c9ea1539b5c92c10499baac02c521e9b 100644 (file)
@@ -292,24 +292,24 @@ static ex beta_series(const ex & arg1,
const ex arg1_pt = arg1.subs(rel);
const ex arg2_pt = arg2.subs(rel);
GINAC_ASSERT(is_ex_exactly_of_type(rel.lhs(),symbol));
-       const symbol *s = static_cast<symbol *>(rel.lhs().bp);
+       const symbol &s = ex_to<symbol>(rel.lhs());
ex arg1_ser, arg2_ser, arg1arg2_ser;
if ((!arg1_pt.info(info_flags::integer) || arg1_pt.info(info_flags::positive)) &&
(!arg2_pt.info(info_flags::integer) || arg2_pt.info(info_flags::positive)))
throw do_taylor();  // caught by function::series()
// trap the case where arg1 is on a pole:
if (arg1.info(info_flags::integer) && !arg1.info(info_flags::positive))
-               arg1_ser = tgamma(arg1+*s).series(rel, order, options);
+               arg1_ser = tgamma(arg1+s).series(rel, order, options);
else
arg1_ser = tgamma(arg1).series(rel,order);
// trap the case where arg2 is on a pole:
if (arg2.info(info_flags::integer) && !arg2.info(info_flags::positive))
-               arg2_ser = tgamma(arg2+*s).series(rel, order, options);
+               arg2_ser = tgamma(arg2+s).series(rel, order, options);
else
arg2_ser = tgamma(arg2).series(rel,order);
// trap the case where arg1+arg2 is on a pole:
if ((arg1+arg2).info(info_flags::integer) && !(arg1+arg2).info(info_flags::positive))
-               arg1arg2_ser = tgamma(arg2+arg1+*s).series(rel, order, options);
+               arg1arg2_ser = tgamma(arg2+arg1+s).series(rel, order, options);
else
arg1arg2_ser = tgamma(arg2+arg1).series(rel,order);
// compose the result (expanding all the terms):
@@ -164,27 +164,35 @@ static ex log_series(const ex &arg,
// This is the branch point: Series expand the argument first, then
// trivially factorize it to isolate that part which has constant
// leading coefficient in this fashion:
-               //   x^n + Order(x^(n+m))  ->  x^n * (1 + Order(x^m)).
+               //   x^n + x^(n+1) +...+ Order(x^(n+m))  ->  x^n * (1 + x +...+ Order(x^m)).
// 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, options));
-               const symbol *s = static_cast<symbol *>(rel.lhs().bp);
+               pseries argser;          // series expansion of log's argument
+               unsigned extra_ord = 0;  // extra expansion order
+               do {
+                       // oops, the argument expanded to a pure Order(x^something)...
+                       argser = ex_to<pseries>(arg.series(rel, order+extra_ord, options));
+                       ++extra_ord;
+               } while (!argser.is_terminating() && argser.nops()==1);
+
+               const symbol &s = ex_to<symbol>(rel.lhs());
const ex point = rel.rhs();
-               const int n = argser.ldegree(*s);
+               const int n = argser.ldegree(s);
epvector seq;
// construct what we carelessly called the n*log(x) term above
-               ex coeff = argser.coeff(*s, n);
+               const ex coeff = argser.coeff(s, n);
// expand the log, but only if coeff is real and > 0, since otherwise
// it would make the branch cut run into the wrong direction
if (coeff.info(info_flags::positive))
-                       seq.push_back(expair(n*log(*s-point)+log(coeff), _ex0()));
+                       seq.push_back(expair(n*log(s-point)+log(coeff), _ex0()));
else
-                       seq.push_back(expair(log(coeff*pow(*s-point, n)), _ex0()));
+                       seq.push_back(expair(log(coeff*pow(s-point, n)), _ex0()));
+
if (!argser.is_terminating() || argser.nops()!=1) {
-                       // in this case n more terms are needed
+                       // in this case n more (or less) terms are needed
// (sadly, to generate them, we have to start from the beginning)
-                       ex newarg = ex_to<pseries>((arg/coeff).series(rel, order+n, options)).shift_exponents(-n).convert_to_poly(true);
+                       const ex newarg = ex_to<pseries>((arg/coeff).series(rel, order+n, options)).shift_exponents(-n).convert_to_poly(true);
} else  // it was a monomial
return pseries(rel, seq);
@@ -194,10 +202,10 @@ static ex log_series(const ex &arg,
// method:
// This is the branch cut: assemble the primitive series manually and
// then add the corresponding complex step function.
-               const symbol *s = static_cast<symbol *>(rel.lhs().bp);
+               const symbol &s = ex_to<symbol>(rel.lhs());
const ex point = rel.rhs();
const symbol foo;
-               const ex replarg = series(log(arg), *s==foo, order).subs(foo==point);
+               const ex replarg = series(log(arg), s==foo, order).subs(foo==point);
epvector seq;
seq.push_back(expair(-I*csgn(arg*I)*Pi, _ex0()));
seq.push_back(expair(Order(_ex1()), order));
@@ -639,10 +647,10 @@ static ex atan_series(const ex &arg,
// method:
// This is the branch cut: assemble the primitive series manually and
// then add the corresponding complex step function.
-               const symbol *s = static_cast<symbol *>(rel.lhs().bp);
+               const symbol &s = ex_to<symbol>(rel.lhs());
const ex point = rel.rhs();
const symbol foo;
-               const ex replarg = series(atan(arg), *s==foo, order).subs(foo==point);
+               const ex replarg = series(atan(arg), s==foo, order).subs(foo==point);
ex Order0correction = replarg.op(0)+csgn(arg)*Pi*_ex_1_2();
if ((I*arg_pt)<_ex0())
Order0correction += log((I*arg_pt+_ex_1())/(I*arg_pt+_ex1()))*I*_ex_1_2();
@@ -1024,10 +1032,10 @@ static ex atanh_series(const ex &arg,
// method:
// This is the branch cut: assemble the primitive series manually and
// then add the corresponding complex step function.
-               const symbol *s = static_cast<symbol *>(rel.lhs().bp);
+               const symbol &s = ex_to<symbol>(rel.lhs());
const ex point = rel.rhs();
const symbol foo;
-               const ex replarg = series(atanh(arg), *s==foo, order).subs(foo==point);
+               const ex replarg = series(atanh(arg), s==foo, order).subs(foo==point);
ex Order0correction = replarg.op(0)+csgn(I*arg)*Pi*I*_ex1_2();
if (arg_pt<_ex0())
Order0correction += log((arg_pt+_ex_1())/(arg_pt+_ex1()))*_ex1_2();
index 5812982447a68816d4a5d956efc0cae8fec124e7..0ec5a61a49bb7e3f4380872f10c2c8eed21d4632 100644 (file)
@@ -77,10 +77,10 @@ DEFAULT_DESTROY(pseries)
pseries::pseries(const ex &rel_, const epvector &ops_) : basic(TINFO_pseries), seq(ops_)
{
debugmsg("pseries ctor from ex,epvector", LOGLEVEL_CONSTRUCT);
-       GINAC_ASSERT(is_ex_exactly_of_type(rel_, relational));
-       GINAC_ASSERT(is_ex_exactly_of_type(rel_.lhs(),symbol));
+       GINAC_ASSERT(is_exactly_a<relational>(rel_));
+       GINAC_ASSERT(is_exactly_a<symbol>(rel_.lhs()));
point = rel_.rhs();
-       var = *static_cast<symbol *>(rel_.lhs().bp);
+       var = rel_.lhs();
}

@@ -330,7 +330,7 @@ ex pseries::coeff(const ex &s, int n) const
int lo = 0, hi = seq.size() - 1;
while (lo <= hi) {
int mid = (lo + hi) / 2;
-                       GINAC_ASSERT(is_ex_exactly_of_type(seq[mid].coeff, numeric));
+                       GINAC_ASSERT(is_exactly_a<numeric>(seq[mid].coeff));
int cmp = ex_to<numeric>(seq[mid].coeff).compare(looking_for);
switch (cmp) {
case -1:
@@ -492,7 +492,7 @@ ex basic::series(const relational & r, int order, unsigned options) const
numeric fac(1);
ex deriv = *this;
ex coeff = deriv.subs(r);
-       const symbol &s = static_cast<symbol &>(*r.lhs().bp);
+       const symbol &s = ex_to<symbol>(r.lhs());

if (!coeff.is_zero())
seq.push_back(expair(coeff, _ex0()));
@@ -524,10 +524,9 @@ ex symbol::series(const relational & r, int order, unsigned options) const
{
epvector seq;
const ex point = r.rhs();
-       GINAC_ASSERT(is_ex_exactly_of_type(r.lhs(),symbol));
-       ex s = r.lhs();
-
-       if (this->is_equal(*s.bp)) {
+       GINAC_ASSERT(is_exactly_a<symbol>(r.lhs()));
+
+       if (this->is_equal_same_type(ex_to<symbol>(r.lhs()))) {
if (order > 0 && !point.is_zero())
seq.push_back(expair(point, _ex0()));
if (order > 1)
@@ -875,8 +874,8 @@ ex power::series(const relational & r, int order, unsigned options) 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));
-       const symbol &s = static_cast<symbol &>(*r.lhs().bp);
+       GINAC_ASSERT(is_exactly_a<symbol>(r.lhs()));
+       const symbol &s = ex_to<symbol>(r.lhs());

if (var.is_equal(s) && point.is_equal(p)) {
if (order > degree(s))