- As advertised: we are calling the Gamma function tgamma() now!

diff --git a/NEWS b/NEWS
index 699284cc2b10ab6d4874348c3b439232a4ef7505..25d89c6a7ddf4c3727f354014c2a19825dd5e029 100644 (file)
--- 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)
   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
   - 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.
 * #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
 
 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 16788c54d2c658f675d7172739b27e3952a6047e..bdb5e7d7c9cfc60350eb41f855447d4dda3ae651 100644 (file)
--- a/check/exam_inifcns.cpp
+++ b/check/exam_inifcns.cpp
@@ -93,44 +93,44 @@ static unsigned inifcns_consist_trans(void)
     return result;
 }
 
     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;
     
  * 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)
     for (int i=2; i<8; ++i)

-        e += Gamma(ex(i));
+        e += tgamma(ex(i));

     if (e != numeric(874)) {
     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 << " instead of 874" << endl;
         ++result;
     }
     

-    e = Gamma(ex(1));
+    e = tgamma(ex(1));

     for (int i=2; i<8; ++i)
     for (int i=2; i<8; ++i)

-        e *= Gamma(ex(i));    
+        e *= tgamma(ex(i));    

     if (e != numeric(24883200)) {
     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 << " 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) {
     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 << " 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)
     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) {
     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;
     }
              << 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
     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)) {
     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;
     }
              << e << " instead of 2*log(2)" << endl;
         ++result;
     }

diff --git a/check/exam_pseries.cpp b/check/exam_pseries.cpp
index ba7b31933b1f17409dea89604bc553fe05157baf..78a55e2da2d865b4f23b79849b229a80320e9563 100644 (file)
--- 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;
     
     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) +
     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 c92f29f7521309e2c5e09febfcf2794615a1a9ff..0bfbe302db0bb5ce2db89bbaa113162b54bc4dea 100644 (file)
--- a/check/time_gammaseries.cpp
+++ b/check/time_gammaseries.cpp
@@ -22,12 +22,12 @@
 
 #include "times.h"
 
 
 #include "times.h"
 

-unsigned Gammaseries(unsigned order)
+unsigned tgammaseries(unsigned order)

 {
     unsigned result = 0;
     symbol x;
     
 {
     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
     // 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
     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;
     }
              << last_coeff << ", violating a simple estimate." << endl;
         ++result;
     }


@@ -61,7 +61,7 @@ unsigned time_gammaseries(void)
     
     for (vector<unsigned>::iterator i=sizes.begin(); i!=sizes.end(); ++i) {
         omega.start();
     
     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;
     }
         times.push_back(omega.read());
         cout << '.' << flush;
     }

diff --git a/check/time_lw_O.cpp b/check/time_lw_O.cpp
index c6af8c8c2753949049d708944df9e27ef62b5836..d0376ad1efa84232f6a667e1885e94e599b2b50c 100644 (file)
--- 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)) {
     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;
         return 1;
     }
     return 0;

diff --git a/ginac/inifcns.cpp b/ginac/inifcns.cpp
index 3e74166d601ed46fdf451c229baea45d20e7a699..96f4c9a596869ab47886ed8eac5ab61dc693fd17 100644 (file)
--- 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). */
 
 /** 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
 unsigned force_include_zeta1 = function_index_zeta1;
 
 #ifndef NO_NAMESPACE_GINAC

diff --git a/ginac/inifcns.h b/ginac/inifcns.h
index cff2bb2484309200a3ff777665e1c5b6dbb90979..75599ea9bfe5e759dedff6ec631bb819b87bfad6 100644 (file)
--- a/ginac/inifcns.h
+++ b/ginac/inifcns.h
@@ -97,10 +97,11 @@ inline function zeta(const ex & p1, const ex & p2) {
 }
 
 /** Gamma-function. */
 }
 
 /** Gamma-function. */

-DECLARE_FUNCTION_1P(Gamma)
+DECLARE_FUNCTION_1P(lgamma)
+DECLARE_FUNCTION_1P(tgamma)

 
 /** Beta-function. */
 
 /** Beta-function. */

-DECLARE_FUNCTION_2P(Beta)
+DECLARE_FUNCTION_2P(beta)

 
 // overloading at work: we cannot use the macros
 /** Psi-function (aka digamma-function). */
 
 // 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 5b35963f078a8ea975b705bc3dc4326a83a2eedc..9e1b48f4d75bf106a184e5a1b2285dbb2d21ec59 100644 (file)
--- a/ginac/inifcns_gamma.cpp
+++ b/ginac/inifcns_gamma.cpp
@@ -39,40 +39,93 @@ namespace GiNaC {
 #endif // ndef NO_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)
 {
     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)) {
 {
     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 {
             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)
             }
         }
         // 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()));
             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)
                 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());
             }
         }
                 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);
     
 {
     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
 {
     // 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
     // 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()
     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;
     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)
 {
     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)) {
 {
     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)
         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
                 else

-                    throw (std::domain_error("Beta_eval(): simple pole"));
+                    throw (std::domain_error("beta_eval(): simple pole"));

             }
             if (ny.is_negative()) {
                 if (ny<=-nx)
             }
             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
                 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() &&
         }
         // 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:
             !(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;
     
 {
     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)
     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)
     if (deriv_param==1)

-        retval = (psi(y)-psi(x+y))*Beta(x,y);
+        retval = (psi(y)-psi(x+y))*beta(x,y);

     return retval;
 }
 
 
     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:
 {
     // 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);
     // 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))
         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
     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))
     // 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
     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))
     // 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
     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);
 }
 
 
     // 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(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()))
     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);
     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 1f4ed18b0a54e0e29ae2fd69aeb02a666add0c0d..37716290f77317e4916b1705073ba3834332b1ce 100644 (file)
--- 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
 {
  *  @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());
         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;
     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) {
     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;
             numeric_flag = false;

+        }

         if ((*r).info(info_flags::rational_function) &&
         if ((*r).info(info_flags::rational_function) &&

-            !(*r).info(info_flags::crational_polynomial))
+            !(*r).info(info_flags::crational_polynomial)) {

             normal_flag = true;
             normal_flag = true;

+        }

     }
     
     if (numeric_flag)
     }
     
     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();

 }
 
 
 }
 
 


@@ -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; 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")
 /** 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")


@@ -856,7 +930,7 @@ ex matrix::determinant_numeric(void) const
  *
  *  @return the determinant as a new expression (in expanded form)
  *  @see matrix::determinant() */
  *
  *  @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)
 {
     // for small matrices the algorithm does not make any sense:
     if (this->row==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 {
             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?
             det = _ex0();
             for (unsigned r=0; r<this->col-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:
             // 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];
             // 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
 /** 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
 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; 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
 /** 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 c0e35dcdaf5e7c59e648f694f6d3756cd5fb4f74..8b0555c8eb102bb937e538e922f5c755f83efca7 100644 (file)
--- 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;
     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);
     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);
     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 fffdc925b09fe5d4e83040f153ed97d258a12c26..a8f64be82a988a330c6cb70997b60b86ae858af2 100644 (file)
--- 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";
 
 {
 //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;
     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
         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();
                }
 
             throw gcdheu_failed();
                }
 

diff --git a/ginac/numeric.cpp b/ginac/numeric.cpp
index c0649f6f70b550987df6177efcea50f44b0e8f9d..191dee5dafcb1a8448f5abd75167b2944e6648a7 100644 (file)
--- 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! */
 
 /** 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);
          << "): 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
 
 
 /** 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
  *
  *  @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();
         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.
     // 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();
     static vector<numeric> 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
 ex gammaEvalf(void)
 { 
     return numeric(::cl_eulerconst(cl_default_float_format));  // -> CLN

diff --git a/ginac/numeric.h b/ginac/numeric.h
index bd05bf663c12d57b4c8a298f2089f4c0712d8fab..3306eeadf7a0819daa2271b20ee266508edb983c 100644 (file)
--- 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 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);
 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 6931d18944f3d64dbada3089b77723bb17db51fe..65a41f5a33a18917e56c0ea6bae303f463ad66f7 100644 (file)
--- 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("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("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");
        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");