+/* 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;
+}
+