+ // d/dx psi(x) -> psi(1,x)
+ return psi(_ex1(), x);
+}
+
+static ex psi1_series(ex const & x, symbol const & s, ex const & point, int order)
+{
+ // method:
+ // Taylor series where there is no pole falls back to polygamma function
+ // evaluation.
+ // On a pole at -m use the recurrence relation
+ // psi(x) == psi(x+1) - 1/z
+ // from which follows
+ // series(psi(x),x,-m,order) ==
+ // series(psi(x+m+1) - 1/x - 1/(x+1) - 1/(x+m)),x,-m,order);
+ ex xpoint = x.subs(s==point);
+ if (!xpoint.info(info_flags::integer) || xpoint.info(info_flags::positive))
+ throw do_taylor(); // caught by function::series()
+ // if we got here we have to care for a simple pole at -m:
+ numeric m = -ex_to_numeric(xpoint);
+ ex recur;
+ for (numeric p; p<=m; ++p)
+ recur += power(x+p,_ex_1());
+ return (psi(x+m+_ex1())-recur).series(s, point, order);
+}
+
+const unsigned function_index_psi1 = function::register_new("psi", psi1_eval, psi1_evalf, psi1_diff, psi1_series);
+
+//////////
+// Psi-functions (aka polygamma-functions) psi(0,x)==psi(x)
+//////////
+
+static ex psi2_evalf(ex const & n, ex const & x)