Hello, On Fri, Dec 08, 2006 at 05:51:24PM +0100, Chris Dams wrote:

Dear all,

I added numerical evaluation for the functions tgamma, lgamma and beta to the CVS. The method employed is the Lanczos approximation as described on Wikepedia: http://en.wikipedia.org/wiki/Lanczos_approximation . (I also corrected some small errors on that page.)

This method uses a fixed precision and tables are used to store coefficients. Currently the tables available are for maximally 20, 50, 100 and 200 digits. If Digits exceeds 200, no numerical evaluation is available.

+class lanczos_coeffs +{ + public: + lanczos_coeffs(); + bool sufficiently_accurate(int digits); + int get_order() const { return current_vector->size(); } + numeric calc_lanczos_A(const numeric &) const; + private: + static std::vector<numeric> coeffs_12; // Use in case Digits <= 20 + static std::vector<numeric> coeffs_30; // Use in case Digits <= 50 + static std::vector<numeric> coeffs_60; // Use in case Digits <= 100 + static std::vector<numeric> coeffs_120; // Use in case Digits <= 200 With such a code one need to break ABI (add extra class members) in order to increase the precision. May be replacing all these with static std::vector<std::vector<numeric> > (see attached patch) would be better solution? +lanczos_coeffs::lanczos_coeffs() +{ if (coeffs_12[0] != 0) + return; I think coeffs_12[0] might be not initialized at this stage (and contain random garbage). const numeric lgamma(const numeric &x) { - throw dunno(); + lanczos_coeffs lc; + if (lc.sufficiently_accurate(Digits)) { + numeric pi_val(cln::pi(cln::default_float_format)); What happens here if Digits > default_float_format? + numeric result + = sqrt(numeric(2).mul(pi_val)) + .mul(temp.power(x.add(numeric(-1,2)))) + .mul(exp(temp.mul(-1))) + .mul(A); + All these foo.mul(bar).add(baz) are plain ugly. Any objections against s/numeric/cl_N/g (so it is possible to use natural syntax)? Best regards, Alexei -- All science is either physics or stamp collecting.