index e51625b2fe2e60b74fd08f9ca588556bd81cdc51..c30ddefec5316266970297022066c4816b52f1cc 100644 (file)
@@ -4,7 +4,7 @@
*  functions. */

/*
*  functions. */

/*
- *  GiNaC Copyright (C) 1999 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2000 Johannes Gutenberg University Mainz, Germany
*
*  This program is free software; you can redistribute it and/or modify
*  it under the terms of the GNU General Public License as published by
*
*  This program is free software; you can redistribute it and/or modify
*  it under the terms of the GNU General Public License as published by
*  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

*  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

-#include <ginac/ginac.h>
+#include "ginac.h"
+
+#ifndef NO_NAMESPACE_GINAC
using namespace GiNaC;
using namespace GiNaC;
+#endif // ndef NO_NAMESPACE_GINAC

-/* Simple tests on the sine trigonometric function. */
+/* Some tests on the sine trigonometric function. */
static unsigned inifcns_consist_sin(void)
{
unsigned result = 0;
static unsigned inifcns_consist_sin(void)
{
unsigned result = 0;
@@ -33,11 +36,12 @@ static unsigned inifcns_consist_sin(void)
// sin(n*Pi) == 0?
errorflag = false;
for (int n=-10; n<=10; ++n) {
// sin(n*Pi) == 0?
errorflag = false;
for (int n=-10; n<=10; ++n) {
-        if (  sin(n*Pi).eval() != numeric(0) ||
-             !sin(n*Pi).eval().info(info_flags::integer) )
+        if (sin(n*Pi).eval() != numeric(0) ||
+            !sin(n*Pi).eval().info(info_flags::integer))
errorflag = true;
}
errorflag = true;
}
-    if ( errorflag ) {
+    if (errorflag) {
+        // we don't count each of those errors
clog << "sin(n*Pi) with integer n does not always return exact 0"
<< endl;
++result;
clog << "sin(n*Pi) with integer n does not always return exact 0"
<< endl;
++result;
@@ -46,17 +50,36 @@ static unsigned inifcns_consist_sin(void)
// sin((n+1/2)*Pi) == {+|-}1?
errorflag = false;
for (int n=-10; n<=10; ++n) {
// sin((n+1/2)*Pi) == {+|-}1?
errorflag = false;
for (int n=-10; n<=10; ++n) {
-        if ( ! sin((n+numeric(1,2))*Pi).eval().info(info_flags::integer) ||
-             !(sin((n+numeric(1,2))*Pi).eval() == numeric(1) ||
-               sin((n+numeric(1,2))*Pi).eval() == numeric(-1)) )
+        if (!sin((n+numeric(1,2))*Pi).eval().info(info_flags::integer) ||
+            !(sin((n+numeric(1,2))*Pi).eval() == numeric(1) ||
+              sin((n+numeric(1,2))*Pi).eval() == numeric(-1)))
errorflag = true;
}
errorflag = true;
}
-    if ( errorflag ) {
+    if (errorflag) {
clog << "sin((n+1/2)*Pi) with integer n does not always return exact {+|-}1"
<< endl;
++result;
}

clog << "sin((n+1/2)*Pi) with integer n does not always return exact {+|-}1"
<< endl;
++result;
}

+    // compare sin((q*Pi).evalf()) with sin(q*Pi).eval().evalf() at various
+    // points.  E.g. if sin(Pi/10) returns something symbolic this should be
+    // equal to sqrt(5)/4-1/4.  This routine will spot programming mistakes
+    // of this kind:
+    errorflag = false;
+    ex argument;
+    numeric epsilon(double(1e-8));
+    for (int n=-240; n<=240; ++n) {
+        argument = n*Pi/60;
+        if (abs(sin(evalf(argument))-evalf(sin(argument)))>epsilon) {
+            clog << "sin(" << argument << ") returns "
+                 << sin(argument) << endl;
+            errorflag = true;
+        }
+    }
+    if (errorflag) {
+        ++result;
+    }
+
return result;
}

return result;
}

@@ -69,11 +92,11 @@ static unsigned inifcns_consist_cos(void)
// cos((n+1/2)*Pi) == 0?
errorflag = false;
for (int n=-10; n<=10; ++n) {
// cos((n+1/2)*Pi) == 0?
errorflag = false;
for (int n=-10; n<=10; ++n) {
-        if (  cos((n+numeric(1,2))*Pi).eval() != numeric(0) ||
-             !cos((n+numeric(1,2))*Pi).eval().info(info_flags::integer) )
+        if (cos((n+numeric(1,2))*Pi).eval() != numeric(0) ||
+            !cos((n+numeric(1,2))*Pi).eval().info(info_flags::integer))
errorflag = true;
}
errorflag = true;
}
-    if ( errorflag ) {
+    if (errorflag) {
clog << "cos((n+1/2)*Pi) with integer n does not always return exact 0"
<< endl;
++result;
clog << "cos((n+1/2)*Pi) with integer n does not always return exact 0"
<< endl;
++result;
@@ -82,17 +105,66 @@ static unsigned inifcns_consist_cos(void)
// cos(n*Pi) == 0?
errorflag = false;
for (int n=-10; n<=10; ++n) {
// cos(n*Pi) == 0?
errorflag = false;
for (int n=-10; n<=10; ++n) {
-        if ( ! cos(n*Pi).eval().info(info_flags::integer) ||
-             !(cos(n*Pi).eval() == numeric(1) ||
-               cos(n*Pi).eval() == numeric(-1)) )
+        if (!cos(n*Pi).eval().info(info_flags::integer) ||
+            !(cos(n*Pi).eval() == numeric(1) ||
+              cos(n*Pi).eval() == numeric(-1)))
errorflag = true;
}
errorflag = true;
}
-    if ( errorflag ) {
+    if (errorflag) {
clog << "cos(n*Pi) with integer n does not always return exact {+|-}1"
<< endl;
++result;
}

clog << "cos(n*Pi) with integer n does not always return exact {+|-}1"
<< endl;
++result;
}

+    // compare cos((q*Pi).evalf()) with cos(q*Pi).eval().evalf() at various
+    // points.  E.g. if cos(Pi/12) returns something symbolic this should be
+    // equal to 1/4*(1+1/3*sqrt(3))*sqrt(6).  This routine will spot
+    // programming mistakes of this kind:
+    errorflag = false;
+    ex argument;
+    numeric epsilon(double(1e-8));
+    for (int n=-240; n<=240; ++n) {
+        argument = n*Pi/60;
+        if (abs(cos(evalf(argument))-evalf(cos(argument)))>epsilon) {
+            clog << "cos(" << argument << ") returns "
+                 << cos(argument) << endl;
+            errorflag = true;
+        }
+    }
+    if (errorflag) {
+        ++result;
+    }
+
+    return result;
+}
+
+/* Simple tests on the tangent trigonometric function. */
+static unsigned inifcns_consist_tan(void)
+{
+    unsigned result = 0;
+    bool errorflag;
+
+    // compare tan((q*Pi).evalf()) with tan(q*Pi).eval().evalf() at various
+    // points.  E.g. if tan(Pi/12) returns something symbolic this should be
+    // equal to 2-sqrt(3).  This routine will spot programming mistakes of
+    // this kind:
+    errorflag = false;
+    ex argument;
+    numeric epsilon(double(1e-8));
+    for (int n=-240; n<=240; ++n) {
+        if (!(n%30) && (n%60))  // skip poles
+            ++n;
+        argument = n*Pi/60;
+        if (abs(tan(evalf(argument))-evalf(tan(argument)))>epsilon) {
+            clog << "tan(" << argument << ") returns "
+                 << tan(argument) << endl;
+            errorflag = true;
+        }
+    }
+    if (errorflag) {
+        ++result;
+    }
+
return result;
}

return result;
}

@@ -166,46 +238,43 @@ static unsigned inifcns_consist_trans(void)
return result;
}

return result;
}

-/* Simple tests on the Gamma combinatorial function.  We stuff in arguments
- * where the result exists in closed form and check if it's ok. */
+/* Simple tests on the Gamma 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));
static unsigned inifcns_consist_gamma(void)
{
unsigned result = 0;
ex e;

e = gamma(ex(1));
-    for (int i=2; i<8; ++i) {
+    for (int i=2; i<8; ++i)
e += gamma(ex(i));
e += gamma(ex(i));
-    }
-    if ( e != numeric(874) ) {
+    if (e != numeric(874)) {
clog << "gamma(1)+...+gamma(7) erroneously returned "
<< e << " instead of 874" << endl;
++result;
}

e = gamma(ex(1));
clog << "gamma(1)+...+gamma(7) erroneously returned "
<< e << " instead of 874" << endl;
++result;
}

e = gamma(ex(1));
-    for (int i=2; i<8; ++i) {
-        e *= gamma(ex(i));
-    }
-    if ( e != numeric(24883200) ) {
+    for (int i=2; i<8; ++i)
+        e *= gamma(ex(i));
+    if (e != numeric(24883200)) {
clog << "gamma(1)*...*gamma(7) erroneously returned "
<< e << " instead of 24883200" << endl;
++result;
}
clog << "gamma(1)*...*gamma(7) erroneously returned "
<< e << " instead of 24883200" << endl;
++result;
}
-
+
e = gamma(ex(numeric(5, 2)))*gamma(ex(numeric(9, 2)))*64;
e = gamma(ex(numeric(5, 2)))*gamma(ex(numeric(9, 2)))*64;
-    if ( e != 315*Pi ) {
+    if (e != 315*Pi) {
clog << "64*gamma(5/2)*gamma(9/2) erroneously returned "
<< e << " instead of 315*Pi" << endl;
++result;
}

e = gamma(ex(numeric(-13, 2)));
clog << "64*gamma(5/2)*gamma(9/2) erroneously returned "
<< e << " instead of 315*Pi" << endl;
++result;
}

e = gamma(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 += gamma(ex(numeric(i, 2)));
-    }
e = (e*gamma(ex(numeric(15, 2)))*numeric(512));
e = (e*gamma(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 "
<< e << " instead of 633935*Pi" << endl;
++result;
clog << "512*(gamma(-13/2)+...+gamma(5/2))*gamma(15/2) erroneously returned "
<< e << " instead of 633935*Pi" << endl;
++result;
@@ -214,6 +283,55 @@ static unsigned inifcns_consist_gamma(void)
return result;
}

return result;
}

+/* Simple tests on the Psi-function (aka polygamma-function).  We stuff in
+   arguments where the result exists in closed form and check if it's ok. */
+static unsigned inifcns_consist_psi(void)
+{
+    unsigned result = 0;
+    symbol x;
+    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));
+    if (e!=2*log(2)) {
+        clog << "gamma(1)'/gamma(1) - gamma(1/2)'/gamma(1/2) erroneously returned "
+             << e << " instead of 2*log(2)" << endl;
+        ++result;
+    }
+
+    return result;
+}
+
+/* Simple tests on the Riemann Zeta function.  We stuff in arguments where the
+ * result exists in closed form and check if it's ok.  Of course, this checks
+ * the Bernoulli numbers as a side effect. */
+static unsigned inifcns_consist_zeta(void)
+{
+    unsigned result = 0;
+    ex e;
+
+    for (int i=0; i<13; i+=2)
+        e += zeta(i)/pow(Pi,i);
+    if (e!=numeric(-204992279,638512875)) {
+        clog << "zeta(0) + zeta(2) + ... + zeta(12) erroneously returned "
+             << e << " instead of -204992279/638512875" << endl;
+        ++result;
+    }
+
+    e = 0;
+    for (int i=-1; i>-16; i--)
+        e += zeta(i);
+    if (e!=numeric(487871,1633632)) {
+        clog << "zeta(-1) + zeta(-2) + ... + zeta(-15) erroneously returned "
+             << e << " instead of 487871/1633632" << endl;
+        ++result;
+    }
+
+    return result;
+}
+
unsigned inifcns_consist(void)
{
unsigned result = 0;
unsigned inifcns_consist(void)
{
unsigned result = 0;
@@ -223,10 +341,13 @@ unsigned inifcns_consist(void)

result += inifcns_consist_sin();
result += inifcns_consist_cos();

result += inifcns_consist_sin();
result += inifcns_consist_cos();
+    result += inifcns_consist_tan();
result += inifcns_consist_trans();
result += inifcns_consist_gamma();
result += inifcns_consist_trans();
result += inifcns_consist_gamma();
+    result += inifcns_consist_psi();
+    result += inifcns_consist_zeta();

-    if ( !result ) {
+    if (!result) {
cout << " passed ";
clog << "(no output)" << endl;
} else {
cout << " passed ";
clog << "(no output)" << endl;
} else {