From:	kunst AT prl DOT philips DOT nl
Subject:	Re: djgpp 1.12 math emulator
To:	s9001200 AT csdvax DOT csd DOT unsw DOT edu DOT au
Date:	Mon, 5 Sep 1994 17:17:32 +0100 (METDST)
Cc:	djgpp AT sun DOT soe DOT clarkson DOT edu (DJGPP users list)

Minxian Yang (s9001200 AT csdvax DOT csd DOT unsw DOT edu DOT au) writes:
> I tried to port a Gamma function routine in Numerical Recipes Fortran to 
> djgpp 1.11. The result was something like Gamma(7.0) = 712.906 (the true
> value should be 720.000). My computer is an NEC 386D20. Though my problem
> is not with djgpp's emulator, I do need some help on the precision of djgpp.

Hello Minxian,

I tried the sample program below (gamma.c) with DJGPP 1.12m1 to check the 
gamma function given in Numerical Recipes. 
Gamma(7) returns 719.999999868 on my PC, but I have a 486DX (which includes 
a coprocessor). I'm sure some of the DJGPP users with a 386SX machine 
(i.e. no coprocessor, just EMU387) can compile the program below to see what
result is obtained with EMU387 (and post the result to you, as well as to the
list).

Regards,

Pieter Kunst (kunst AT prl DOT philips DOT nl)

================================ gamma.c ======================================
/*
 * gamma(7) = 719.999999868
 */

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

double gammln (double xx)   /* from Numerical Recipes */
{
  double x,tmp,ser;
  static double cof[6]={76.18009173,-86.50532033,24.01409822,
	  -1.231739516,0.120858003e-2,-0.536382e-5};
  int j;

  x=xx-1.0;
  tmp=x+5.5;
  tmp -= (x+0.5)*log(tmp);
  ser=1.0;
  for (j=0;j<=5;j++) 
  {
    x += 1.0;
    ser += cof[j]/x;
  }
  return -tmp+log(2.50662827465*ser);
}

double gamma (double x)
{
  return exp(gammln(x));
}

int main (int argc, char *argv[])
{
  if (argc == 1)
  {
    fprintf (stderr, "Usage: gamma x\n");
    exit (1);
  }
  printf ("gamma(%s) = %12.9f\n", argv[1], gamma (atof(argv[1])));
  return 0;
}

- Raw text -