Provide Gamma function for Apfloat and Apcomplex

[axkr]

Please provide Gamma(z) function for Apfloat and Apcomplex number argument.

See:

• Wikipedia - Gamma function
• Abramowitz and Stegun - Gamma function

[mtommila]

This has been requested a lot. Unfortunately, no "fast" algorithm is known for the calculation of the gamma function. Generally, the algorithms in the apfloat library are "fast" algorithms.

Of course, it's a valid question if no algorithm is better than a "slow" algorithm.

See: Wikipedia - Computational complexity of mathematical operations

[axkr]

What about Spouges approximation.
See http://fredrikj.net/blog

[mtommila]

Well, Spouge's approximation is not a "fast" algorithm either.

[axkr]

Ok what does this mean for the future of the apfloat project?
Will there be any new functions implemented?

For example a port of the mpf_gamma function from the mpmath project starting at line 1571?
https://github.com/fredrik-johansson/mpmath/blob/master/mpmath/libmp/gammazeta.py

[mtommila]

Here's a simple implementation using Spouge's approximation. It's not complete, e.g. it only works with positive input values. I suppose the same algorithm works with an Apcomplex as the argument as well. Also the the definition of the working precision isn't really based on anything.

public static Apfloat gamma(Apfloat x) {
  long precision = x.precision();
  if (precision == Apfloat.INFINITE) {
    throw new InfiniteExpansionException("Cannot calculate gamma function to infinite precision");
  }
  int radix = x.radix();
  if (x.equals(new Apint(1, radix))) {
    return x;
  }
  Apint one = new Apint(1, radix);
  long a1 = (long) (precision / Math.log(2 * Math.PI) * Math.log(radix));
  long workingPrecision = (long) (precision * 1.4); // increase intermediate precision - ck are large and alternating in sign, lots of precision loss there
  x = x.precision(workingPrecision).subtract(one);
  Apint a = new Apint(a1 + 1, radix);
  Apint two = new Apint(2, radix);
  Apfloat c0 = ApfloatMath.sqrt(ApfloatMath.pi(workingPrecision, radix).multiply(two));
  Apfloat sum = c0;
  Apfloat e = ApfloatMath.exp(one.precision(workingPrecision));
  Apfloat divisor = ApfloatMath.exp(new Apfloat(-a1, workingPrecision, radix));
  for (long k = 1; k <= a1; k++) {
    Apint kk = new Apint(k, radix);
    Apfloat ak = a.subtract(kk).precision(workingPrecision);
    Apfloat ck = ApfloatMath.inverseRoot(ak, 2).multiply(ApfloatMath.pow(ak, k)).divide(divisor);
    sum = sum.add(ck.divide(x.add(kk)));
    if (k < a1) {
      divisor = divisor.multiply(e).multiply(kk).negate();
    }
  }
  Aprational half = new Aprational(one, two);
  return ApfloatMath.pow(x.add(a), x.add(half)).multiply(ApfloatMath.exp(x.negate().subtract(a))).multiply(sum).precision(precision);
}

The complexity of this algorithm is something like O(n² log n), which means it's slow.

If I can calculate a result to 4000 digits of precision in half a minute then it means I can calculate:

• 40000 digits in 50 minutes
• 400000 digits in three and half days
• 4000000 digits in a year

Calculating a square root to 4000000 digits takes about one second so compared to that one year is pretty slow.

The mpmath library seems to be using Stirling's approximation but it isn't fundamentally different. I think the asymptotic complexity is the same.

Another consideration for implementation is the caching of the series coefficients (Spouge's or Stirling's approximation), I suppose it would take on the order of O(n²) space so for a million digits of precision the cache would take on the order of a terabyte of storage space.

What comes to new functions to be added to the apfloat library, do you have any other suggestions than the gamma function?

[axkr]

Although high-precision calculations can be slow, I think they are nethertheless useful additions for a limited range of precision.

If the Gamma function is implemented you have a base to create other related functions like: Gamma Functions and Inverses

[mtommila]

I've been working on an implementation of the gamma function. Among other changes, it looks like the next release will require Java 9, though.

[axkr]

I'm using Java 8 in the near future and don't upgrade to Java 9.
Maybe it's possible to use some "Cross-Compilation Options" under Java 9 for the javac compiler?

[mtommila]

Apparently this problem could be solved by building a multi-release jar (JEP 238) so that it both works with Java 8 and also utilizes Java 9 where available.

[axkr]

Thanks for creating Gamma() function.

What about other related functions like: Gamma Functions and Inverses

Should I open a new issue or implement them in my own project?
I think from the "topical focus" they should belong to apfloat.

[mtommila]

Only the plain gamma function is implemented as of now, not the incomplete gamma function.

It would seem like most of the functions listed on that page require at least the incomplete gamma function, they cannot be represented just by the gamma function.

Surely enough, implementing e.g. the generalized hypergeometric function and the Riemann zeta function would allow representing a whole lot of new special functions. This would probably not be easy. Also neither has a known fast algorithm for efficient numeric computation.

Please feel free to either implement your own, or create a new issue, or both.