- Banned exZERO(), exONE(), exMINUSHALF() and all this from the interface.
authorRichard Kreckel <Richard.Kreckel@uni-mainz.de>
Fri, 17 Dec 1999 19:58:25 +0000 (19:58 +0000)
committerRichard Kreckel <Richard.Kreckel@uni-mainz.de>
Fri, 17 Dec 1999 19:58:25 +0000 (19:58 +0000)
  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

36 files changed:
check/differentiation.cpp
check/normalization.cpp
check/numeric_consist.cpp
check/paranoia_check.cpp
check/poly_gcd.cpp
check/series_expansion.cpp
ginac/add.cpp
ginac/basic.cpp
ginac/color.cpp
ginac/diff.cpp
ginac/ex.cpp
ginac/ex.h
ginac/expair.h
ginac/expairseq.cpp
ginac/inifcns.cpp
ginac/inifcns.h
ginac/inifcns_gamma.cpp
ginac/inifcns_trans.cpp
ginac/inifcns_zeta.cpp
ginac/lortensor.cpp
ginac/matrix.cpp
ginac/mul.cpp
ginac/ncmul.cpp
ginac/normal.cpp
ginac/numeric.cpp
ginac/numeric.h
ginac/operators.cpp
ginac/power.cpp
ginac/power.h
ginac/relational.cpp
ginac/series.cpp
ginac/simp_lor.cpp
ginac/symbol.cpp
ginac/utils.cpp
ginac/utils.h
ginsh/ginsh_parser.yy

index 348b5fd..ece8bc3 100644 (file)
@@ -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<series *>(ed.bp)->convert_to_poly();
     d = static_cast<series *>(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;
index 7ce1e5d..78ae631 100644 (file)
@@ -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
index d7edb3b..494cc57 100644 (file)
@@ -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;
     }
index 2e0c40d..f337060 100644 (file)
@@ -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;
index 1448ae3..12de267 100644 (file)
@@ -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;
                        }
index b67705d..b8e8020 100644 (file)
@@ -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<series *>(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;
 }
index 804f07a..ad0ac1f 100644 (file)
@@ -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<mul *>(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<mul *>(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;
index acacd23..7b6a6fb 100644 (file)
@@ -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
index 696d70f..0f5f604 100644 (file)
@@ -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; k<j; ++k) {
                             S.push_back(Tvecs[rl][k]);
                         }
-                        t1=exONE();
-                        t2=exONE();
+                        t1=_ex1();
+                        t2=_ex1();
                         ex term1=numeric(-1)/numeric(6)*nonsimplified_ncmul(recombine_color_string(
                                  delta8vec,fvec,dvec,Tvecs,ONEvecs,unknownvec));
                         for (unsigned k=i+1; k<j; ++k) {
-                            S.push_back(exONE());
+                            S.push_back(_ex1());
                         }
                         t1=color_trace_of_one_representation_label(S);
                         ex term2=numeric(1)/numeric(2)*nonsimplified_ncmul(recombine_color_string(
@@ -875,7 +876,7 @@ ex simplify_color(ex const & e)
 
     // simplification of sum=sum of simplifications
     if (is_ex_exactly_of_type(e_expanded,add)) {
-        ex sum=exZERO();
+        ex sum=_ex0();
         for (int i=0; i<e_expanded.nops(); ++i) {
             sum += simplify_color(e_expanded.op(i));
         }
@@ -884,7 +885,7 @@ ex simplify_color(ex const & e)
 
     // simplification of commutative product=commutative product of simplifications
     if (is_ex_exactly_of_type(e_expanded,mul)) {
-        ex prod=exONE();
+        ex prod=_ex1();
         for (int i=0; i<e_expanded.nops(); ++i) {
             prod *= simplify_color(e_expanded.op(i));
         }
@@ -935,7 +936,7 @@ ex brute_force_sum_color_indices(ex const & e)
         counter[l]=1;
     }
 
-    ex sum=exZERO();
+    ex sum=_ex0();
     
     while (1) {
         ex term=e;
index 65ee609..ae85b55 100644 (file)
@@ -36,6 +36,7 @@
 #include "relational.h"
 #include "series.h"
 #include "symbol.h"
+#include "utils.h"
 
 #ifndef NO_GINAC_NAMESPACE
 namespace GiNaC {
@@ -54,7 +55,7 @@ ex basic::diff(symbol const & s) const
  *  @see ex::diff */
 ex numeric::diff(symbol const & s) const
 {
-    return exZERO();
+    return _ex0();
 }
 
 
@@ -65,9 +66,9 @@ ex numeric::diff(symbol const & s) const
 ex symbol::diff(symbol const & s) const
 {
     if (compare_same_type(s)) {
-        return exZERO();
+        return _ex0();
     } else {
-        return exONE();
+        return _ex1();
     }
 }
 
@@ -76,7 +77,7 @@ ex symbol::diff(symbol const & s) const
  *  @see ex::diff */
 ex constant::diff(symbol const & s) const
 {
-    return exZERO();
+    return _ex0();
 }
 
 /** Implementation of ex::diff() for multiple differentiation of a symbol.
@@ -92,13 +93,13 @@ ex symbol::diff(symbol const & s, unsigned nth) const
             return s;
             break;
         case 1:
-            return exONE();
+            return _ex1();
             break;
         default:
-            return exZERO();
+            return _ex0();
         }
     } else {
-        return exONE();
+        return _ex1();
     }
 }
 
@@ -107,7 +108,7 @@ ex symbol::diff(symbol const & s, unsigned nth) const
  *  @see ex::diff */
 ex indexed::diff(symbol const & s) const
 {
-        return exZERO();
+        return _ex0();
 }
 
 
@@ -151,7 +152,7 @@ ex mul::diff(symbol const & s) const
  *  @see ex::diff */
 ex ncmul::diff(symbol const & s) const
 {
-    return exZERO();
+    return _ex0();
 }
 
 
@@ -161,7 +162,7 @@ ex power::diff(symbol const & s) const
 {
     if (exponent.info(info_flags::real)) {
         // D(b^r) = r * b^(r-1) * D(b) (faster than the formula below)
-        return mul(mul(exponent, power(basis, exponent - exONE())), basis.diff(s));
+        return mul(mul(exponent, power(basis, exponent - _ex1())), basis.diff(s));
     } else {
         // D(b^e) = b^e * (D(e)*ln(b) + e*D(b)/b)
         return mul(power(basis, exponent),
index cdef625..6a7a0aa 100644 (file)
@@ -29,6 +29,7 @@
 #include "numeric.h"
 #include "power.h"
 #include "debugmsg.h"
+#include "utils.h"
 
 #ifndef NO_GINAC_NAMESPACE
 namespace GiNaC {
@@ -42,11 +43,11 @@ namespace GiNaC {
 
 #ifndef INLINE_EX_CONSTRUCTORS
 
-ex::ex() : bp(exZERO().bp)
+ex::ex() : bp(ex0().bp)
 {
     debugmsg("ex default constructor",LOGLEVEL_CONSTRUCT);
-    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;
 }
@@ -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<n.nops(); i++) {
         if (!is_ex_exactly_of_type(n.op(i), power) || !n.op(i).op(1).info(info_flags::negint))
             res *= n.op(i);
@@ -336,7 +337,7 @@ ex ex::denom(bool normalize) const
 
     // polynomial
     if (n.info(info_flags::polynomial))
-        return exONE();
+        return _ex1();
 
     // something^(-int)
     if (is_ex_exactly_of_type(n, power) && n.op(1).info(info_flags::negint))
@@ -344,8 +345,8 @@ ex ex::denom(bool normalize) const
 
     // something^(int) * something^(int) * ...
     if (!is_ex_exactly_of_type(n, mul))
-        return exONE();
-    ex res = exONE();
+        return _ex1();
+    ex res = _ex1();
     for (int i=0; i<n.nops(); i++) {
         if (is_ex_exactly_of_type(n.op(i), power) && n.op(i).op(1).info(info_flags::negint))
             res *= power(n.op(i), -1);
@@ -543,47 +544,8 @@ void ex::construct_from_basic(basic const & other)
 // global functions
 //////////
 
-ex const & exZERO(void)
-{
-    static ex * eZERO=new ex(numZERO());
-    return *eZERO;
-}
-
-ex const & exONE(void)
-{
-    static ex * eONE=new ex(numONE());
-    return *eONE;
-}
-
-ex const & exTWO(void)
-{
-    static ex * eTWO=new ex(numTWO());
-    return *eTWO;
-}
-
-ex const & exTHREE(void)
-{
-    static ex * eTHREE=new ex(numTHREE());
-    return *eTHREE;
-}
-
-ex const & exMINUSONE(void)
-{
-    static ex * eMINUSONE=new ex(numMINUSONE());
-    return *eMINUSONE;
-}
-
-ex const & exHALF(void)
-{
-    static ex * eHALF=new ex(ex(1)/ex(2));
-    return *eHALF;
-}
+// none
 
-ex const & exMINUSHALF(void)
-{
-    static ex * eMINUSHALF=new ex(numeric(-1,2));
-    return *eMINUSHALF;
-}
 
 #ifndef NO_GINAC_NAMESPACE
 } // namespace GiNaC
index f2c0a03..f8d2970 100644 (file)
@@ -38,17 +38,9 @@ class status_flags;
 class symbol;
 class lst;
 
-// typedef vector<ex> 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<ex> 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;
index b130465..8580598 100644 (file)
@@ -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;
             }
index 209c725..1fa592e 100644 (file)
@@ -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;
         }
     }
index 8b7f327..0a38c4b 100644 (file)
@@ -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;
index d7642ba..fdd3310 100644 (file)
@@ -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);
index ecec37a..3f7fc67 100644 (file)
@@ -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
index 8838360..4473fe0 100644 (file)
@@ -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);
index 6829d56..9fdce1e 100644 (file)
@@ -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);
index 7cbd2f4..1854af0 100644 (file)
 #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; i<n; ++i) {
         ex f=m.op(i);
-        if (is_ex_exactly_of_type(f,power)&&f.op(1).is_equal(exTWO())) {
+        if (is_ex_exactly_of_type(f,power)&&f.op(1).is_equal(_ex2())) {
             v_contracted.push_back(f.op(0));
             v_contracted.push_back(f.op(0));
         } else {
             v_contracted.push_back(f);
-       }
+        }
     }
 
     unsigned replacements;
@@ -426,7 +426,7 @@ ex simplify_lortensor_mul(ex const & m)
                 } else {
                     // a contracted index should occur exactly once
                     GINAC_ASSERT(replacements==1);
-                    *it=exONE();
+                    *it=_ex1();
                     something_changed=true;
                 }
             }
@@ -442,7 +442,7 @@ ex simplify_lortensor_mul(ex const & m)
                 } else {
                     // a contracted index should occur exactly once
                     GINAC_ASSERT(replacements==1);
-                    *it=exONE();
+                    *it=_ex1();
                     something_changed=true;
                 }
             }
@@ -462,7 +462,7 @@ ex simplify_lortensor(ex const & e)
 
     // simplification of sum=sum of simplifications
     if (is_ex_exactly_of_type(e_expanded,add)) {
-        ex sum=exZERO();
+        ex sum=_ex0();
         for (int i=0; i<e_expanded.nops(); ++i) {
             sum += simplify_lortensor(e_expanded.op(i));
         }
index dd0caa5..eddc620 100644 (file)
@@ -25,6 +25,7 @@
 
 #include "matrix.h"
 #include "debugmsg.h"
+#include "utils.h"
 
 #ifndef NO_GINAC_NAMESPACE
 namespace GiNaC {
@@ -42,7 +43,7 @@ matrix::matrix()
     : basic(TINFO_matrix), row(1), col(1)
 {
     debugmsg("matrix default constructor",LOGLEVEL_CONSTRUCT);
-    m.push_back(exZERO());
+    m.push_back(_ex0());
 }
 
 matrix::~matrix()
@@ -95,7 +96,7 @@ matrix::matrix(int r, int c)
     : basic(TINFO_matrix), row(r), col(c)
 {
     debugmsg("matrix constructor from int,int",LOGLEVEL_CONSTRUCT);
-    m.resize(r*c, exZERO());
+    m.resize(r*c, _ex0());
 }
 
 // protected
@@ -390,16 +391,16 @@ ex determinant_numeric(const matrix & M)
 {
     GINAC_ASSERT(M.rows()==M.cols());  // cannot happen, just in case...
     matrix tmp(M);
-    ex det=exONE();
+    ex det=_ex1();
     ex piv;
     
     for (int r1=0; r1<M.rows(); ++r1) {
         int indx = tmp.pivot(r1);
         if (indx == -1) {
-            return exZERO();
+            return _ex0();
         }
         if (indx != 0) {
-            det *= exMINUSONE();
+            det *= _ex_1();
         }
         det = det * tmp.m[r1*M.cols()+r1];
         for (int r2=r1+1; r2<M.rows(); ++r2) {
@@ -609,7 +610,7 @@ matrix matrix::inverse(void) const
     matrix tmp(row,col);
     // set tmp to the unit matrix
     for (int i=0; i<col; ++i) {
-        tmp.m[i*col+i] = exONE();
+        tmp.m[i*col+i] = _ex1();
     }
     // create a copy of this matrix
     matrix cpy(*this);
@@ -689,7 +690,7 @@ matrix matrix::fraction_free_elim(matrix const & vars,
     for (int k=1; (k<=n)&&(r<=m); ++k) {
         // find a nonzero pivot
         int p;
-        for (p=r; (p<=m)&&(a.ffe_get(p,k).is_equal(exZERO())); ++p) {}
+        for (p=r; (p<=m)&&(a.ffe_get(p,k).is_equal(_ex0())); ++p) {}
         // pivot is in row p
         if (p<=m) {
             if (p!=r) {
@@ -734,7 +735,7 @@ matrix matrix::fraction_free_elim(matrix const & vars,
     for (int r=1; r<=m; ++r) {
         int zero_in_this_row=0;
         for (int c=1; c<=n; ++c) {
-            if (a.ffe_get(r,c).is_equal(exZERO())) {
+            if (a.ffe_get(r,c).is_equal(_ex0())) {
                zero_in_this_row++;
             } else {
                 break;
index e8595b6..9bee46f 100644 (file)
@@ -27,6 +27,7 @@
 #include "add.h"
 #include "power.h"
 #include "debugmsg.h"
+#include "utils.h"
 
 #ifndef NO_GINAC_NAMESPACE
 namespace GiNaC {
@@ -88,7 +89,7 @@ mul::mul(ex const & lh, ex const & rh)
 {
     debugmsg("mul constructor from ex,ex",LOGLEVEL_CONSTRUCT);
     tinfo_key = TINFO_mul;
-    overall_coeff=exONE();
+    overall_coeff=_ex1();
     construct_from_2_ex(lh,rh);
     GINAC_ASSERT(is_canonical());
 }
@@ -97,7 +98,7 @@ mul::mul(exvector const & v)
 {
     debugmsg("mul constructor from exvector",LOGLEVEL_CONSTRUCT);
     tinfo_key = TINFO_mul;
-    overall_coeff=exONE();
+    overall_coeff=_ex1();
     construct_from_exvector(v);
     GINAC_ASSERT(is_canonical());
 }
@@ -123,7 +124,7 @@ mul::mul(epvector const & v)
 {
     debugmsg("mul constructor from epvector",LOGLEVEL_CONSTRUCT);
     tinfo_key = TINFO_mul;
-    overall_coeff=exONE();
+    overall_coeff=_ex1();
     construct_from_epvector(v);
     GINAC_ASSERT(is_canonical());
 }
@@ -157,7 +158,7 @@ mul::mul(ex const & lh, ex const & mh, ex const & rh)
     factors.push_back(lh);
     factors.push_back(mh);
     factors.push_back(rh);
-    overall_coeff=exONE();
+    overall_coeff=_ex1();
     construct_from_exvector(factors);
     GINAC_ASSERT(is_canonical());
 }
@@ -181,10 +182,10 @@ void mul::print(ostream & os, unsigned upper_precedence) const
     bool first=true;
     // first print the overall numeric coefficient:
     if (ex_to_numeric(overall_coeff).csgn()==-1) os << '-';
-    if (!overall_coeff.is_equal(exONE()) &&
-        !overall_coeff.is_equal(exMINUSONE())) {
+    if (!overall_coeff.is_equal(_ex1()) &&
+        !overall_coeff.is_equal(_ex_1())) {
         if (ex_to_numeric(overall_coeff).csgn()==-1)
-            (numMINUSONE()*overall_coeff).print(os, precedence);
+            (_num_1()*overall_coeff).print(os, precedence);
         else
             overall_coeff.print(os, precedence);
         os << '*';
@@ -223,7 +224,7 @@ void mul::printcsrc(ostream & os, unsigned type, unsigned upper_precedence) cons
     if (precedence <= upper_precedence)
         os << "(";
 
-    if (!overall_coeff.is_equal(exONE())) {
+    if (!overall_coeff.is_equal(_ex1())) {
         overall_coeff.bp->printcsrc(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/<expr>"
-        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; l<number_of_adds; l++) {
             add const & addref=ex_to_add(expanded_seq[positions_of_adds[l]].rest);
-            GINAC_ASSERT(term[positions_of_adds[l]].coeff.compare(exONE())==0);
+            GINAC_ASSERT(term[positions_of_adds[l]].coeff.compare(_ex1())==0);
             term[positions_of_adds[l]]=split_ex_to_pair(addref.op(k[l]));
         }
         /*
index d94e61c..5dce23f 100644 (file)
@@ -29,6 +29,7 @@
 #include "add.h"
 #include "mul.h"
 #include "debugmsg.h"
+#include "utils.h"
 
 #ifndef NO_GINAC_NAMESPACE
 namespace GiNaC {
@@ -302,7 +303,7 @@ ex ncmul::coeff(symbol const & s, int const n) const
 
     if (coeff_found) return (new ncmul(coeffseq,1))->setflag(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];
     }
index 9a3cd15..d6661d9 100644 (file)
@@ -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; i<e.nops(); i++) {
             c = lcmcoeff(e.op(i), c);
         }
@@ -210,7 +211,7 @@ static numeric lcmcoeff(const ex &e, const numeric &l)
 
 static numeric lcm_of_coefficients_denominators(const ex &e)
 {
-    return lcmcoeff(e.expand(), numONE());
+    return lcmcoeff(e.expand(), _num1());
 }
 
 
@@ -228,7 +229,7 @@ numeric ex::integer_content(void) const
 
 numeric basic::integer_content(void) const
 {
-    return numONE();
+    return _num1();
 }
 
 numeric numeric::integer_content(void) const
@@ -240,7 +241,7 @@ numeric add::integer_content(void) const
 {
     epvector::const_iterator it = seq.begin();
     epvector::const_iterator itend = seq.end();
-    numeric c = numZERO();
+    numeric c = _num0();
     while (it != itend) {
         GINAC_ASSERT(!is_ex_exactly_of_type(it->rest,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, ex2_less> 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<numeric> 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;
index b5cc40b..472d9c3 100644 (file)
@@ -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<numeric const &>((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<const numeric *>(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<const numeric *>(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<const numeric *>(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<const numeric *>(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<const numeric *>(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<const numeric *>(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<numeric> 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.
index 68c227d..27153df 100644 (file)
@@ -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);
index c21e0ea..0fc42e7 100644 (file)
@@ -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;
 }
 
index 1c9a342..51e5bfc 100644 (file)
@@ -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 << ")";
 
     // <expr>^-1 is printed as "1.0/<expr>" 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<numeric const &>(*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
index b1136c3..17c06c7 100644 (file)
@@ -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
index 19a2e3a..614c162 100644 (file)
@@ -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;
index 0cc2be3..4189b09 100644 (file)
@@ -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; i<deg; i++) {
-        ex sum = exZERO();
+        ex sum = _ex0();
         for (int j=1; j<=i; j++) {
             ex c = coeff(*s, j + ldeg);
             if (is_order_function(c)) {
-                co.push_back(Order(exONE()));
+                co.push_back(Order(_ex1()));
                 break;
             } else
                 sum += (p * j - (i - j)) * co[i - j] * c;
@@ -602,7 +603,7 @@ ex series::power_const(const numeric &p, int deg) const
         }
     }
     if (!higher_order && !all_sums_zero)
-        new_seq.push_back(expair(Order(exONE()), numeric(deg) + p * ldeg));
+        new_seq.push_back(expair(Order(_ex1()), numeric(deg) + p * ldeg));
     return series::series(var, point, new_seq);
 }
 
index a102b74..ae07019 100644 (file)
@@ -33,6 +33,7 @@
 #include "mul.h"
 #include "symbol.h"
 #include "debugmsg.h"
+#include "utils.h"
 
 #ifndef NO_GINAC_NAMESPACE
 namespace GiNaC {
@@ -203,7 +204,7 @@ ex simp_lor::eval(int level) const
         int sig=canonicalize_indices(iv,false); // symmetric
         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)*simp_lor(type,name,iv);
         }
         lorentzidx const & idx1=ex_to_lorentzidx(seq[0]);
@@ -214,19 +215,19 @@ ex simp_lor::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)) {
@@ -339,7 +340,7 @@ ex simplify_simp_lor_mul(ex const & m, scalar_products const & sp)
     v_contracted.reserve(2*n);
     for (int i=0; i<n; ++i) {
         ex f=m.op(i);
-        if (is_ex_exactly_of_type(f,power)&&f.op(1).is_equal(exTWO())) {
+        if (is_ex_exactly_of_type(f,power)&&f.op(1).is_equal(_ex2())) {
             v_contracted.push_back(f.op(0));
             v_contracted.push_back(f.op(0));
         } else {
@@ -374,7 +375,7 @@ ex simplify_simp_lor_mul(ex const & m, scalar_products const & sp)
                 } else {
                     // a contracted index should occur exactly once
                     GINAC_ASSERT(replacements==1);
-                    *it=exONE();
+                    *it=_ex1();
                     something_changed=true;
                 }
             }
@@ -390,7 +391,7 @@ ex simplify_simp_lor_mul(ex const & m, scalar_products const & sp)
                 } else {
                     // a contracted index should occur exactly once
                     GINAC_ASSERT(replacements==1);
-                    *it=exONE();
+                    *it=_ex1();
                     something_changed=true;
                 }
             }
@@ -418,7 +419,7 @@ ex simplify_simp_lor_mul(ex const & m, scalar_products const & sp)
                         idx1.is_co_contra_pair(idx2) &&
                         sp.is_defined(vec1,vec2)) {
                         *it1=sp.evaluate(vec1,vec2);
-                        *it2=exONE();
+                        *it2=_ex1();
                         something_changed=true;
                         jump_to_next=true;
                     }
@@ -442,7 +443,7 @@ ex simplify_simp_lor(ex const & e, scalar_products const & sp)
 
     // simplification of sum=sum of simplifications
     if (is_ex_exactly_of_type(e_expanded,add)) {
-        ex sum=exZERO();
+        ex sum=_ex0();
         for (int i=0; i<e_expanded.nops(); ++i) {
             sum += simplify_simp_lor(e_expanded.op(i),sp);
         }
@@ -458,11 +459,11 @@ ex simplify_simp_lor(ex const & e, scalar_products const & sp)
     return e_expanded;
 }
 
-ex Dim(void)
-{
-    static symbol * d=new symbol("dim");
-    return *d;
-}
+//ex Dim(void)   // FIXME: what's going on here?
+//{
+//    static symbol * d=new symbol("dim");
+//    return *d;
+//}
 
 //////////
 // helper classes
index 2ea8ed5..94c008c 100644 (file)
@@ -177,9 +177,9 @@ int symbol::ldegree(symbol const & s) const
 ex symbol::coeff(symbol const & s, int const n) const
 {
     if (compare_same_type(s)==0) {
-        return n==1 ? exONE() : exZERO();
+        return n==1 ? _ex1() : _ex0();
     } else {
-        return n==0 ? *this : exZERO();
+        return n==0 ? *this : _ex0();
     }
 }
 
@@ -277,7 +277,7 @@ void symbol::unassign(void)
 {
     if (asexinfop->is_assigned) {
         asexinfop->is_assigned=0;
-        asexinfop->assigned_expression=exZERO();
+        asexinfop->assigned_expression=_ex0();
     }
     setflag(status_flags::evaluated);
 }
index 7592e40..e113c04 100644 (file)
@@ -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<const numeric *>(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<const numeric *>(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<const numeric *>(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<const numeric *>(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<const numeric *>(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<const numeric *>(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<const numeric *>(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<const numeric *>(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<const numeric *>(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<const numeric *>(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<const numeric *>(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<const numeric *>(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<const numeric *>(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<const numeric *>(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<const numeric *>(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<const numeric *>(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<const numeric *>(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<const numeric *>(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<const numeric *>(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<const numeric *>(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<const numeric *>(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<const numeric *>(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<const numeric *>(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<const numeric *>(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<const numeric *>(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<const numeric *>(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<const numeric *>(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<const numeric *>(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<const numeric *>(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
index 7be44f7..0f6082a 100644 (file)
@@ -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
index 2e02b6c..0134ab1 100644 (file)
@@ -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;
 }