• Ali Riggs <aliriggs live.com> May 14 2010
• Walter Bright <newshound1 digitalmars.com> May 15 2010

Ali Riggs <aliriggs live.com> writes:

cgammal1.c:
#include <stdio.h>
#include <math.h>
#include <complex.h>

#ifdef LD128BITS
#define NGITER 50.0L
#define NLGITER 50.0L
#define LGMXINT 44.4L
#define GSMALL 1.e-17L
#else
#define NGITER 20.0L
#define NLGITER 16.0L
#define LGMXINT 78.3L
#define GSMALL 1.e-9L
#endif

#define MAXGAM  1755.455L
static long double LOGPIL =  1.1447298858494001741434273513530587116473L;

/* Stirling's formula for the gamma function */
#define NSTIR 18
static long double STIR[NSTIR] = {
  1.50561130400264244123842218771311273E-2L,
  1.79540117061234856107699407722226331E-1L,
-2.48174360026499773091565836874346432E-3L,
-2.95278809456991205054406510546938244E-2L,
  5.40164767892604515180467508570241736E-4L,
  6.40336283380806979482363809026579583E-3L,
-1.62516262783915816898635123980270998E-4L,
-1.91443849856547752650089885832852254E-3L,
  7.20489541602001055908571930225015052E-5L,
  8.39498720672087279993357516764983445E-4L,
-5.17179090826059219337057843002058823E-5L,
-5.92166437353693882864836225604401187E-4L,
  6.97281375836585777429398828575783308E-5L,
  7.84039221720066627474034881442288850E-4L,
-2.29472093621399176954732510288065844E-4L,
-2.68132716049382716049382716049382716E-3L,
  3.47222222222222222222222222222222222E-3L,
  8.33333333333333333333333333333333333E-2L,
};
#define MAXSTIR 1024.0L
static long double SQTPIL = 2.50662827463100050241576528481104525L;

//extern long double MAXLOGL, MAXNUML, PIL;
#define MAXNUML	10E200l
#define PIL		3.1415926535897932384626433832795028841971693993751l

inline long double complex conjl(long double complex z)
{
  return creall(z) - I*cimagl(z);
}

/* Gamma function computed by Stirling's formula.  */

/* static double complex cstirf(x) */
long double complex cstirfl(long double complex x)
{
long double complex y, w;
int i;

w = 1.0L/x;

y = STIR[0];
for (i = 1; i < NSTIR; i++)
   {
     y = y * w + STIR[i];
   }

w = 1.0L + w * y;
#if 1
y = cpowl( x, x - 0.5L ) * cexpl(-x);
#else
y = (x - 0.5L) * clogl(x) - x;
y = cexpl(y);
#endif
y = SQTPIL * y * w;
return( y );
}



long double complex cgammal(long double complex x)
{
long double complex c, u;
long double p, q;
int cj, k;

cj = 0;
if (cimagl(x) < 0.0L)
   {
     cj = 1;
     x = conjl(x);
   }

if( fabsl(creall(x)) > NGITER )
	{
	if( creall(x) < 0.0L )
		{
		q = creall(x);
		p = floorl(q);
		if(( p == q ) && (cimagl(x) == 0.0L))
			{
//			mtherr( "cgammal", OVERFLOW );
			c = MAXNUML + I * MAXNUML;
			goto gamdone;
			}
		/*	c = csinl( PIL * x );*/
		/* Compute sin(pi x)  */
		k = q - 2.0L * floorl (0.5L * q);
		q = PIL * (q - p);
		p = PIL * cimagl(x);
		c = sinl(q) * coshl(p) + cosl(q) * sinhl(p) * I;
		if (k & 1)
		  c = -c;
		c = PIL/(c * cgammal(1.0L - x) );
		goto gamdone;
		}
	else
		{
		  c = cstirfl(x);
		  goto gamdone;
		}
	}
c = 1.0L;
p = 0.0L;
u = x;
while( creall(u) < NGITER )
	{
	if((fabsl (creall(u)) < GSMALL) && (fabsl (cimagl(u)) < GSMALL))
		goto small;
	c *= u;
	p += 1.0L;
	u = x + p;
	}
u = cstirfl(u);
c = u / c;
goto gamdone;


small:
if((creall(x) == 0.0L) && (cimagl(x) == 0.0L))
	{
//	mtherr( "cgammal", SING );
	c = MAXNUML + MAXNUML * I;
	goto gamdone;
	}
else
	c = 1.0L/(((1.0L + 0.57721566490153286060651209008240243L * u) * u)*c);

gamdone:

if (cj)
   c = conjl(c);
return( c );
}

int main(int argc, char *argv[]){
  long double complex w, z = 1.3l + I*0.7l;
  w = cgammal(z);

  printf("Gamma(%lg + i*(%lg)) = %lg + i*(%lg)", (double) creall(z), 
(double) cimagl(z),
                                                 (double) creall(w), 
(double) cimagl(w));
  return 0;
}

dmc cgammal1.c -o
Internal error: cg87 1364
--- errorlevel 1

dmc cgammal1.c & cgammal1
link cgammal1,,,user32+kernel32/noi;

Gamma(1.3 + i*(0.7)) = 0.62366 + i*(0.188284)

This is not correct.

With MinGW gcc the result is correct:
gcc -O3 cgammal1.c -o cgammal1.exe & cgammal1.exe
Gamma(1.3 + i*(0.7)) = 0.692657 + i*(-0.0402559)

Why DMC does not support conj and conjl?

Ali

Walter Bright <newshound1 digitalmars.com> writes:

Ali Riggs wrote:
 dmc cgammal1.c -o
 Internal error: cg87 1364
 --- errorlevel 1


Thanks for the report, I've taken the liberty to submit it to the DMC++ bug
list:

http://bugzilla.digitalmars.com/issues/show_bug.cgi?id=57



 Why DMC does not support conj and conjl?


Nobody has asked for them before!