#include "f2c.h" #undef abs #undef min #undef max #include "stdio.h" static integer memfailure = 3; #ifdef KR_headers extern char *malloc(); extern void exit_(); char * F77_aloc(Len, whence) integer Len; char *whence; #else #include "stdlib.h" extern void exit_(integer*); char * F77_aloc(integer Len, char *whence) #endif { char *rv; unsigned int uLen = (unsigned int) Len; /* for K&R C */ if (!(rv = (char*)malloc(uLen))) { fprintf(stderr, "malloc(%u) failure in %s\n", uLen, whence); exit_(&memfailure); } return rv; } static char junk[] = "\n@(#)LIBF77 VERSION 19960619\n"; /* 2.00 11 June 1980. File version.c added to library. 2.01 31 May 1988. s_paus() flushes stderr; names of hl_* fixed [ d]erf[c ] added 8 Aug. 1989: #ifdefs for f2c -i2 added to s_cat.c 29 Nov. 1989: s_cmp returns long (for f2c) 30 Nov. 1989: arg types from f2c.h 12 Dec. 1989: s_rnge allows long names 19 Dec. 1989: getenv_ allows unsorted environment 28 Mar. 1990: add exit(0) to end of main() 2 Oct. 1990: test signal(...) == SIG_IGN rather than & 01 in main 17 Oct. 1990: abort() calls changed to sig_die(...,1) 22 Oct. 1990: separate sig_die from main 25 Apr. 1991: minor, theoretically invisible tweaks to s_cat, sig_die 31 May 1991: make system_ return status 18 Dec. 1991: change long to ftnlen (for -i2) many places 28 Feb. 1992: repair z_sqrt.c (scribbled on input, gave wrong answer) 18 July 1992: for n < 0, repair handling of 0**n in pow_[dr]i.c and m**n in pow_hh.c and pow_ii.c; catch SIGTRAP in main() for error msg before abort 23 July 1992: switch to ANSI prototypes unless KR_headers is #defined 23 Oct. 1992: fix botch in signal_.c (erroneous deref of 2nd arg); change Cabs to f__cabs. 12 March 1993: various tweaks for C++ 2 June 1994: adjust so abnormal terminations invoke f_exit just once 16 Sept. 1994: s_cmp: treat characters as unsigned in comparisons. 19 Sept. 1994: s_paus: flush after end of PAUSE; add -DMSDOS 12 Jan. 1995: pow_[dhiqrz][hiq]: adjust x**i to work on machines that sign-extend right shifts when i is the most negative integer. 26 Jan. 1995: adjust s_cat.c, s_copy.c to permit the left-hand side of character assignments to appear on the right-hand side (unless compiled with -DNO_OVERWRITE). 27 Jan. 1995: minor tweak to s_copy.c: copy forward whenever possible (for better cache behavior). 30 May 1995: added subroutine exit(rc) integer rc. Version not changed. 29 Aug. 1995: add F77_aloc.c; use it in s_cat.c and system_.c. 6 Sept. 1995: fix return type of system_ under -DKR_headers. 19 Dec. 1995: s_cat.c: fix bug when 2nd or later arg overlaps lhs. 19 Mar. 1996: s_cat.c: supply missing break after overlap detection. 13 May 1996: add [lq]bitbits.c and [lq]bitshft.c (f90 bit intrinsics). 19 June 1996: add casts to unsigned in [lq]bitshft.c. */ #include "stdio.h" #include "f2c.h" #ifdef KR_headers extern VOID sig_die(); int abort_() #else extern void sig_die(char*,int); int abort_(void) #endif { sig_die("Fortran abort routine called", 1); #ifdef __cplusplus return 0; #endif } #include "f2c.h" #ifdef KR_headers extern double f__cabs(); double c_abs(z) complex *z; #else extern double f__cabs(double, double); double c_abs(complex *z) #endif { return( f__cabs( z->r, z->i ) ); } #include "f2c.h" #ifdef KR_headers extern double sin(), cos(), sinh(), cosh(); VOID c_cos(r, z) complex *r, *z; #else #undef abs #include "math.h" void c_cos(complex *r, complex *z) #endif { r->r = cos(z->r) * cosh(z->i); r->i = - sin(z->r) * sinh(z->i); } #include "f2c.h" #ifdef KR_headers extern VOID sig_die(); VOID c_div(c, a, b) complex *a, *b, *c; #else extern void sig_die(char*,int); void c_div(complex *c, complex *a, complex *b) #endif { double ratio, den; double abr, abi; if( (abr = b->r) < 0.) abr = - abr; if( (abi = b->i) < 0.) abi = - abi; if( abr <= abi ) { if(abi == 0) sig_die("complex division by zero", 1); ratio = (double)b->r / b->i ; den = b->i * (1 + ratio*ratio); c->r = (a->r*ratio + a->i) / den; c->i = (a->i*ratio - a->r) / den; } else { ratio = (double)b->i / b->r ; den = b->r * (1 + ratio*ratio); c->r = (a->r + a->i*ratio) / den; c->i = (a->i - a->r*ratio) / den; } } #include "f2c.h" #ifdef KR_headers extern double exp(), cos(), sin(); VOID c_exp(r, z) complex *r, *z; #else #undef abs #include "math.h" void c_exp(complex *r, complex *z) #endif { double expx; expx = exp(z->r); r->r = expx * cos(z->i); r->i = expx * sin(z->i); } #include "f2c.h" #ifdef KR_headers extern double log(), f__cabs(), atan2(); VOID c_log(r, z) complex *r, *z; #else #undef abs #include "math.h" extern double f__cabs(double, double); void c_log(complex *r, complex *z) #endif { r->i = atan2(z->i, z->r); r->r = log( f__cabs(z->r, z->i) ); } #include "f2c.h" #ifdef KR_headers extern double sin(), cos(), sinh(), cosh(); VOID c_sin(r, z) complex *r, *z; #else #undef abs #include "math.h" void c_sin(complex *r, complex *z) #endif { r->r = sin(z->r) * cosh(z->i); r->i = cos(z->r) * sinh(z->i); } #include "f2c.h" #ifdef KR_headers extern double sqrt(), f__cabs(); VOID c_sqrt(r, z) complex *r, *z; #else #undef abs #include "math.h" extern double f__cabs(double, double); void c_sqrt(complex *r, complex *z) #endif { double mag, t; if( (mag = f__cabs(z->r, z->i)) == 0.) r->r = r->i = 0.; else if(z->r > 0) { r->r = t = sqrt(0.5 * (mag + z->r) ); t = z->i / t; r->i = 0.5 * t; } else { t = sqrt(0.5 * (mag - z->r) ); if(z->i < 0) t = -t; r->i = t; t = z->i / t; r->r = 0.5 * t; } } #ifdef KR_headers extern double sqrt(); double f__cabs(real, imag) double real, imag; #else #undef abs #include "math.h" double f__cabs(double real, double imag) #endif { double temp; if(real < 0) real = -real; if(imag < 0) imag = -imag; if(imag > real){ temp = real; real = imag; imag = temp; } if((real+imag) == real) return(real); temp = imag/real; temp = real*sqrt(1.0 + temp*temp); /*overflow!!*/ return(temp); } #include "f2c.h" #ifdef KR_headers double d_abs(x) doublereal *x; #else double d_abs(doublereal *x) #endif { if(*x >= 0) return(*x); return(- *x); } #include "f2c.h" #ifdef KR_headers double acos(); double d_acos(x) doublereal *x; #else #undef abs #include "math.h" double d_acos(doublereal *x) #endif { return( acos(*x) ); } #include "f2c.h" #ifdef KR_headers double asin(); double d_asin(x) doublereal *x; #else #undef abs #include "math.h" double d_asin(doublereal *x) #endif { return( asin(*x) ); } #include "f2c.h" #ifdef KR_headers double atan(); double d_atan(x) doublereal *x; #else #undef abs #include "math.h" double d_atan(doublereal *x) #endif { return( atan(*x) ); } #include "f2c.h" #ifdef KR_headers double atan2(); double d_atn2(x,y) doublereal *x, *y; #else #undef abs #include "math.h" double d_atn2(doublereal *x, doublereal *y) #endif { return( atan2(*x,*y) ); } #include "f2c.h" VOID #ifdef KR_headers d_cnjg(r, z) doublecomplex *r, *z; #else d_cnjg(doublecomplex *r, doublecomplex *z) #endif { r->r = z->r; r->i = - z->i; } #include "f2c.h" #ifdef KR_headers double cos(); double d_cos(x) doublereal *x; #else #undef abs #include "math.h" double d_cos(doublereal *x) #endif { return( cos(*x) ); } #include "f2c.h" #ifdef KR_headers double cosh(); double d_cosh(x) doublereal *x; #else #undef abs #include "math.h" double d_cosh(doublereal *x) #endif { return( cosh(*x) ); } #include "f2c.h" #ifdef KR_headers double d_dim(a,b) doublereal *a, *b; #else double d_dim(doublereal *a, doublereal *b) #endif { return( *a > *b ? *a - *b : 0); } #include "f2c.h" #ifdef KR_headers double exp(); double d_exp(x) doublereal *x; #else #undef abs #include "math.h" double d_exp(doublereal *x) #endif { return( exp(*x) ); } #include "f2c.h" #ifdef KR_headers double d_imag(z) doublecomplex *z; #else double d_imag(doublecomplex *z) #endif { return(z->i); } #include "f2c.h" #ifdef KR_headers double floor(); double d_int(x) doublereal *x; #else #undef abs #include "math.h" double d_int(doublereal *x) #endif { return( (*x>0) ? floor(*x) : -floor(- *x) ); } #include "f2c.h" #define log10e 0.43429448190325182765 #ifdef KR_headers double log(); double d_lg10(x) doublereal *x; #else #undef abs #include "math.h" double d_lg10(doublereal *x) #endif { return( log10e * log(*x) ); } #include "f2c.h" #ifdef KR_headers double log(); double d_log(x) doublereal *x; #else #undef abs #include "math.h" double d_log(doublereal *x) #endif { return( log(*x) ); } #include "f2c.h" #ifdef KR_headers #ifdef IEEE_drem double drem(); #else double floor(); #endif double d_mod(x,y) doublereal *x, *y; #else #ifdef IEEE_drem double drem(double, double); #else #undef abs #include "math.h" #endif double d_mod(doublereal *x, doublereal *y) #endif { #ifdef IEEE_drem double xa, ya, z; if ((ya = *y) < 0.) ya = -ya; z = drem(xa = *x, ya); if (xa > 0) { if (z < 0) z += ya; } else if (z > 0) z -= ya; return z; #else double quotient; if( (quotient = *x / *y) >= 0) quotient = floor(quotient); else quotient = -floor(-quotient); return(*x - (*y) * quotient ); #endif } #include "f2c.h" #ifdef KR_headers double floor(); double d_nint(x) doublereal *x; #else #undef abs #include "math.h" double d_nint(doublereal *x) #endif { return( (*x)>=0 ? floor(*x + .5) : -floor(.5 - *x) ); } #include "f2c.h" #ifdef KR_headers double d_prod(x,y) real *x, *y; #else double d_prod(real *x, real *y) #endif { return( (*x) * (*y) ); } #include "f2c.h" #ifdef KR_headers double d_sign(a,b) doublereal *a, *b; #else double d_sign(doublereal *a, doublereal *b) #endif { double x; x = (*a >= 0 ? *a : - *a); return( *b >= 0 ? x : -x); } #include "f2c.h" #ifdef KR_headers double sin(); double d_sin(x) doublereal *x; #else #undef abs #include "math.h" double d_sin(doublereal *x) #endif { return( sin(*x) ); } #include "f2c.h" #ifdef KR_headers double sinh(); double d_sinh(x) doublereal *x; #else #undef abs #include "math.h" double d_sinh(doublereal *x) #endif { return( sinh(*x) ); } #include "f2c.h" #ifdef KR_headers double sqrt(); double d_sqrt(x) doublereal *x; #else #undef abs #include "math.h" double d_sqrt(doublereal *x) #endif { return( sqrt(*x) ); } #include "f2c.h" #ifdef KR_headers double tan(); double d_tan(x) doublereal *x; #else #undef abs #include "math.h" double d_tan(doublereal *x) #endif { return( tan(*x) ); } #include "f2c.h" #ifdef KR_headers double tanh(); double d_tanh(x) doublereal *x; #else #undef abs #include "math.h" double d_tanh(doublereal *x) #endif { return( tanh(*x) ); } #include "f2c.h" #ifdef KR_headers double erf(); double derf_(x) doublereal *x; #else extern double erf(double); double derf_(doublereal *x) #endif { return( erf(*x) ); } #include "f2c.h" #ifdef KR_headers extern double erfc(); double derfc_(x) doublereal *x; #else extern double erfc(double); double derfc_(doublereal *x) #endif { return( erfc(*x) ); } /* EFL support routine to copy string b to string a */ #include "f2c.h" #define M ( (long) (sizeof(long) - 1) ) #define EVEN(x) ( ( (x)+ M) & (~M) ) #ifdef KR_headers extern VOID s_copy(); ef1asc_(a, la, b, lb) ftnint *a, *b; ftnlen *la, *lb; #else extern void s_copy(char*,char*,ftnlen,ftnlen); int ef1asc_(ftnint *a, ftnlen *la, ftnint *b, ftnlen *lb) #endif { s_copy( (char *)a, (char *)b, EVEN(*la), *lb ); #ifdef __cplusplus return 0; #endif } /* EFL support routine to compare two character strings */ #include "f2c.h" #ifdef KR_headers extern integer s_cmp(); integer ef1cmc_(a, la, b, lb) ftnint *a, *b; ftnlen *la, *lb; #else extern integer s_cmp(char*,char*,ftnlen,ftnlen); integer ef1cmc_(ftnint *a, ftnlen *la, ftnint *b, ftnlen *lb) #endif { return( s_cmp( (char *)a, (char *)b, *la, *lb) ); } #include "f2c.h" #ifdef KR_headers double erf(); double erf_(x) real *x; #else extern double erf(double); double erf_(real *x) #endif { return( erf(*x) ); } #include "f2c.h" #ifdef KR_headers double erfc(); double erfc_(x) real *x; #else extern double erfc(double); double erfc_(real *x) #endif { return( erfc(*x) ); } /* This gives the effect of subroutine exit(rc) integer*4 rc stop end * with the added side effect of supplying rc as the program's exit code. */ #include "f2c.h" #undef abs #undef min #undef max #ifndef KR_headers #include "stdlib.h" #ifdef __cplusplus extern "C" { #endif extern void f_exit(void); #endif void #ifdef KR_headers exit_(rc) integer *rc; #else exit_(integer *rc) #endif { #ifdef NO_ONEXIT f_exit(); #endif exit(*rc); } #ifdef __cplusplus } #endif #include "f2c.h" /* * subroutine getarg(k, c, &strlen, max_len) * returns the kth unix command argument in fortran character * variable argument c * * SYNTAX ADJUSTED TO WHAT IS USED IN PMXAB.C!! */ #ifdef KR_headers VOID getarg_(n, s, len, ls) ftnint *n; register char *s; integer *len; ftnlen ls; #else void getarg_(ftnint *n, register char *s, integer *len, ftnlen ls) #endif { extern int xargc; extern char **xargv; register char *t; register int i; if(*n>=0 && *n= 0) (*len)++; /* correct from -1 = zero-length to ordinary counting ... */ for( ; i= 0) return(*x); return(- *x); } #include "f2c.h" #ifdef KR_headers shortint h_dim(a,b) shortint *a, *b; #else shortint h_dim(shortint *a, shortint *b) #endif { return( *a > *b ? *a - *b : 0); } #include "f2c.h" #ifdef KR_headers double floor(); shortint h_dnnt(x) doublereal *x; #else #undef abs #include "math.h" shortint h_dnnt(doublereal *x) #endif { return( (*x)>=0 ? floor(*x + .5) : -floor(.5 - *x) ); } #include "f2c.h" #ifdef KR_headers shortint h_indx(a, b, la, lb) char *a, *b; ftnlen la, lb; #else shortint h_indx(char *a, char *b, ftnlen la, ftnlen lb) #endif { ftnlen i, n; char *s, *t, *bend; n = la - lb + 1; bend = b + lb; for(i = 0 ; i < n ; ++i) { s = a + i; t = b; while(t < bend) if(*s++ != *t++) goto no; return((shortint)i+1); no: ; } return(0); } #include "f2c.h" #ifdef KR_headers shortint h_len(s, n) char *s; ftnlen n; #else shortint h_len(char *s, ftnlen n) #endif { return(n); } #include "f2c.h" #ifdef KR_headers shortint h_mod(a,b) short *a, *b; #else shortint h_mod(short *a, short *b) #endif { return( *a % *b); } #include "f2c.h" #ifdef KR_headers double floor(); shortint h_nint(x) real *x; #else #undef abs #include "math.h" shortint h_nint(real *x) #endif { return( (*x)>=0 ? floor(*x + .5) : -floor(.5 - *x) ); } #include "f2c.h" #ifdef KR_headers shortint h_sign(a,b) shortint *a, *b; #else shortint h_sign(shortint *a, shortint *b) #endif { shortint x; x = (*a >= 0 ? *a : - *a); return( *b >= 0 ? x : -x); } #include "f2c.h" #ifdef KR_headers extern integer s_cmp(); shortlogical hl_ge(a,b,la,lb) char *a, *b; ftnlen la, lb; #else extern integer s_cmp(char *, char *, ftnlen, ftnlen); shortlogical hl_ge(char *a, char *b, ftnlen la, ftnlen lb) #endif { return(s_cmp(a,b,la,lb) >= 0); } #include "f2c.h" #ifdef KR_headers extern integer s_cmp(); shortlogical hl_gt(a,b,la,lb) char *a, *b; ftnlen la, lb; #else extern integer s_cmp(char *, char *, ftnlen, ftnlen); shortlogical hl_gt(char *a, char *b, ftnlen la, ftnlen lb) #endif { return(s_cmp(a,b,la,lb) > 0); } #include "f2c.h" #ifdef KR_headers extern integer s_cmp(); shortlogical hl_le(a,b,la,lb) char *a, *b; ftnlen la, lb; #else extern integer s_cmp(char *, char *, ftnlen, ftnlen); shortlogical hl_le(char *a, char *b, ftnlen la, ftnlen lb) #endif { return(s_cmp(a,b,la,lb) <= 0); } #include "f2c.h" #ifdef KR_headers extern integer s_cmp(); shortlogical hl_lt(a,b,la,lb) char *a, *b; ftnlen la, lb; #else extern integer s_cmp(char *, char *, ftnlen, ftnlen); shortlogical hl_lt(char *a, char *b, ftnlen la, ftnlen lb) #endif { return(s_cmp(a,b,la,lb) < 0); } #include "f2c.h" #ifdef KR_headers integer i_abs(x) integer *x; #else integer i_abs(integer *x) #endif { if(*x >= 0) return(*x); return(- *x); } #include "f2c.h" #ifdef KR_headers integer i_dim(a,b) integer *a, *b; #else integer i_dim(integer *a, integer *b) #endif { return( *a > *b ? *a - *b : 0); } #include "f2c.h" #ifdef KR_headers double floor(); integer i_dnnt(x) doublereal *x; #else #undef abs #include "math.h" integer i_dnnt(doublereal *x) #endif { return( (*x)>=0 ? floor(*x + .5) : -floor(.5 - *x) ); } #include "f2c.h" #ifdef KR_headers integer i_indx(a, b, la, lb) char *a, *b; ftnlen la, lb; #else integer i_indx(char *a, char *b, ftnlen la, ftnlen lb) #endif { ftnlen i, n; char *s, *t, *bend; n = la - lb + 1; bend = b + lb; for(i = 0 ; i < n ; ++i) { s = a + i; t = b; while(t < bend) if(*s++ != *t++) goto no; return(i+1); no: ; } return(0); } #include "f2c.h" #ifdef KR_headers integer i_len(s, n) char *s; ftnlen n; #else integer i_len(char *s, ftnlen n) #endif { return(n); } #include "f2c.h" #ifdef KR_headers integer i_mod(a,b) integer *a, *b; #else integer i_mod(integer *a, integer *b) #endif { return( *a % *b); } #include "f2c.h" #ifdef KR_headers double floor(); integer i_nint(x) real *x; #else #undef abs #include "math.h" integer i_nint(real *x) #endif { return( (*x)>=0 ? floor(*x + .5) : -floor(.5 - *x) ); } #include "f2c.h" #ifdef KR_headers integer i_sign(a,b) integer *a, *b; #else integer i_sign(integer *a, integer *b) #endif { integer x; x = (*a >= 0 ? *a : - *a); return( *b >= 0 ? x : -x); } #include "f2c.h" #ifdef KR_headers ftnint iargc_() #else ftnint iargc_(void) #endif { extern int xargc; return ( xargc - 1 ); } #include "f2c.h" #ifdef KR_headers extern integer s_cmp(); logical l_ge(a,b,la,lb) char *a, *b; ftnlen la, lb; #else extern integer s_cmp(char *, char *, ftnlen, ftnlen); logical l_ge(char *a, char *b, ftnlen la, ftnlen lb) #endif { return(s_cmp(a,b,la,lb) >= 0); } #include "f2c.h" #ifdef KR_headers extern integer s_cmp(); logical l_gt(a,b,la,lb) char *a, *b; ftnlen la, lb; #else extern integer s_cmp(char *, char *, ftnlen, ftnlen); logical l_gt(char *a, char *b, ftnlen la, ftnlen lb) #endif { return(s_cmp(a,b,la,lb) > 0); } #include "f2c.h" #ifdef KR_headers extern integer s_cmp(); logical l_le(a,b,la,lb) char *a, *b; ftnlen la, lb; #else extern integer s_cmp(char *, char *, ftnlen, ftnlen); logical l_le(char *a, char *b, ftnlen la, ftnlen lb) #endif { return(s_cmp(a,b,la,lb) <= 0); } #include "f2c.h" #ifdef KR_headers extern integer s_cmp(); logical l_lt(a,b,la,lb) char *a, *b; ftnlen la, lb; #else extern integer s_cmp(char *, char *, ftnlen, ftnlen); logical l_lt(char *a, char *b, ftnlen la, ftnlen lb) #endif { return(s_cmp(a,b,la,lb) < 0); } #include "f2c.h" #ifndef LONGBITS #define LONGBITS 32 #endif integer #ifdef KR_headers lbit_bits(a, b, len) integer a, b, len; #else lbit_bits(integer a, integer b, integer len) #endif { /* Assume 2's complement arithmetic */ unsigned long x, y; x = (unsigned long) a; y = (unsigned long)-1L; x >>= b; y <<= len; return (integer)(x & ~y); } integer #ifdef KR_headers lbit_cshift(a, b, len) integer a, b, len; #else lbit_cshift(integer a, integer b, integer len) #endif { unsigned long x, y, z; x = (unsigned long)a; if (len <= 0) { if (len == 0) return 0; goto full_len; } if (len >= LONGBITS) { full_len: if (b >= 0) { b %= LONGBITS; return (integer)(x << b | x >> LONGBITS -b ); } b = -b; b %= LONGBITS; return (integer)(x << LONGBITS - b | x >> b); } y = z = (unsigned long)-1; y <<= len; z &= ~y; y &= x; x &= z; if (b >= 0) { b %= len; return (integer)(y | z & (x << b | x >> len - b)); } b = -b; b %= len; return (integer)(y | z & (x >> b | x << len - b)); } #include "f2c.h" integer #ifdef KR_headers lbit_shift(a, b) integer a; integer b; #else lbit_shift(integer a, integer b) #endif { return b >= 0 ? a << b : (integer)((uinteger)a >> -b); } /* STARTUP PROCEDURE FOR UNIX FORTRAN PROGRAMS */ #include "stdio.h" #include "signal.h" #ifndef SIGIOT #ifdef SIGABRT #define SIGIOT SIGABRT #endif #endif #ifndef KR_headers #undef VOID #include "stdlib.h" #endif #ifndef VOID #define VOID void #endif #ifdef __cplusplus extern "C" { #endif #ifdef NO__STDC #define ONEXIT onexit extern VOID f_exit(); #else #ifndef KR_headers extern void f_exit(void); #ifndef NO_ONEXIT #define ONEXIT atexit extern int atexit(void (*)(void)); #endif #else #ifndef NO_ONEXIT #define ONEXIT onexit extern VOID f_exit(); #endif #endif #endif #ifdef KR_headers extern VOID f_init(), sig_die(); extern int MAIN__(); #define Int /* int */ #else extern void f_init(void), sig_die(char*, int); extern int MAIN__(void); #define Int int #endif static VOID sigfdie(Int n) { sig_die("Floating Exception", 1); } static VOID sigidie(Int n) { sig_die("IOT Trap", 1); } #ifdef SIGQUIT static VOID sigqdie(Int n) { sig_die("Quit signal", 1); } #endif static VOID sigindie(Int n) { sig_die("Interrupt", 0); } static VOID sigtdie(Int n) { sig_die("Killed", 0); } #ifdef SIGTRAP static VOID sigtrdie(Int n) { sig_die("Trace trap", 1); } #endif int xargc; char **xargv; #ifdef KR_headers main(argc, argv) int argc; char **argv; #else main(int argc, char **argv) #endif { xargc = argc; xargv = argv; signal(SIGFPE, sigfdie); /* ignore underflow, enable overflow */ #ifdef SIGIOT signal(SIGIOT, sigidie); #endif #ifdef SIGTRAP signal(SIGTRAP, sigtrdie); #endif #ifdef SIGQUIT if(signal(SIGQUIT,sigqdie) == SIG_IGN) signal(SIGQUIT, SIG_IGN); #endif if(signal(SIGINT, sigindie) == SIG_IGN) signal(SIGINT, SIG_IGN); signal(SIGTERM,sigtdie); #ifdef pdp11 ldfps(01200); /* detect overflow as an exception */ #endif f_init(); #ifndef NO_ONEXIT ONEXIT(f_exit); #endif MAIN__(); #ifdef NO_ONEXIT f_exit(); #endif exit(0); /* exit(0) rather than return(0) to bypass Cray bug */ return 0; /* For compilers that complain of missing return values; */ /* others will complain that this is unreachable code. */ } #ifdef __cplusplus } #endif #include "f2c.h" #ifdef KR_headers VOID pow_ci(p, a, b) /* p = a**b */ complex *p, *a; integer *b; #else extern void pow_zi(doublecomplex*, doublecomplex*, integer*); void pow_ci(complex *p, complex *a, integer *b) /* p = a**b */ #endif { doublecomplex p1, a1; a1.r = a->r; a1.i = a->i; pow_zi(&p1, &a1, b); p->r = p1.r; p->i = p1.i; } #include "f2c.h" #ifdef KR_headers double pow(); double pow_dd(ap, bp) doublereal *ap, *bp; #else #undef abs #include "math.h" double pow_dd(doublereal *ap, doublereal *bp) #endif { return(pow(*ap, *bp) ); } #include "f2c.h" #ifdef KR_headers double pow_di(ap, bp) doublereal *ap; integer *bp; #else double pow_di(doublereal *ap, integer *bp) #endif { double pow, x; integer n; unsigned long u; pow = 1; x = *ap; n = *bp; if(n != 0) { if(n < 0) { n = -n; x = 1/x; } for(u = n; ; ) { if(u & 01) pow *= x; if(u >>= 1) x *= x; else break; } } return(pow); } #include "f2c.h" #ifdef KR_headers shortint pow_hh(ap, bp) shortint *ap, *bp; #else shortint pow_hh(shortint *ap, shortint *bp) #endif { shortint pow, x, n; unsigned u; x = *ap; n = *bp; if (n <= 0) { if (n == 0 || x == 1) return 1; if (x != -1) return x == 0 ? 1/x : 0; n = -n; } u = n; for(pow = 1; ; ) { if(u & 01) pow *= x; if(u >>= 1) x *= x; else break; } return(pow); } #include "f2c.h" #ifdef KR_headers integer pow_ii(ap, bp) integer *ap, *bp; #else integer pow_ii(integer *ap, integer *bp) #endif { integer pow, x, n; unsigned long u; x = *ap; n = *bp; if (n <= 0) { if (n == 0 || x == 1) return 1; if (x != -1) return x == 0 ? 1/x : 0; n = -n; } u = n; for(pow = 1; ; ) { if(u & 01) pow *= x; if(u >>= 1) x *= x; else break; } return(pow); } #include "f2c.h" #ifdef KR_headers double pow_ri(ap, bp) real *ap; integer *bp; #else double pow_ri(real *ap, integer *bp) #endif { double pow, x; integer n; unsigned long u; pow = 1; x = *ap; n = *bp; if(n != 0) { if(n < 0) { n = -n; x = 1/x; } for(u = n; ; ) { if(u & 01) pow *= x; if(u >>= 1) x *= x; else break; } } return(pow); } #include "f2c.h" #ifdef KR_headers VOID pow_zi(p, a, b) /* p = a**b */ doublecomplex *p, *a; integer *b; #else extern void z_div(doublecomplex*, doublecomplex*, doublecomplex*); void pow_zi(doublecomplex *p, doublecomplex *a, integer *b) /* p = a**b */ #endif { integer n; unsigned long u; double t; doublecomplex x; static doublecomplex one = {1.0, 0.0}; n = *b; p->r = 1; p->i = 0; if(n == 0) return; if(n < 0) { n = -n; z_div(&x, &one, a); } else { x.r = a->r; x.i = a->i; } for(u = n; ; ) { if(u & 01) { t = p->r * x.r - p->i * x.i; p->i = p->r * x.i + p->i * x.r; p->r = t; } if(u >>= 1) { t = x.r * x.r - x.i * x.i; x.i = 2 * x.r * x.i; x.r = t; } else break; } } #include "f2c.h" #ifdef KR_headers double log(), exp(), cos(), sin(), atan2(), f__cabs(); VOID pow_zz(r,a,b) doublecomplex *r, *a, *b; #else #undef abs #include "math.h" extern double f__cabs(double,double); void pow_zz(doublecomplex *r, doublecomplex *a, doublecomplex *b) #endif { double logr, logi, x, y; logr = log( f__cabs(a->r, a->i) ); logi = atan2(a->i, a->r); x = exp( logr * b->r - logi * b->i ); y = logr * b->i + logi * b->r; r->r = x * cos(y); r->i = x * sin(y); } #include "f2c.h" #ifdef KR_headers double r_abs(x) real *x; #else double r_abs(real *x) #endif { if(*x >= 0) return(*x); return(- *x); } #include "f2c.h" #ifdef KR_headers double acos(); double r_acos(x) real *x; #else #undef abs #include "math.h" double r_acos(real *x) #endif { return( acos(*x) ); } #include "f2c.h" #ifdef KR_headers double asin(); double r_asin(x) real *x; #else #undef abs #include "math.h" double r_asin(real *x) #endif { return( asin(*x) ); } #include "f2c.h" #ifdef KR_headers double atan(); double r_atan(x) real *x; #else #undef abs #include "math.h" double r_atan(real *x) #endif { return( atan(*x) ); } #include "f2c.h" #ifdef KR_headers double atan2(); double r_atn2(x,y) real *x, *y; #else #undef abs #include "math.h" double r_atn2(real *x, real *y) #endif { return( atan2(*x,*y) ); } #include "f2c.h" #ifdef KR_headers VOID r_cnjg(r, z) complex *r, *z; #else VOID r_cnjg(complex *r, complex *z) #endif { r->r = z->r; r->i = - z->i; } #include "f2c.h" #ifdef KR_headers double cos(); double r_cos(x) real *x; #else #undef abs #include "math.h" double r_cos(real *x) #endif { return( cos(*x) ); } #include "f2c.h" #ifdef KR_headers double cosh(); double r_cosh(x) real *x; #else #undef abs #include "math.h" double r_cosh(real *x) #endif { return( cosh(*x) ); } #include "f2c.h" #ifdef KR_headers double r_dim(a,b) real *a, *b; #else double r_dim(real *a, real *b) #endif { return( *a > *b ? *a - *b : 0); } #include "f2c.h" #ifdef KR_headers double exp(); double r_exp(x) real *x; #else #undef abs #include "math.h" double r_exp(real *x) #endif { return( exp(*x) ); } #include "f2c.h" #ifdef KR_headers double r_imag(z) complex *z; #else double r_imag(complex *z) #endif { return(z->i); } #include "f2c.h" #ifdef KR_headers double floor(); double r_int(x) real *x; #else #undef abs #include "math.h" double r_int(real *x) #endif { return( (*x>0) ? floor(*x) : -floor(- *x) ); } #include "f2c.h" #define log10e 0.43429448190325182765 #ifdef KR_headers double log(); double r_lg10(x) real *x; #else #undef abs #include "math.h" double r_lg10(real *x) #endif { return( log10e * log(*x) ); } #include "f2c.h" #ifdef KR_headers double log(); double r_log(x) real *x; #else #undef abs #include "math.h" double r_log(real *x) #endif { return( log(*x) ); } #include "f2c.h" #ifdef KR_headers #ifdef IEEE_drem double drem(); #else double floor(); #endif double r_mod(x,y) real *x, *y; #else #ifdef IEEE_drem double drem(double, double); #else #undef abs #include "math.h" #endif double r_mod(real *x, real *y) #endif { #ifdef IEEE_drem double xa, ya, z; if ((ya = *y) < 0.) ya = -ya; z = drem(xa = *x, ya); if (xa > 0) { if (z < 0) z += ya; } else if (z > 0) z -= ya; return z; #else double quotient; if( (quotient = (double)*x / *y) >= 0) quotient = floor(quotient); else quotient = -floor(-quotient); return(*x - (*y) * quotient ); #endif } #include "f2c.h" #ifdef KR_headers double floor(); double r_nint(x) real *x; #else #undef abs #include "math.h" double r_nint(real *x) #endif { return( (*x)>=0 ? floor(*x + .5) : -floor(.5 - *x) ); } #include "f2c.h" #ifdef KR_headers double r_sign(a,b) real *a, *b; #else double r_sign(real *a, real *b) #endif { double x; x = (*a >= 0 ? *a : - *a); return( *b >= 0 ? x : -x); } #include "f2c.h" #ifdef KR_headers double sin(); double r_sin(x) real *x; #else #undef abs #include "math.h" double r_sin(real *x) #endif { return( sin(*x) ); } #include "f2c.h" #ifdef KR_headers double sinh(); double r_sinh(x) real *x; #else #undef abs #include "math.h" double r_sinh(real *x) #endif { return( sinh(*x) ); } #include "f2c.h" #ifdef KR_headers double sqrt(); double r_sqrt(x) real *x; #else #undef abs #include "math.h" double r_sqrt(real *x) #endif { return( sqrt(*x) ); } #include "f2c.h" #ifdef KR_headers double tan(); double r_tan(x) real *x; #else #undef abs #include "math.h" double r_tan(real *x) #endif { return( tan(*x) ); } #include "f2c.h" #ifdef KR_headers double tanh(); double r_tanh(x) real *x; #else #undef abs #include "math.h" double r_tanh(real *x) #endif { return( tanh(*x) ); } /* Unless compiled with -DNO_OVERWRITE, this variant of s_cat allows the * target of a concatenation to appear on its right-hand side (contrary * to the Fortran 77 Standard, but in accordance with Fortran 90). */ #include "f2c.h" #ifndef NO_OVERWRITE #include "stdio.h" #undef abs #ifdef KR_headers extern char *F77_aloc(); extern void free(); extern void exit_(); #else #include "stdlib.h" extern char *F77_aloc(ftnlen, char*); #endif #include "string.h" #endif /* NO_OVERWRITE */ VOID #ifdef KR_headers s_cat(lp, rpp, rnp, np, ll) char *lp, *rpp[]; ftnlen rnp[], *np, ll; #else s_cat(char *lp, char *rpp[], ftnlen rnp[], ftnlen *np, ftnlen ll) #endif { ftnlen i, nc; char *rp; ftnlen n = *np; #ifndef NO_OVERWRITE ftnlen L, m; char *lp0, *lp1; lp0 = 0; lp1 = lp; L = ll; i = 0; while(i < n) { rp = rpp[i]; m = rnp[i++]; if (rp >= lp1 || rp + m <= lp) { if ((L -= m) <= 0) { n = i; break; } lp1 += m; continue; } lp0 = lp; lp = lp1 = F77_aloc(L = ll, "s_cat"); break; } lp1 = lp; #endif /* NO_OVERWRITE */ for(i = 0 ; i < n ; ++i) { nc = ll; if(rnp[i] < nc) nc = rnp[i]; ll -= nc; rp = rpp[i]; while(--nc >= 0) *lp++ = *rp++; } while(--ll >= 0) *lp++ = ' '; #ifndef NO_OVERWRITE if (lp0) { memcpy(lp0, lp1, L); free(lp1); } #endif } #include "f2c.h" /* compare two strings */ #ifdef KR_headers integer s_cmp(a0, b0, la, lb) char *a0, *b0; ftnlen la, lb; #else integer s_cmp(char *a0, char *b0, ftnlen la, ftnlen lb) #endif { register unsigned char *a, *aend, *b, *bend; a = (unsigned char *)a0; b = (unsigned char *)b0; aend = a + la; bend = b + lb; if(la <= lb) { while(a < aend) if(*a != *b) return( *a - *b ); else { ++a; ++b; } while(b < bend) if(*b != ' ') return( ' ' - *b ); else ++b; } else { while(b < bend) if(*a == *b) { ++a; ++b; } else return( *a - *b ); while(a < aend) if(*a != ' ') return(*a - ' '); else ++a; } return(0); } /* Unless compiled with -DNO_OVERWRITE, this variant of s_copy allows the * target of an assignment to appear on its right-hand side (contrary * to the Fortran 77 Standard, but in accordance with Fortran 90), * as in a(2:5) = a(4:7) . */ #include "f2c.h" /* assign strings: a = b */ #ifdef KR_headers VOID s_copy(a, b, la, lb) register char *a, *b; ftnlen la, lb; #else void s_copy(register char *a, register char *b, ftnlen la, ftnlen lb) #endif { register char *aend, *bend; aend = a + la; if(la <= lb) #ifndef NO_OVERWRITE if (a <= b || a >= b + la) #endif while(a < aend) *a++ = *b++; #ifndef NO_OVERWRITE else for(b += la; a < aend; ) *--aend = *--b; #endif else { bend = b + lb; #ifndef NO_OVERWRITE if (a <= b || a >= bend) #endif while(b < bend) *a++ = *b++; #ifndef NO_OVERWRITE else { a += lb; while(b < bend) *--a = *--bend; a += lb; } #endif while(a < aend) *a++ = ' '; } } #include "stdio.h" #include "f2c.h" #define PAUSESIG 15 #ifdef KR_headers #define Void /* void */ #define Int /* int */ #else #define Void void #define Int int #undef abs #undef min #undef max #include "stdlib.h" #include "signal.h" #ifdef __cplusplus extern "C" { #endif extern int getpid(void), isatty(int), pause(void); #endif extern VOID f_exit(Void); static VOID waitpause(Int n) { n = n; /* shut up compiler warning */ return; } static VOID #ifdef KR_headers s_1paus(fin) FILE *fin; #else s_1paus(FILE *fin) #endif { fprintf(stderr, "To resume execution, type go. Other input will terminate the job.\n"); fflush(stderr); if( getc(fin)!='g' || getc(fin)!='o' || getc(fin)!='\n' ) { fprintf(stderr, "STOP\n"); #ifdef NO_ONEXIT f_exit(); #endif exit(0); } } int #ifdef KR_headers s_paus(s, n) char *s; ftnlen n; #else s_paus(char *s, ftnlen n) #endif { fprintf(stderr, "PAUSE "); if(n > 0) fprintf(stderr, " %.*s", (int)n, s); fprintf(stderr, " statement executed\n"); if( isatty(fileno(stdin)) ) s_1paus(stdin); else { #ifdef MSDOS FILE *fin; fin = fopen("con", "r"); if (!fin) { fprintf(stderr, "s_paus: can't open con!\n"); fflush(stderr); exit(1); } s_1paus(fin); fclose(fin); #else fprintf(stderr, "To resume execution, execute a kill -%d %d command\n", PAUSESIG, getpid() ); signal(PAUSESIG, waitpause); fflush(stderr); pause(); #endif } fprintf(stderr, "Execution resumes after PAUSE.\n"); fflush(stderr); return 0; /* NOT REACHED */ #ifdef __cplusplus } #endif } #include "stdio.h" #include "f2c.h" /* called when a subscript is out of range */ #ifdef KR_headers extern VOID sig_die(); integer s_rnge(varn, offset, procn, line) char *varn, *procn; ftnint offset, line; #else extern VOID sig_die(char*,int); integer s_rnge(char *varn, ftnint offset, char *procn, ftnint line) #endif { #ifdef GIVE_WARNINGS fprintf(stderr,"WARNING: array index out of range\n"); #else printf("."); #endif #ifdef __cplusplus return 0; #endif } #include "stdio.h" #include "f2c.h" #ifdef KR_headers extern void f_exit(); VOID s_stop(s, n) char *s; ftnlen n; #else #undef abs #undef min #undef max #include "stdlib.h" #ifdef __cplusplus extern "C" { #endif void f_exit(void); int s_stop(char *s, ftnlen n) #endif { int i; if(n > 0) { fprintf(stderr, "STOP "); for(i = 0; ir, z->i ) ); } #include "f2c.h" #ifdef KR_headers double sin(), cos(), sinh(), cosh(); VOID z_cos(r, z) doublecomplex *r, *z; #else #undef abs #include "math.h" void z_cos(doublecomplex *r, doublecomplex *z) #endif { r->r = cos(z->r) * cosh(z->i); r->i = - sin(z->r) * sinh(z->i); } #include "f2c.h" #ifdef KR_headers extern void sig_die(); VOID z_div(c, a, b) doublecomplex *a, *b, *c; #else extern void sig_die(char*, int); void z_div(doublecomplex *c, doublecomplex *a, doublecomplex *b) #endif { double ratio, den; double abr, abi; if( (abr = b->r) < 0.) abr = - abr; if( (abi = b->i) < 0.) abi = - abi; if( abr <= abi ) { if(abi == 0) sig_die("complex division by zero", 1); ratio = b->r / b->i ; den = b->i * (1 + ratio*ratio); c->r = (a->r*ratio + a->i) / den; c->i = (a->i*ratio - a->r) / den; } else { ratio = b->i / b->r ; den = b->r * (1 + ratio*ratio); c->r = (a->r + a->i*ratio) / den; c->i = (a->i - a->r*ratio) / den; } } #include "f2c.h" #ifdef KR_headers double exp(), cos(), sin(); VOID z_exp(r, z) doublecomplex *r, *z; #else #undef abs #include "math.h" void z_exp(doublecomplex *r, doublecomplex *z) #endif { double expx; expx = exp(z->r); r->r = expx * cos(z->i); r->i = expx * sin(z->i); } #include "f2c.h" #ifdef KR_headers double log(), f__cabs(), atan2(); VOID z_log(r, z) doublecomplex *r, *z; #else #undef abs #include "math.h" extern double f__cabs(double, double); void z_log(doublecomplex *r, doublecomplex *z) #endif { r->i = atan2(z->i, z->r); r->r = log( f__cabs( z->r, z->i ) ); } #include "f2c.h" #ifdef KR_headers double sin(), cos(), sinh(), cosh(); VOID z_sin(r, z) doublecomplex *r, *z; #else #undef abs #include "math.h" void z_sin(doublecomplex *r, doublecomplex *z) #endif { r->r = sin(z->r) * cosh(z->i); r->i = cos(z->r) * sinh(z->i); } #include "f2c.h" #ifdef KR_headers double sqrt(), f__cabs(); VOID z_sqrt(r, z) doublecomplex *r, *z; #else #undef abs #include "math.h" extern double f__cabs(double, double); void z_sqrt(doublecomplex *r, doublecomplex *z) #endif { double mag; if( (mag = f__cabs(z->r, z->i)) == 0.) r->r = r->i = 0.; else if(z->r > 0) { r->r = sqrt(0.5 * (mag + z->r) ); r->i = z->i / r->r / 2; } else { r->i = sqrt(0.5 * (mag - z->r) ); if(z->i < 0) r->i = - r->i; r->r = z->i / r->i / 2; } }