From f30eb24864efd57bd7271cfccd8e526321682266 Mon Sep 17 00:00:00 2001 From: Richard Kreckel Date: Wed, 27 Jun 2001 23:07:30 +0000 Subject: [PATCH] - eta(x,y) was broken all along since ages: eta(I,I) for instance ought to return 0 but... that formula that reduces eta(x,y) to theta functions is basically to naive. --- check/exam_inifcns.cpp | 25 ++++++++++++- ginac/inifcns.cpp | 85 ++++++++++++++++++++++++++---------------- 2 files changed, 76 insertions(+), 34 deletions(-) diff --git a/check/exam_inifcns.cpp b/check/exam_inifcns.cpp index 47cc1040..96ea9cc7 100644 --- a/check/exam_inifcns.cpp +++ b/check/exam_inifcns.cpp @@ -90,6 +90,27 @@ static unsigned inifcns_consist_trans(void) ++result; } + // check consistency of log and eta phases: + for (int r1=-1; r1<=1; ++r1) { + for (int i1=-1; i1<=1; ++i1) { + ex x1 = r1+I*i1; + if (x1.is_zero()) + continue; + for (int r2=-1; r2<=1; ++r2) { + for (int i2=-1; i2<=1; ++i2) { + ex x2 = r2+I*i2; + if (x2.is_zero()) + continue; + if (abs(evalf(eta(x1,x2)-log(x1*x2)+log(x1)+log(x2)))>.1e-12) { + clog << "either eta(x,y), log(x), log(y) or log(x*y) is wrong" + << " at x==" << x1 << ", y==" << x2 << endl; + ++result; + } + } + } + } + } + return result; } @@ -100,7 +121,7 @@ static unsigned inifcns_consist_gamma(void) unsigned result = 0; ex e; - e = tgamma(ex(1)); + e = tgamma(1); for (int i=2; i<8; ++i) e += tgamma(ex(i)); if (e != numeric(874)) { @@ -109,7 +130,7 @@ static unsigned inifcns_consist_gamma(void) ++result; } - e = tgamma(ex(1)); + e = tgamma(1); for (int i=2; i<8; ++i) e *= tgamma(ex(i)); if (e != numeric(24883200)) { diff --git a/ginac/inifcns.cpp b/ginac/inifcns.cpp index f64597e5..316eea6c 100644 --- a/ginac/inifcns.cpp +++ b/ginac/inifcns.cpp @@ -126,57 +126,78 @@ REGISTER_FUNCTION(csgn, eval_func(csgn_eval). ////////// -// Eta function: log(x*y) == log(x) + log(y) + eta(x,y). +// Eta function: eta(x,y) == log(x*y) - log(x) - log(y). ////////// -static ex eta_evalf(const ex & x, const ex & y) +static ex eta_evalf(const ex &x, const ex &y) { - BEGIN_TYPECHECK - TYPECHECK(x,numeric) - TYPECHECK(y,numeric) - END_TYPECHECK(eta(x,y)) - - numeric xim = imag(ex_to(x)); - numeric yim = imag(ex_to(y)); - numeric xyim = imag(ex_to(x*y)); - return evalf(I/4*Pi)*((csgn(-xim)+1)*(csgn(-yim)+1)*(csgn(xyim)+1)-(csgn(xim)+1)*(csgn(yim)+1)*(csgn(-xyim)+1)); + // It seems like we basically have to replicate the eval function here, + // since the expression might not be fully evaluated yet. + if (x.info(info_flags::positive) || y.info(info_flags::positive)) + return _ex0(); + + if (x.info(info_flags::numeric) && y.info(info_flags::numeric)) { + const numeric nx = ex_to(x); + const numeric ny = ex_to(y); + const numeric nxy = ex_to(x*y); + int cut = 0; + if (nx.is_real() && nx.is_negative()) + cut -= 4; + if (ny.is_real() && ny.is_negative()) + cut -= 4; + if (nxy.is_real() && nxy.is_negative()) + cut += 4; + return evalf(I/4*Pi)*((csgn(-imag(nx))+1)*(csgn(-imag(ny))+1)*(csgn(imag(nxy))+1)- + (csgn(imag(nx))+1)*(csgn(imag(ny))+1)*(csgn(-imag(nxy))+1)+cut); + } + + return eta(x,y).hold(); } -static ex eta_eval(const ex & x, const ex & y) +static ex eta_eval(const ex &x, const ex &y) { - if (is_ex_exactly_of_type(x, numeric) && - is_ex_exactly_of_type(y, numeric)) { + // trivial: eta(x,c) -> 0 if c is real and positive + if (x.info(info_flags::positive) || y.info(info_flags::positive)) + return _ex0(); + + if (x.info(info_flags::numeric) && y.info(info_flags::numeric)) { // don't call eta_evalf here because it would call Pi.evalf()! - numeric xim = imag(ex_to(x)); - numeric yim = imag(ex_to(y)); - numeric xyim = imag(ex_to(x*y)); - return (I/4)*Pi*((csgn(-xim)+1)*(csgn(-yim)+1)*(csgn(xyim)+1)-(csgn(xim)+1)*(csgn(yim)+1)*(csgn(-xyim)+1)); + const numeric nx = ex_to(x); + const numeric ny = ex_to(y); + const numeric nxy = ex_to(x*y); + int cut = 0; + if (nx.is_real() && nx.is_negative()) + cut -= 4; + if (ny.is_real() && ny.is_negative()) + cut -= 4; + if (nxy.is_real() && nxy.is_negative()) + cut += 4; + return (I/4)*Pi*((csgn(-imag(nx))+1)*(csgn(-imag(ny))+1)*(csgn(imag(nxy))+1)- + (csgn(imag(nx))+1)*(csgn(imag(ny))+1)*(csgn(-imag(nxy))+1)+cut); } return eta(x,y).hold(); } -static ex eta_series(const ex & arg1, - const ex & arg2, +static ex eta_series(const ex & x, const ex & y, const relational & rel, int order, unsigned options) { - const ex arg1_pt = arg1.subs(rel); - const ex arg2_pt = arg2.subs(rel); - if (ex_to(arg1_pt).imag().is_zero() || - ex_to(arg2_pt).imag().is_zero() || - ex_to(arg1_pt*arg2_pt).imag().is_zero()) { - throw (std::domain_error("eta_series(): on discontinuity")); - } + const ex x_pt = x.subs(rel); + const ex y_pt = y.subs(rel); + if ((x_pt.info(info_flags::numeric) && x_pt.info(info_flags::negative)) || + (y_pt.info(info_flags::numeric) && y_pt.info(info_flags::negative)) || + ((x_pt*y_pt).info(info_flags::numeric) && (x_pt*y_pt).info(info_flags::negative))) + throw (std::domain_error("eta_series(): on discontinuity")); epvector seq; - seq.push_back(expair(eta(arg1_pt,arg2_pt), _ex0())); + seq.push_back(expair(eta(x_pt,y_pt), _ex0())); return pseries(rel,seq); } REGISTER_FUNCTION(eta, eval_func(eta_eval). evalf_func(eta_evalf). - series_func(eta_series). + series_func(eta_series). latex_name("\\eta")); @@ -419,7 +440,7 @@ ex lsolve(const ex &eqns, const ex &symbols) if (eqns.info(info_flags::relation_equal)) { if (!symbols.info(info_flags::symbol)) throw(std::invalid_argument("lsolve(): 2nd argument must be a symbol")); - ex sol=lsolve(lst(eqns),lst(symbols)); + const ex sol = lsolve(lst(eqns),lst(symbols)); GINAC_ASSERT(sol.nops()==1); GINAC_ASSERT(is_ex_exactly_of_type(sol.op(0),relational)); @@ -451,10 +472,10 @@ ex lsolve(const ex &eqns, const ex &symbols) matrix vars(symbols.nops(),1); for (unsigned r=0; r(symbols.op(c)),1); + const ex co = eq.coeff(ex_to(symbols.op(c)),1); linpart -= co*symbols.op(c); sys(r,c) = co; } -- 2.44.0