+ if ((_num2()*nx).is_integer()) {
+ // half integer case
+ if (nx.is_equal(_num1_2()))
+ // use psi(n,1/2) == (-)^(n+1) * n! * (2^(n+1)-1) * zeta(n+1)
+ return pow(_num_1(),nn+_num1())*factorial(nn)*(pow(_num2(),nn+_num1()) + _num_1())*zeta(ex(nn+_num1()));
+ if (nx.is_positive()) {
+ numeric m = nx - _num1_2();
+ // use the multiplication formula
+ // psi(n,2*m) == (psi(n,m) + psi(n,m+1/2)) / 2^(n+1)
+ // to revert to positive integer case
+ return psi(n,_num2()*m)*pow(_num2(),nn+_num1())-psi(n,m);
+ } else {
+ // use the recurrence relation
+ // psi(n,-m-1/2) == psi(n,-m-1/2+1) - (-)^n * n! / (-m-1/2)^(n+1)
+ // to relate psi(n,-m-1/2) to psi(n,1/2):
+ // psi(n,-m-1/2) == psi(n,1/2) + r
+ // where r == (-)^(n+1) * n! * ((-1/2)^(-n-1) + ... + (-m-1/2)^(-n-1))
+ numeric recur(0);
+ for (numeric p(nx); p<0; ++p)
+ recur += pow(p, -nn+_num_1());
+ recur *= factorial(nn)*pow(_num_1(), nn+_num_1());
+ return recur+psi(n,_ex1_2());
+ }
+ }
+ // psi2_evalf should be called here once it becomes available