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
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;
}
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;
}
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) +
#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
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;
}
for (vector<unsigned>::iterator i=sizes.begin(); i!=sizes.end(); ++i) {
omega.start();
- result += Gammaseries(*i);
+ result += tgammaseries(*i);
times.push_back(omega.read());
cout << '.' << flush;
}
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;
/** 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
}
/** 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). */
#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()));
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()
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() &&
!(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);
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));
//////////
// 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);
* @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<unsigned,unsigned> uintpair; // # of zeros, column
+ vector<uintpair> c_zeros; // number of zeros in column
+ for (unsigned c=0; c<col; ++c) {
+ unsigned acc = 0;
+ for (unsigned r=0; r<row; ++r)
+ if (m[r*col+c].is_zero())
+ ++acc;
+ c_zeros.push_back(uintpair(acc,c));
+ }
+ sort(c_zeros.begin(),c_zeros.end());
+ vector<unsigned> pre_sort; // unfortunately vector<uintpair> can't be used
+ // for permutation_sign.
+ for (vector<uintpair>::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<unsigned>::iterator i=pre_sort.begin();
+ i!=pre_sort.end();
+ ++i,++c) {
+ for (unsigned r=0; r<row; ++r)
+ result[r*col+c] = m[r*col+(*i)];
+ }
+
+ if (0*sparse_count>row*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();
}
*}*/
+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; r1<this->row; ++r1) {
+ // shortcut if element(r1,0) vanishes
+ if (m[r1*col].is_zero())
+ continue;
+ // assemble the minor matrix
+ for (unsigned r=0; r<minorM.rows(); ++r) {
+ for (unsigned c=0; c<minorM.cols(); ++c) {
+ if (r<r1)
+ minorM.set(r,c,m[r*col+c+1]);
+ else
+ minorM.set(r,c,m[(r+1)*col+c+1]);
+ }
+ }
+ // recurse down and care for sign:
+ if (r1%2)
+ det -= m[r1*col] * minorM.determinant_minor_sparse();
+ else
+ det += m[r1*col] * minorM.determinant_minor_sparse();
+ }
+ return det.expand();
+}
+
/** Recursive determinant for small matrices having at least one symbolic
* entry. The basic algorithm, known as Laplace-expansion, is enhanced by
* some bookkeeping to avoid calculation of the same submatrices ("minors")
*
* @return the determinant as a new expression (in expanded form)
* @see matrix::determinant() */
-ex matrix::determinant_minor(void) const
+ex matrix::determinant_minor_dense(void) const
{
// for small matrices the algorithm does not make any sense:
if (this->row==1)
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; r<this->col-c; ++r) {
// maybe there is nothing to do?
// 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];
}
+/* 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
}
+/** 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; r1<row-1; ++r1) {
+ int indx = pivot(r1);
+ if (indx==-1)
+ return 0; // Note: leaves *this in a messy state.
+ if (indx>0)
+ sign = -sign;
+ for (unsigned r2=r1+1; r2<row; ++r2) {
+ for (unsigned c=r1+1; c<col; ++c)
+ this->m[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; r1<row-1; ++r1) {
+ int indx = pivot(r1);
+ if (indx==-1)
+ return 0; // Note: leaves *this in a messy state.
+ if (indx>0)
+ sign = -sign;
+ if (r1>0)
+ divisor = this->m[(r1-1)*col + (r1-1)];
+ for (unsigned r2=r1+1; r2<row; ++r2) {
+ for (unsigned c=r1+1; c<col; ++c)
+ this->m[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
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);
{
//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;
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();
}
/** 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);
/** 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
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<numeric> results;
static int highest_result = -1;
int n = nn.sub(_num2()).div(_num2()).to_int();
}
-/** 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
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);
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");