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 ";
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;
// Expansion
e = pow(x, 2) - (x+1)*(x-1) - 1;
- d = exZERO();
+ d = ex(0);
result += check_normal(e, d);
// Expansion inside functions
// 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;
}
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;
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;
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;
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;
}
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;
}
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;
}
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;
}
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;
}
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;
}
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;
}
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;
}
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;
}
{
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;
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;
}
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;
}
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;
}
#include "add.h"
#include "mul.h"
#include "debugmsg.h"
+#include "utils.h"
#ifndef NO_GINAC_NAMESPACE
namespace GiNaC {
{
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());
}
{
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());
}
{
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());
}
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 << '*';
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 << "/";
// 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);
}
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;
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;
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;
}
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()));
}
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,
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);
}
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)));
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;
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
#include "numeric.h"
#include "relational.h"
#include "debugmsg.h"
+#include "utils.h"
#ifndef NO_GINAC_NAMESPACE
namespace GiNaC {
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);
}
}
// 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();
}
}
}
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
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)) {
} 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;
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)) {
} 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;
}
} else {
// a contracted index should occur exactly twice
GINAC_ASSERT(replacements==2);
- *it=exONE();
+ *it=_ex1();
something_changed=true;
}
}
} else {
// a contracted index should occur exactly twice
GINAC_ASSERT(replacements==2);
- *it=exONE();
+ *it=_ex1();
something_changed=true;
}
}
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();
}
}
}
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));
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));
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();
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(
// 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));
}
// 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));
}
counter[l]=1;
}
- ex sum=exZERO();
+ ex sum=_ex0();
while (1) {
ex term=e;
#include "relational.h"
#include "series.h"
#include "symbol.h"
+#include "utils.h"
#ifndef NO_GINAC_NAMESPACE
namespace GiNaC {
* @see ex::diff */
ex numeric::diff(symbol const & s) const
{
- return exZERO();
+ return _ex0();
}
ex symbol::diff(symbol const & s) const
{
if (compare_same_type(s)) {
- return exZERO();
+ return _ex0();
} else {
- return exONE();
+ return _ex1();
}
}
* @see ex::diff */
ex constant::diff(symbol const & s) const
{
- return exZERO();
+ return _ex0();
}
/** Implementation of ex::diff() for multiple differentiation of a symbol.
return s;
break;
case 1:
- return exONE();
+ return _ex1();
break;
default:
- return exZERO();
+ return _ex0();
}
} else {
- return exONE();
+ return _ex1();
}
}
* @see ex::diff */
ex indexed::diff(symbol const & s) const
{
- return exZERO();
+ return _ex0();
}
* @see ex::diff */
ex ncmul::diff(symbol const & s) const
{
- return exZERO();
+ return _ex0();
}
{
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),
#include "numeric.h"
#include "power.h"
#include "debugmsg.h"
+#include "utils.h"
#ifndef NO_GINAC_NAMESPACE
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;
}
// 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);
// 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))
// 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);
// 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
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
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
#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;
{
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
*/
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;
}
{
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;
}
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,
ex expairseq::default_overall_coeff(void) const
{
- return exZERO();
+ return _ex0();
}
void expairseq::combine_overall_coeff(ex const & c)
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 {
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;
}
}
#include "relational.h"
#include "series.h"
#include "symbol.h"
+#include "utils.h"
#ifndef NO_GINAC_NAMESPACE
namespace GiNaC {
{
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();
}
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)) {
{
// 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);
}
ex ncpower(ex const &basis, unsigned exponent)
{
if (exponent==0) {
- return exONE();
+ return _ex1();
}
exvector v;
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);
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"));
}
// 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);
}
}
{
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);
}
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"));
}
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();
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;
}
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)
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)
{
// 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
}
}
{
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)
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();
}
// 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
{
// 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
{
GINAC_ASSERT(diff_param==0);
+ // d/dx exp(x) -> exp(x)
return exp(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))
{
GINAC_ASSERT(diff_param==0);
+ // d/dx log(x) -> 1/x
return power(x, -1);
}
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)) {
// 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
{
GINAC_ASSERT(diff_param==0);
+ // d/dx sin(x) -> cos(x)
return cos(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)) {
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
{
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);
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
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
{
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)
// 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);
}
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);
{
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);
{
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))
{
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);
{
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);
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);
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();
{
GINAC_ASSERT(diff_param==0);
+ // d/dx sinh(x) -> cosh(x)
return cosh(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);
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();
{
GINAC_ASSERT(diff_param==0);
+ // d/dx cosh(x) -> sinh(x)
return sinh(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
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();
{
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)
// 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);
}
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);
{
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);
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))
{
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);
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))
{
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);
#include "numeric.h"
#include "power.h"
#include "symbol.h"
+#include "utils.h"
#ifndef NO_GINAC_NAMESPACE
namespace GiNaC {
// 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();
}
}
}
{
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);
#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 {
//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();
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;
} else {
// a contracted index should occur exactly once
GINAC_ASSERT(replacements==1);
- *it=exONE();
+ *it=_ex1();
something_changed=true;
}
}
} else {
// a contracted index should occur exactly once
GINAC_ASSERT(replacements==1);
- *it=exONE();
+ *it=_ex1();
something_changed=true;
}
}
// 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));
}
#include "matrix.h"
#include "debugmsg.h"
+#include "utils.h"
#ifndef NO_GINAC_NAMESPACE
namespace GiNaC {
: basic(TINFO_matrix), row(1), col(1)
{
debugmsg("matrix default constructor",LOGLEVEL_CONSTRUCT);
- m.push_back(exZERO());
+ m.push_back(_ex0());
}
matrix::~matrix()
: 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
{
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) {
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);
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) {
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;
#include "add.h"
#include "power.h"
#include "debugmsg.h"
+#include "utils.h"
#ifndef NO_GINAC_NAMESPACE
namespace GiNaC {
{
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());
}
{
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());
}
{
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());
}
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());
}
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 << '*';
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 << "*";
}
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
}
// 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:
// 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 << "*";
return (new mul(coeffseq))->setflag(status_flags::dynallocated);
}
- return exZERO();
+ return _ex0();
}
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;
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,
// 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));
// 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));
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);
*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;
}
ex mul::default_overall_coeff(void) const
{
- return exONE();
+ return _ex1();
}
void mul::combine_overall_coeff(ex const & c)
// 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
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();
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]));
}
/*
#include "add.h"
#include "mul.h"
#include "debugmsg.h"
+#include "utils.h"
#ifndef NO_GINAC_NAMESPACE
namespace GiNaC {
if (coeff_found) return (new ncmul(coeffseq,1))->setflag(status_flags::dynallocated);
- return exZERO();
+ return _ex0();
}
unsigned ncmul::count_factors(ex const & e) 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;
ex simplified_ncmul(exvector const & v)
{
if (v.size()==0) {
- return exONE();
+ return _ex1();
} else if (v.size()==1) {
return v[0];
}
#include "relational.h"
#include "series.h"
#include "symbol.h"
+#include "utils.h"
#ifndef NO_GINAC_NAMESPACE
namespace GiNaC {
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);
}
static numeric lcm_of_coefficients_denominators(const ex &e)
{
- return lcmcoeff(e.expand(), numONE());
+ return lcmcoeff(e.expand(), _num1());
}
numeric basic::integer_content(void) const
{
- return numONE();
+ return _num1();
}
numeric numeric::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));
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;
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"));
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;
}
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;
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)) {
return false;
#if FAST_COMPARE
if (a.is_equal(b)) {
- q = exONE();
+ q = _ex1();
return true;
}
#endif
* @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;
}
}
#if FAST_COMPARE
if (a.is_equal(b)) {
- q = exONE();
+ q = _ex1();
return true;
}
#endif
// 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
{
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))
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();
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;
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);
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))
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 (;;) {
numeric basic::max_coefficient(void) const
{
- return numONE();
+ return _num1();
}
numeric numeric::max_coefficient(void) 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++) {
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);
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;
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
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);
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))
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);
{
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"));
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
// 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;
}
// 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;
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);
}
// Add fractions
- if (den.is_equal(exONE()))
+ if (den.is_equal(_ex1()))
return (new add(o))->setflag(status_flags::dynallocated);
else {
exvector num_seq;
#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:
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() << "\")";
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:
* 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));
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"));
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)))->
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"));
numeric numeric::denom(void) const
{
if (is_integer()) {
- return numONE();
+ return _num1();
}
#ifdef SANE_LINKER
if (instanceof(*value, cl_RA_ring)) {
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))
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))
}
#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
* 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). */
{
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
}
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];
}
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++) {
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];
}
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++) {
{
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);
}
}
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
// 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))
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).
* @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.
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.
return rem_quo.remainder;
}
else {
- q = numZERO();
- return numZERO(); // Throw?
+ q = _num0();
+ return _num0(); // Throw?
}
}
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.
r = rem_quo.remainder;
return rem_quo.quotient;
} else {
- r = numZERO();
- return numZERO(); // Throw?
+ r = _num0();
+ return _num0(); // Throw?
}
}
::isqrt(The(cl_I)(*x.value), &root); // -> CLN
return root;
} else
- return numZERO(); // Throw?
+ return _num0(); // Throw?
}
/** Greatest Common Divisor.
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.
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
// 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);
#include "power.h"
#include "relational.h"
#include "debugmsg.h"
+#include "utils.h"
#ifndef NO_GINAC_NAMESPACE
namespace GiNaC {
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)
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)
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
return lh.div(rh);
}
+
// binary arithmetic assignment operators with ex
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
ex operator-(ex const & lh)
{
- return exMINUSONE()*lh;
+ return lh.exmul(_ex_1());
}
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;
}
numeric operator++(numeric & lh, int)
{
numeric tmp = lh;
- lh = lh+numONE();
+ lh = lh+_num1();
return tmp;
}
numeric operator--(numeric & lh, int)
{
numeric tmp = lh;
- lh = lh-numONE();
+ lh = lh-_num1();
return tmp;
}
#include "relational.h"
#include "symbol.h"
#include "debugmsg.h"
+#include "utils.h"
#ifndef NO_GINAC_NAMESPACE
namespace GiNaC {
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
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();
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();
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
// ^(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)
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(),
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();
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),
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
!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())));
}
}
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);
{
// 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;
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).
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
/** 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
#include "relational.h"
#include "numeric.h"
#include "debugmsg.h"
+#include "utils.h"
#ifndef NO_GINAC_NAMESPACE
namespace GiNaC {
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;
#include "relational.h"
#include "symbol.h"
#include "debugmsg.h"
+#include "utils.h"
#ifndef NO_GINAC_NAMESPACE
namespace GiNaC {
/** 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)
}
// 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();
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);
}
// 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);
}
// 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);
}
} 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;
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++;
}
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
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
// 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);
}
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);
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);
}
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++;
}
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
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
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;
}
}
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);
}
#include "mul.h"
#include "symbol.h"
#include "debugmsg.h"
+#include "utils.h"
#ifndef NO_GINAC_NAMESPACE
namespace GiNaC {
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]);
// 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)) {
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 {
} else {
// a contracted index should occur exactly once
GINAC_ASSERT(replacements==1);
- *it=exONE();
+ *it=_ex1();
something_changed=true;
}
}
} else {
// a contracted index should occur exactly once
GINAC_ASSERT(replacements==1);
- *it=exONE();
+ *it=_ex1();
something_changed=true;
}
}
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;
}
// 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);
}
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
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();
}
}
{
if (asexinfop->is_assigned) {
asexinfop->is_assigned=0;
- asexinfop->assigned_expression=exZERO();
+ asexinfop->assigned_expression=_ex0();
}
setflag(status_flags::evaluated);
}
/** @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
* 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
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
// private
// none
+
#ifndef NO_GINAC_NAMESPACE
} // namespace GiNaC
#endif // ndef NO_GINAC_NAMESPACE
/** @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
}
}
+// 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
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)
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)
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)
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;
}