[GiNaC-devel] Numerical evaluation for tgamma, lgamma and beta.

[Alexei Sheplyakov]

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.

-------------- next part --------------
A non-text attachment was scrubbed...
Name: 0001-lanczos_coeffs-fix-several-implementation-deficiencies.txt.gz
Type: application/octet-stream
Size: 40445 bytes
Desc: not available
Url : http://www.cebix.net/pipermail/ginac-devel/attachments/20061212/6a55ae88/0001-lanczos_coeffs-fix-several-implementation-deficiencies.txt-0001.obj
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 827 bytes
Desc: Digital signature
Url : http://www.cebix.net/pipermail/ginac-devel/attachments/20061212/6a55ae88/attachment-0001.pgp