From: Richard Kreckel Date: Sun, 26 Mar 2000 22:57:01 +0000 (+0000) Subject: - As advertised: we are calling the Gamma function tgamma() now! X-Git-Tag: release_0-6-0~48 X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=commitdiff_plain;h=5b090bb7e4951b48a28c20ba21bf9810c86eb0ca - As advertised: we are calling the Gamma function tgamma() now! --- diff --git a/NEWS b/NEWS index 699284cc..25d89c6a 100644 --- a/NEWS +++ b/NEWS @@ -5,12 +5,14 @@ This file records noteworthy changes. much clearer but break compatibility with older versions: - f(x).series(x,p[,o]) -> f(x).series(x==p,o) - series(f(x),x,p[,o]) -> series(f(x),x==p,o) - - gamma() -> Gamma() + - gamma() -> tgamma() (The true Gamma function, there is now also + log(tgamma()), called lgamma(), in accord with ISO/IEC 9899:1999.) - EulerGamma -> gamma - - beta() -> Beta() * #include'ing ginac.h defines the preprocessor symbols GINACLIB_MAJOR_VERSION, GINACLIB_MINOR_VERSION, and GINACLIB_MICRO_VERSION with the respective GiNaC library version numbers. +* Several new timings in the check target. Some of them may be rather rude + at your machine. 0.5.4 (15 March 2000) * Some algorithms in class matrix (notably determinant) were replaced by diff --git a/check/exam_inifcns.cpp b/check/exam_inifcns.cpp index 16788c54..bdb5e7d7 100644 --- a/check/exam_inifcns.cpp +++ b/check/exam_inifcns.cpp @@ -93,44 +93,44 @@ static unsigned inifcns_consist_trans(void) return result; } -/* Simple tests on the Gamma function. We stuff in arguments where the results +/* Simple tests on the tgamma function. We stuff in arguments where the results * exists in closed form and check if it's ok. */ static unsigned inifcns_consist_gamma(void) { unsigned result = 0; ex e; - e = Gamma(ex(1)); + e = tgamma(ex(1)); for (int i=2; i<8; ++i) - e += Gamma(ex(i)); + e += tgamma(ex(i)); if (e != numeric(874)) { - clog << "Gamma(1)+...+Gamma(7) erroneously returned " + clog << "tgamma(1)+...+tgamma(7) erroneously returned " << e << " instead of 874" << endl; ++result; } - e = Gamma(ex(1)); + e = tgamma(ex(1)); for (int i=2; i<8; ++i) - e *= Gamma(ex(i)); + e *= tgamma(ex(i)); if (e != numeric(24883200)) { - clog << "Gamma(1)*...*Gamma(7) erroneously returned " + clog << "tgamma(1)*...*tgamma(7) erroneously returned " << e << " instead of 24883200" << endl; ++result; } - e = Gamma(ex(numeric(5, 2)))*Gamma(ex(numeric(9, 2)))*64; + e = tgamma(ex(numeric(5, 2)))*tgamma(ex(numeric(9, 2)))*64; if (e != 315*Pi) { - clog << "64*Gamma(5/2)*Gamma(9/2) erroneously returned " + clog << "64*tgamma(5/2)*tgamma(9/2) erroneously returned " << e << " instead of 315*Pi" << endl; ++result; } - e = Gamma(ex(numeric(-13, 2))); + e = tgamma(ex(numeric(-13, 2))); for (int i=-13; i<7; i=i+2) - e += Gamma(ex(numeric(i, 2))); - e = (e*Gamma(ex(numeric(15, 2)))*numeric(512)); + e += tgamma(ex(numeric(i, 2))); + e = (e*tgamma(ex(numeric(15, 2)))*numeric(512)); if (e != numeric(633935)*Pi) { - clog << "512*(Gamma(-13/2)+...+Gamma(5/2))*Gamma(15/2) erroneously returned " + clog << "512*(tgamma(-13/2)+...+tgamma(5/2))*tgamma(15/2) erroneously returned " << e << " instead of 633935*Pi" << endl; ++result; } @@ -147,11 +147,11 @@ static unsigned inifcns_consist_psi(void) ex e, f; // We check psi(1) and psi(1/2) implicitly by calculating the curious - // little identity Gamma(1)'/Gamma(1) - Gamma(1/2)'/Gamma(1/2) == 2*log(2). - e += (Gamma(x).diff(x)/Gamma(x)).subs(x==numeric(1)); - e -= (Gamma(x).diff(x)/Gamma(x)).subs(x==numeric(1,2)); + // little identity tgamma(1)'/tgamma(1) - tgamma(1/2)'/tgamma(1/2) == 2*log(2). + e += (tgamma(x).diff(x)/tgamma(x)).subs(x==numeric(1)); + e -= (tgamma(x).diff(x)/tgamma(x)).subs(x==numeric(1,2)); if (e!=2*log(2)) { - clog << "Gamma(1)'/Gamma(1) - Gamma(1/2)'/Gamma(1/2) erroneously returned " + clog << "tgamma(1)'/tgamma(1) - tgamma(1/2)'/tgamma(1/2) erroneously returned " << e << " instead of 2*log(2)" << endl; ++result; } diff --git a/check/exam_pseries.cpp b/check/exam_pseries.cpp index ba7b3193..78a55e2d 100644 --- a/check/exam_pseries.cpp +++ b/check/exam_pseries.cpp @@ -153,8 +153,8 @@ static unsigned exam_series5(void) unsigned result = 0; ex e, d; - // Gamma(-1): - e = Gamma(2*x); + // tgamma(-1): + e = tgamma(2*x); d = pow(x+1,-1)*numeric(1,4) + pow(x+1,0)*(numeric(3,4) - numeric(1,2)*gamma) + diff --git a/check/time_gammaseries.cpp b/check/time_gammaseries.cpp index c92f29f7..0bfbe302 100644 --- a/check/time_gammaseries.cpp +++ b/check/time_gammaseries.cpp @@ -22,12 +22,12 @@ #include "times.h" -unsigned Gammaseries(unsigned order) +unsigned tgammaseries(unsigned order) { unsigned result = 0; symbol x; - ex myseries = series(Gamma(x),x==0,order); + ex myseries = series(tgamma(x),x==0,order); // compute the last coefficient numerically: ex last_coeff = myseries.coeff(x,order-1).evalf(); // compute a bound for that coefficient using a variation of the leading @@ -35,7 +35,7 @@ unsigned Gammaseries(unsigned order) ex bound = evalf(exp(ex(-.57721566490153286*(order-1)))/(order-1)); if (evalf(abs((last_coeff-pow(-1,order))/bound)) > numeric(1)) { clog << "The " << order-1 - << "th order coefficient in the power series expansion of Gamma(0) was erroneously found to be " + << "th order coefficient in the power series expansion of tgamma(0) was erroneously found to be " << last_coeff << ", violating a simple estimate." << endl; ++result; } @@ -61,7 +61,7 @@ unsigned time_gammaseries(void) for (vector::iterator i=sizes.begin(); i!=sizes.end(); ++i) { omega.start(); - result += Gammaseries(*i); + result += tgammaseries(*i); times.push_back(omega.read()); cout << '.' << flush; } diff --git a/check/time_lw_O.cpp b/check/time_lw_O.cpp index c6af8c8c..d0376ad1 100644 --- a/check/time_lw_O.cpp +++ b/check/time_lw_O.cpp @@ -92,7 +92,7 @@ static unsigned test1(void) unsigned nops3 = nops(d3.determinant()); cout << '.' << flush; if ((nops1 != 37490) || (nops2 != 37490) || (nops3 != 37490)) { - clog << "Determinant from van der Waerden's example were miscalculated" << endl; + clog << "Determinants were miscalculated" << endl; return 1; } return 0; diff --git a/ginac/inifcns.cpp b/ginac/inifcns.cpp index 3e74166d..96f4c9a5 100644 --- a/ginac/inifcns.cpp +++ b/ginac/inifcns.cpp @@ -301,7 +301,7 @@ ex ncpower(const ex &basis, unsigned exponent) /** Force inclusion of functions from initcns_gamma and inifcns_zeta * for static lib (so ginsh will see them). */ -unsigned force_include_gamma = function_index_Gamma; +unsigned force_include_tgamma = function_index_tgamma; unsigned force_include_zeta1 = function_index_zeta1; #ifndef NO_NAMESPACE_GINAC diff --git a/ginac/inifcns.h b/ginac/inifcns.h index cff2bb24..75599ea9 100644 --- a/ginac/inifcns.h +++ b/ginac/inifcns.h @@ -97,10 +97,11 @@ inline function zeta(const ex & p1, const ex & p2) { } /** Gamma-function. */ -DECLARE_FUNCTION_1P(Gamma) +DECLARE_FUNCTION_1P(lgamma) +DECLARE_FUNCTION_1P(tgamma) /** Beta-function. */ -DECLARE_FUNCTION_2P(Beta) +DECLARE_FUNCTION_2P(beta) // overloading at work: we cannot use the macros /** Psi-function (aka digamma-function). */ diff --git a/ginac/inifcns_gamma.cpp b/ginac/inifcns_gamma.cpp index 5b35963f..9e1b48f4 100644 --- a/ginac/inifcns_gamma.cpp +++ b/ginac/inifcns_gamma.cpp @@ -39,40 +39,93 @@ namespace GiNaC { #endif // ndef NO_NAMESPACE_GINAC ////////// -// Gamma-function +// Logarithm of Gamma function ////////// -static ex Gamma_evalf(const ex & x) +static ex lgamma_evalf(const ex & x) { BEGIN_TYPECHECK TYPECHECK(x,numeric) - END_TYPECHECK(Gamma(x)) + END_TYPECHECK(lgamma(x)) - return Gamma(ex_to_numeric(x)); + return lgamma(ex_to_numeric(x)); } -/** Evaluation of Gamma(x). Knows about integer arguments, half-integer - * arguments and that's it. Somebody ought to provide some good numerical - * evaluation some day... +/** Evaluation of lgamma(x), the natural logarithm of the Gamma function. + * Knows about integer arguments and that's it. Somebody ought to provide + * some good numerical evaluation some day... * - * @exception std::domain_error("Gamma_eval(): simple pole") */ -static ex Gamma_eval(const ex & x) + * @exception std::domain_error("lgamma_eval(): simple pole") */ +static ex lgamma_eval(const ex & x) { if (x.info(info_flags::numeric)) { // trap integer arguments: if (x.info(info_flags::integer)) { - // Gamma(n+1) -> n! for postitive n + // lgamma(n) -> log((n-1)!) for postitive n + if (x.info(info_flags::posint)) { + return log(factorial(x.exadd(_ex_1()))); + } else { + throw (std::domain_error("lgamma_eval(): logarithmic pole")); + } + } + // lgamma_evalf should be called here once it becomes available + } + + return lgamma(x).hold(); +} + + +static ex lgamma_deriv(const ex & x, unsigned deriv_param) +{ + GINAC_ASSERT(deriv_param==0); + + // d/dx lgamma(x) -> psi(x) + return psi(x); +} + +// need to implement lgamma_series. + +REGISTER_FUNCTION(lgamma, eval_func(lgamma_eval). + evalf_func(lgamma_evalf). + derivative_func(lgamma_deriv)); + + +////////// +// true Gamma function +////////// + +static ex tgamma_evalf(const ex & x) +{ + BEGIN_TYPECHECK + TYPECHECK(x,numeric) + END_TYPECHECK(tgamma(x)) + + return tgamma(ex_to_numeric(x)); +} + + +/** Evaluation of tgamma(x), the true Gamma function. Knows about integer + * arguments, half-integer arguments and that's it. Somebody ought to provide + * some good numerical evaluation some day... + * + * @exception std::domain_error("tgamma_eval(): simple pole") */ +static ex tgamma_eval(const ex & x) +{ + if (x.info(info_flags::numeric)) { + // trap integer arguments: + if (x.info(info_flags::integer)) { + // tgamma(n) -> (n-1)! for postitive n if (x.info(info_flags::posint)) { return factorial(ex_to_numeric(x).sub(_num1())); } else { - throw (std::domain_error("Gamma_eval(): simple pole")); + throw (std::domain_error("tgamma_eval(): simple pole")); } } // trap half integer arguments: if ((x*2).info(info_flags::integer)) { // trap positive x==(n+1/2) - // Gamma(n+1/2) -> Pi^(1/2)*(1*3*..*(2*n-1))/(2^n) + // tgamma(n+1/2) -> Pi^(1/2)*(1*3*..*(2*n-1))/(2^n) if ((x*_ex2()).info(info_flags::posint)) { numeric n = ex_to_numeric(x).sub(_num1_2()); numeric coefficient = doublefactorial(n.mul(_num2()).sub(_num1())); @@ -80,40 +133,39 @@ static ex Gamma_eval(const ex & x) return coefficient * pow(Pi,_ex1_2()); } else { // trap negative x==(-n+1/2) - // Gamma(-n+1/2) -> Pi^(1/2)*(-2)^n/(1*3*..*(2*n-1)) + // tgamma(-n+1/2) -> Pi^(1/2)*(-2)^n/(1*3*..*(2*n-1)) numeric n = abs(ex_to_numeric(x).sub(_num1_2())); numeric coefficient = pow(_num_2(), n); coefficient = coefficient.div(doublefactorial(n.mul(_num2()).sub(_num1())));; return coefficient*power(Pi,_ex1_2()); } } - // Gamma_evalf should be called here once it becomes available + // tgamma_evalf should be called here once it becomes available } - return Gamma(x).hold(); -} + return tgamma(x).hold(); +} -static ex Gamma_deriv(const ex & x, unsigned deriv_param) +static ex tgamma_deriv(const ex & x, unsigned deriv_param) { GINAC_ASSERT(deriv_param==0); - // d/dx log(Gamma(x)) -> psi(x) - // d/dx Gamma(x) -> psi(x)*Gamma(x) - return psi(x)*Gamma(x); + // d/dx tgamma(x) -> psi(x)*tgamma(x) + return psi(x)*tgamma(x); } -static ex Gamma_series(const ex & x, const relational & r, int order) +static ex tgamma_series(const ex & x, const relational & r, int order) { // method: // 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 + // tgamma(x) == tgamma(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); + // series(tgamma(x),x,-m,order) == + // series(tgamma(x+m+1)/(x*(x+1)*...*(x+m)),x,-m,order+1); const ex x_pt = x.subs(r); if (!x_pt.info(info_flags::integer) || x_pt.info(info_flags::positive)) throw do_taylor(); // caught by function::series() @@ -122,54 +174,54 @@ static ex Gamma_series(const ex & x, const relational & r, int order) ex ser_denom = _ex1(); for (numeric p; p<=m; ++p) ser_denom *= x+p; - return (Gamma(x+m+_ex1())/ser_denom).series(r, order+1); + return (tgamma(x+m+_ex1())/ser_denom).series(r, order+1); } -REGISTER_FUNCTION(Gamma, eval_func(Gamma_eval). - evalf_func(Gamma_evalf). - derivative_func(Gamma_deriv). - series_func(Gamma_series)); +REGISTER_FUNCTION(tgamma, eval_func(tgamma_eval). + evalf_func(tgamma_evalf). + derivative_func(tgamma_deriv). + series_func(tgamma_series)); ////////// -// Beta-function +// beta-function ////////// -static ex Beta_evalf(const ex & x, const ex & y) +static ex beta_evalf(const ex & x, const ex & y) { BEGIN_TYPECHECK TYPECHECK(x,numeric) TYPECHECK(y,numeric) - END_TYPECHECK(Beta(x,y)) + END_TYPECHECK(beta(x,y)) - return Gamma(ex_to_numeric(x))*Gamma(ex_to_numeric(y))/Gamma(ex_to_numeric(x+y)); + return tgamma(ex_to_numeric(x))*tgamma(ex_to_numeric(y))/tgamma(ex_to_numeric(x+y)); } -static ex Beta_eval(const ex & x, const ex & y) +static ex beta_eval(const ex & x, const ex & y) { if (x.info(info_flags::numeric) && y.info(info_flags::numeric)) { - // treat all problematic x and y that may not be passed into Gamma, - // because they would throw there although Beta(x,y) is well-defined - // using the formula Beta(x,y) == (-1)^y * Beta(1-x-y, y) + // treat all problematic x and y that may not be passed into tgamma, + // because they would throw there although beta(x,y) is well-defined + // using the formula beta(x,y) == (-1)^y * beta(1-x-y, y) numeric nx(ex_to_numeric(x)); numeric ny(ex_to_numeric(y)); if (nx.is_real() && nx.is_integer() && ny.is_real() && ny.is_integer()) { if (nx.is_negative()) { if (nx<=-ny) - return pow(_num_1(), ny)*Beta(1-x-y, y); + return pow(_num_1(), ny)*beta(1-x-y, y); else - throw (std::domain_error("Beta_eval(): simple pole")); + throw (std::domain_error("beta_eval(): simple pole")); } if (ny.is_negative()) { if (ny<=-nx) - return pow(_num_1(), nx)*Beta(1-y-x, x); + return pow(_num_1(), nx)*beta(1-y-x, x); else - throw (std::domain_error("Beta_eval(): simple pole")); + throw (std::domain_error("beta_eval(): simple pole")); } - return Gamma(x)*Gamma(y)/Gamma(x+y); + return tgamma(x)*tgamma(y)/tgamma(x+y); } // no problem in numerator, but denominator has pole: if ((nx+ny).is_real() && @@ -177,34 +229,34 @@ static ex Beta_eval(const ex & x, const ex & y) !(nx+ny).is_positive()) return _ex0(); // everything is ok: - return Gamma(x)*Gamma(y)/Gamma(x+y); + return tgamma(x)*tgamma(y)/tgamma(x+y); } - return Beta(x,y).hold(); + return beta(x,y).hold(); } -static ex Beta_deriv(const ex & x, const ex & y, unsigned deriv_param) +static ex beta_deriv(const ex & x, const ex & y, unsigned deriv_param) { GINAC_ASSERT(deriv_param<2); ex retval; - // d/dx Beta(x,y) -> (psi(x)-psi(x+y)) * Beta(x,y) + // d/dx beta(x,y) -> (psi(x)-psi(x+y)) * beta(x,y) if (deriv_param==0) - retval = (psi(x)-psi(x+y))*Beta(x,y); - // d/dy Beta(x,y) -> (psi(y)-psi(x+y)) * Beta(x,y) + retval = (psi(x)-psi(x+y))*beta(x,y); + // d/dy beta(x,y) -> (psi(y)-psi(x+y)) * beta(x,y) if (deriv_param==1) - retval = (psi(y)-psi(x+y))*Beta(x,y); + retval = (psi(y)-psi(x+y))*beta(x,y); return retval; } -static ex Beta_series(const ex & x, const ex & y, const relational & r, int order) +static ex beta_series(const ex & x, const ex & y, const relational & r, int order) { // method: - // Taylor series where there is no pole of one of the Gamma functions - // falls back to Beta function evaluation. Otherwise, fall back to - // Gamma series directly. + // Taylor series where there is no pole of one of the tgamma functions + // falls back to beta function evaluation. Otherwise, fall back to + // tgamma series directly. // FIXME: this could need some testing, maybe it's wrong in some cases? const ex x_pt = x.subs(r); const ex y_pt = y.subs(r); @@ -216,28 +268,28 @@ static ex Beta_series(const ex & x, const ex & y, const relational & r, int orde throw do_taylor(); // caught by function::series() // trap the case where x is on a pole directly: if (x.info(info_flags::integer) && !x.info(info_flags::positive)) - x_ser = Gamma(x+*s).series(r,order); + x_ser = tgamma(x+*s).series(r,order); else - x_ser = Gamma(x).series(r,order); + x_ser = tgamma(x).series(r,order); // trap the case where y is on a pole directly: if (y.info(info_flags::integer) && !y.info(info_flags::positive)) - y_ser = Gamma(y+*s).series(r,order); + y_ser = tgamma(y+*s).series(r,order); else - y_ser = Gamma(y).series(r,order); + y_ser = tgamma(y).series(r,order); // trap the case where y is on a pole directly: if ((x+y).info(info_flags::integer) && !(x+y).info(info_flags::positive)) - xy_ser = Gamma(y+x+*s).series(r,order); + xy_ser = tgamma(y+x+*s).series(r,order); else - xy_ser = Gamma(y+x).series(r,order); + xy_ser = tgamma(y+x).series(r,order); // compose the result: return (x_ser*y_ser/xy_ser).series(r,order); } -REGISTER_FUNCTION(Beta, eval_func(Beta_eval). - evalf_func(Beta_evalf). - derivative_func(Beta_deriv). - series_func(Beta_series)); +REGISTER_FUNCTION(beta, eval_func(beta_eval). + evalf_func(beta_evalf). + derivative_func(beta_deriv). + series_func(beta_series)); ////////// @@ -356,9 +408,9 @@ static ex psi2_eval(const ex & n, const ex & x) // psi(0,x) -> psi(x) if (n.is_zero()) return psi(x); - // psi(-1,x) -> log(Gamma(x)) + // psi(-1,x) -> log(tgamma(x)) if (n.is_equal(_ex_1())) - return log(Gamma(x)); + return log(tgamma(x)); if (n.info(info_flags::numeric) && n.info(info_flags::posint) && x.info(info_flags::numeric)) { numeric nn = ex_to_numeric(n); diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp index 1f4ed18b..37716290 100644 --- a/ginac/matrix.cpp +++ b/ginac/matrix.cpp @@ -444,30 +444,68 @@ matrix matrix::transpose(void) const * @exception logic_error (matrix not square) */ ex matrix::determinant(void) const { - if (row != col) + if (row!=col) throw (std::logic_error("matrix::determinant(): matrix not square")); GINAC_ASSERT(row*col==m.capacity()); + if (this->row==1) // continuation would be pointless + return m[0]; - ex det; bool numeric_flag = true; bool normal_flag = false; + unsigned sparse_count = 0; // count non-zero elements for (exvector::const_iterator r=m.begin(); r!=m.end(); ++r) { - if (!(*r).info(info_flags::numeric)) + if (!(*r).is_zero()) { + ++sparse_count; + } + if (!(*r).info(info_flags::numeric)) { numeric_flag = false; + } if ((*r).info(info_flags::rational_function) && - !(*r).info(info_flags::crational_polynomial)) + !(*r).info(info_flags::crational_polynomial)) { normal_flag = true; + } } if (numeric_flag) - det = determinant_numeric(); - else - if (normal_flag) - det = determinant_minor().normal(); - else - det = determinant_minor(); // is already expanded! + return determinant_numeric(); - return det; + // Now come the minor expansion schemes. We always develop such that the + // smallest minors (i.e, the trivial 1x1 ones) are on the rightmost column. + // For this to be efficient it turns out that the emptiest columns (i.e. + // the ones with most zeros) should be the ones on the right hand side. + // Therefore we presort the columns of the matrix: + typedef pair uintpair; // # of zeros, column + vector c_zeros; // number of zeros in column + for (unsigned c=0; c pre_sort; // unfortunately vector can't be used + // for permutation_sign. + for (vector::iterator i=c_zeros.begin(); i!=c_zeros.end(); ++i) + pre_sort.push_back(i->second); + int sign = permutation_sign(pre_sort); + exvector result(row*col); // represents sorted matrix + unsigned c = 0; + for (vector::iterator i=pre_sort.begin(); + i!=pre_sort.end(); + ++i,++c) { + for (unsigned r=0; rrow*col) { // MAGIC, maybe 10 some day? + if (normal_flag) + return sign*matrix(row,col,result).determinant_minor_sparse().normal(); + return sign*matrix(row,col,result).determinant_minor_sparse(); + } + if (normal_flag) + return sign*matrix(row,col,result).determinant_minor_dense().normal(); + return sign*matrix(row,col,result).determinant_minor_dense(); } @@ -846,6 +884,42 @@ ex matrix::determinant_numeric(void) const *}*/ +ex matrix::determinant_minor_sparse(void) const +{ + // for small matrices the algorithm does not make any sense: + if (this->row==1) + return m[0]; + if (this->row==2) + return (m[0]*m[3]-m[2]*m[1]).expand(); + if (this->row==3) + return (m[0]*m[4]*m[8]-m[0]*m[5]*m[7]- + m[1]*m[3]*m[8]+m[2]*m[3]*m[7]+ + m[1]*m[5]*m[6]-m[2]*m[4]*m[6]).expand(); + + ex det; + matrix minorM(this->row-1,this->col-1); + for (unsigned r1=0; r1row; ++r1) { + // shortcut if element(r1,0) vanishes + if (m[r1*col].is_zero()) + continue; + // assemble the minor matrix + for (unsigned r=0; rrow==1) @@ -926,7 +1000,6 @@ ex matrix::determinant_minor(void) const Pkey.push_back(i); unsigned fc = 0; // controls logic for our strange flipper counter do { - A.insert(Rmap_value(Pkey,_ex0())); det = _ex0(); for (unsigned r=0; rcol-c; ++r) { // maybe there is nothing to do? @@ -946,7 +1019,8 @@ ex matrix::determinant_minor(void) const // prevent build-up of deep nesting of expressions saves time: det = det.expand(); // store the new determinant at its place in B: - B.insert(Rmap_value(Pkey,det)); + if (!det.is_zero()) + B.insert(Rmap_value(Pkey,det)); // increment our strange flipper counter for (fc=this->col-c; fc>0; --fc) { ++Pkey[fc-1]; @@ -966,8 +1040,25 @@ ex matrix::determinant_minor(void) const } +/* Determinant using a simple Bareiss elimination scheme. Suited for + * sparse matrices. + * + * @return the determinant as a new expression (in expanded form) + * @see matrix::determinant() */ +ex matrix::determinant_bareiss(void) const +{ + matrix M(*this); + int sign = M.fraction_free_elimination(); + if (sign) + return sign*M(row-1,col-1); + else + return _ex0(); +} + + /** Determinant built by application of the full permutation group. This - * routine is only called internally by matrix::determinant(). */ + * routine is only called internally by matrix::determinant(). + * NOTE: it is probably inefficient in all cases and may be eliminated. */ ex matrix::determinant_perm(void) const { if (rows()==1) // speed things up @@ -1017,6 +1108,64 @@ int matrix::gauss_elimination(void) } +/** Perform the steps of division free elimination to bring the matrix + * into an upper echelon form. + * + * @return sign is 1 if an even number of rows was swapped, -1 if an odd + * number of rows was swapped and 0 if the matrix is singular. */ +int matrix::division_free_elimination(void) +{ + int sign = 1; + ensure_if_modifiable(); + for (unsigned r1=0; r10) + sign = -sign; + for (unsigned r2=r1+1; r2m[r2*col+c] = this->m[r1*col+r1]*this->m[r2*col+c] - this->m[r2*col+r1]*this->m[r1*col+c]; + for (unsigned c=0; c<=r1; ++c) + this->m[r2*col+c] = _ex0(); + } + } + + return sign; +} + + +/** Perform the steps of Bareiss' one-step fraction free elimination to bring + * the matrix into an upper echelon form. + * + * @return sign is 1 if an even number of rows was swapped, -1 if an odd + * number of rows was swapped and 0 if the matrix is singular. */ +int matrix::fraction_free_elimination(void) +{ + int sign = 1; + ex divisor = 1; + + ensure_if_modifiable(); + for (unsigned r1=0; r10) + sign = -sign; + if (r1>0) + divisor = this->m[(r1-1)*col + (r1-1)]; + for (unsigned r2=r1+1; r2m[r2*col+c] = ((this->m[r1*col+r1]*this->m[r2*col+c] - this->m[r2*col+r1]*this->m[r1*col+c])/divisor).normal(); + for (unsigned c=0; c<=r1; ++c) + this->m[r2*col+c] = _ex0(); + } + } + + return sign; +} + + /** Partial pivoting method. * Usual pivoting (symbolic==false) returns the index to the element with the * largest absolute value in column ro and swaps the current row with the one diff --git a/ginac/matrix.h b/ginac/matrix.h index c0e35dcd..8b0555c8 100644 --- a/ginac/matrix.h +++ b/ginac/matrix.h @@ -96,10 +96,13 @@ public: matrix old_solve(const matrix & v) const; // FIXME: may be removed protected: ex determinant_numeric(void) const; - ex determinant_minor(void) const; + ex determinant_minor_sparse(void) const; + ex determinant_minor_dense(void) const; + ex determinant_bareiss(void) const; ex determinant_perm(void) const; int gauss_elimination(void); int fraction_free_elimination(void); + int division_free_elimination(void); int pivot(unsigned ro, bool symbolic=true); private: // FIXME: these should be obsoleted void ffe_swap(unsigned r1, unsigned c1, unsigned r2 ,unsigned c2); diff --git a/ginac/normal.cpp b/ginac/normal.cpp index fffdc925..a8f64be8 100644 --- a/ginac/normal.cpp +++ b/ginac/normal.cpp @@ -1025,7 +1025,7 @@ static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const { //clog << "heur_gcd(" << a << "," << b << ")\n"; - // Trivial cases + // GCD of two numeric values -> CLN if (is_ex_exactly_of_type(a, numeric) && is_ex_exactly_of_type(b, numeric)) { numeric g = gcd(ex_to_numeric(a), ex_to_numeric(b)); numeric rg; @@ -1057,9 +1057,9 @@ static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const xi = mp * _num2() + _num2(); // 6 tries maximum - for (int t=0; t<6; t++) { - if (xi.int_length() * maxdeg > 100000) { -//clog << "giving up heur_gcd, xi.int_length = " << xi.int_length() << ", maxdeg = " << maxdeg << endl; + for (int t=0; t<6; t++) { // MAGIC + if (xi.int_length() * maxdeg > 100000) { // MAGIC +// clog << "giving up heur_gcd, xi.int_length = " << xi.int_length() << ", maxdeg = " << maxdeg << endl; throw gcdheu_failed(); } diff --git a/ginac/numeric.cpp b/ginac/numeric.cpp index c0649f6f..191dee5d 100644 --- a/ginac/numeric.cpp +++ b/ginac/numeric.cpp @@ -1384,9 +1384,16 @@ const numeric zeta(const numeric & x) /** The Gamma function. * This is only a stub! */ -const numeric Gamma(const numeric & x) +const numeric lgamma(const numeric & x) { - clog << "Gamma(" << x + clog << "lgamma(" << x + << "): Does anybody know good way to calculate this numerically?" + << endl; + return numeric(0); +} +const numeric tgamma(const numeric & x) +{ + clog << "tgamma(" << x << "): Does anybody know good way to calculate this numerically?" << endl; return numeric(0); @@ -1428,7 +1435,7 @@ const numeric factorial(const numeric & n) /** The double factorial combinatorial function. (Scarcely used, but still - * useful in cases, like for exact results of Gamma(n+1/2) for instance.) + * useful in cases, like for exact results of tgamma(n+1/2) for instance.) * * @param n integer argument >= -1 * @return n!! == n * (n-2) * (n-4) * ... * ({1|2}) with 0!! == (-1)!! == 1 @@ -1482,13 +1489,13 @@ const numeric bernoulli(const numeric & nn) return numeric(-1,2); if (nn.is_odd()) return _num0(); - // Until somebody has the Blues and comes up with a much better idea and + // 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 defining formula // B(nn) == - 1/(nn+1) * sum_{k=0}^{nn-1}(binomial(nn+1,k)*B(k)) // whith B(0) == 1. - // Be warned, though: the Bernoulli numbers are probably computationally - // very expensive anyhow and you shouldn't expect miracles to happen. + // Be warned, though: the Bernoulli numbers are computationally very + // expensive anyhow and you shouldn't expect miracles to happen. static vector results; static int highest_result = -1; int n = nn.sub(_num2()).div(_num2()).to_int(); @@ -1731,7 +1738,7 @@ ex PiEvalf(void) } -/** Floating point evaluation of Euler's constant Gamma. */ +/** Floating point evaluation of Euler's constant gamma. */ ex gammaEvalf(void) { return numeric(::cl_eulerconst(cl_default_float_format)); // -> CLN diff --git a/ginac/numeric.h b/ginac/numeric.h index bd05bf66..3306eead 100644 --- a/ginac/numeric.h +++ b/ginac/numeric.h @@ -232,7 +232,8 @@ const numeric asinh(const numeric & x); const numeric acosh(const numeric & x); const numeric atanh(const numeric & x); const numeric zeta(const numeric & x); -const numeric Gamma(const numeric & x); +const numeric lgamma(const numeric & x); +const numeric tgamma(const numeric & x); const numeric psi(const numeric & x); const numeric psi(const numeric & n, const numeric & x); const numeric factorial(const numeric & n); diff --git a/ginsh/ginsh_parser.yy b/ginsh/ginsh_parser.yy index 6931d189..65a41f5a 100644 --- a/ginsh/ginsh_parser.yy +++ b/ginsh/ginsh_parser.yy @@ -758,13 +758,14 @@ int main(int argc, char **argv) insert_fcn_help("atan", "inverse tangent function"); insert_fcn_help("atan2", "inverse tangent function with two arguments"); insert_fcn_help("atanh", "inverse hyperbolic tangent function"); - insert_fcn_help("Beta", "Beta function"); + insert_fcn_help("beta", "Beta function"); insert_fcn_help("binomial", "binomial function"); insert_fcn_help("cos", "cosine function"); insert_fcn_help("cosh", "hyperbolic cosine function"); insert_fcn_help("exp", "exponential function"); insert_fcn_help("factorial", "factorial function"); - insert_fcn_help("Gamma", "Gamma function"); + insert_fcn_help("lgamma", "natural logarithm of Gamma function"); + insert_fcn_help("tgamma", "Gamma function"); insert_fcn_help("log", "natural logarithm"); insert_fcn_help("psi", "psi function\npsi(x) is the digamma function, psi(n,x) the nth polygamma function"); insert_fcn_help("sin", "sine function");