Conference rusure::math

i figured this would be answered immediately with full details, but 
here is a little help.  the next term is also negative, and with a
much larger denominator. most of the series i saw were for log(gamma(1+x)).

at least one promised 7 significant digits from 9(i think) terms.

to get the best advice, tell us what accuracy you need, and for what
range of argument. also, do you have constraints on the error?
for instance, min value of max error, or min average error over what range.

do you have a math library, or are you restricted to an algebraic expansion?

> However I can't understand where the 4th term comes from (the one with
> 139 as the numerator), nor can I work out what *should* be the next
> term.
							571
    Me neither, but the next term is     -	---------------------
    						   2488320 * x^4
    
    Doncha just hate it when somebody gives you an infinite series without
    bothering to explain where the coefficients come from?
    
      John

   There is a better way to calculate gamma - that's to calculate
   log(gamma) instead, and then exponentiate to get gamma itself.

   One reason is that the the asymptotic series is more accurate for
   a given number of terms, and only has odd powers of 1/x^k in the
   expansion, and also that the coefficients have a very simple form
   in terms of Bernoulli numbers so you can in principle have as many
   terms as you want.

   The mysterious coefficients in the series you give come from
   exponentiating the series for log(gamma), (known in terms of
   Bernoulli numbers) by formal power series manipulation.  There
   is no simple closed form for them.

   Also, the series is asymptotic - this means that accuracy for a
   given number of terms in the series increases as x increases, but
   does *not* mean that for a given x, you adding more terms to the
   series increases accuracy without limit.  In fact, accuracy gets worse
   after a point in the series!

   A simple "fix" for this last problem is to use the relation

	gamma(z+1) = z*gamma(z)

   and work backwards to get an accurate gamma for a small number.

   For example, suppose you want gamma(4.5).  Since the series is
   more accurate for larger x, instead calculate gamma(10.5) (say)
   and then use

	gamma(10.5) = 9.5*8.5*7.5*6.5*5.5*4.5*gamma(4.5)

   to find gamma(4.5).

   This is now the multiprecision package "MP" calculates gamma to
   arbitrary precision.

   - Jim

>	gamma(10.5) = 9.5*8.5*7.5*6.5*5.5*4.5*gamma(4.5)
>
>   to find gamma(4.5).
    
    An excellent idea! But it raises a precision question: for oddball
    multipliers (like, say, 10 1/3 * 9 1/3 * ��� * 2 1/3) isn't there
    a danger of loss of precision from so many multiplications?
    
      John

>    An excellent idea! But it raises a precision question: for oddball
>    multipliers (like, say, 10 1/3 * 9 1/3 * ��� * 2 1/3) isn't there
>    a danger of loss of precision from so many multiplications?
>    
>      John

   I think it's true that you won't get accuracy to the LSB this way,
   but in practice the error accumulation isn't too serious (compared
   to the relatively low accuracy of Stirling's approximation.)

   In the package MP, you can just add on a few more bits of precision,
   so the error accumulation can be kept under control without too much
   cost.  I don't know how Maple does it, but possibly in a similar way.

   The book Numerical Recipes mentions another way to approximate
   Gamma, invented by Lanczos.  He starts with the basic Stirlings
   approximation, and then instead of an asymptotic series, he makes
   a "correction" series in terms like this

	c1    c2
       --- + --- + ...
       x+1   x+2

   Since the gamma function has simple poles at the negative integers,
   these terms sort of "fit" the effect of these poles (in an ad-hoc way)
   and give a good approximation.  I don't know how the constants are
   chosen since I never read the origional paper.

   - Jim   

   (Simple log(gamma) routine using asymptotic approximation)


#include <math.h>
#include <stdio.h>

/*
 *  Return ln(gamma(x)) using asymptotoc expansion
 */

double gamma_ln_ser(double x, int n)
{
    static double bn[26] = { 1.0, -1.0, 1.0, -1.0, 5.0, -691.0, 7.0, -3617.0,
                            43867.0, -174611.0, 854513.0, -236364091.0,
			    8553103.0, -23749461029.0, 8615841276005.0,
			    -7709321041217.0,
			    2577687858367.0,
			    -26315271553053477373.0,
			    2929993913841559.0,
			    -261082718496449122051.0,
			    1520097643918070802691.0,
			    -27833269579301024235023.0,
			    596451111593912163277961.0,
			    -5609403368997817686249127547.0,
			    495057205241079648212477525.0,
			    -801165718135489957347924991853.0 };

    static double bd[26] = { 6.0, 30.0, 42.0, 30.0, 66.0, 2730.0, 6.0, 510.0,
			     798.0, 330.0, 138.0, 2730.0, 6.0, 870.0, 14322.0,
			     510.0, 6.0, 1919190.0, 6.0, 13530.0, 1806.0,
			     690.0, 282.0, 46410.0, 66.0, 1590.0};
    int
	i;
    double
	s, p, x2 = x*x;

    if (n > 25)
	n = 25;

    for (s = 0.0, p = x, i = n;
	 i >= 0;
	 i--) {
	s = bn[i]/(bd[i]*(2*(i+1)*(2*(i+1)-1)))+s/x2;
    }
    s /= x;

    return (x-0.5)*log(x)-x+0.5*log(2.0*3.14159265358979323846264338)+s;
}

double ask_d(query)
    char *query;
{
    static response[132];
    double answer;

    while(1) {
	printf(query);
	gets(response);
	if (sscanf(response, "%lf", &answer) > 0) return answer;
    }
}

main()
{
    double z, g;
    int i, n;

    z = ask_d("enter z: ");

    for (i = 1; i < 26; i++) {
	g = gamma_ln_ser(z, i);
	printf("%5d %22.18f %22.18f\n",
	       i, z, exp(g));
    }
}