Function digamma_reflection_step 

fn digamma_reflection_step(frac: f64, frac_cot: f64) -> f64

ψ(2 - f) - ψ(1 + f). For f <= 1/2 it is π cot(πf) - 1/f + 1/(1 - f) (DLMF 5.5.2, 5.5.4) from frac_cot = f π cot(πf), with an absolute error of about ε/f that the callers’ factor f removes; otherwise two digamma values.