c++ - Internal error: cg87 1364
- Ali Riggs (159/159) May 14 2010 cgammal1.c:
- Walter Bright (4/8) May 15 2010 Thanks for the report, I've taken the liberty to submit it to the DMC++ ...
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
May 14 2010
Ali Riggs wrote:dmc cgammal1.c -o Internal error: cg87 1364 --- errorlevel 1Thanks 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=57Why DMC does not support conj and conjl?Nobody has asked for them before!
May 15 2010