From 955cb185a85535ab328ffedbfccdc508ce80fa91 Mon Sep 17 00:00:00 2001 From: Richard Kreckel Date: Fri, 17 Dec 1999 19:58:25 +0000 Subject: [PATCH] - Banned exZERO(), exONE(), exMINUSHALF() and all this from the interface. There is now _ex1() meaning 1, _ex_1() meaning -1, _ex1_2() meaning 1/2 and so on defined in utils.h and implemented in utils.cpp. Feel free to extend them as it pleases you but use them inside the library only. - Added more evaluations of trigonometric functions - Added series expansions for psi(x) and psi(n,x) at all their poles --- check/differentiation.cpp | 8 +- check/normalization.cpp | 2 +- check/numeric_consist.cpp | 2 +- check/paranoia_check.cpp | 14 +- check/poly_gcd.cpp | 21 +- check/series_expansion.cpp | 36 ++-- ginac/add.cpp | 59 +++--- ginac/basic.cpp | 2 +- ginac/color.cpp | 53 ++--- ginac/diff.cpp | 21 +- ginac/ex.cpp | 58 +----- ginac/ex.h | 20 +- ginac/expair.h | 14 +- ginac/expairseq.cpp | 8 +- ginac/inifcns.cpp | 11 +- ginac/inifcns.h | 2 +- ginac/inifcns_gamma.cpp | 127 ++++++++---- ginac/inifcns_trans.cpp | 311 ++++++++++++++------------- ginac/inifcns_zeta.cpp | 13 +- ginac/lortensor.cpp | 20 +- ginac/matrix.cpp | 17 +- ginac/mul.cpp | 57 ++--- ginac/ncmul.cpp | 7 +- ginac/normal.cpp | 117 +++++------ ginac/numeric.cpp | 128 ++++-------- ginac/numeric.h | 12 -- ginac/operators.cpp | 120 +---------- ginac/power.cpp | 109 ++++------ ginac/power.h | 3 +- ginac/relational.cpp | 3 +- ginac/series.cpp | 39 ++-- ginac/simp_lor.cpp | 31 +-- ginac/symbol.cpp | 6 +- ginac/utils.cpp | 415 ++++++++++++++++++++++++++++++++++++- ginac/utils.h | 70 ++++++- ginsh/ginsh_parser.yy | 8 +- 36 files changed, 1137 insertions(+), 807 deletions(-) diff --git a/check/differentiation.cpp b/check/differentiation.cpp index 348b5fd5..ece8bc37 100644 --- a/check/differentiation.cpp +++ b/check/differentiation.cpp @@ -30,7 +30,7 @@ static unsigned check_diff(const ex &e, const symbol &x, const ex &d, unsigned nth=1) { ex ed = e.diff(x, nth); - if ((ed - d).compare(exZERO()) != 0) { + if ((ed - d).compare(ex(0)) != 0) { switch (nth) { case 0: clog << "zeroth "; @@ -258,13 +258,13 @@ static unsigned differentiation6(void) symbol x("x"); ex e, d, ed; - e = sin(x).series(x, exZERO(), 8); - d = cos(x).series(x, exZERO(), 7); + e = sin(x).series(x, 0, 8); + d = cos(x).series(x, 0, 7); ed = e.diff(x); ed = static_cast(ed.bp)->convert_to_poly(); d = static_cast(d.bp)->convert_to_poly(); - if ((ed - d).compare(exZERO()) != 0) { + if ((ed - d).compare(ex(0)) != 0) { clog << "derivative of " << e << " by " << x << " returned " << ed << " instead of " << d << ")" << endl; return 1; diff --git a/check/normalization.cpp b/check/normalization.cpp index 7ce1e5dc..78ae6317 100644 --- a/check/normalization.cpp +++ b/check/normalization.cpp @@ -46,7 +46,7 @@ static unsigned normal1(void) // Expansion e = pow(x, 2) - (x+1)*(x-1) - 1; - d = exZERO(); + d = ex(0); result += check_normal(e, d); // Expansion inside functions diff --git a/check/numeric_consist.cpp b/check/numeric_consist.cpp index d7edb3b6..494cc579 100644 --- a/check/numeric_consist.cpp +++ b/check/numeric_consist.cpp @@ -161,7 +161,7 @@ static unsigned numeric_consist2(void) // The fix in the workaround left a whole which was fixed hours later... ex another_zero = pow(zero,numeric(1)/numeric(2)); - if (another_zero.compare(exZERO())) { + if (!another_zero.is_zero()) { clog << "pow(0,1/2) erroneously returned" << another_zero << endl; ++result; } diff --git a/check/paranoia_check.cpp b/check/paranoia_check.cpp index 2e0c40df..f3370602 100644 --- a/check/paranoia_check.cpp +++ b/check/paranoia_check.cpp @@ -74,19 +74,19 @@ static unsigned paranoia_check2(void) g = f - e*y; // After .expand(), g should be zero: - if (!g.expand().is_equal(exZERO())) { + if (!g.expand().is_zero()) { clog << "e = (x + z*x); f = e*y; expand(f - e*y) erroneously returned " << g.expand() << endl; ++result; } // After .eval(), g should be zero: - if (!g.eval().is_equal(exZERO())) { + if (!g.eval().is_zero()) { clog << "e = (x + z*x); f = e*y; eval(f - e*y) erroneously returned " << g.eval() << endl; ++result; } - // This actually worked already back in April. But we are very paranoic! - if (!g.expand().eval().is_equal(exZERO())) { + // This actually worked already back in April 1999. But we are very paranoic! + if (!g.expand().eval().is_zero()) { clog << "e = (x + z*x); f = e*y; eval(expand(f - e*y)) erroneously returned " << g.expand().eval() << endl; ++result; @@ -136,12 +136,12 @@ static unsigned paranoia_check4(void) f = pow(x, 2) + x + 1; g = e - f; - if (!g.is_equal(exZERO())) { + if (!g.is_zero()) { clog << "e = pow(x,2) + x + 1; f = pow(x,2) + x + 1; g = e-f; g erroneously returned " << g << endl; ++result; } - if (!g.is_equal(exZERO())) { + if (!g.is_zero()) { clog << "e = pow(x,2) + x + 1; f = pow(x,2) + x + 1; g = e-f; g.eval() erroneously returned " << g.eval() << endl; ++result; @@ -160,7 +160,7 @@ static unsigned paranoia_check5(void) e = pow(x*y + 1, 2); f = pow(x, 2) * pow(y, 2) + 2*x*y + 1; - if (!(e-f).expand().is_equal(exZERO())) { + if (!(e-f).expand().is_zero()) { clog << "e = pow(x*y+1,2); f = pow(x,2)*pow(y,2) + 2*x*y + 1; (e-f).expand() erroneously returned " << (e-f).expand() << endl; ++result; diff --git a/check/poly_gcd.cpp b/check/poly_gcd.cpp index 1448ae3c..12de267c 100644 --- a/check/poly_gcd.cpp +++ b/check/poly_gcd.cpp @@ -46,7 +46,7 @@ static unsigned poly_gcd1(void) ex f = (e1 + 1) * (e1 + 2); ex g = e2 * (-pow(x, 2) * y[0] * 3 + pow(y[0], 2) - 1); ex r = gcd(f, g); - if (r != exONE()) { + if (r != 1) { clog << "case 1, gcd(" << f << "," << g << ") = " << r << " (should be 1)" << endl; return 1; } @@ -69,10 +69,7 @@ static unsigned poly_gcd2(void) ex f = d * pow(e2 - 2, 2); ex g = d * pow(e1 + 2, 2); ex r = gcd(f, g); - ex re=r.expand(); - ex df1=r-d; - ex df=(r-d).expand(); - if ((r - d).expand().compare(exZERO()) != 0) { + if (!(r - d).expand().is_zero()) { clog << "case 2, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl; return 1; } @@ -92,7 +89,7 @@ static unsigned poly_gcd3(void) ex f = d * (e1 - 2); ex g = d * (e1 + 2); ex r = gcd(f, g); - if ((r - d).expand().compare(exZERO()) != 0) { + if (!(r - d).expand().is_zero()) { clog << "case 3, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl; return 1; } @@ -115,7 +112,7 @@ static unsigned poly_gcd3p(void) ex f = d * (e1 - 2); ex g = d * (e2 + 2); ex r = gcd(f, g); - if ((r - d).expand().compare(exZERO()) != 0) { + if (!(r - d).expand().is_zero()) { clog << "case 3p, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl; return 1; } @@ -140,7 +137,7 @@ static unsigned poly_gcd4(void) ex f = d * (e2 - 1); ex g = d * pow(e3 + 2, 2); ex r = gcd(f, g); - if ((r - d).expand().compare(exZERO()) != 0) { + if (!(r - d).expand().is_zero()) { clog << "case 4, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl; return 1; } @@ -165,7 +162,7 @@ static unsigned poly_gcd5(void) ex f = d * (e2 + 3); ex g = d * (e3 - 3); ex r = gcd(f, g); - if ((r - d).expand().compare(exZERO()) != 0) { + if (!(r - d).expand().is_zero()) { clog << "case 5, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl; return 1; } @@ -185,7 +182,7 @@ static unsigned poly_gcd5p(void) ex f = d * (e1 + 3); ex g = d * (e1 - 3); ex r = gcd(f, g); - if ((r - d).expand().compare(exZERO()) != 0) { + if (!(r - d).expand().is_zero()) { clog << "case 5p, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl; return 1; } @@ -203,7 +200,7 @@ static unsigned poly_gcd6(void) ex f = d * (pow(x, j) + pow(y, j + 1) * pow(z, j) + 1); ex g = d * (pow(x, j + 1) + pow(y, j) * pow(z, j + 1) - 7); ex r = gcd(f, g); - if ((r - d).expand().compare(exZERO()) != 0) { + if (!(r - d).expand().is_zero()) { clog << "case 6, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl; return 1; } @@ -224,7 +221,7 @@ static unsigned poly_gcd7(void) ex f = pow(p, j) * pow(q, k); ex g = pow(p, k) * pow(q, j); ex r = gcd(f, g); - if ((r - d).expand().compare(exZERO()) != 0 && (r + d).expand().compare(exZERO()) != 0) { + if (!(r - d).expand().is_zero() && !(r + d).expand().is_zero()) { clog << "case 7, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl; return 1; } diff --git a/check/series_expansion.cpp b/check/series_expansion.cpp index b67705d7..b8e80204 100644 --- a/check/series_expansion.cpp +++ b/check/series_expansion.cpp @@ -32,7 +32,7 @@ static unsigned check_series(const ex &e, const ex &point, const ex &d, int orde { ex es = e.series(x, point, order); ex ep = static_cast(es.bp)->convert_to_poly(); - if ((ep - d).compare(exZERO()) != 0) { + if (!(ep - d).is_zero()) { clog << "series expansion of " << e << " at " << point << " erroneously returned " << ep << " (instead of " << d << ")" << endl; @@ -50,57 +50,57 @@ static unsigned series1(void) e = sin(x); d = x - pow(x, 3) / 6 + pow(x, 5) / 120 - pow(x, 7) / 5040 + Order(pow(x, 8)); - result += check_series(e, exZERO(), d); + result += check_series(e, 0, d); e = cos(x); d = 1 - pow(x, 2) / 2 + pow(x, 4) / 24 - pow(x, 6) / 720 + Order(pow(x, 8)); - result += check_series(e, exZERO(), d); + result += check_series(e, 0, d); e = exp(x); d = 1 + x + pow(x, 2) / 2 + pow(x, 3) / 6 + pow(x, 4) / 24 + pow(x, 5) / 120 + pow(x, 6) / 720 + pow(x, 7) / 5040 + Order(pow(x, 8)); - result += check_series(e, exZERO(), d); + result += check_series(e, 0, d); e = pow(1 - x, -1); d = 1 + x + pow(x, 2) + pow(x, 3) + pow(x, 4) + pow(x, 5) + pow(x, 6) + pow(x, 7) + Order(pow(x, 8)); - result += check_series(e, exZERO(), d); + result += check_series(e, 0, d); e = x + pow(x, -1); d = x + pow(x, -1); - result += check_series(e, exZERO(), d); + result += check_series(e, 0, d); e = x + pow(x, -1); d = 2 + pow(x-1, 2) - pow(x-1, 3) + pow(x-1, 4) - pow(x-1, 5) + pow(x-1, 6) - pow(x-1, 7) + Order(pow(x-1, 8)); - result += check_series(e, exONE(), d); + result += check_series(e, 1, d); e = pow(x + pow(x, 3), -1); d = pow(x, -1) - x + pow(x, 3) - pow(x, 5) + Order(pow(x, 7)); - result += check_series(e, exZERO(), d); + result += check_series(e, 0, d); e = pow(pow(x, 2) + pow(x, 4), -1); d = pow(x, -2) - 1 + pow(x, 2) - pow(x, 4) + Order(pow(x, 6)); - result += check_series(e, exZERO(), d); + result += check_series(e, 0, d); e = pow(sin(x), -2); d = pow(x, -2) + numeric(1,3) + pow(x, 2) / 15 + pow(x, 4) * 2/189 + Order(pow(x, 5)); - result += check_series(e, exZERO(), d); + result += check_series(e, 0, d); e = sin(x) / cos(x); d = x + pow(x, 3) / 3 + pow(x, 5) * 2/15 + pow(x, 7) * 17/315 + Order(pow(x, 8)); - result += check_series(e, exZERO(), d); + result += check_series(e, 0, d); e = cos(x) / sin(x); d = pow(x, -1) - x / 3 - pow(x, 3) / 45 - pow(x, 5) * 2/945 + Order(pow(x, 6)); - result += check_series(e, exZERO(), d); + result += check_series(e, 0, d); e = pow(numeric(2), x); ex t = log(ex(2)) * x; d = 1 + t + pow(t, 2) / 2 + pow(t, 3) / 6 + pow(t, 4) / 24 + pow(t, 5) / 120 + pow(t, 6) / 720 + pow(t, 7) / 5040 + Order(pow(x, 8)); - result += check_series(e, exZERO(), d.expand()); + result += check_series(e, 0, d.expand()); e = pow(Pi, x); t = log(Pi) * x; d = 1 + t + pow(t, 2) / 2 + pow(t, 3) / 6 + pow(t, 4) / 24 + pow(t, 5) / 120 + pow(t, 6) / 720 + pow(t, 7) / 5040 + Order(pow(x, 8)); - result += check_series(e, exZERO(), d.expand()); + result += check_series(e, 0, d.expand()); return result; } @@ -111,9 +111,9 @@ static unsigned series2(void) unsigned result = 0; ex e, d; - e = pow(sin(x), -1).series(x, exZERO(), 8) + pow(sin(-x), -1).series(x, exZERO(), 12); + e = pow(sin(x), -1).series(x, 0, 8) + pow(sin(-x), -1).series(x, 0, 12); d = Order(pow(x, 6)); - result += check_series(e, exZERO(), d); + result += check_series(e, 0, d); return result; } @@ -124,9 +124,9 @@ static unsigned series3(void) unsigned result = 0; ex e, d; - e = sin(x).series(x, exZERO(), 8) * pow(sin(x), -1).series(x, exZERO(), 12); + e = sin(x).series(x, 0, 8) * pow(sin(x), -1).series(x, 0, 12); d = 1 + Order(pow(x, 7)); - result += check_series(e, exZERO(), d); + result += check_series(e, 0, d); return result; } diff --git a/ginac/add.cpp b/ginac/add.cpp index 804f07ae..ad0ac1fc 100644 --- a/ginac/add.cpp +++ b/ginac/add.cpp @@ -26,6 +26,7 @@ #include "add.h" #include "mul.h" #include "debugmsg.h" +#include "utils.h" #ifndef NO_GINAC_NAMESPACE namespace GiNaC { @@ -87,7 +88,7 @@ add::add(ex const & lh, ex const & rh) { debugmsg("add constructor from ex,ex",LOGLEVEL_CONSTRUCT); tinfo_key = TINFO_add; - overall_coeff=exZERO(); + overall_coeff=_ex0(); construct_from_2_ex(lh,rh); GINAC_ASSERT(is_canonical()); } @@ -96,7 +97,7 @@ add::add(exvector const & v) { debugmsg("add constructor from exvector",LOGLEVEL_CONSTRUCT); tinfo_key = TINFO_add; - overall_coeff=exZERO(); + overall_coeff=_ex0(); construct_from_exvector(v); GINAC_ASSERT(is_canonical()); } @@ -122,7 +123,7 @@ add::add(epvector const & v) { debugmsg("add constructor from epvector",LOGLEVEL_CONSTRUCT); tinfo_key = TINFO_add; - overall_coeff=exZERO(); + overall_coeff=_ex0(); construct_from_epvector(v); GINAC_ASSERT(is_canonical()); } @@ -173,10 +174,10 @@ void add::print(ostream & os, unsigned upper_precedence) const if (coeff.csgn()==-1) os << '-'; first=false; } - if (!coeff.is_equal(numONE()) && - !coeff.is_equal(numMINUSONE())) { + if (!coeff.is_equal(_num1()) && + !coeff.is_equal(_num_1())) { if (coeff.csgn()==-1) - (numMINUSONE()*coeff).print(os, precedence); + (_num_1()*coeff).print(os, precedence); else coeff.print(os, precedence); os << '*'; @@ -219,16 +220,16 @@ void add::printcsrc(ostream & os, unsigned type, unsigned upper_precedence) cons while (it != itend) { // If the coefficient is -1, it is replaced by a single minus sign - if (it->coeff.compare(numONE()) == 0) { + if (it->coeff.compare(_num1()) == 0) { it->rest.bp->printcsrc(os, type, precedence); - } else if (it->coeff.compare(numMINUSONE()) == 0) { + } else if (it->coeff.compare(_num_1()) == 0) { os << "-"; it->rest.bp->printcsrc(os, type, precedence); - } else if (ex_to_numeric(it->coeff).numer().compare(numONE()) == 0) { + } else if (ex_to_numeric(it->coeff).numer().compare(_num1()) == 0) { it->rest.bp->printcsrc(os, type, precedence); os << "/"; ex_to_numeric(it->coeff).denom().printcsrc(os, type, precedence); - } else if (ex_to_numeric(it->coeff).numer().compare(numMINUSONE()) == 0) { + } else if (ex_to_numeric(it->coeff).numer().compare(_num_1()) == 0) { os << "-"; it->rest.bp->printcsrc(os, type, precedence); os << "/"; @@ -241,11 +242,11 @@ void add::printcsrc(ostream & os, unsigned type, unsigned upper_precedence) cons // Separator is "+", except if the following expression would have a leading minus sign it++; - if (it != itend && !(it->coeff.compare(numZERO()) < 0 || (it->coeff.compare(numONE()) == 0 && is_ex_exactly_of_type(it->rest, numeric) && it->rest.compare(numZERO()) < 0))) + if (it != itend && !(it->coeff.compare(_num0()) < 0 || (it->coeff.compare(_num1()) == 0 && is_ex_exactly_of_type(it->rest, numeric) && it->rest.compare(_num0()) < 0))) os << "+"; } - if (!overall_coeff.is_equal(exZERO())) { + if (!overall_coeff.is_equal(_ex0())) { if (overall_coeff.info(info_flags::positive)) os << '+'; overall_coeff.bp->printcsrc(os,type,precedence); } @@ -276,7 +277,7 @@ bool add::info(unsigned inf) const int add::degree(symbol const & s) const { int deg=INT_MIN; - if (!overall_coeff.is_equal(exZERO())) { + if (!overall_coeff.is_equal(_ex0())) { deg=0; } int cur_deg; @@ -290,7 +291,7 @@ int add::degree(symbol const & s) const int add::ldegree(symbol const & s) const { int deg=INT_MAX; - if (!overall_coeff.is_equal(exZERO())) { + if (!overall_coeff.is_equal(_ex0())) { deg=0; } int cur_deg; @@ -344,7 +345,7 @@ ex add::eval(int level) const if (flags & status_flags::evaluated) { GINAC_ASSERT(seq.size()>0); - GINAC_ASSERT((seq.size()>1)||!overall_coeff.is_equal(exZERO())); + GINAC_ASSERT((seq.size()>1)||!overall_coeff.is_equal(_ex0())); return *this; } @@ -352,7 +353,7 @@ ex add::eval(int level) const if (seq_size==0) { // +(;c) -> c return overall_coeff; - } else if ((seq_size==1)&&overall_coeff.is_equal(exZERO())) { + } else if ((seq_size==1)&&overall_coeff.is_equal(_ex0())) { // +(x;0) -> x return recombine_pair_to_ex(*(seq.begin())); } @@ -423,12 +424,12 @@ expair add::split_ex_to_pair(ex const & e) const ex numfactor=mulref.overall_coeff; // mul * mulcopyp=static_cast(mulref.duplicate()); mul * mulcopyp=new mul(mulref); - mulcopyp->overall_coeff=exONE(); + mulcopyp->overall_coeff=_ex1(); mulcopyp->clearflag(status_flags::evaluated); mulcopyp->clearflag(status_flags::hash_calculated); return expair(mulcopyp->setflag(status_flags::dynallocated),numfactor); } - return expair(e,exONE()); + return expair(e,_ex1()); } expair add::combine_ex_with_coeff_to_pair(ex const & e, @@ -440,21 +441,21 @@ expair add::combine_ex_with_coeff_to_pair(ex const & e, ex numfactor=mulref.overall_coeff; //mul * mulcopyp=static_cast(mulref.duplicate()); mul * mulcopyp=new mul(mulref); - mulcopyp->overall_coeff=exONE(); + mulcopyp->overall_coeff=_ex1(); mulcopyp->clearflag(status_flags::evaluated); mulcopyp->clearflag(status_flags::hash_calculated); - if (are_ex_trivially_equal(c,exONE())) { + if (are_ex_trivially_equal(c,_ex1())) { return expair(mulcopyp->setflag(status_flags::dynallocated),numfactor); - } else if (are_ex_trivially_equal(numfactor,exONE())) { + } else if (are_ex_trivially_equal(numfactor,_ex1())) { return expair(mulcopyp->setflag(status_flags::dynallocated),c); } return expair(mulcopyp->setflag(status_flags::dynallocated), ex_to_numeric(numfactor).mul_dyn(ex_to_numeric(c))); } else if (is_ex_exactly_of_type(e,numeric)) { - if (are_ex_trivially_equal(c,exONE())) { - return expair(e,exONE()); + if (are_ex_trivially_equal(c,_ex1())) { + return expair(e,_ex1()); } - return expair(ex_to_numeric(e).mul_dyn(ex_to_numeric(c)),exONE()); + return expair(ex_to_numeric(e).mul_dyn(ex_to_numeric(c)),_ex1()); } return expair(e,c); } @@ -466,8 +467,8 @@ expair add::combine_pair_with_coeff_to_pair(expair const & p, GINAC_ASSERT(is_ex_exactly_of_type(c,numeric)); if (is_ex_exactly_of_type(p.rest,numeric)) { - GINAC_ASSERT(ex_to_numeric(p.coeff).is_equal(numONE())); // should be normalized - return expair(ex_to_numeric(p.rest).mul_dyn(ex_to_numeric(c)),exONE()); + GINAC_ASSERT(ex_to_numeric(p.coeff).is_equal(_num1())); // should be normalized + return expair(ex_to_numeric(p.rest).mul_dyn(ex_to_numeric(c)),_ex1()); } return expair(p.rest,ex_to_numeric(p.coeff).mul_dyn(ex_to_numeric(c))); @@ -475,9 +476,9 @@ expair add::combine_pair_with_coeff_to_pair(expair const & p, ex add::recombine_pair_to_ex(expair const & p) const { - //if (p.coeff.compare(exONE())==0) { - //if (are_ex_trivially_equal(p.coeff,exONE())) { - if (ex_to_numeric(p.coeff).is_equal(numONE())) { + //if (p.coeff.compare(_ex1())==0) { + //if (are_ex_trivially_equal(p.coeff,_ex1())) { + if (ex_to_numeric(p.coeff).is_equal(_num1())) { return p.rest; } else { return p.rest*p.coeff; diff --git a/ginac/basic.cpp b/ginac/basic.cpp index acacd234..7b6a6fbe 100644 --- a/ginac/basic.cpp +++ b/ginac/basic.cpp @@ -226,7 +226,7 @@ int basic::ldegree(symbol const & s) const ex basic::coeff(symbol const & s, int const n) const { - return n==0 ? *this : exZERO(); + return n==0 ? *this : _ex0(); } ex basic::collect(symbol const & s) const diff --git a/ginac/color.cpp b/ginac/color.cpp index 696d70f7..0f5f6047 100644 --- a/ginac/color.cpp +++ b/ginac/color.cpp @@ -34,6 +34,7 @@ #include "numeric.h" #include "relational.h" #include "debugmsg.h" +#include "utils.h" #ifndef NO_GINAC_NAMESPACE namespace GiNaC { @@ -242,7 +243,7 @@ ex color::eval(int level) const int sig=canonicalize_indices(iv,antisymmetric); if (sig!=INT_MAX) { // something has changed while sorting indices, more evaluations later - if (sig==0) return exZERO(); + if (sig==0) return _ex0(); return ex(sig)*color(type,iv,representation_label); } } @@ -267,9 +268,9 @@ ex color::eval(int level) const // check for delta8_{a,b} where a and b are numeric indices, replace by 0 or 1 if ((!idx1.is_symbolic())&&(!idx2.is_symbolic())) { if ((idx1.get_value()!=idx2.get_value())) { - return exONE(); + return _ex1(); } else { - return exZERO(); + return _ex0(); } } } @@ -283,9 +284,9 @@ ex color::eval(int level) const coloridx const & idx3=ex_to_coloridx(seq[2]); if (idx1.is_equal_same_type(idx2) && idx1.is_symbolic()) { - return exZERO(); + return _ex0(); } else if (idx2.is_equal_same_type(idx3) && idx2.is_symbolic()) { - return exZERO(); + return _ex0(); } // check for three numeric indices @@ -294,9 +295,9 @@ ex color::eval(int level) const GINAC_ASSERT(idx2.get_value()<=idx3.get_value()); if (CMPINDICES(1,4,6)||CMPINDICES(1,5,7)||CMPINDICES(2,5,6)|| CMPINDICES(3,4,4)||CMPINDICES(3,5,5)) { - return exHALF(); + return _ex1_2(); } else if (CMPINDICES(2,4,7)||CMPINDICES(3,6,6)||CMPINDICES(3,7,7)) { - return -exHALF(); + return -_ex1_2(); } else if (CMPINDICES(1,1,8)||CMPINDICES(2,2,8)||CMPINDICES(3,3,8)) { return 1/sqrt(numeric(3)); } else if (CMPINDICES(8,8,8)) { @@ -304,7 +305,7 @@ ex color::eval(int level) const } else if (CMPINDICES(4,4,8)||CMPINDICES(5,5,8)||CMPINDICES(6,6,8)||CMPINDICES(7,7,8)) { return -1/(2*sqrt(numeric(3))); } - return exZERO(); + return _ex0(); } } break; @@ -320,12 +321,12 @@ ex color::eval(int level) const GINAC_ASSERT(idx1.get_value()<=idx2.get_value()); GINAC_ASSERT(idx2.get_value()<=idx3.get_value()); if (CMPINDICES(1,2,3)) { - return exONE(); + return _ex1(); } else if (CMPINDICES(1,4,7)||CMPINDICES(2,4,6)|| CMPINDICES(2,5,7)||CMPINDICES(3,4,5)) { - return exHALF(); + return _ex1_2(); } else if (CMPINDICES(1,5,6)||CMPINDICES(3,6,7)) { - return -exHALF(); + return -_ex1_2(); } else if (CMPINDICES(4,5,8)||CMPINDICES(6,7,8)) { return sqrt(numeric(3))/2; } else if (CMPINDICES(8,8,8)) { @@ -333,7 +334,7 @@ ex color::eval(int level) const } else if (CMPINDICES(4,4,8)||CMPINDICES(5,5,8)||CMPINDICES(6,6,8)||CMPINDICES(7,7,8)) { return -1/(2*sqrt(numeric(3))); } - return exZERO(); + return _ex0(); } break; } @@ -405,7 +406,7 @@ ex color::simplify_ncmul(exvector const & v) const } else { // a contracted index should occur exactly twice GINAC_ASSERT(replacements==2); - *it=exONE(); + *it=_ex1(); something_changed=true; } } @@ -420,7 +421,7 @@ ex color::simplify_ncmul(exvector const & v) const } else { // a contracted index should occur exactly twice GINAC_ASSERT(replacements==2); - *it=exONE(); + *it=_ex1(); something_changed=true; } } @@ -457,7 +458,7 @@ ex color::simplify_ncmul(exvector const & v) const color const & col1=ex_to_color(*it1); color const & col2=ex_to_color(*it2); exvector iv_intersect=idx_intersect(col1.seq,col2.seq); - if (iv_intersect.size()>=2) return exZERO(); + if (iv_intersect.size()>=2) return _ex0(); } } } @@ -474,13 +475,13 @@ ex color::simplify_ncmul(exvector const & v) const if (iv_intersect.size()>=2) { if (iv_intersect.size()==3) { *it1=numeric(40)/numeric(3); - *it2=exONE(); + *it2=_ex1(); } else { int sig1, sig2; // unimportant, since symmetric ex idx1=permute_free_index_to_front(col1.seq,iv_intersect,false,&sig1); ex idx2=permute_free_index_to_front(col2.seq,iv_intersect,false,&sig2); *it1=numeric(5)/numeric(3)*color(color_delta8,idx1,idx2); - *it2=exONE(); + *it2=_ex1(); } return nonsimplified_ncmul(recombine_color_string( delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec)); @@ -501,13 +502,13 @@ ex color::simplify_ncmul(exvector const & v) const if (iv_intersect.size()>=2) { if (iv_intersect.size()==3) { *it1=numeric(24); - *it2=exONE(); + *it2=_ex1(); } else { int sig1, sig2; ex idx1=permute_free_index_to_front(col1.seq,iv_intersect,true,&sig1); ex idx2=permute_free_index_to_front(col2.seq,iv_intersect,true,&sig2); *it1=numeric(sig1*sig2*5)/numeric(3)*color(color_delta8,idx1,idx2); - *it2=exONE(); + *it2=_ex1(); } return nonsimplified_ncmul(recombine_color_string( delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec)); @@ -732,7 +733,7 @@ ex color_trace_of_one_representation_label(exvector const & v) return numeric(COLOR_THREE); } else if (v.size()==1) { GINAC_ASSERT(is_ex_exactly_of_type(*(v.begin()),color)); - return exZERO(); + return _ex0(); } exvector v1=v; ex last_element=v1.back(); @@ -846,12 +847,12 @@ ex simplify_pure_color_string(ex const & e) for (unsigned k=i+1; kflags & status_flags::dynallocated); + GINAC_ASSERT(ex0().bp!=0); + GINAC_ASSERT(ex0().bp->flags & status_flags::dynallocated); GINAC_ASSERT(bp!=0); ++bp->refcount; } @@ -313,12 +314,12 @@ ex ex::numer(bool normalize) const // something^(-int) if (is_ex_exactly_of_type(n, power) && n.op(1).info(info_flags::negint)) - return exONE(); + return _ex1(); // something^(int) * something^(int) * ... if (!is_ex_exactly_of_type(n, mul)) return n; - ex res = exONE(); + ex res = _ex1(); for (int i=0; i exvector; - -// enum definitions +extern ex const & _ex0(void); /* FIXME: should this pollute headers? */ -ex const & exZERO(void); -ex const & exONE(void); -ex const & exTWO(void); -ex const & exTHREE(void); -ex const & exMINUSONE(void); -ex const & exHALF(void); -ex const & exMINUSHALF(void); +// typedef vector exvector; #define INLINE_EX_CONSTRUCTORS @@ -65,10 +57,10 @@ class ex public: ex() #ifdef INLINE_EX_CONSTRUCTORS - : bp(exZERO().bp) + : bp(_ex0().bp) { - GINAC_ASSERT(exZERO().bp!=0); - GINAC_ASSERT(exZERO().bp->flags & status_flags::dynallocated); + GINAC_ASSERT(_ex0().bp!=0); + GINAC_ASSERT(_ex0().bp->flags & status_flags::dynallocated); GINAC_ASSERT(bp!=0); ++bp->refcount; #ifdef OBSCURE_CINT_HACK @@ -224,7 +216,7 @@ public: #else ; #endif // def INLINE_EX_CONSTRUCTORS - bool is_zero(void) const {return compare(exZERO()) == 0;}; + bool is_zero(void) const {return compare(_ex0())==0;}; unsigned return_type(void) const; unsigned return_type_tinfo(void) const; diff --git a/ginac/expair.h b/ginac/expair.h index b130465b..8580598d 100644 --- a/ginac/expair.h +++ b/ginac/expair.h @@ -59,7 +59,7 @@ public: { GINAC_ASSERT(is_ex_exactly_of_type(coeff,numeric)); return is_ex_exactly_of_type(rest,numeric) && - (ex_to_numeric(coeff).compare(numONE())==0); + (coeff.is_equal(ex(1))); } bool is_equal(expair const & other) const @@ -101,14 +101,14 @@ public: */ if (is_ex_exactly_of_type(rest,numeric) && is_ex_exactly_of_type(other.rest,numeric)) { - if (ex_to_numeric(coeff).compare(numONE())==0) { - if (ex_to_numeric(other.coeff).compare(numONE())==0) { + if (coeff.is_equal(ex(1))) { + if ((other.coeff).is_equal(ex(1))) { // both have coeff 1: compare rests return rest.compare(other.rest)<0; } // only this has coeff 1: > return false; - } else if (ex_to_numeric(other.coeff).compare(numONE())==0) { + } else if ((other.coeff).is_equal(ex(1))) { // only other has coeff 1: < return true; } @@ -122,14 +122,14 @@ public: { if (is_ex_exactly_of_type(rest,numeric) && is_ex_exactly_of_type(other.rest,numeric)) { - if (ex_to_numeric(coeff).compare(numONE())==0) { - if (ex_to_numeric(other.coeff).compare(numONE())==0) { + if ((coeff).is_equal(ex(1))) { + if ((other.coeff).is_equal(ex(1))) { // both have coeff 1: compare rests return rest.compare(other.rest); } // only this has coeff 1: > return 1; - } else if (ex_to_numeric(other.coeff).compare(numONE())==0) { + } else if ((other.coeff).is_equal(ex(1))) { // only other has coeff 1: < return -1; } diff --git a/ginac/expairseq.cpp b/ginac/expairseq.cpp index 209c7252..1fa592e7 100644 --- a/ginac/expairseq.cpp +++ b/ginac/expairseq.cpp @@ -539,7 +539,7 @@ void expairseq::printseq(ostream & os, char delim, unsigned this_precedence, expair expairseq::split_ex_to_pair(ex const & e) const { - return expair(e,exONE()); + return expair(e,_ex1()); } expair expairseq::combine_ex_with_coeff_to_pair(ex const & e, @@ -571,7 +571,7 @@ bool expairseq::expair_needs_further_processing(epp it) ex expairseq::default_overall_coeff(void) const { - return exZERO(); + return _ex0(); } void expairseq::combine_overall_coeff(ex const & c) @@ -1314,7 +1314,7 @@ void expairseq::drop_coeff_0_terms(epvector::iterator & first_numeric, if (!touched[i]) { ++current; ++i; - } else if (!ex_to_numeric((*current).coeff).is_equal(numZERO())) { + } else if (!ex_to_numeric((*current).coeff).is_equal(_num0())) { ++current; ++i; } else { @@ -1356,7 +1356,7 @@ void expairseq::drop_coeff_0_terms(epvector::iterator & first_numeric, bool expairseq::has_coeff_0(void) const { for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) { - if ((*cit).coeff.is_equal(exZERO())) { + if ((*cit).coeff.is_equal(_ex0())) { return true; } } diff --git a/ginac/inifcns.cpp b/ginac/inifcns.cpp index 8b7f327a..0a38c4b0 100644 --- a/ginac/inifcns.cpp +++ b/ginac/inifcns.cpp @@ -35,6 +35,7 @@ #include "relational.h" #include "series.h" #include "symbol.h" +#include "utils.h" #ifndef NO_GINAC_NAMESPACE namespace GiNaC { @@ -48,9 +49,9 @@ static ex Li2_eval(ex const & x) { if (x.is_zero()) return x; - if (x.is_equal(exONE())) + if (x.is_equal(_ex1())) return power(Pi, 2) / 6; - if (x.is_equal(exMINUSONE())) + if (x.is_equal(_ex_1())) return -power(Pi, 2) / 12; return Li2(x).hold(); } @@ -117,7 +118,7 @@ static ex Order_eval(ex const & x) if (is_ex_exactly_of_type(x, numeric)) { // O(c)=O(1) - return Order(exONE()).hold(); + return Order(_ex1()).hold(); } else if (is_ex_exactly_of_type(x, mul)) { @@ -135,7 +136,7 @@ static ex Order_series(ex const & x, symbol const & s, ex const & point, int ord { // Just wrap the function into a series object epvector new_seq; - new_seq.push_back(expair(Order(exONE()), numeric(min(x.ldegree(s), order)))); + new_seq.push_back(expair(Order(_ex1()), numeric(min(x.ldegree(s), order)))); return series(s, point, new_seq); } @@ -237,7 +238,7 @@ ex lsolve(ex const &eqns, ex const &symbols) ex ncpower(ex const &basis, unsigned exponent) { if (exponent==0) { - return exONE(); + return _ex1(); } exvector v; diff --git a/ginac/inifcns.h b/ginac/inifcns.h index d7642ba2..fdd33108 100644 --- a/ginac/inifcns.h +++ b/ginac/inifcns.h @@ -100,7 +100,7 @@ DECLARE_FUNCTION_1P(gamma) DECLARE_FUNCTION_2P(beta) // overloading at work: we cannot use the macros -/** Psi-function (aka polygamma-function). */ +/** Psi-function (aka digamma-function). */ extern const unsigned function_index_psi1; inline function psi(ex const & p1) { return function(function_index_psi1, p1); diff --git a/ginac/inifcns_gamma.cpp b/ginac/inifcns_gamma.cpp index ecec37a6..3f7fc677 100644 --- a/ginac/inifcns_gamma.cpp +++ b/ginac/inifcns_gamma.cpp @@ -63,7 +63,7 @@ static ex gamma_eval(ex const & x) if (x.info(info_flags::integer)) { // gamma(n+1) -> n! for postitive n if (x.info(info_flags::posint)) { - return factorial(ex_to_numeric(x).sub(numONE())); + return factorial(ex_to_numeric(x).sub(_num1())); } else { throw (std::domain_error("gamma_eval(): simple pole")); } @@ -73,16 +73,16 @@ static ex gamma_eval(ex const & x) // trap positive x==(n+1/2) // gamma(n+1/2) -> Pi^(1/2)*(1*3*..*(2*n-1))/(2^n) if ((x*2).info(info_flags::posint)) { - numeric n = ex_to_numeric(x).sub(numHALF()); - numeric coefficient = doublefactorial(n.mul(numTWO()).sub(numONE())); - coefficient = coefficient.div(numTWO().power(n)); - return coefficient * pow(Pi,numHALF()); + numeric n = ex_to_numeric(x).sub(_num1_2()); + numeric coefficient = doublefactorial(n.mul(_num2()).sub(_num1())); + coefficient = coefficient.div(_num2().power(n)); + return coefficient * pow(Pi,_num1_2()); } else { // trap negative x==(-n+1/2) // gamma(-n+1/2) -> Pi^(1/2)*(-2)^n/(1*3*..*(2*n-1)) - numeric n = abs(ex_to_numeric(x).sub(numHALF())); + numeric n = abs(ex_to_numeric(x).sub(_num1_2())); numeric coefficient = numeric(-2).power(n); - coefficient = coefficient.div(doublefactorial(n.mul(numTWO()).sub(numONE())));; + coefficient = coefficient.div(doublefactorial(n.mul(_num2()).sub(_num1())));; return coefficient*sqrt(Pi); } } @@ -94,24 +94,28 @@ static ex gamma_diff(ex const & x, unsigned diff_param) { GINAC_ASSERT(diff_param==0); - return psi(x)*gamma(x); // diff(log(gamma(x)),x)==psi(x) + // d/dx log(gamma(x)) -> psi(x) + // d/dx gamma(x) -> psi(x)*gamma(x) + return psi(x)*gamma(x); } static ex gamma_series(ex const & x, symbol const & s, ex const & point, int order) { // method: - // Taylor series where there is no pole falls back to psi functions. - // On a pole at -n use the identity - // series(GAMMA(x),x=-n,order) == - // series(GAMMA(x+n+1)/(x*(x+1)...*(x+n)),x=-n,order+1); + // Taylor series where there is no pole falls back to psi function evaluation. + // On a pole at -m use the recurrence relation + // gamma(x) == gamma(x+1) / x + // from which follows + // series(gamma(x),x,-m,order) == + // series(gamma(x+m+1)/(x*(x+1)...*(x+m)),x,-m,order+1); ex xpoint = x.subs(s==point); if (!xpoint.info(info_flags::integer) || xpoint.info(info_flags::positive)) - throw do_taylor(); - // if we got here we have to care for a simple pole at -n: - numeric n = -ex_to_numeric(xpoint); - ex ser_numer = gamma(x+n+exONE()); - ex ser_denom = exONE(); - for (numeric p; p<=n; ++p) + throw do_taylor(); // caught by function::series() + // if we got here we have to care for a simple pole at -m: + numeric m = -ex_to_numeric(xpoint); + ex ser_numer = gamma(x+m+_ex1()); + ex ser_denom = _ex1(); + for (numeric p; p<=m; ++p) ser_denom *= x+p; return (ser_numer/ser_denom).series(s, point, order+1); } @@ -144,13 +148,13 @@ static ex beta_eval(ex const & x, ex const & y) ny.is_real() && ny.is_integer()) { if (nx.is_negative()) { if (nx<=-ny) - return numMINUSONE().power(ny)*beta(1-x-y, y); + return _num_1().power(ny)*beta(1-x-y, y); else throw (std::domain_error("beta_eval(): simple pole")); } if (ny.is_negative()) { if (ny<=-nx) - return numMINUSONE().power(nx)*beta(1-y-x, x); + return _num_1().power(nx)*beta(1-y-x, x); else throw (std::domain_error("beta_eval(): simple pole")); } @@ -160,7 +164,7 @@ static ex beta_eval(ex const & x, ex const & y) if ((nx+ny).is_real() && (nx+ny).is_integer() && !(nx+ny).is_positive()) - return exZERO(); + return _ex0(); return gamma(x)*gamma(y)/gamma(x+y); } return beta(x,y).hold(); @@ -171,9 +175,11 @@ static ex beta_diff(ex const & x, ex const & y, unsigned diff_param) GINAC_ASSERT(diff_param<2); ex retval; - if (diff_param==0) // d/dx beta(x,y) + // d/dx beta(x,y) -> (psi(x)-psi(x+y)) * beta(x,y) + if (diff_param==0) retval = (psi(x)-psi(x+y))*beta(x,y); - if (diff_param==1) // d/dy beta(x,y) + // d/dy beta(x,y) -> (psi(y)-psi(x+y)) * beta(x,y) + if (diff_param==1) retval = (psi(y)-psi(x+y))*beta(x,y); return retval; } @@ -181,7 +187,7 @@ static ex beta_diff(ex const & x, ex const & y, unsigned diff_param) REGISTER_FUNCTION(beta, beta_eval, beta_evalf, beta_diff, NULL); ////////// -// Psi-function (aka polygamma-function) +// Psi-function (aka digamma-function) ////////// static ex psi1_evalf(ex const & x) @@ -193,7 +199,7 @@ static ex psi1_evalf(ex const & x) return psi(ex_to_numeric(x)); } -/** Evaluation of polygamma-function psi(x). +/** Evaluation of digamma-function psi(x). * Somebody ought to provide some good numerical evaluation some day... */ static ex psi1_eval(ex const & x) { @@ -204,18 +210,18 @@ static ex psi1_eval(ex const & x) // psi(n) -> 1 + 1/2 +...+ 1/(n-1) - EulerGamma if (x.info(info_flags::integer)) { numeric rat(0); - for (numeric i(ex_to_numeric(x)-numONE()); i.is_positive(); --i) + for (numeric i(ex_to_numeric(x)-_num1()); i.is_positive(); --i) rat += i.inverse(); return rat-EulerGamma; } // psi((2m+1)/2) -> 2/(2m+1) + 2/2m +...+ 2/1 - EulerGamma - 2log(2) - if ((exTWO()*x).info(info_flags::integer)) { + if ((_ex2()*x).info(info_flags::integer)) { numeric rat(0); - for (numeric i((ex_to_numeric(x)-numONE())*numTWO()); i.is_positive(); i-=numTWO()) - rat += numTWO()*i.inverse(); - return rat-EulerGamma-exTWO()*log(exTWO()); + for (numeric i((ex_to_numeric(x)-_num1())*_num2()); i.is_positive(); i-=_num2()) + rat += _num2()*i.inverse(); + return rat-EulerGamma-_ex2()*log(_ex2()); } - if (x.compare(exONE())==1) { + if (x.compare(_ex1())==1) { // should call numeric, since >1 } } @@ -227,10 +233,32 @@ static ex psi1_diff(ex const & x, unsigned diff_param) { GINAC_ASSERT(diff_param==0); - return psi(exONE(), x); + // d/dx psi(x) -> psi(1,x) + return psi(_ex1(), x); } -const unsigned function_index_psi1 = function::register_new("psi", psi1_eval, psi1_evalf, psi1_diff, NULL); +static ex psi1_series(ex const & x, symbol const & s, ex const & point, int order) +{ + // method: + // Taylor series where there is no pole falls back to polygamma function + // evaluation. + // On a pole at -m use the recurrence relation + // psi(x) == psi(x+1) - 1/z + // from which follows + // series(psi(x),x,-m,order) == + // series(psi(x+m+1) - 1/x - 1/(x+1) - 1/(x+m)),x,-m,order); + ex xpoint = x.subs(s==point); + if (!xpoint.info(info_flags::integer) || xpoint.info(info_flags::positive)) + throw do_taylor(); // caught by function::series() + // if we got here we have to care for a simple pole at -m: + numeric m = -ex_to_numeric(xpoint); + ex recur; + for (numeric p; p<=m; ++p) + recur += power(x+p,_ex_1()); + return (psi(x+m+_ex1())-recur).series(s, point, order); +} + +const unsigned function_index_psi1 = function::register_new("psi", psi1_eval, psi1_evalf, psi1_diff, psi1_series); ////////// // Psi-functions (aka polygamma-functions) psi(0,x)==psi(x) @@ -254,14 +282,14 @@ static ex psi2_eval(ex const & n, ex const & x) if (n.is_zero()) return psi(x); // psi(-1,x) -> log(gamma(x)) - if (n.is_equal(exMINUSONE())) + if (n.is_equal(_ex_1())) return log(gamma(x)); if (n.info(info_flags::numeric) && n.info(info_flags::posint) && x.info(info_flags::numeric)) { numeric nn = ex_to_numeric(n); numeric nx = ex_to_numeric(x); - if (x.is_equal(exONE())) - return numMINUSONE().power(nn+numONE())*factorial(nn)*zeta(ex(nn+numONE())); + if (x.is_equal(_ex1())) + return _num_1().power(nn+_num1())*factorial(nn)*zeta(ex(nn+_num1())); } return psi(n, x).hold(); } @@ -274,11 +302,34 @@ static ex psi2_diff(ex const & n, ex const & x, unsigned diff_param) // d/dn psi(n,x) throw(std::logic_error("cannot diff psi(n,x) with respect to n")); } - // d/dx psi(n,x) + // d/dx psi(n,x) -> psi(n+1,x) return psi(n+1, x); } -const unsigned function_index_psi2 = function::register_new("psi", psi2_eval, psi2_evalf, psi2_diff, NULL); +static ex psi2_series(ex const & n, ex const & x, symbol const & s, ex const & point, int order) +{ + // method: + // Taylor series where there is no pole falls back to polygamma function + // evaluation. + // On a pole at -m use the recurrence relation + // psi(n,x) == psi(n,x+1) - (-)^n * n! / z^(n+1) + // from which follows + // series(psi(x),x,-m,order) == + // series(psi(x+m+1) - (-1)^n * n! + // * ((x)^(-n-1) + (x+1)^(-n-1) + (x+m)^(-n-1))),x,-m,order); + ex xpoint = x.subs(s==point); + if (!xpoint.info(info_flags::integer) || xpoint.info(info_flags::positive)) + throw do_taylor(); // caught by function::series() + // if we got here we have to care for a pole of order n+1 at -m: + numeric m = -ex_to_numeric(xpoint); + ex recur; + for (numeric p; p<=m; ++p) + recur += power(x+p,-n+_ex_1()); + recur *= factorial(n)*power(_ex_1(),n); + return (psi(n, x+m+_ex1())-recur).series(s, point, order); +} + +const unsigned function_index_psi2 = function::register_new("psi", psi2_eval, psi2_evalf, psi2_diff, psi2_series); #ifndef NO_GINAC_NAMESPACE } // namespace GiNaC diff --git a/ginac/inifcns_trans.cpp b/ginac/inifcns_trans.cpp index 8838360a..4473fe02 100644 --- a/ginac/inifcns_trans.cpp +++ b/ginac/inifcns_trans.cpp @@ -54,19 +54,19 @@ static ex exp_eval(ex const & x) { // exp(0) -> 1 if (x.is_zero()) { - return exONE(); + return _ex1(); } // exp(n*Pi*I/2) -> {+1|+I|-1|-I} - ex TwoExOverPiI=(2*x)/(Pi*I); + ex TwoExOverPiI=(_ex2()*x)/(Pi*I); if (TwoExOverPiI.info(info_flags::integer)) { - numeric z=mod(ex_to_numeric(TwoExOverPiI),numeric(4)); - if (z.is_equal(numZERO())) - return exONE(); - if (z.is_equal(numONE())) + numeric z=mod(ex_to_numeric(TwoExOverPiI),_num4()); + if (z.is_equal(_num0())) + return _ex1(); + if (z.is_equal(_num1())) return ex(I); - if (z.is_equal(numTWO())) - return exMINUSONE(); - if (z.is_equal(numTHREE())) + if (z.is_equal(_num2())) + return _ex_1(); + if (z.is_equal(_num3())) return ex(-I); } // exp(log(x)) -> x @@ -84,6 +84,7 @@ static ex exp_diff(ex const & x, unsigned diff_param) { GINAC_ASSERT(diff_param==0); + // d/dx exp(x) -> exp(x) return exp(x); } @@ -105,20 +106,15 @@ static ex log_evalf(ex const & x) static ex log_eval(ex const & x) { if (x.info(info_flags::numeric)) { - // log(1) -> 0 - if (x.is_equal(exONE())) - return exZERO(); - // log(-1) -> I*Pi - if (x.is_equal(exMINUSONE())) - return (I*Pi); - // log(I) -> Pi*I/2 - if (x.is_equal(I)) - return (I*Pi*numeric(1,2)); - // log(-I) -> -Pi*I/2 - if (x.is_equal(-I)) - return (I*Pi*numeric(-1,2)); - // log(0) -> throw singularity - if (x.is_equal(exZERO())) + if (x.is_equal(_ex1())) // log(1) -> 0 + return _ex0(); + if (x.is_equal(_ex_1())) // log(-1) -> I*Pi + return (I*Pi); + if (x.is_equal(I)) // log(I) -> Pi*I/2 + return (Pi*I*_num1_2()); + if (x.is_equal(-I)) // log(-I) -> -Pi*I/2 + return (Pi*I*_num_1_2()); + if (x.is_equal(_ex0())) // log(0) -> infinity throw(std::domain_error("log_eval(): log(0)")); // log(float) if (!x.info(info_flags::crational)) @@ -138,6 +134,7 @@ static ex log_diff(ex const & x, unsigned diff_param) { GINAC_ASSERT(diff_param==0); + // d/dx log(x) -> 1/x return power(x, -1); } @@ -158,18 +155,34 @@ static ex sin_evalf(ex const & x) static ex sin_eval(ex const & x) { - ex xOverPi=x/Pi; - if (xOverPi.info(info_flags::numeric)) { - // sin(n*Pi) -> 0 - if (xOverPi.info(info_flags::integer)) - return exZERO(); - - // sin((2n+1)*Pi/2) -> {+|-}1 - ex xOverPiMinusHalf=xOverPi-exHALF(); - if (xOverPiMinusHalf.info(info_flags::even)) - return exONE(); - else if (xOverPiMinusHalf.info(info_flags::odd)) - return exMINUSONE(); + // sin(n*Pi/6) -> { 0 | +/-1/2 | +/-sqrt(3)/2 | +/-1 } + ex SixExOverPi = _ex6()*x/Pi; + if (SixExOverPi.info(info_flags::integer)) { + numeric z = smod(ex_to_numeric(SixExOverPi),_num12()); + if (z.is_equal(_num_5())) // sin(7*Pi/6) -> -1/2 + return _ex_1_2(); + if (z.is_equal(_num_4())) // sin(8*Pi/6) -> -sqrt(3)/2 + return _ex_1_2()*power(_ex3(),_ex1_2()); + if (z.is_equal(_num_3())) // sin(9*Pi/6) -> -1 + return _ex_1(); + if (z.is_equal(_num_2())) // sin(10*Pi/6) -> -sqrt(3)/2 + return _ex_1_2()*power(_ex3(),_ex1_2()); + if (z.is_equal(_num_1())) // sin(11*Pi/6) -> -1/2 + return _ex_1_2(); + if (z.is_equal(_num0())) // sin(0) -> 0 + return _ex0(); + if (z.is_equal(_num1())) // sin(1*Pi/6) -> 1/2 + return _ex1_2(); + if (z.is_equal(_num2())) // sin(2*Pi/6) -> sqrt(3)/2 + return _ex1_2()*power(_ex3(),_ex1_2()); + if (z.is_equal(_num3())) // sin(3*Pi/6) -> 1 + return _ex1(); + if (z.is_equal(_num4())) // sin(4*Pi/6) -> sqrt(3)/2 + return _ex1_2()*power(_ex3(),_ex1_2()); + if (z.is_equal(_num5())) // sin(5*Pi/6) -> 1/2 + return _ex1_2(); + if (z.is_equal(_num6())) // sin(6*Pi/6) -> 0 + return _ex0(); } if (is_ex_exactly_of_type(x, function)) { @@ -177,12 +190,12 @@ static ex sin_eval(ex const & x) // sin(asin(x)) -> x if (is_ex_the_function(x, asin)) return t; - // sin(acos(x)) -> (1-x^2)^(1/2) + // sin(acos(x)) -> sqrt(1-x^2) if (is_ex_the_function(x, acos)) - return power(exONE()-power(t,exTWO()),exHALF()); + return power(_ex1()-power(t,_ex2()),_ex1_2()); // sin(atan(x)) -> x*(1+x^2)^(-1/2) if (is_ex_the_function(x, atan)) - return t*power(exONE()+power(t,exTWO()),exMINUSHALF()); + return t*power(_ex1()+power(t,_ex2()),_ex_1_2()); } // sin(float) -> float @@ -196,6 +209,7 @@ static ex sin_diff(ex const & x, unsigned diff_param) { GINAC_ASSERT(diff_param==0); + // d/dx sin(x) -> cos(x) return cos(x); } @@ -216,18 +230,34 @@ static ex cos_evalf(ex const & x) static ex cos_eval(ex const & x) { - ex xOverPi=x/Pi; - if (xOverPi.info(info_flags::numeric)) { - // cos(n*Pi) -> {+|-}1 - if (xOverPi.info(info_flags::even)) - return exONE(); - else if (xOverPi.info(info_flags::odd)) - return exMINUSONE(); - - // cos((2n+1)*Pi/2) -> 0 - ex xOverPiMinusHalf=xOverPi-exHALF(); - if (xOverPiMinusHalf.info(info_flags::integer)) - return exZERO(); + // cos(n*Pi/6) -> { 0 | +/-1/2 | +/-sqrt(3)/2 | +/-1 } + ex SixExOverPi = _ex6()*x/Pi; + if (SixExOverPi.info(info_flags::integer)) { + numeric z = smod(ex_to_numeric(SixExOverPi),_num12()); + if (z.is_equal(_num_5())) // cos(7*Pi/6) -> -sqrt(3)/2 + return _ex_1_2()*power(_ex3(),_ex1_2()); + if (z.is_equal(_num_4())) // cos(8*Pi/6) -> -1/2 + return _ex_1_2(); + if (z.is_equal(_num_3())) // cos(9*Pi/6) -> 0 + return _ex0(); + if (z.is_equal(_num_2())) // cos(10*Pi/6) -> 1/2 + return _ex1_2(); + if (z.is_equal(_num_1())) // cos(11*Pi/6) -> sqrt(3)/2 + return _ex1_2()*power(_ex3(),_ex1_2()); + if (z.is_equal(_num0())) // cos(0) -> 1 + return _ex1(); + if (z.is_equal(_num1())) // cos(1*Pi/6) -> sqrt(3)/2 + return _ex1_2()*power(_ex3(),_ex1_2()); + if (z.is_equal(_num2())) // cos(2*Pi/6) -> 1/2 + return _ex1_2(); + if (z.is_equal(_num3())) // cos(3*Pi/6) -> 0 + return _ex0(); + if (z.is_equal(_num4())) // cos(4*Pi/6) -> -1/2 + return _ex_1_2(); + if (z.is_equal(_num5())) // cos(5*Pi/6) -> -sqrt(3)/2 + return _ex_1_2()*power(_ex3(),_ex1_2()); + if (z.is_equal(_num6())) // cos(6*Pi/6) -> -1 + return _ex_1(); } if (is_ex_exactly_of_type(x, function)) { @@ -237,10 +267,10 @@ static ex cos_eval(ex const & x) return t; // cos(asin(x)) -> (1-x^2)^(1/2) if (is_ex_the_function(x, asin)) - return power(exONE()-power(t,exTWO()),exHALF()); + return power(_ex1()-power(t,_ex2()),_ex1_2()); // cos(atan(x)) -> (1+x^2)^(-1/2) if (is_ex_the_function(x, atan)) - return power(exONE()+power(t,exTWO()),exMINUSHALF()); + return power(_ex1()+power(t,_ex2()),_ex_1_2()); } // cos(float) -> float @@ -254,7 +284,8 @@ static ex cos_diff(ex const & x, unsigned diff_param) { GINAC_ASSERT(diff_param==0); - return numMINUSONE()*sin(x); + // d/dx cos(x) -> -sin(x) + return _ex_1()*sin(x); } REGISTER_FUNCTION(cos, cos_eval, cos_evalf, cos_diff, NULL); @@ -274,23 +305,24 @@ static ex tan_evalf(ex const & x) static ex tan_eval(ex const & x) { - // tan(n*Pi/3) -> {0|3^(1/2)|-(3^(1/2))} - ex ThreeExOverPi=numTHREE()*x/Pi; - if (ThreeExOverPi.info(info_flags::integer)) { - numeric z=mod(ex_to_numeric(ThreeExOverPi),numeric(3)); - if (z.is_equal(numZERO())) - return exZERO(); - if (z.is_equal(numONE())) - return power(exTHREE(),exHALF()); - if (z.is_equal(numTWO())) - return -power(exTHREE(),exHALF()); + // tan(n*Pi/6) -> { 0 | +/-sqrt(3) | +/-sqrt(3)/2 } + ex SixExOverPi = _ex6()*x/Pi; + if (SixExOverPi.info(info_flags::integer)) { + numeric z = smod(ex_to_numeric(SixExOverPi),_num6()); + if (z.is_equal(_num_2())) // tan(4*Pi/6) -> -sqrt(3) + return _ex_1()*power(_ex3(),_ex1_2()); + if (z.is_equal(_num_1())) // tan(5*Pi/6) -> -sqrt(3)/3 + return _ex_1_3()*power(_ex3(),_ex1_2()); + if (z.is_equal(_num0())) // tan(0) -> 0 + return _ex0(); + if (z.is_equal(_num1())) // tan(1*Pi/6) -> sqrt(3)/3 + return _ex1_3()*power(_ex3(),_ex1_2()); + if (z.is_equal(_num2())) // tan(2*Pi/6) -> sqrt(3) + return power(_ex3(),_ex1_2()); + if (z.is_equal(_num3())) // tan(3*Pi/6) -> infinity + throw (std::domain_error("tan_eval(): infinity")); } - - // tan((2n+1)*Pi/2) -> throw - ex ExOverPiMinusHalf=x/Pi-exHALF(); - if (ExOverPiMinusHalf.info(info_flags::integer)) - throw (std::domain_error("tan_eval(): infinity")); - + if (is_ex_exactly_of_type(x, function)) { ex t=x.op(0); // tan(atan(x)) -> x @@ -298,10 +330,10 @@ static ex tan_eval(ex const & x) return t; // tan(asin(x)) -> x*(1+x^2)^(-1/2) if (is_ex_the_function(x, asin)) - return t*power(exONE()-power(t,exTWO()),exMINUSHALF()); + return t*power(_ex1()-power(t,_ex2()),_ex_1_2()); // tan(acos(x)) -> (1-x^2)^(1/2)/x if (is_ex_the_function(x, acos)) - return power(t,exMINUSONE())*power(exONE()-power(t,exTWO()),exHALF()); + return power(t,_ex_1())*power(_ex1()-power(t,_ex2()),_ex1_2()); } // tan(float) -> float @@ -316,7 +348,8 @@ static ex tan_diff(ex const & x, unsigned diff_param) { GINAC_ASSERT(diff_param==0); - return (1+power(tan(x),exTWO())); + // d/dx tan(x) -> 1+tan(x)^2; + return (1+power(tan(x),_ex2())); } static ex tan_series(ex const & x, symbol const & s, ex const & point, int order) @@ -326,7 +359,7 @@ static ex tan_series(ex const & x, symbol const & s, ex const & point, int order // On a pole simply expand sin(x)/cos(x). ex xpoint = x.subs(s==point); if (!(2*xpoint/Pi).info(info_flags::odd)) - throw do_taylor(); + throw do_taylor(); // caught by function::series() // if we got here we have to care for a simple pole return (sin(x)/cos(x)).series(s, point, order+2); } @@ -353,17 +386,17 @@ static ex asin_eval(ex const & x) if (x.is_zero()) return x; // asin(1/2) -> Pi/6 - if (x.is_equal(exHALF())) + if (x.is_equal(_ex1_2())) return numeric(1,6)*Pi; // asin(1) -> Pi/2 - if (x.is_equal(exONE())) - return numeric(1,2)*Pi; + if (x.is_equal(_ex1())) + return _num1_2()*Pi; // asin(-1/2) -> -Pi/6 - if (x.is_equal(exMINUSHALF())) + if (x.is_equal(_ex_1_2())) return numeric(-1,6)*Pi; // asin(-1) -> -Pi/2 - if (x.is_equal(exMINUSONE())) - return numeric(-1,2)*Pi; + if (x.is_equal(_ex_1())) + return _num_1_2()*Pi; // asin(float) -> float if (!x.info(info_flags::crational)) return asin_evalf(x); @@ -376,7 +409,8 @@ static ex asin_diff(ex const & x, unsigned diff_param) { GINAC_ASSERT(diff_param==0); - return power(1-power(x,exTWO()),exMINUSHALF()); + // d/dx asin(x) -> 1/sqrt(1-x^2) + return power(1-power(x,_ex2()),_ex_1_2()); } REGISTER_FUNCTION(asin, asin_eval, asin_evalf, asin_diff, NULL); @@ -398,19 +432,19 @@ static ex acos_eval(ex const & x) { if (x.info(info_flags::numeric)) { // acos(1) -> 0 - if (x.is_equal(exONE())) - return exZERO(); + if (x.is_equal(_ex1())) + return _ex0(); // acos(1/2) -> Pi/3 - if (x.is_equal(exHALF())) + if (x.is_equal(_ex1_2())) return numeric(1,3)*Pi; // acos(0) -> Pi/2 if (x.is_zero()) return numeric(1,2)*Pi; // acos(-1/2) -> 2/3*Pi - if (x.is_equal(exMINUSHALF())) + if (x.is_equal(_ex_1_2())) return numeric(2,3)*Pi; // acos(-1) -> Pi - if (x.is_equal(exMINUSONE())) + if (x.is_equal(_ex_1())) return Pi; // acos(float) -> float if (!x.info(info_flags::crational)) @@ -424,7 +458,8 @@ static ex acos_diff(ex const & x, unsigned diff_param) { GINAC_ASSERT(diff_param==0); - return numMINUSONE()*power(1-power(x,exTWO()),exMINUSHALF()); + // d/dx acos(x) -> -1/sqrt(1-x^2) + return _ex_1()*power(1-power(x,_ex2()),_ex_1_2()); } REGISTER_FUNCTION(acos, acos_eval, acos_evalf, acos_diff, NULL); @@ -446,8 +481,8 @@ static ex atan_eval(ex const & x) { if (x.info(info_flags::numeric)) { // atan(0) -> 0 - if (x.is_equal(exZERO())) - return exZERO(); + if (x.is_equal(_ex0())) + return _ex0(); // atan(float) -> float if (!x.info(info_flags::crational)) return atan_evalf(x); @@ -519,27 +554,15 @@ static ex sinh_evalf(ex const & x) static ex sinh_eval(ex const & x) { if (x.info(info_flags::numeric)) { - // sinh(0) -> 0 - if (x.is_zero()) - return exZERO(); - // sinh(float) -> float - if (!x.info(info_flags::crational)) + if (x.is_zero()) // sinh(0) -> 0 + return _ex0(); + if (!x.info(info_flags::crational)) // sinh(float) -> float return sinh_evalf(x); } - ex xOverPiI=x/Pi/I; - if (xOverPiI.info(info_flags::numeric)) { - // sinh(n*Pi*I) -> 0 - if (xOverPiI.info(info_flags::integer)) - return exZERO(); - - // sinh((2n+1)*Pi*I/2) -> {+|-}I - ex xOverPiIMinusHalf=xOverPiI-exHALF(); - if (xOverPiIMinusHalf.info(info_flags::even)) - return I; - else if (xOverPiIMinusHalf.info(info_flags::odd)) - return -I; - } + if ((x/Pi).info(info_flags::numeric) && + ex_to_numeric(x/Pi).real().is_zero()) // sinh(I*x) -> I*sin(x) + return I*sin(x/I); if (is_ex_exactly_of_type(x, function)) { ex t=x.op(0); @@ -548,10 +571,10 @@ static ex sinh_eval(ex const & x) return t; // sinh(acosh(x)) -> (x-1)^(1/2) * (x+1)^(1/2) if (is_ex_the_function(x, acosh)) - return power(t-exONE(),exHALF())*power(t+exONE(),exHALF()); + return power(t-_ex1(),_ex1_2())*power(t+_ex1(),_ex1_2()); // sinh(atanh(x)) -> x*(1-x^2)^(-1/2) if (is_ex_the_function(x, atanh)) - return t*power(exONE()-power(t,exTWO()),exMINUSHALF()); + return t*power(_ex1()-power(t,_ex2()),_ex_1_2()); } return sinh(x).hold(); @@ -561,6 +584,7 @@ static ex sinh_diff(ex const & x, unsigned diff_param) { GINAC_ASSERT(diff_param==0); + // d/dx sinh(x) -> cosh(x) return cosh(x); } @@ -582,27 +606,15 @@ static ex cosh_evalf(ex const & x) static ex cosh_eval(ex const & x) { if (x.info(info_flags::numeric)) { - // cosh(0) -> 1 - if (x.is_zero()) - return exONE(); - // cosh(float) -> float - if (!x.info(info_flags::crational)) + if (x.is_zero()) // cosh(0) -> 1 + return _ex1(); + if (!x.info(info_flags::crational)) // cosh(float) -> float return cosh_evalf(x); } - ex xOverPiI=x/Pi/I; - if (xOverPiI.info(info_flags::numeric)) { - // cosh(n*Pi*I) -> {+|-}1 - if (xOverPiI.info(info_flags::even)) - return exONE(); - else if (xOverPiI.info(info_flags::odd)) - return exMINUSONE(); - - // cosh((2n+1)*Pi*I/2) -> 0 - ex xOverPiIMinusHalf=xOverPiI-exHALF(); - if (xOverPiIMinusHalf.info(info_flags::integer)) - return exZERO(); - } + if ((x/Pi).info(info_flags::numeric) && + ex_to_numeric(x/Pi).real().is_zero()) // cosh(I*x) -> cos(x) + return cos(x/I); if (is_ex_exactly_of_type(x, function)) { ex t=x.op(0); @@ -611,10 +623,10 @@ static ex cosh_eval(ex const & x) return t; // cosh(asinh(x)) -> (1+x^2)^(1/2) if (is_ex_the_function(x, asinh)) - return power(exONE()+power(t,exTWO()),exHALF()); + return power(_ex1()+power(t,_ex2()),_ex1_2()); // cosh(atanh(x)) -> (1-x^2)^(-1/2) if (is_ex_the_function(x, atanh)) - return power(exONE()-power(t,exTWO()),exMINUSHALF()); + return power(_ex1()-power(t,_ex2()),_ex_1_2()); } return cosh(x).hold(); @@ -624,6 +636,7 @@ static ex cosh_diff(ex const & x, unsigned diff_param) { GINAC_ASSERT(diff_param==0); + // d/dx cosh(x) -> sinh(x) return sinh(x); } @@ -645,14 +658,16 @@ static ex tanh_evalf(ex const & x) static ex tanh_eval(ex const & x) { if (x.info(info_flags::numeric)) { - // tanh(0) -> 0 - if (x.is_zero()) - return exZERO(); - // tanh(float) -> float - if (!x.info(info_flags::crational)) + if (x.is_zero()) // tanh(0) -> 0 + return _ex0(); + if (!x.info(info_flags::crational)) // tanh(float) -> float return tanh_evalf(x); } + if ((x/Pi).info(info_flags::numeric) && + ex_to_numeric(x/Pi).real().is_zero()) // tanh(I*x) -> I*tan(x); + return I*tan(x/I); + if (is_ex_exactly_of_type(x, function)) { ex t=x.op(0); // tanh(atanh(x)) -> x @@ -660,10 +675,10 @@ static ex tanh_eval(ex const & x) return t; // tanh(asinh(x)) -> x*(1+x^2)^(-1/2) if (is_ex_the_function(x, asinh)) - return t*power(exONE()+power(t,exTWO()),exMINUSHALF()); + return t*power(_ex1()+power(t,_ex2()),_ex_1_2()); // tanh(acosh(x)) -> (x-1)^(1/2)*(x+1)^(1/2)/x if (is_ex_the_function(x, acosh)) - return power(t-exONE(),exHALF())*power(t+exONE(),exHALF())*power(t,exMINUSONE()); + return power(t-_ex1(),_ex1_2())*power(t+_ex1(),_ex1_2())*power(t,_ex_1()); } return tanh(x).hold(); @@ -673,7 +688,8 @@ static ex tanh_diff(ex const & x, unsigned diff_param) { GINAC_ASSERT(diff_param==0); - return exONE()-power(tanh(x),exTWO()); + // d/dx tanh(x) -> 1-tanh(x)^2 + return _ex1()-power(tanh(x),_ex2()); } static ex tanh_series(ex const & x, symbol const & s, ex const & point, int order) @@ -683,7 +699,7 @@ static ex tanh_series(ex const & x, symbol const & s, ex const & point, int orde // On a pole simply expand sinh(x)/cosh(x). ex xpoint = x.subs(s==point); if (!(2*I*xpoint/Pi).info(info_flags::odd)) - throw do_taylor(); + throw do_taylor(); // caught by function::series() // if we got here we have to care for a simple pole return (sinh(x)/cosh(x)).series(s, point, order+2); } @@ -708,7 +724,7 @@ static ex asinh_eval(ex const & x) if (x.info(info_flags::numeric)) { // asinh(0) -> 0 if (x.is_zero()) - return exZERO(); + return _ex0(); // asinh(float) -> float if (!x.info(info_flags::crational)) return asinh_evalf(x); @@ -721,7 +737,8 @@ static ex asinh_diff(ex const & x, unsigned diff_param) { GINAC_ASSERT(diff_param==0); - return power(1+power(x,exTWO()),exMINUSHALF()); + // d/dx asinh(x) -> 1/sqrt(1+x^2) + return power(1+power(x,_ex2()),_ex_1_2()); } REGISTER_FUNCTION(asinh, asinh_eval, asinh_evalf, asinh_diff, NULL); @@ -746,10 +763,10 @@ static ex acosh_eval(ex const & x) if (x.is_zero()) return Pi*I*numeric(1,2); // acosh(1) -> 0 - if (x.is_equal(exONE())) - return exZERO(); + if (x.is_equal(_ex1())) + return _ex0(); // acosh(-1) -> Pi*I - if (x.is_equal(exMINUSONE())) + if (x.is_equal(_ex_1())) return Pi*I; // acosh(float) -> float if (!x.info(info_flags::crational)) @@ -763,7 +780,8 @@ static ex acosh_diff(ex const & x, unsigned diff_param) { GINAC_ASSERT(diff_param==0); - return power(x-1,exMINUSHALF())*power(x+1,exMINUSHALF()); + // d/dx acosh(x) -> 1/(sqrt(x-1)*sqrt(x+1)) + return power(x+_ex_1(),_ex_1_2())*power(x+_ex1(),_ex_1_2()); } REGISTER_FUNCTION(acosh, acosh_eval, acosh_evalf, acosh_diff, NULL); @@ -786,9 +804,9 @@ static ex atanh_eval(ex const & x) if (x.info(info_flags::numeric)) { // atanh(0) -> 0 if (x.is_zero()) - return exZERO(); + return _ex0(); // atanh({+|-}1) -> throw - if (x.is_equal(exONE()) || x.is_equal(exONE())) + if (x.is_equal(_ex1()) || x.is_equal(_ex1())) throw (std::domain_error("atanh_eval(): infinity")); // atanh(float) -> float if (!x.info(info_flags::crational)) @@ -802,7 +820,8 @@ static ex atanh_diff(ex const & x, unsigned diff_param) { GINAC_ASSERT(diff_param==0); - return power(exONE()-power(x,exTWO()),exMINUSONE()); + // d/dx atanh(x) -> 1/(1-x^2) + return power(_ex1()-power(x,_ex2()),_ex_1()); } REGISTER_FUNCTION(atanh, atanh_eval, atanh_evalf, atanh_diff, NULL); diff --git a/ginac/inifcns_zeta.cpp b/ginac/inifcns_zeta.cpp index 6829d565..9fdce1e8 100644 --- a/ginac/inifcns_zeta.cpp +++ b/ginac/inifcns_zeta.cpp @@ -29,6 +29,7 @@ #include "numeric.h" #include "power.h" #include "symbol.h" +#include "utils.h" #ifndef NO_GINAC_NAMESPACE namespace GiNaC { @@ -54,19 +55,19 @@ static ex zeta1_eval(ex const & x) // trap integer arguments: if (y.is_integer()) { if (y.is_zero()) - return -exHALF(); - if (x.is_equal(exONE())) + return -_ex1_2(); + if (x.is_equal(_ex1())) throw(std::domain_error("zeta(1): infinity")); if (x.info(info_flags::posint)) { if (x.info(info_flags::odd)) return zeta(x).hold(); else - return abs(bernoulli(y))*pow(Pi,x)*numTWO().power(y-numONE())/factorial(y); + return abs(bernoulli(y))*pow(Pi,x)*_num2().power(y-_num1())/factorial(y); } else { if (x.info(info_flags::odd)) - return -bernoulli(numONE()-y)/(numONE()-y); + return -bernoulli(_num1()-y)/(_num1()-y); else - return numZERO(); + return _num0(); } } } @@ -77,7 +78,7 @@ static ex zeta1_diff(ex const & x, unsigned diff_param) { GINAC_ASSERT(diff_param==0); - return zeta(exONE(), x); + return zeta(_ex1(), x); } const unsigned function_index_zeta1 = function::register_new("zeta", zeta1_eval, zeta1_evalf, zeta1_diff, NULL); diff --git a/ginac/lortensor.cpp b/ginac/lortensor.cpp index 7cbd2f43..1854af0b 100644 --- a/ginac/lortensor.cpp +++ b/ginac/lortensor.cpp @@ -35,11 +35,11 @@ #include "flags.h" #include "lst.h" #include "lortensor.h" -#include "utils.h" #include "operators.h" #include "tinfos.h" #include "power.h" #include "symbol.h" +#include "utils.h" #ifndef NO_GINAC_NAMESPACE namespace GiNaC { @@ -254,19 +254,19 @@ ex lortensor::eval(int level) const //both on diagonal if (idx1.get_value()==0){ // (0,0) - return exONE(); + return _ex1(); } else { if (idx1.is_covariant() != idx2.is_covariant()) { // (_i,~i) or (~i,_i), i = 1...3 - return exONE(); + return _ex1(); } else { // (_i,_i) or (~i,~i), i= 1...3 - return exMINUSONE(); + return _ex_1(); } } } else { // at least one off-diagonal - return exZERO(); + return _ex0(); } } else if (idx1.is_symbolic() && idx1.is_co_contra_pair(idx2)) { return Dim()-idx1.get_dim_parallel_space(); @@ -391,12 +391,12 @@ ex simplify_lortensor_mul(ex const & m) v_contracted.reserve(2*n); for (int i=0; iprintcsrc(os,type,precedence); os << "*"; } @@ -234,7 +235,7 @@ void mul::printcsrc(ostream & os, unsigned type, unsigned upper_precedence) cons while (it != itend) { // If the first argument is a negative integer power, it gets printed as "1.0/" - if (it == seq.begin() && ex_to_numeric(it->coeff).is_integer() && it->coeff.compare(numZERO()) < 0) { + if (it == seq.begin() && ex_to_numeric(it->coeff).is_integer() && it->coeff.compare(_num0()) < 0) { if (type == csrc_types::ctype_cl_N) os << "recip("; else @@ -242,7 +243,7 @@ void mul::printcsrc(ostream & os, unsigned type, unsigned upper_precedence) cons } // If the exponent is 1 or -1, it is left out - if (it->coeff.compare(exONE()) == 0 || it->coeff.compare(numMINUSONE()) == 0) + if (it->coeff.compare(_ex1()) == 0 || it->coeff.compare(_num_1()) == 0) it->rest.bp->printcsrc(os, type, precedence); else // outer parens around ex needed for broken gcc-2.95 parser: @@ -251,7 +252,7 @@ void mul::printcsrc(ostream & os, unsigned type, unsigned upper_precedence) cons // Separator is "/" for negative integer powers, "*" otherwise it++; if (it != itend) { - if (ex_to_numeric(it->coeff).is_integer() && it->coeff.compare(numZERO()) < 0) + if (ex_to_numeric(it->coeff).is_integer() && it->coeff.compare(_num0()) < 0) os << "/"; else os << "*"; @@ -335,7 +336,7 @@ ex mul::coeff(symbol const & s, int const n) const return (new mul(coeffseq))->setflag(status_flags::dynallocated); } - return exZERO(); + return _ex0(); } ex mul::eval(int level) const @@ -373,23 +374,23 @@ ex mul::eval(int level) const if (flags & status_flags::evaluated) { GINAC_ASSERT(seq.size()>0); - GINAC_ASSERT((seq.size()>1)||!overall_coeff.is_equal(exONE())); + GINAC_ASSERT((seq.size()>1)||!overall_coeff.is_equal(_ex1())); return *this; } int seq_size=seq.size(); - if (overall_coeff.is_equal(exZERO())) { + if (overall_coeff.is_equal(_ex0())) { // *(...,x;0) -> 0 - return exZERO(); + return _ex0(); } else if (seq_size==0) { // *(;c) -> c return overall_coeff; - } else if ((seq_size==1)&&overall_coeff.is_equal(exONE())) { + } else if ((seq_size==1)&&overall_coeff.is_equal(_ex1())) { // *(x;1) -> x return recombine_pair_to_ex(*(seq.begin())); } else if ((seq_size==1) && is_ex_exactly_of_type((*seq.begin()).rest,add) && - ex_to_numeric((*seq.begin()).coeff).is_equal(numONE())) { + ex_to_numeric((*seq.begin()).coeff).is_equal(_num1())) { // *(+(x,y,...);c) -> +(*(x,c),*(y,c),...) (c numeric(), no powers of +()) add const & addref=ex_to_add((*seq.begin()).rest); epvector distrseq; @@ -503,7 +504,7 @@ expair mul::split_ex_to_pair(ex const & e) const return expair(powerref.basis,powerref.exponent); } } - return expair(e,exONE()); + return expair(e,_ex1()); } expair mul::combine_ex_with_coeff_to_pair(ex const & e, @@ -513,7 +514,7 @@ expair mul::combine_ex_with_coeff_to_pair(ex const & e, // we create a temporary power object // otherwise it would be hard to correctly simplify // expression like (4^(1/3))^(3/2) - if (are_ex_trivially_equal(c,exONE())) { + if (are_ex_trivially_equal(c,_ex1())) { return split_ex_to_pair(e); } return split_ex_to_pair(power(e,c)); @@ -526,7 +527,7 @@ expair mul::combine_pair_with_coeff_to_pair(expair const & p, // we create a temporary power object // otherwise it would be hard to correctly simplify // expression like (4^(1/3))^(3/2) - if (are_ex_trivially_equal(c,exONE())) { + if (are_ex_trivially_equal(c,_ex1())) { return p; } return split_ex_to_pair(power(recombine_pair_to_ex(p),c)); @@ -534,9 +535,9 @@ expair mul::combine_pair_with_coeff_to_pair(expair const & p, ex mul::recombine_pair_to_ex(expair const & p) const { - // if (p.coeff.compare(exONE())==0) { - // if (are_ex_trivially_equal(p.coeff,exONE())) { - if (ex_to_numeric(p.coeff).is_equal(numONE())) { + // if (p.coeff.compare(_ex1())==0) { + // if (are_ex_trivially_equal(p.coeff,_ex1())) { + if (ex_to_numeric(p.coeff).is_equal(_num1())) { return p.rest; } else { return power(p.rest,p.coeff); @@ -558,7 +559,7 @@ bool mul::expair_needs_further_processing(epp it) *it=ep; return true; } - if (ex_to_numeric((*it).coeff).is_equal(numONE())) { + if (ex_to_numeric((*it).coeff).is_equal(_num1())) { // combined pair has coeff 1 and must be moved to the end return true; } @@ -568,7 +569,7 @@ bool mul::expair_needs_further_processing(epp it) ex mul::default_overall_coeff(void) const { - return exONE(); + return _ex1(); } void mul::combine_overall_coeff(ex const & c) @@ -593,7 +594,7 @@ bool mul::can_make_flat(expair const & p) const // this assertion will probably fail somewhere // it would require a more careful make_flat, obeying the power laws // probably should return true only if p.coeff is integer - return ex_to_numeric(p.coeff).is_equal(numONE()); + return ex_to_numeric(p.coeff).is_equal(_num1()); } ex mul::expand(unsigned options) const @@ -616,7 +617,7 @@ ex mul::expand(unsigned options) const epvector::const_iterator last=expanded_seq.end(); for (epvector::const_iterator cit=expanded_seq.begin(); cit!=last; ++cit) { if (is_ex_exactly_of_type((*cit).rest,add)&& - (ex_to_numeric((*cit).coeff).is_equal(numONE()))) { + (ex_to_numeric((*cit).coeff).is_equal(_num1()))) { positions_of_adds[number_of_adds]=current_position; add const & expanded_addref=ex_to_add((*cit).rest); int addref_nops=expanded_addref.nops(); @@ -652,7 +653,7 @@ ex mul::expand(unsigned options) const term=expanded_seq; for (l=0; lsetflag(status_flags::dynallocated); - return exZERO(); + return _ex0(); } unsigned ncmul::count_factors(ex const & e) const @@ -374,7 +375,7 @@ ex ncmul::eval(int level) const if (assocseq.size()==1) return *(seq.begin()); // ncmul() -> 1 - if (assocseq.size()==0) return exONE(); + if (assocseq.size()==0) return _ex1(); // determine return types unsignedvector rettypes; @@ -618,7 +619,7 @@ ex nonsimplified_ncmul(exvector const & v) ex simplified_ncmul(exvector const & v) { if (v.size()==0) { - return exONE(); + return _ex1(); } else if (v.size()==1) { return v[0]; } diff --git a/ginac/normal.cpp b/ginac/normal.cpp index 9a3cd159..d6661d99 100644 --- a/ginac/normal.cpp +++ b/ginac/normal.cpp @@ -45,6 +45,7 @@ #include "relational.h" #include "series.h" #include "symbol.h" +#include "utils.h" #ifndef NO_GINAC_NAMESPACE namespace GiNaC { @@ -190,7 +191,7 @@ static numeric lcmcoeff(const ex &e, const numeric &l) if (e.info(info_flags::rational)) return lcm(ex_to_numeric(e).denom(), l); else if (is_ex_exactly_of_type(e, add) || is_ex_exactly_of_type(e, mul)) { - numeric c = numONE(); + numeric c = _num1(); for (int i=0; irest,numeric)); GINAC_ASSERT(is_ex_exactly_of_type(it->coeff,numeric)); @@ -289,13 +290,13 @@ ex quo(const ex &a, const ex &b, const symbol &x, bool check_args) return a / b; #if FAST_COMPARE if (a.is_equal(b)) - return exONE(); + return _ex1(); #endif if (check_args && (!a.info(info_flags::rational_polynomial) || !b.info(info_flags::rational_polynomial))) throw(std::invalid_argument("quo: arguments must be polynomials over the rationals")); // Polynomial long division - ex q = exZERO(); + ex q = _ex0(); ex r = a.expand(); if (r.is_zero()) return r; @@ -338,13 +339,13 @@ ex rem(const ex &a, const ex &b, const symbol &x, bool check_args) throw(std::overflow_error("rem: division by zero")); if (is_ex_exactly_of_type(a, numeric)) { if (is_ex_exactly_of_type(b, numeric)) - return exZERO(); + return _ex0(); else return b; } #if FAST_COMPARE if (a.is_equal(b)) - return exZERO(); + return _ex0(); #endif if (check_args && (!a.info(info_flags::rational_polynomial) || !b.info(info_flags::rational_polynomial))) throw(std::invalid_argument("rem: arguments must be polynomials over the rationals")); @@ -390,7 +391,7 @@ ex prem(const ex &a, const ex &b, const symbol &x, bool check_args) throw(std::overflow_error("prem: division by zero")); if (is_ex_exactly_of_type(a, numeric)) { if (is_ex_exactly_of_type(b, numeric)) - return exZERO(); + return _ex0(); else return b; } @@ -406,18 +407,18 @@ ex prem(const ex &a, const ex &b, const symbol &x, bool check_args) if (bdeg <= rdeg) { blcoeff = eb.coeff(x, bdeg); if (bdeg == 0) - eb = exZERO(); + eb = _ex0(); else eb -= blcoeff * power(x, bdeg); } else - blcoeff = exONE(); + blcoeff = _ex1(); int delta = rdeg - bdeg + 1, i = 0; while (rdeg >= bdeg && !r.is_zero()) { ex rlcoeff = r.coeff(x, rdeg); ex term = (power(x, rdeg - bdeg) * eb * rlcoeff).expand(); if (rdeg == 0) - r = exZERO(); + r = _ex0(); else r -= rlcoeff * power(x, rdeg); r = (blcoeff * r).expand() - term; @@ -440,7 +441,7 @@ ex prem(const ex &a, const ex &b, const symbol &x, bool check_args) bool divide(const ex &a, const ex &b, ex &q, bool check_args) { - q = exZERO(); + q = _ex0(); if (b.is_zero()) throw(std::overflow_error("divide: division by zero")); if (is_ex_exactly_of_type(b, numeric)) { @@ -450,7 +451,7 @@ bool divide(const ex &a, const ex &b, ex &q, bool check_args) return false; #if FAST_COMPARE if (a.is_equal(b)) { - q = exONE(); + q = _ex1(); return true; } #endif @@ -525,10 +526,10 @@ typedef map ex2_exbool_remember; * @see get_symbol_stats, heur_gcd */ static bool divide_in_z(const ex &a, const ex &b, ex &q, sym_desc_vec::const_iterator var) { - q = exZERO(); + q = _ex0(); if (b.is_zero()) throw(std::overflow_error("divide_in_z: division by zero")); - if (b.is_equal(exONE())) { + if (b.is_equal(_ex1())) { q = a; return true; } @@ -541,7 +542,7 @@ static bool divide_in_z(const ex &a, const ex &b, ex &q, sym_desc_vec::const_ite } #if FAST_COMPARE if (a.is_equal(b)) { - q = exONE(); + q = _ex1(); return true; } #endif @@ -601,19 +602,19 @@ static bool divide_in_z(const ex &a, const ex &b, ex &q, sym_desc_vec::const_ite // Compute values at evaluation points 0..adeg vector alpha; alpha.reserve(adeg + 1); exvector u; u.reserve(adeg + 1); - numeric point = numZERO(); + numeric point = _num0(); ex c; for (i=0; i<=adeg; i++) { ex bs = b.subs(*x == point); while (bs.is_zero()) { - point += numONE(); + point += _num1(); bs = b.subs(*x == point); } if (!divide_in_z(a.subs(*x == point), bs, c, var+1)) return false; alpha.push_back(point); u.push_back(c); - point += numONE(); + point += _num1(); } // Compute inverses @@ -665,7 +666,7 @@ ex ex::unit(const symbol &x) const { ex c = expand().lcoeff(x); if (is_ex_exactly_of_type(c, numeric)) - return c < exZERO() ? exMINUSONE() : exONE(); + return c < _ex0() ? _ex_1() : _ex1(); else { const symbol *y; if (get_first_symbol(c, y)) @@ -686,12 +687,12 @@ ex ex::unit(const symbol &x) const ex ex::content(const symbol &x) const { if (is_zero()) - return exZERO(); + return _ex0(); if (is_ex_exactly_of_type(*this, numeric)) return info(info_flags::negative) ? -*this : *this; ex e = expand(); if (e.is_zero()) - return exZERO(); + return _ex0(); // First, try the integer content ex c = e.integer_content(); @@ -705,7 +706,7 @@ ex ex::content(const symbol &x) const int ldeg = e.ldegree(x); if (deg == ldeg) return e.lcoeff(x) / e.unit(x); - c = exZERO(); + c = _ex0(); for (int i=ldeg; i<=deg; i++) c = gcd(e.coeff(x, i), c, NULL, NULL, false); return c; @@ -722,13 +723,13 @@ ex ex::content(const symbol &x) const ex ex::primpart(const symbol &x) const { if (is_zero()) - return exZERO(); + return _ex0(); if (is_ex_exactly_of_type(*this, numeric)) - return exONE(); + return _ex1(); ex c = content(x); if (c.is_zero()) - return exZERO(); + return _ex0(); ex u = unit(x); if (is_ex_exactly_of_type(c, numeric)) return *this / (c * u); @@ -748,11 +749,11 @@ ex ex::primpart(const symbol &x) const ex ex::primpart(const symbol &x, const ex &c) const { if (is_zero()) - return exZERO(); + return _ex0(); if (c.is_zero()) - return exZERO(); + return _ex0(); if (is_ex_exactly_of_type(*this, numeric)) - return exONE(); + return _ex1(); ex u = unit(x); if (is_ex_exactly_of_type(c, numeric)) @@ -803,7 +804,7 @@ static ex sr_gcd(const ex &a, const ex &b, const symbol *x) d = d.primpart(*x, cont_d); // First element of subresultant sequence - ex r = exZERO(), ri = exONE(), psi = exONE(); + ex r = _ex0(), ri = _ex1(), psi = _ex1(); int delta = cdeg - ddeg; for (;;) { @@ -849,7 +850,7 @@ numeric ex::max_coefficient(void) const numeric basic::max_coefficient(void) const { - return numONE(); + return _num1(); } numeric numeric::max_coefficient(void) const @@ -1013,9 +1014,9 @@ static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const numeric mp = p.max_coefficient(), mq = q.max_coefficient(); numeric xi; if (mp > mq) - xi = mq * numTWO() + numTWO(); + xi = mq * _num2() + _num2(); else - xi = mp * numTWO() + numTWO(); + xi = mp * _num2() + _num2(); // 6 tries maximum for (int t=0; t<6; t++) { @@ -1027,7 +1028,7 @@ static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const if (!is_ex_exactly_of_type(gamma, fail)) { // Reconstruct polynomial from GCD of mapped polynomials - ex g = exZERO(); + ex g = _ex0(); numeric rxi = xi.inverse(); for (int i=0; !gamma.is_zero(); i++) { ex gi = gamma.smod(xi); @@ -1042,7 +1043,7 @@ static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const if (divide_in_z(p, g, ca ? *ca : dummy, var) && divide_in_z(q, g, cb ? *cb : dummy, var)) { g *= gc; ex lc = g.lcoeff(*x); - if (is_ex_exactly_of_type(lc, numeric) && lc.compare(exZERO()) < 0) + if (is_ex_exactly_of_type(lc, numeric) && lc.compare(_ex0()) < 0) return -g; else return g; @@ -1071,31 +1072,31 @@ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args) ex aex = a.expand(), bex = b.expand(); if (aex.is_zero()) { if (ca) - *ca = exZERO(); + *ca = _ex0(); if (cb) - *cb = exONE(); + *cb = _ex1(); return b; } if (bex.is_zero()) { if (ca) - *ca = exONE(); + *ca = _ex1(); if (cb) - *cb = exZERO(); + *cb = _ex0(); return a; } - if (aex.is_equal(exONE()) || bex.is_equal(exONE())) { + if (aex.is_equal(_ex1()) || bex.is_equal(_ex1())) { if (ca) *ca = a; if (cb) *cb = b; - return exONE(); + return _ex1(); } #if FAST_COMPARE if (a.is_equal(b)) { if (ca) - *ca = exONE(); + *ca = _ex1(); if (cb) - *cb = exONE(); + *cb = _ex1(); return a; } #endif @@ -1154,7 +1155,7 @@ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args) g = *new ex(fail()); } if (is_ex_exactly_of_type(g, fail)) { -//clog << "heuristics failed\n"; +// clog << "heuristics failed" << endl; g = sr_gcd(aex, bex, x); if (ca) divide(aex, g, *ca, false); @@ -1197,8 +1198,8 @@ static ex univariate_gcd(const ex &a, const ex &b, const symbol &x) return b; if (b.is_zero()) return a; - if (a.is_equal(exONE()) || b.is_equal(exONE())) - return exONE(); + if (a.is_equal(_ex1()) || b.is_equal(_ex1())) + return _ex1(); if (is_ex_of_type(a, numeric) && is_ex_of_type(b, numeric)) return gcd(ex_to_numeric(a), ex_to_numeric(b)); if (!a.info(info_flags::rational_polynomial) || !b.info(info_flags::rational_polynomial)) @@ -1233,11 +1234,11 @@ static ex univariate_gcd(const ex &a, const ex &b, const symbol &x) ex sqrfree(const ex &a, const symbol &x) { int i = 1; - ex res = exONE(); + ex res = _ex1(); ex b = a.diff(x); ex c = univariate_gcd(a, b, x); ex w; - if (c.is_equal(exONE())) { + if (c.is_equal(_ex1())) { w = a; } else { w = quo(a, c, x); @@ -1326,11 +1327,11 @@ static ex frac_cancel(const ex &n, const ex &d) { ex num = n; ex den = d; - ex pre_factor = exONE(); + ex pre_factor = _ex1(); // Handle special cases where numerator or denominator is 0 if (num.is_zero()) - return exZERO(); + return _ex0(); if (den.expand().is_zero()) throw(std::overflow_error("frac_cancel: division by zero in frac_cancel")); @@ -1338,7 +1339,7 @@ static ex frac_cancel(const ex &n, const ex &d) if (is_ex_exactly_of_type(den, numeric)) return num / den; if (num.is_zero()) - return exZERO(); + return _ex0(); // Bring numerator and denominator to Z[X] by multiplying with // LCM of all coefficients' denominators @@ -1350,7 +1351,7 @@ static ex frac_cancel(const ex &n, const ex &d) // Cancel GCD from numerator and denominator ex cnum, cden; - if (gcd(num, den, &cnum, &cden, false) != exONE()) { + if (gcd(num, den, &cnum, &cden, false) != _ex1()) { num = cnum; den = cden; } @@ -1359,9 +1360,9 @@ static ex frac_cancel(const ex &n, const ex &d) // as defined by get_first_symbol() is made positive) const symbol *x; if (get_first_symbol(den, x)) { - if (den.unit(*x).compare(exZERO()) < 0) { - num *= exMINUSONE(); - den *= exMINUSONE(); + if (den.unit(*x).compare(_ex0()) < 0) { + num *= _ex_1(); + den *= _ex_1(); } } return pre_factor * num / den; @@ -1393,7 +1394,7 @@ ex add::normal(lst &sym_lst, lst &repl_lst, int level) const o.push_back(overall_coeff.bp->normal(sym_lst, repl_lst, level-1)); // Determine common denominator - ex den = exONE(); + ex den = _ex1(); exvector::const_iterator ait = o.begin(), aitend = o.end(); while (ait != aitend) { den = lcm((*ait).denom(false), den, false); @@ -1401,7 +1402,7 @@ ex add::normal(lst &sym_lst, lst &repl_lst, int level) const } // Add fractions - if (den.is_equal(exONE())) + if (den.is_equal(_ex1())) return (new add(o))->setflag(status_flags::dynallocated); else { exvector num_seq; diff --git a/ginac/numeric.cpp b/ginac/numeric.cpp index b5cc40b3..472d9c31 100644 --- a/ginac/numeric.cpp +++ b/ginac/numeric.cpp @@ -31,6 +31,7 @@ #include "ex.h" #include "config.h" #include "debugmsg.h" +#include "utils.h" // CLN should not pollute the global namespace, hence we include it here // instead of in some header file where it would propagate to other parts: @@ -290,7 +291,7 @@ void numeric::printcsrc(ostream & os, unsigned type, unsigned upper_precedence) ios::fmtflags oldflags = os.flags(); os.setf(ios::scientific); if (is_rational() && !is_integer()) { - if (compare(numZERO()) > 0) { + if (compare(_num0()) > 0) { os << "("; if (type == csrc_types::ctype_cl_N) os << "cl_F(\"" << numer().evalf() << "\")"; @@ -344,11 +345,11 @@ bool numeric::info(unsigned inf) const case info_flags::negative: return is_negative(); case info_flags::nonnegative: - return compare(numZERO())>=0; + return compare(_num0())>=0; case info_flags::posint: return is_pos_integer(); case info_flags::negint: - return is_integer() && (compare(numZERO())<0); + return is_integer() && (compare(_num0())<0); case info_flags::nonnegint: return is_nonneg_integer(); case info_flags::even: @@ -439,10 +440,10 @@ numeric numeric::sub(numeric const & other) const * result as a new numeric object. */ numeric numeric::mul(numeric const & other) const { - static const numeric * numONEp=&numONE(); - if (this==numONEp) { + static const numeric * _num1p=&_num1(); + if (this==_num1p) { return other; - } else if (&other==numONEp) { + } else if (&other==_num1p) { return *this; } return numeric((*value)*(*other.value)); @@ -461,8 +462,8 @@ numeric numeric::div(numeric const & other) const numeric numeric::power(numeric const & other) const { - static const numeric * numONEp=&numONE(); - if (&other==numONEp) + static const numeric * _num1p=&_num1(); + if (&other==_num1p) return *this; if (::zerop(*value) && other.is_real() && ::minusp(realpart(*other.value))) throw (std::overflow_error("division by zero")); @@ -489,10 +490,10 @@ numeric const & numeric::sub_dyn(numeric const & other) const numeric const & numeric::mul_dyn(numeric const & other) const { - static const numeric * numONEp=&numONE(); - if (this==numONEp) { + static const numeric * _num1p=&_num1(); + if (this==_num1p) { return other; - } else if (&other==numONEp) { + } else if (&other==_num1p) { return *this; } return static_cast((new numeric((*value)*(*other.value)))-> @@ -509,8 +510,8 @@ numeric const & numeric::div_dyn(numeric const & other) const numeric const & numeric::power_dyn(numeric const & other) const { - static const numeric * numONEp=&numONE(); - if (&other==numONEp) + static const numeric * _num1p=&_num1(); + if (&other==_num1p) return *this; if (::zerop(*value) && other.is_real() && ::minusp(realpart(*other.value))) throw (std::overflow_error("division by zero")); @@ -852,7 +853,7 @@ numeric numeric::numer(void) const numeric numeric::denom(void) const { if (is_integer()) { - return numONE(); + return _num1(); } #ifdef SANE_LINKER if (instanceof(*value, cl_RA_ring)) { @@ -862,7 +863,7 @@ numeric numeric::denom(void) const cl_R r = realpart(*value); cl_R i = imagpart(*value); if (::instanceof(r, cl_I_ring) && ::instanceof(i, cl_I_ring)) - return numONE(); + return _num1(); if (::instanceof(r, cl_I_ring) && ::instanceof(i, cl_RA_ring)) return numeric(::denominator(The(cl_RA)(i))); if (::instanceof(r, cl_RA_ring) && ::instanceof(i, cl_I_ring)) @@ -878,7 +879,7 @@ numeric numeric::denom(void) const cl_R r = realpart(*value); cl_R i = imagpart(*value); if (instanceof(r, cl_I_ring) && instanceof(i, cl_I_ring)) - return numONE(); + return _num1(); if (instanceof(r, cl_I_ring) && instanceof(i, cl_RA_ring)) return numeric(TheRatio(i)->denominator); if (instanceof(r, cl_RA_ring) && instanceof(i, cl_I_ring)) @@ -888,7 +889,7 @@ numeric numeric::denom(void) const } #endif // def SANE_LINKER // at least one float encountered - return numONE(); + return _num1(); } /** Size in binary notation. For integers, this is the smallest n >= 0 such @@ -924,52 +925,6 @@ type_info const & typeid_numeric=typeid(some_numeric); * natively handing complex numbers anyways. */ const numeric I = numeric(complex(cl_I(0),cl_I(1))); -////////// -// global functions -////////// - -numeric const & numZERO(void) -{ - const static ex eZERO = ex((new numeric(0))->setflag(status_flags::dynallocated)); - const static numeric * nZERO = static_cast(eZERO.bp); - return *nZERO; -} - -numeric const & numONE(void) -{ - const static ex eONE = ex((new numeric(1))->setflag(status_flags::dynallocated)); - const static numeric * nONE = static_cast(eONE.bp); - return *nONE; -} - -numeric const & numTWO(void) -{ - const static ex eTWO = ex((new numeric(2))->setflag(status_flags::dynallocated)); - const static numeric * nTWO = static_cast(eTWO.bp); - return *nTWO; -} - -numeric const & numTHREE(void) -{ - const static ex eTHREE = ex((new numeric(3))->setflag(status_flags::dynallocated)); - const static numeric * nTHREE = static_cast(eTHREE.bp); - return *nTHREE; -} - -numeric const & numMINUSONE(void) -{ - const static ex eMINUSONE = ex((new numeric(-1))->setflag(status_flags::dynallocated)); - const static numeric * nMINUSONE = static_cast(eMINUSONE.bp); - return *nMINUSONE; -} - -numeric const & numHALF(void) -{ - const static ex eHALF = ex((new numeric(1, 2))->setflag(status_flags::dynallocated)); - const static numeric * nHALF = static_cast(eHALF.bp); - return *nHALF; -} - /** Exponential function. * * @return arbitrary precision numerical exp(x). */ @@ -1039,7 +994,7 @@ numeric atan(numeric const & x) { if (!x.is_real() && x.real().is_zero() && - !abs(x.imag()).is_equal(numONE())) + !abs(x.imag()).is_equal(_num1())) throw (std::overflow_error("atan(): logarithmic singularity")); return ::atan(*x.value); // -> CLN } @@ -1189,13 +1144,13 @@ numeric doublefactorial(numeric const & nn) static int highest_oddresult = -1; if (nn == numeric(-1)) { - return numONE(); + return _num1(); } if (!nn.is_nonneg_integer()) { throw (std::range_error("numeric::doublefactorial(): argument must be integer >= -1")); } if (nn.is_even()) { - int n = nn.div(numTWO()).to_int(); + int n = nn.div(_num2()).to_int(); if (n <= highest_evenresult) { return evenresults[n]; } @@ -1203,7 +1158,7 @@ numeric doublefactorial(numeric const & nn) evenresults.reserve(n+1); } if (highest_evenresult < 0) { - evenresults.push_back(numONE()); + evenresults.push_back(_num1()); highest_evenresult=0; } for (int i=highest_evenresult+1; i<=n; i++) { @@ -1212,7 +1167,7 @@ numeric doublefactorial(numeric const & nn) highest_evenresult=n; return evenresults[n]; } else { - int n = nn.sub(numONE()).div(numTWO()).to_int(); + int n = nn.sub(_num1()).div(_num2()).to_int(); if (n <= highest_oddresult) { return oddresults[n]; } @@ -1220,7 +1175,7 @@ numeric doublefactorial(numeric const & nn) oddresults.reserve(n+1); } if (highest_oddresult < 0) { - oddresults.push_back(numONE()); + oddresults.push_back(_num1()); highest_oddresult=0; } for (int i=highest_oddresult+1; i<=n; i++) { @@ -1239,12 +1194,12 @@ numeric binomial(numeric const & n, numeric const & k) { if (n.is_integer() && k.is_integer()) { if (n.is_nonneg_integer()) { - if (k.compare(n)!=1 && k.compare(numZERO())!=-1) + if (k.compare(n)!=1 && k.compare(_num0())!=-1) return numeric(::binomial(n.to_int(),k.to_int())); // -> CLN else - return numZERO(); + return _num0(); } else { - return numMINUSONE().power(k)*binomial(k-n-numONE(),k); + return _num_1().power(k)*binomial(k-n-_num1(),k); } } @@ -1262,11 +1217,11 @@ numeric bernoulli(numeric const & nn) if (!nn.is_integer() || nn.is_negative()) throw (std::range_error("numeric::bernoulli(): argument must be integer >= 0")); if (nn.is_zero()) - return numONE(); - if (!nn.compare(numONE())) + return _num1(); + if (!nn.compare(_num1())) return numeric(-1,2); if (nn.is_odd()) - return numZERO(); + return _num0(); // Until somebody has the Blues and comes up with a much better idea and // codes it (preferably in CLN) we make this a remembering function which // computes its results using the formula @@ -1274,7 +1229,7 @@ numeric bernoulli(numeric const & nn) // whith B(0) == 1. static vector results; static int highest_result = -1; - int n = nn.sub(numTWO()).div(numTWO()).to_int(); + int n = nn.sub(_num2()).div(_num2()).to_int(); if (n <= highest_result) return results[n]; if (results.capacity() < (unsigned)(n+1)) @@ -1312,7 +1267,7 @@ numeric mod(numeric const & a, numeric const & b) if (a.is_integer() && b.is_integer()) return ::mod(The(cl_I)(*a.value), The(cl_I)(*b.value)); // -> CLN else - return numZERO(); // Throw? + return _num0(); // Throw? } /** Modulus (in symmetric representation). @@ -1321,11 +1276,12 @@ numeric mod(numeric const & a, numeric const & b) * @return a mod b in the range [-iquo(abs(m)-1,2), iquo(abs(m),2)]. */ numeric smod(numeric const & a, numeric const & b) { + // FIXME: Should this become a member function? if (a.is_integer() && b.is_integer()) { cl_I b2 = The(cl_I)(ceiling1(The(cl_I)(*b.value) / 2)) - 1; return ::mod(The(cl_I)(*a.value) + b2, The(cl_I)(*b.value)) - b2; } else - return numZERO(); // Throw? + return _num0(); // Throw? } /** Numeric integer remainder. @@ -1339,7 +1295,7 @@ numeric irem(numeric const & a, numeric const & b) if (a.is_integer() && b.is_integer()) return ::rem(The(cl_I)(*a.value), The(cl_I)(*b.value)); // -> CLN else - return numZERO(); // Throw? + return _num0(); // Throw? } /** Numeric integer remainder. @@ -1357,8 +1313,8 @@ numeric irem(numeric const & a, numeric const & b, numeric & q) return rem_quo.remainder; } else { - q = numZERO(); - return numZERO(); // Throw? + q = _num0(); + return _num0(); // Throw? } } @@ -1371,7 +1327,7 @@ numeric iquo(numeric const & a, numeric const & b) if (a.is_integer() && b.is_integer()) return truncate1(The(cl_I)(*a.value), The(cl_I)(*b.value)); // -> CLN else - return numZERO(); // Throw? + return _num0(); // Throw? } /** Numeric integer quotient. @@ -1387,8 +1343,8 @@ numeric iquo(numeric const & a, numeric const & b, numeric & r) r = rem_quo.remainder; return rem_quo.quotient; } else { - r = numZERO(); - return numZERO(); // Throw? + r = _num0(); + return _num0(); // Throw? } } @@ -1413,7 +1369,7 @@ numeric isqrt(numeric const & x) ::isqrt(The(cl_I)(*x.value), &root); // -> CLN return root; } else - return numZERO(); // Throw? + return _num0(); // Throw? } /** Greatest Common Divisor. @@ -1425,7 +1381,7 @@ numeric gcd(numeric const & a, numeric const & b) if (a.is_integer() && b.is_integer()) return ::gcd(The(cl_I)(*a.value), The(cl_I)(*b.value)); // -> CLN else - return numONE(); + return _num1(); } /** Least Common Multiple. diff --git a/ginac/numeric.h b/ginac/numeric.h index 68c227d1..27153df3 100644 --- a/ginac/numeric.h +++ b/ginac/numeric.h @@ -92,12 +92,6 @@ class numeric : public basic friend numeric isqrt(numeric const & x); friend numeric gcd(numeric const & a, numeric const & b); friend numeric lcm(numeric const & a, numeric const & b); - friend numeric const & numZERO(void); - friend numeric const & numONE(void); - friend numeric const & numTWO(void); - friend numeric const & numTHREE(void); - friend numeric const & numMINUSONE(void); - friend numeric const & numHALF(void); // member functions @@ -223,12 +217,6 @@ extern _numeric_digits Digits; // global functions -numeric const & numZERO(void); -numeric const & numONE(void); -numeric const & numTWO(void); -numeric const & numMINUSONE(void); -numeric const & numHALF(void); - numeric exp(numeric const & x); numeric log(numeric const & x); numeric sin(numeric const & x); diff --git a/ginac/operators.cpp b/ginac/operators.cpp index c21e0ea6..0fc42e78 100644 --- a/ginac/operators.cpp +++ b/ginac/operators.cpp @@ -30,6 +30,7 @@ #include "power.h" #include "relational.h" #include "debugmsg.h" +#include "utils.h" #ifndef NO_GINAC_NAMESPACE namespace GiNaC { @@ -46,7 +47,7 @@ ex operator+(ex const & lh, ex const & rh) ex operator-(ex const & lh, ex const & rh) { debugmsg("operator-(ex,ex)",LOGLEVEL_OPERATOR); - return lh.exadd(rh.exmul(exMINUSONE())); + return lh.exadd(rh.exmul(_ex_1())); } ex operator*(ex const & lh, ex const & rh) @@ -58,7 +59,7 @@ ex operator*(ex const & lh, ex const & rh) ex operator/(ex const & lh, ex const & rh) { debugmsg("operator*(ex,ex)",LOGLEVEL_OPERATOR); - return lh.exmul(power(rh,exMINUSONE())); + return lh.exmul(power(rh,_ex_1())); } ex operator%(ex const & lh, ex const & rh) @@ -67,73 +68,6 @@ ex operator%(ex const & lh, ex const & rh) return lh.exncmul(rh); } -/* - -// binary arithmetic operators ex with numeric - -ex operator+(ex const & lh, numeric const & rh) -{ - debugmsg("operator+(ex,numeric)",LOGLEVEL_OPERATOR); - return lh+ex(rh); -} - -ex operator-(ex const & lh, numeric const & rh) -{ - debugmsg("operator-(ex,numeric)",LOGLEVEL_OPERATOR); - return lh-ex(rh); -} - -ex operator*(ex const & lh, numeric const & rh) -{ - debugmsg("operator*(ex,numeric)",LOGLEVEL_OPERATOR); - return lh*ex(rh); -} - -ex operator/(ex const & lh, numeric const & rh) -{ - debugmsg("operator/(ex,numeric)",LOGLEVEL_OPERATOR); - return lh/ex(rh); -} - -ex operator%(ex const & lh, numeric const & rh) -{ - debugmsg("operator%(ex,numeric)",LOGLEVEL_OPERATOR); - return lh%ex(rh); -} - -// binary arithmetic operators numeric with ex - -ex operator+(numeric const & lh, ex const & rh) -{ - debugmsg("operator+(numeric,ex)",LOGLEVEL_OPERATOR); - return ex(lh)+rh; -} - -ex operator-(numeric const & lh, ex const & rh) -{ - debugmsg("operator-(numeric,ex)",LOGLEVEL_OPERATOR); - return ex(lh)-rh; -} - -ex operator*(numeric const & lh, ex const & rh) -{ - debugmsg("operator*(numeric,ex)",LOGLEVEL_OPERATOR); - return ex(lh)*rh; -} - -ex operator/(numeric const & lh, ex const & rh) -{ - debugmsg("operator/(numeric,ex)",LOGLEVEL_OPERATOR); - return ex(lh)/rh; -} - -ex operator%(numeric const & lh, ex const & rh) -{ - debugmsg("operator%(numeric,ex)",LOGLEVEL_OPERATOR); - return ex(lh)%rh; -} - -*/ // binary arithmetic operators numeric with numeric @@ -161,6 +95,7 @@ numeric operator/(numeric const & lh, numeric const & rh) return lh.div(rh); } + // binary arithmetic assignment operators with ex ex const & operator+=(ex & lh, ex const & rh) @@ -193,41 +128,6 @@ ex const & operator%=(ex & lh, ex const & rh) return (lh=lh%rh); } -/* - -// binary arithmetic assignment operators with numeric - -ex const & operator+=(ex & lh, numeric const & rh) -{ - debugmsg("operator+=(ex,numeric)",LOGLEVEL_OPERATOR); - return (lh=lh+ex(rh)); -} - -ex const & operator-=(ex & lh, numeric const & rh) -{ - debugmsg("operator-=(ex,numeric)",LOGLEVEL_OPERATOR); - return (lh=lh-ex(rh)); -} - -ex const & operator*=(ex & lh, numeric const & rh) -{ - debugmsg("operator*=(ex,numeric)",LOGLEVEL_OPERATOR); - return (lh=lh*ex(rh)); -} - -ex const & operator/=(ex & lh, numeric const & rh) -{ - debugmsg("operator/=(ex,numeric)",LOGLEVEL_OPERATOR); - return (lh=lh/ex(rh)); -} - -ex const & operator%=(ex & lh, numeric const & rh) -{ - debugmsg("operator%=(ex,numeric)",LOGLEVEL_OPERATOR); - return (lh=lh%ex(rh)); -} - -*/ // binary arithmetic assignment operators with numeric @@ -264,7 +164,7 @@ ex operator+(ex const & lh) ex operator-(ex const & lh) { - return exMINUSONE()*lh; + return lh.exmul(_ex_1()); } numeric operator+(numeric const & lh) @@ -274,20 +174,20 @@ numeric operator+(numeric const & lh) numeric operator-(numeric const & lh) { - return numMINUSONE()*lh; + return _num_1()*lh; } /** Numeric prefix increment. Adds 1 and returns incremented number. */ numeric& operator++(numeric & rh) { - rh = rh+numONE(); + rh = rh+_num1(); return rh; } /** Numeric prefix decrement. Subtracts 1 and returns decremented number. */ numeric& operator--(numeric & rh) { - rh = rh-numONE(); + rh = rh-_num1(); return rh; } @@ -296,7 +196,7 @@ numeric& operator--(numeric & rh) numeric operator++(numeric & lh, int) { numeric tmp = lh; - lh = lh+numONE(); + lh = lh+_num1(); return tmp; } @@ -305,7 +205,7 @@ numeric operator++(numeric & lh, int) numeric operator--(numeric & lh, int) { numeric tmp = lh; - lh = lh-numONE(); + lh = lh-_num1(); return tmp; } diff --git a/ginac/power.cpp b/ginac/power.cpp index 1c9a342d..51e5bfc5 100644 --- a/ginac/power.cpp +++ b/ginac/power.cpp @@ -32,6 +32,7 @@ #include "relational.h" #include "symbol.h" #include "debugmsg.h" +#include "utils.h" #ifndef NO_GINAC_NAMESPACE namespace GiNaC { @@ -192,7 +193,7 @@ void power::printcsrc(ostream & os, unsigned type, unsigned upper_precedence) co os << ")"; // ^-1 is printed as "1.0/" or with the recip() function of CLN - } else if (exponent.compare(numMINUSONE()) == 0) { + } else if (exponent.compare(_num_1()) == 0) { if (type == csrc_types::ctype_cl_N) os << "recip("; else @@ -244,7 +245,7 @@ ex & power::let_op(int const i) int power::degree(symbol const & s) const { if (is_exactly_of_type(*exponent.bp,numeric)) { - if ((*basis.bp).compare(s)==0) + if ((*basis.bp).compare(s)==0) return ex_to_numeric(exponent).to_int(); else return basis.degree(s) * ex_to_numeric(exponent).to_int(); @@ -255,7 +256,7 @@ int power::degree(symbol const & s) const int power::ldegree(symbol const & s) const { if (is_exactly_of_type(*exponent.bp,numeric)) { - if ((*basis.bp).compare(s)==0) + if ((*basis.bp).compare(s)==0) return ex_to_numeric(exponent).to_int(); else return basis.ldegree(s) * ex_to_numeric(exponent).to_int(); @@ -270,14 +271,14 @@ ex power::coeff(symbol const & s, int const n) const if (n==0) { return *this; } else { - return exZERO(); + return _ex0(); } } else if (is_exactly_of_type(*exponent.bp,numeric)&& (static_cast(*exponent.bp).compare(numeric(n))==0)) { - return exONE(); + return _ex1(); } - return exZERO(); + return _ex0(); } ex power::eval(int level) const @@ -319,10 +320,10 @@ ex power::eval(int level) const // ^(x,0) -> 1 (0^0 also handled here) if (eexponent.is_zero()) - return exONE(); + return _ex1(); // ^(x,1) -> x - if (eexponent.is_equal(exONE())) + if (eexponent.is_equal(_ex1())) return ebasis; // ^(0,x) -> 0 (except if x is real and negative) @@ -330,12 +331,12 @@ ex power::eval(int level) const if (exponent_is_numerical && num_exponent->is_negative()) { throw(std::overflow_error("power::eval(): division by zero")); } else - return exZERO(); + return _ex0(); } // ^(1,x) -> 1 - if (ebasis.is_equal(exONE())) - return exONE(); + if (ebasis.is_equal(_ex1())) + return _ex1(); if (basis_is_numerical && exponent_is_numerical) { // ^(c1,c2) -> c1^c2 (c1, c2 numeric(), @@ -359,14 +360,14 @@ ex power::eval(int level) const q = iquo(n, m, r); if (r.is_negative()) { r = r.add(m); - q = q.sub(numONE()); + q = q.sub(_num1()); } if (q.is_zero()) // the exponent was in the allowed range 0<(n/m)<1 return this->hold(); else { epvector res(2); res.push_back(expair(ebasis,r.div(m))); - res.push_back(expair(ex(num_basis->power(q)),exONE())); + res.push_back(expair(ex(num_basis->power(q)),_ex1())); return (new mul(res))->setflag(status_flags::dynallocated | status_flags::evaluated); /*return mul(num_basis->power(q), power(ex(*num_basis),ex(r.div(m)))).hold(); @@ -405,22 +406,22 @@ ex power::eval(int level) const if (exponent_is_numerical && is_ex_exactly_of_type(ebasis,mul)) { GINAC_ASSERT(!num_exponent->is_integer()); // should have been handled above mul const & mulref=ex_to_mul(ebasis); - if (!mulref.overall_coeff.is_equal(exONE())) { + if (!mulref.overall_coeff.is_equal(_ex1())) { numeric const & num_coeff=ex_to_numeric(mulref.overall_coeff); if (num_coeff.is_real()) { if (num_coeff.is_positive()>0) { mul * mulp=new mul(mulref); - mulp->overall_coeff=exONE(); + mulp->overall_coeff=_ex1(); mulp->clearflag(status_flags::evaluated); mulp->clearflag(status_flags::hash_calculated); return (new mul(power(*mulp,exponent), power(num_coeff,*num_exponent)))-> setflag(status_flags::dynallocated); } else { - GINAC_ASSERT(num_coeff.compare(numZERO())<0); - if (num_coeff.compare(numMINUSONE())!=0) { + GINAC_ASSERT(num_coeff.compare(_num0())<0); + if (num_coeff.compare(_num_1())!=0) { mul * mulp=new mul(mulref); - mulp->overall_coeff=exMINUSONE(); + mulp->overall_coeff=_ex_1(); mulp->clearflag(status_flags::evaluated); mulp->clearflag(status_flags::hash_calculated); return (new mul(power(*mulp,exponent), @@ -644,41 +645,6 @@ ex power::expand_add(add const & a, int const n) const return (new add(sum))->setflag(status_flags::dynallocated); } -/* -ex power::expand_add_2(add const & a) const -{ - // special case: expand a^2 where a is an add - - epvector sum; - sum.reserve((a.seq.size()*(a.seq.size()+1))/2); - epvector::const_iterator last=a.seq.end(); - - for (epvector::const_iterator cit0=a.seq.begin(); cit0!=last; ++cit0) { - ex const & b=a.recombine_pair_to_ex(*cit0); - GINAC_ASSERT(!is_ex_exactly_of_type(b,add)); - GINAC_ASSERT(!is_ex_exactly_of_type(b,power)|| - !is_ex_exactly_of_type(ex_to_power(b).exponent,numeric)|| - !ex_to_numeric(ex_to_power(b).exponent).is_pos_integer()); - if (is_ex_exactly_of_type(b,mul)) { - sum.push_back(a.split_ex_to_pair(expand_mul(ex_to_mul(b),numTWO()))); - } else { - sum.push_back(a.split_ex_to_pair((new power(b,exTWO()))-> - setflag(status_flags::dynallocated))); - } - for (epvector::const_iterator cit1=cit0+1; cit1!=last; ++cit1) { - sum.push_back(a.split_ex_to_pair((new mul(a.recombine_pair_to_ex(*cit0), - a.recombine_pair_to_ex(*cit1)))-> - setflag(status_flags::dynallocated), - exTWO())); - } - } - - GINAC_ASSERT(sum.size()==(a.seq.size()*(a.seq.size()+1))/2); - - return (new add(sum))->setflag(status_flags::dynallocated); -} -*/ - ex power::expand_add_2(add const & a) const { // special case: expand a^2 where a is an add @@ -702,20 +668,20 @@ ex power::expand_add_2(add const & a) const !is_ex_exactly_of_type(ex_to_power(r).basis,mul)|| !is_ex_exactly_of_type(ex_to_power(r).basis,power)); - if (are_ex_trivially_equal(c,exONE())) { + if (are_ex_trivially_equal(c,_ex1())) { if (is_ex_exactly_of_type(r,mul)) { - sum.push_back(expair(expand_mul(ex_to_mul(r),numTWO()),exONE())); + sum.push_back(expair(expand_mul(ex_to_mul(r),_num2()),_ex1())); } else { - sum.push_back(expair((new power(r,exTWO()))->setflag(status_flags::dynallocated), - exONE())); + sum.push_back(expair((new power(r,_ex2()))->setflag(status_flags::dynallocated), + _ex1())); } } else { if (is_ex_exactly_of_type(r,mul)) { - sum.push_back(expair(expand_mul(ex_to_mul(r),numTWO()), - ex_to_numeric(c).power_dyn(numTWO()))); + sum.push_back(expair(expand_mul(ex_to_mul(r),_num2()), + ex_to_numeric(c).power_dyn(_num2()))); } else { - sum.push_back(expair((new power(r,exTWO()))->setflag(status_flags::dynallocated), - ex_to_numeric(c).power_dyn(numTWO()))); + sum.push_back(expair((new power(r,_ex2()))->setflag(status_flags::dynallocated), + ex_to_numeric(c).power_dyn(_num2()))); } } @@ -723,18 +689,18 @@ ex power::expand_add_2(add const & a) const ex const & r1=(*cit1).rest; ex const & c1=(*cit1).coeff; sum.push_back(a.combine_ex_with_coeff_to_pair((new mul(r,r1))->setflag(status_flags::dynallocated), - numTWO().mul(ex_to_numeric(c)).mul_dyn(ex_to_numeric(c1)))); + _num2().mul(ex_to_numeric(c)).mul_dyn(ex_to_numeric(c1)))); } } GINAC_ASSERT(sum.size()==(a.seq.size()*(a.seq.size()+1))/2); // second part: add terms coming from overall_factor (if != 0) - if (!a.overall_coeff.is_equal(exZERO())) { + if (!a.overall_coeff.is_equal(_ex0())) { for (epvector::const_iterator cit=a.seq.begin(); cit!=a.seq.end(); ++cit) { - sum.push_back(a.combine_pair_with_coeff_to_pair(*cit,ex_to_numeric(a.overall_coeff).mul_dyn(numTWO()))); + sum.push_back(a.combine_pair_with_coeff_to_pair(*cit,ex_to_numeric(a.overall_coeff).mul_dyn(_num2()))); } - sum.push_back(expair(ex_to_numeric(a.overall_coeff).power_dyn(numTWO()),exONE())); + sum.push_back(expair(ex_to_numeric(a.overall_coeff).power_dyn(_num2()),_ex1())); } GINAC_ASSERT(sum.size()==(a_nops*(a_nops+1))/2); @@ -746,8 +712,8 @@ ex power::expand_mul(mul const & m, numeric const & n) const { // expand m^n where m is a mul and n is and integer - if (n.is_equal(numZERO())) { - return exONE(); + if (n.is_equal(_num0())) { + return _ex1(); } epvector distrseq; @@ -801,7 +767,7 @@ ex power::expand_commutative_3(ex const & basis, numeric const & exponent, ex power::expand_noncommutative(ex const & basis, numeric const & exponent, unsigned options) const { - ex rest_power=ex(power(basis,exponent.add(numMINUSONE()))). + ex rest_power=ex(power(basis,exponent.add(_num_1()))). expand(options | expand_options::internal_do_not_expand_power_operands); return ex(mul(rest_power,basis),0). @@ -824,6 +790,13 @@ unsigned power::precedence=60; const power some_power; type_info const & typeid_power=typeid(some_power); +// helper function + +ex sqrt(ex const & a) +{ + return power(a,_ex1_2()); +} + #ifndef NO_GINAC_NAMESPACE } // namespace GiNaC #endif // ndef NO_GINAC_NAMESPACE diff --git a/ginac/power.h b/ginac/power.h index b1136c35..17c06c77 100644 --- a/ginac/power.h +++ b/ginac/power.h @@ -124,8 +124,7 @@ inline ex pow(ex const & b, ex const & e) /** Square root expression. Returns a power-object with exponent 1/2 as a new * expression. */ -inline ex sqrt(ex const & a) -{ return power(a,exHALF()); } +ex sqrt(ex const & a); #ifndef NO_GINAC_NAMESPACE } // namespace GiNaC diff --git a/ginac/relational.cpp b/ginac/relational.cpp index 19a2e3a2..614c1629 100644 --- a/ginac/relational.cpp +++ b/ginac/relational.cpp @@ -25,6 +25,7 @@ #include "relational.h" #include "numeric.h" #include "debugmsg.h" +#include "utils.h" #ifndef NO_GINAC_NAMESPACE namespace GiNaC { @@ -330,7 +331,7 @@ relational::operator bool() const if (!is_ex_exactly_of_type(df,numeric)) { return o==not_equal ? true : false; // cannot decide on non-numerical results } - int cmpval=ex_to_numeric(df).compare(numZERO()); + int cmpval=ex_to_numeric(df).compare(_num0()); switch (o) { case equal: return cmpval==0; diff --git a/ginac/series.cpp b/ginac/series.cpp index 0cc2be34..4189b098 100644 --- a/ginac/series.cpp +++ b/ginac/series.cpp @@ -29,6 +29,7 @@ #include "relational.h" #include "symbol.h" #include "debugmsg.h" +#include "utils.h" #ifndef NO_GINAC_NAMESPACE namespace GiNaC { @@ -87,7 +88,7 @@ void series::destroy(bool call_parent) /** Construct series from a vector of coefficients and powers. * expair.rest holds the coefficient, expair.coeff holds the power. * The powers must be integers (positive or negative) and in ascending order; - * the last coefficient can be Order(exONE()) to represent a truncated, + * the last coefficient can be Order(_ex1()) to represent a truncated, * non-terminating series. * * @param var_ series variable (must hold a symbol) @@ -177,7 +178,7 @@ int series::ldegree(symbol const &s) const } // Coefficient of variable -ex series::coeff(symbol const &s, int n) const +ex series::coeff(symbol const &s, int const n) const { if (var.is_equal(s)) { epvector::const_iterator it = seq.begin(), itend = seq.end(); @@ -186,10 +187,10 @@ ex series::coeff(symbol const &s, int n) const if (pow == n) return it->rest; if (pow > n) - return exZERO(); + return _ex0(); it++; } - return exZERO(); + return _ex0(); } else return convert_to_poly().coeff(s, n); } @@ -271,7 +272,7 @@ ex basic::series(symbol const & s, ex const & point, int order) const // Higher-order terms, if present deriv = deriv.diff(s); if (!deriv.is_zero()) - seq.push_back(expair(Order(exONE()), numeric(n))); + seq.push_back(expair(Order(_ex1()), numeric(n))); return series::series(s, point, seq); } @@ -287,7 +288,7 @@ ex series::add_series(const series &other) const // results in an empty (constant) series if (!is_compatible_to(other)) { epvector nul; - nul.push_back(expair(Order(exONE()), exZERO())); + nul.push_back(expair(Order(_ex1()), _ex0())); return series(var, point, nul); } @@ -335,7 +336,7 @@ ex series::add_series(const series &other) const } else { // Add coefficient of a and b if (is_order_function((*a).rest) || is_order_function((*b).rest)) { - new_seq.push_back(expair(Order(exONE()), (*a).coeff)); + new_seq.push_back(expair(Order(_ex1()), (*a).coeff)); break; // Order term ends the sequence } else { ex sum = (*a).rest + (*b).rest; @@ -366,7 +367,7 @@ ex add::series(symbol const & s, ex const & point, int order) const acc = it->rest; else acc = it->rest.series(s, point, order); - if (!it->coeff.is_equal(exONE())) + if (!it->coeff.is_equal(_ex1())) acc = ex_to_series(acc).mul_const(ex_to_numeric(it->coeff)); it++; } @@ -378,7 +379,7 @@ ex add::series(symbol const & s, ex const & point, int order) const op = it->rest; else op = it->rest.series(s, point, order); - if (!it->coeff.is_equal(exONE())) + if (!it->coeff.is_equal(_ex1())) op = ex_to_series(op).mul_const(ex_to_numeric(it->coeff)); // Series addition @@ -403,7 +404,7 @@ ex add::series(symbol const & s, ex const & point, int order) const op = it->rest; else op = it->rest.series(s, point, order); - if (!it->coeff.is_equal(exONE())) + if (!it->coeff.is_equal(_ex1())) op = ex_to_series(op).mul_const(ex_to_numeric(it->coeff)); // Series addition @@ -446,7 +447,7 @@ ex series::mul_series(const series &other) const // results in an empty (constant) series if (!is_compatible_to(other)) { epvector nul; - nul.push_back(expair(Order(exONE()), exZERO())); + nul.push_back(expair(Order(_ex1()), _ex0())); return series(var, point, nul); } @@ -472,7 +473,7 @@ ex series::mul_series(const series &other) const cdeg_max = higher_order_c - 1; for (int cdeg=cdeg_min; cdeg<=cdeg_max; cdeg++) { - ex co = exZERO(); + ex co = _ex0(); // c(i)=a(0)b(i)+...+a(i)b(0) for (int i=a_min; cdeg-i>=b_min; i++) { ex a_coeff = coeff(*s, i); @@ -484,7 +485,7 @@ ex series::mul_series(const series &other) const new_seq.push_back(expair(co, numeric(cdeg))); } if (higher_order_c < INT_MAX) - new_seq.push_back(expair(Order(exONE()), numeric(higher_order_c))); + new_seq.push_back(expair(Order(_ex1()), numeric(higher_order_c))); return series::series(var, point, new_seq); } @@ -502,7 +503,7 @@ ex mul::series(symbol const & s, ex const & point, int order) const acc = it->rest; else acc = it->rest.series(s, point, order); - if (!it->coeff.is_equal(exONE())) + if (!it->coeff.is_equal(_ex1())) acc = ex_to_series(acc).power_const(ex_to_numeric(it->coeff), order); it++; } @@ -517,7 +518,7 @@ ex mul::series(symbol const & s, ex const & point, int order) const continue; } else if (!is_ex_exactly_of_type(op, series)) op = op.series(s, point, order); - if (!it->coeff.is_equal(exONE())) + if (!it->coeff.is_equal(_ex1())) op = ex_to_series(op).power_const(ex_to_numeric(it->coeff), order); // Series multiplication @@ -549,7 +550,7 @@ ex mul::series(symbol const & s, ex const & point, int order) const continue; } else if (!is_ex_exactly_of_type(op, series)) op = op.series(s, point, order); - if (!it->coeff.is_equal(exONE())) + if (!it->coeff.is_equal(_ex1())) op = ex_to_series(op).power_const(ex_to_numeric(it->coeff), order); // Series multiplication @@ -576,11 +577,11 @@ ex series::power_const(const numeric &p, int deg) const co.push_back(co0 = power(coeff(*s, ldeg), p)); bool all_sums_zero = true; for (i=1; iis_assigned) { asexinfop->is_assigned=0; - asexinfop->assigned_expression=exZERO(); + asexinfop->assigned_expression=_ex0(); } setflag(status_flags::evaluated); } diff --git a/ginac/utils.cpp b/ginac/utils.cpp index 7592e403..e113c04b 100644 --- a/ginac/utils.cpp +++ b/ginac/utils.cpp @@ -1,6 +1,7 @@ /** @file utils.cpp * - * Implementation of several small and furry utilities. */ + * Implementation of several small and furry utilities needed within GiNaC + * but not of any interest to the user of the library. */ /* * GiNaC Copyright (C) 1999 Johannes Gutenberg University Mainz, Germany @@ -20,6 +21,8 @@ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ +#include "ex.h" +#include "numeric.h" #include "utils.h" #ifndef NO_GINAC_NAMESPACE @@ -46,6 +49,415 @@ int compare_pointers(void const * a, void const * b) return 0; } +////////// +// `construct on first use' chest of numbers +////////// + +// numeric -12 +numeric const & _num_12(void) +{ + const static ex e = ex(numeric(-12)); + const static numeric * n = static_cast(e.bp); + return *n; +} + +ex const & _ex_12(void) +{ + static ex * e = new ex(_num_12()); + return *e; +} + +// numeric -11 +numeric const & _num_11(void) +{ + const static ex e = ex(numeric(-11)); + const static numeric * n = static_cast(e.bp); + return *n; +} + +ex const & _ex_11(void) +{ + static ex * e = new ex(_num_11()); + return *e; +} + +// numeric -10 +numeric const & _num_10(void) +{ + const static ex e = ex(numeric(-10)); + const static numeric * n = static_cast(e.bp); + return *n; +} + +ex const & _ex_10(void) +{ + static ex * e = new ex(_num_10()); + return *e; +} + +// numeric -9 +numeric const & _num_9(void) +{ + const static ex e = ex(numeric(-9)); + const static numeric * n = static_cast(e.bp); + return *n; +} + +ex const & _ex_9(void) +{ + static ex * e = new ex(_num_9()); + return *e; +} + +// numeric -8 +numeric const & _num_8(void) +{ + const static ex e = ex(numeric(-8)); + const static numeric * n = static_cast(e.bp); + return *n; +} + +ex const & _ex_8(void) +{ + static ex * e = new ex(_num_8()); + return *e; +} + +// numeric -7 +numeric const & _num_7(void) +{ + const static ex e = ex(numeric(-7)); + const static numeric * n = static_cast(e.bp); + return *n; +} + +ex const & _ex_7(void) +{ + static ex * e = new ex(_num_7()); + return *e; +} + +// numeric -6 +numeric const & _num_6(void) +{ + const static ex e = ex(numeric(-6)); + const static numeric * n = static_cast(e.bp); + return *n; +} + +ex const & _ex_6(void) +{ + static ex * e = new ex(_num_6()); + return *e; +} + +// numeric -5 +numeric const & _num_5(void) +{ + const static ex e = ex(numeric(-5)); + const static numeric * n = static_cast(e.bp); + return *n; +} + +ex const & _ex_5(void) +{ + static ex * e = new ex(_num_5()); + return *e; +} + +// numeric -4 +numeric const & _num_4(void) +{ + const static ex e = ex(numeric(-4)); + const static numeric * n = static_cast(e.bp); + return *n; +} + +ex const & _ex_4(void) +{ + static ex * e = new ex(_num_4()); + return *e; +} + +// numeric -3 +numeric const & _num_3(void) +{ + const static ex e = ex(numeric(-3)); + const static numeric * n = static_cast(e.bp); + return *n; +} + +ex const & _ex_3(void) +{ + static ex * e = new ex(_num_3()); + return *e; +} + +// numeric -2 +numeric const & _num_2(void) +{ + const static ex e = ex(numeric(-2)); + const static numeric * n = static_cast(e.bp); + return *n; +} + +ex const & _ex_2(void) +{ + static ex * e = new ex(_num_2()); + return *e; +} + +// numeric -1 +numeric const & _num_1(void) +{ + const static ex e = ex(numeric(-1)); + const static numeric * n = static_cast(e.bp); + return *n; +} + +ex const & _ex_1(void) +{ + static ex * e = new ex(_num_1()); + return *e; +} + +// numeric -1/2 +numeric const & _num_1_2(void) +{ + const static ex e = ex(numeric(-1,2)); + const static numeric * n = static_cast(e.bp); + return *n; +} + +ex const & _ex_1_2(void) +{ + static ex * e = new ex(_num_1_2()); + return *e; +} + +// numeric -1/3 +numeric const & _num_1_3(void) +{ + const static ex e = ex(numeric(-1,3)); + const static numeric * n = static_cast(e.bp); + return *n; +} + +ex const & _ex_1_3(void) +{ + static ex * e = new ex(_num_1_3()); + return *e; +} + +// numeric 0 +numeric const & _num0(void) +{ + const static ex e = ex(numeric(0)); + const static numeric * n = static_cast(e.bp); + return *n; +} + +ex const & _ex0(void) +{ + static ex * e = new ex(_num0()); + return *e; +} + +// numeric 1/3 +numeric const & _num1_3(void) +{ + const static ex e = ex(numeric(1,3)); + const static numeric * n = static_cast(e.bp); + return *n; +} + +ex const & _ex1_3(void) +{ + static ex * e = new ex(_num1_3()); + return *e; +} + +// numeric 1/2 +numeric const & _num1_2(void) +{ + const static ex e = ex(numeric(1,2)); + const static numeric * n = static_cast(e.bp); + return *n; +} + +ex const & _ex1_2(void) +{ + static ex * e = new ex(_num1_2()); + return *e; +} + +// numeric 1 +numeric const & _num1(void) +{ + const static ex e = ex(numeric(1)); + const static numeric * n = static_cast(e.bp); + return *n; +} + +ex const & _ex1(void) +{ + static ex * e = new ex(_num1()); + return *e; +} + +// numeric 2 +numeric const & _num2(void) +{ + const static ex e = ex(numeric(2)); + const static numeric * n = static_cast(e.bp); + return *n; +} + +ex const & _ex2(void) +{ + static ex * e = new ex(_num2()); + return *e; +} + +// numeric 3 +numeric const & _num3(void) +{ + const static ex e = ex(numeric(3)); + const static numeric * n = static_cast(e.bp); + return *n; +} + +ex const & _ex3(void) +{ + static ex * e = new ex(_num3()); + return *e; +} + +// numeric 4 +numeric const & _num4(void) +{ + const static ex e = ex(numeric(4)); + const static numeric * n = static_cast(e.bp); + return *n; +} + +ex const & _ex4(void) +{ + static ex * e = new ex(_num4()); + return *e; +} + +// numeric 5 +numeric const & _num5(void) +{ + const static ex e = ex(numeric(5)); + const static numeric * n = static_cast(e.bp); + return *n; +} + +ex const & _ex5(void) +{ + static ex * e = new ex(_num5()); + return *e; +} + +// numeric 6 +numeric const & _num6(void) +{ + const static ex e = ex(numeric(6)); + const static numeric * n = static_cast(e.bp); + return *n; +} + +ex const & _ex6(void) +{ + static ex * e = new ex(_num6()); + return *e; +} + +// numeric 7 +numeric const & _num7(void) +{ + const static ex e = ex(numeric(7)); + const static numeric * n = static_cast(e.bp); + return *n; +} + +ex const & _ex7(void) +{ + static ex * e = new ex(_num7()); + return *e; +} + +// numeric 8 +numeric const & _num8(void) +{ + const static ex e = ex(numeric(8)); + const static numeric * n = static_cast(e.bp); + return *n; +} + +ex const & _ex8(void) +{ + static ex * e = new ex(_num8()); + return *e; +} + +// numeric 9 +numeric const & _num9(void) +{ + const static ex e = ex(numeric(9)); + const static numeric * n = static_cast(e.bp); + return *n; +} + +ex const & _ex9(void) +{ + static ex * e = new ex(_num9()); + return *e; +} + +// numeric 10 +numeric const & _num10(void) +{ + const static ex e = ex(numeric(10)); + const static numeric * n = static_cast(e.bp); + return *n; +} + +ex const & _ex10(void) +{ + static ex * e = new ex(_num10()); + return *e; +} + +// numeric 11 +numeric const & _num11(void) +{ + const static ex e = ex(numeric(11)); + const static numeric * n = static_cast(e.bp); + return *n; +} + +ex const & _ex11(void) +{ + static ex * e = new ex(_num11()); + return *e; +} + +// numeric 12 +numeric const & _num12(void) +{ + const static ex e = ex(numeric(12)); + const static numeric * n = static_cast(e.bp); + return *n; +} + +ex const & _ex12(void) +{ + static ex * e = new ex(_num12()); + return *e; +} // comment skeleton for header files @@ -121,6 +533,7 @@ int compare_pointers(void const * a, void const * b) // private // none + #ifndef NO_GINAC_NAMESPACE } // namespace GiNaC #endif // ndef NO_GINAC_NAMESPACE diff --git a/ginac/utils.h b/ginac/utils.h index 7be44f73..0f6082a6 100644 --- a/ginac/utils.h +++ b/ginac/utils.h @@ -1,7 +1,7 @@ /** @file utils.h * * Interface to several small and furry utilities needed within GiNaC but not - * of interest to the user of the library. */ + * of any interest to the user of the library. */ /* * GiNaC Copyright (C) 1999 Johannes Gutenberg University Mainz, Germany @@ -138,6 +138,74 @@ OutputIterator mymerge3(InputIterator1 first1, InputIterator1 last1, } } +// Collection of `construct on first use' wrappers for safely avoiding +// internal object replication without running into the `static +// initialization order fiasco'. This chest of numbers helps speed up +// the library but should not be used outside it since it is +// potentially confusing. + +class numeric; +class ex; + +numeric const & _num_12(void); // -12 +numeric const & _num_11(void); // -11 +numeric const & _num_10(void); // -10 +numeric const & _num_9(void); // -9 +numeric const & _num_8(void); // -8 +numeric const & _num_7(void); // -7 +numeric const & _num_6(void); // -6 +numeric const & _num_5(void); // -5 +numeric const & _num_4(void); // -4 +numeric const & _num_3(void); // -3 +numeric const & _num_2(void); // -2 +numeric const & _num_1(void); // -1 +numeric const & _num_1_2(void); // -1/2 +numeric const & _num_1_3(void); // -1/3 +numeric const & _num0(void); // 0 +numeric const & _num1_3(void); // 1/3 +numeric const & _num1_2(void); // 1/2 +numeric const & _num1(void); // 1 +numeric const & _num2(void); // 2 +numeric const & _num3(void); // 3 +numeric const & _num4(void); // 4 +numeric const & _num5(void); // 5 +numeric const & _num6(void); // 6 +numeric const & _num7(void); // 7 +numeric const & _num8(void); // 8 +numeric const & _num9(void); // 9 +numeric const & _num10(void); // 10 +numeric const & _num11(void); // 11 +numeric const & _num12(void); // 12 +ex const & _ex_12(void); // -12 +ex const & _ex_11(void); // -11 +ex const & _ex_10(void); // -10 +ex const & _ex_9(void); // -9 +ex const & _ex_8(void); // -8 +ex const & _ex_7(void); // -7 +ex const & _ex_6(void); // -6 +ex const & _ex_5(void); // -5 +ex const & _ex_4(void); // -4 +ex const & _ex_3(void); // -3 +ex const & _ex_2(void); // -2 +ex const & _ex_1(void); // -1 +ex const & _ex_1_2(void); // -1/2 +ex const & _ex_1_3(void); // -1/3 +ex const & _ex0(void); // 0 +ex const & _ex1_3(void); // 1/3 +ex const & _ex1_2(void); // 1/2 +ex const & _ex1(void); // 1 +ex const & _ex2(void); // 2 +ex const & _ex3(void); // 3 +ex const & _ex4(void); // 4 +ex const & _ex5(void); // 5 +ex const & _ex6(void); // 6 +ex const & _ex7(void); // 7 +ex const & _ex8(void); // 8 +ex const & _ex9(void); // 9 +ex const & _ex10(void); // 10 +ex const & _ex11(void); // 11 +ex const & _ex12(void); // 12 + #ifndef NO_GINAC_NAMESPACE } // namespace GiNaC #endif // ndef NO_GINAC_NAMESPACE diff --git a/ginsh/ginsh_parser.yy b/ginsh/ginsh_parser.yy index 2e02b6cc..0134ab15 100644 --- a/ginsh/ginsh_parser.yy +++ b/ginsh/ginsh_parser.yy @@ -341,7 +341,7 @@ static ex f_evalf2(const exprseq &e) static ex f_has(const exprseq &e) { - return e[0].has(e[1]) ? exONE() : exZERO(); + return e[0].has(e[1]) ? ex(1) : ex(0); } static ex f_inverse(const exprseq &e) @@ -353,7 +353,7 @@ static ex f_inverse(const exprseq &e) static ex f_is(const exprseq &e) { CHECK_ARG(0, relational, is); - return (bool)ex_to_relational(e[0]) ? exONE() : exZERO(); + return (bool)ex_to_relational(e[0]) ? ex(1) : ex(0); } static ex f_lcoeff(const exprseq &e) @@ -410,7 +410,7 @@ static ex f_rem(const exprseq &e) static ex f_series2(const exprseq &e) { CHECK_ARG(1, symbol, series); - return e[0].series(ex_to_symbol(e[1]), exZERO()); + return e[0].series(ex_to_symbol(e[1]), ex(0)); } static ex f_series3(const exprseq &e) @@ -677,7 +677,7 @@ static ex lst2matrix(const ex &l) if (l.op(i).nops() > j) m.set(i, j, l.op(i).op(j)); else - m.set(i, j, exZERO()); + m.set(i, j, ex(0)); return m; } -- 2.44.0