/* Output from p2c, the Pascal-to-C translator */ /* From input file "meteorJ.p" */ /* to simplify portability, the files p2c.h and p2clib.c are included here explicitly, replacing #include */ /* the following is p2c.h, copied here for convenience */ #ifndef P2C_H #define P2C_H /* Header file for code generated by "p2c", the Pascal-to-C translator */ /* "p2c" Copyright (C) 1989 Dave Gillespie, version 1.16. * This file may be copied, modified, etc. in any way. It is not restricted * by the licence agreement accompanying p2c itself. */ #include /* If the following heuristic fails, compile -DBSD=0 for non-BSD systems, or -DBSD=1 for BSD systems. */ #ifdef M_XENIX # define BSD 0 #endif #ifdef FILE /* a #define in BSD, a typedef in SYSV (hp-ux, at least) */ # ifndef BSD /* (a convenient, but horrible kludge!) */ # define BSD 1 # endif #endif #ifdef BSD # if !BSD # undef BSD # endif #endif #ifdef __STDC__ # include # include # define HAS_STDLIB # define __CAT__(a,b)a##b #else # ifndef BSD # include # endif # include # define __ID__(a)a # define __CAT__(a,b)__ID__(a)b #endif #ifdef BSD # include # define memcpy(a,b,n) (bcopy(b,a,n),a) # define memcmp(a,b,n) bcmp(a,b,n) # define strchr(s,c) index(s,c) # define strrchr(s,c) rindex(s,c) #else # include #endif #include #include #include #include typedef struct __p2c_jmp_buf { struct __p2c_jmp_buf *next; jmp_buf jbuf; } __p2c_jmp_buf; /* Warning: The following will not work if setjmp is used simultaneously. This also violates the ANSI restriction about using vars after longjmp, but a typical implementation of longjmp will get it right anyway. */ #ifndef FAKE_TRY # define TRY(x) do { __p2c_jmp_buf __try_jb; \ __try_jb.next = __top_jb; \ if (!setjmp((__top_jb = &__try_jb)->jbuf)) { # define RECOVER(x) __top_jb = __try_jb.next; } else { # define RECOVER2(x,L) __top_jb = __try_jb.next; } else { \ if (0) { L: __top_jb = __try_jb.next; } # define ENDTRY(x) } } while (0) #else # define TRY(x) if (1) { # define RECOVER(x) } else do { # define RECOVER2(x,L) } else do { L: ; # define ENDTRY(x) } while (0) #endif #ifdef M_XENIX /* avoid compiler bug */ # define SHORT_MAX (32767) # define SHORT_MIN (-32768) #endif /* The following definitions work only on twos-complement machines */ #ifndef SHORT_MAX # define SHORT_MAX (((unsigned short) -1) >> 1) # define SHORT_MIN (~SHORT_MAX) #endif #ifndef INT_MAX # define INT_MAX (((unsigned int) -1) >> 1) # define INT_MIN (~INT_MAX) #endif #ifndef LONG_MAX # define LONG_MAX (((unsigned long) -1) >> 1) # define LONG_MIN (~LONG_MAX) #endif #ifndef SEEK_SET # define SEEK_SET 0 # define SEEK_CUR 1 # define SEEK_END 2 #endif #ifndef EXIT_SUCCESS # define EXIT_SUCCESS 0 # define EXIT_FAILURE 1 #endif #define SETBITS 32 #ifdef __STDC__ # define Signed signed # define Void void /* Void f() = procedure */ # ifndef Const # define Const const # endif # ifndef Volatile # define Volatile volatile # endif # define PP(x) x /* function prototype */ # define PV() (void) /* null function prototype */ typedef void *Anyptr; #else # define Signed # define Void void # ifndef Const # define Const # endif # ifndef Volatile # define Volatile # endif # define PP(x) () # define PV() () typedef char *Anyptr; #endif #ifdef __GNUC__ # define Inline inline #else # define Inline #endif #define Register register /* Register variables */ #define Char char /* Characters (not bytes) */ #ifndef Static # define Static static /* Private global funcs and vars */ #endif #ifndef Local # define Local static /* Nested functions */ #endif typedef Signed char schar; typedef unsigned char uchar; typedef unsigned char boolean; #ifndef true # define true 1 # define false 0 #endif typedef struct { Anyptr proc, link; } _PROCEDURE; #ifndef _FNSIZE # define _FNSIZE 120 #endif extern Void PASCAL_MAIN PP( (int, Char **) ); extern Char **P_argv; extern int P_argc; extern short P_escapecode; extern int P_ioresult; extern __p2c_jmp_buf *__top_jb; #ifdef P2C_H_PROTO /* if you have Ansi C but non-prototyped header files */ extern Char *strcat PP( (Char *, Const Char *) ); extern Char *strchr PP( (Const Char *, int) ); extern int strcmp PP( (Const Char *, Const Char *) ); extern Char *strcpy PP( (Char *, Const Char *) ); extern size_t strlen PP( (Const Char *) ); extern Char *strncat PP( (Char *, Const Char *, size_t) ); extern int strncmp PP( (Const Char *, Const Char *, size_t) ); extern Char *strncpy PP( (Char *, Const Char *, size_t) ); extern Char *strrchr PP( (Const Char *, int) ); extern Anyptr memchr PP( (Const Anyptr, int, size_t) ); extern Anyptr memmove PP( (Anyptr, Const Anyptr, size_t) ); extern Anyptr memset PP( (Anyptr, int, size_t) ); #ifndef memcpy extern Anyptr memcpy PP( (Anyptr, Const Anyptr, size_t) ); extern int memcmp PP( (Const Anyptr, Const Anyptr, size_t) ); #endif extern int atoi PP( (Const Char *) ); extern double atof PP( (Const Char *) ); extern long atol PP( (Const Char *) ); extern double strtod PP( (Const Char *, Char **) ); extern long strtol PP( (Const Char *, Char **, int) ); #endif /*P2C_H_PROTO*/ #ifndef HAS_STDLIB extern Anyptr malloc PP( (size_t) ); extern Void free PP( (Anyptr) ); #endif extern int _OutMem PV(); extern int _CaseCheck PV(); extern int _NilCheck PV(); extern int _Escape PP( (int) ); extern int _EscIO PP( (int) ); extern long ipow PP( (long, long) ); extern Char *strsub PP( (Char *, Char *, int, int) ); extern Char *strltrim PP( (Char *) ); extern Char *strrtrim PP( (Char *) ); extern Char *strrpt PP( (Char *, Char *, int) ); extern Char *strpad PP( (Char *, Char *, int, int) ); extern int strpos2 PP( (Char *, Char *, int) ); extern long memavail PV(); extern int P_peek PP( (FILE *) ); extern int P_eof PP( (FILE *) ); extern int P_eoln PP( (FILE *) ); extern Void P_readpaoc PP( (FILE *, Char *, int) ); extern Void P_readlnpaoc PP( (FILE *, Char *, int) ); extern long P_maxpos PP( (FILE *) ); extern Char *P_trimname PP( (Char *, int) ); extern long *P_setunion PP( (long *, long *, long *) ); extern long *P_setint PP( (long *, long *, long *) ); extern long *P_setdiff PP( (long *, long *, long *) ); extern long *P_setxor PP( (long *, long *, long *) ); extern int P_inset PP( (unsigned, long *) ); extern int P_setequal PP( (long *, long *) ); extern int P_subset PP( (long *, long *) ); extern long *P_addset PP( (long *, unsigned) ); extern long *P_addsetr PP( (long *, unsigned, unsigned) ); extern long *P_remset PP( (long *, unsigned) ); extern long *P_setcpy PP( (long *, long *) ); extern long *P_expset PP( (long *, long) ); extern long P_packset PP( (long *) ); extern int P_getcmdline PP( (int l, int h, Char *line) ); extern Void TimeStamp PP( (int *Day, int *Month, int *Year, int *Hour, int *Min, int *Sec) ); extern Void P_sun_argv PP( (char *, int, int) ); /* I/O error handling */ #define _CHKIO(cond,ior,val,def) ((cond) ? P_ioresult=0,(val) \ : P_ioresult=(ior),(def)) #define _SETIO(cond,ior) (P_ioresult = (cond) ? 0 : (ior)) /* Following defines are suitable for the HP Pascal operating system */ #define FileNotFound 10 #define FileNotOpen 13 #define FileWriteError 38 #define BadInputFormat 14 #define EndOfFile 30 /* Creating temporary files */ #if (defined(BSD) || defined(NO_TMPFILE)) && !defined(HAVE_TMPFILE) # define tmpfile() (fopen(tmpnam(NULL), "w+")) #endif /* File buffers */ #define FILEBUF(f,sc,type) sc int __CAT__(f,_BFLAGS); \ sc type __CAT__(f,_BUFFER) #define RESETBUF(f,type) (__CAT__(f,_BFLAGS) = 1) #define SETUPBUF(f,type) (__CAT__(f,_BFLAGS) = 0) #define GETFBUF(f,type) (*((__CAT__(f,_BFLAGS) == 1 && \ ((__CAT__(f,_BFLAGS) = 2), \ fread(&__CAT__(f,_BUFFER), \ sizeof(type),1,(f)))),\ &__CAT__(f,_BUFFER))) #define AGETFBUF(f,type) ((__CAT__(f,_BFLAGS) == 1 && \ ((__CAT__(f,_BFLAGS) = 2), \ fread(&__CAT__(f,_BUFFER), \ sizeof(type),1,(f)))),\ __CAT__(f,_BUFFER)) #define PUTFBUF(f,type,v) (GETFBUF(f,type) = (v)) #define CPUTFBUF(f,v) (PUTFBUF(f,char,v)) #define APUTFBUF(f,type,v) (memcpy(GETFBUF(f,type), (v), \ sizeof(__CAT__(f,_BUFFER)))) #define GET(f,type) (__CAT__(f,_BFLAGS) == 1 ? \ fread(&__CAT__(f,_BUFFER),sizeof(type),1,(f)) : \ (__CAT__(f,_BFLAGS) = 1)) #define PUT(f,type) (fwrite(&__CAT__(f,_BUFFER),sizeof(type),1,(f)), \ (__CAT__(f,_BFLAGS) = 0)) #define CPUT(f) (PUT(f,char)) #define BUFEOF(f) (__CAT__(f,_BFLAGS) != 2 && P_eof(f)) #define BUFFPOS(f) (ftell(f) - (__CAT__(f,_BFLAGS) == 2)) /* Memory allocation */ #ifdef __GCC__ # define Malloc(n) (malloc(n) ?: (Anyptr)_OutMem()) #else extern Anyptr __MallocTemp__; # define Malloc(n) ((__MallocTemp__ = malloc(n)) ? __MallocTemp__ : (Anyptr)_OutMem()) #endif #define FreeR(p) (free((Anyptr)(p))) /* used if arg is an rvalue */ #define Free(p) (free((Anyptr)(p)), (p)=NULL) /* sign extension */ #define SEXT(x,n) ((x) | -(((x) & (1L<<((n)-1))) << 1)) /* packed arrays */ /* BEWARE: these are untested! */ #define P_getbits_UB(a,i,n,L) ((int)((a)[(i)>>(L)-(n)] >> \ (((~(i))&((1<<(L)-(n))-1)) << (n)) & \ (1<<(1<<(n)))-1)) #define P_getbits_SB(a,i,n,L) ((int)((a)[(i)>>(L)-(n)] << \ (16 - ((((~(i))&((1<<(L)-(n))-1))+1) <<\ (n)) >> (16-(1<<(n)))))) #define P_putbits_UB(a,i,x,n,L) ((a)[(i)>>(L)-(n)] |= \ (x) << (((~(i))&((1<<(L)-(n))-1)) << (n))) #define P_putbits_SB(a,i,x,n,L) ((a)[(i)>>(L)-(n)] |= \ ((x) & (1<<(1<<(n)))-1) << \ (((~(i))&((1<<(L)-(n))-1)) << (n))) #define P_clrbits_B(a,i,n,L) ((a)[(i)>>(L)-(n)] &= \ ~( ((1<<(1<<(n)))-1) << \ (((~(i))&((1<<(L)-(n))-1)) << (n))) ) /* small packed arrays */ #define P_getbits_US(v,i,n) ((int)((v) >> ((i)<<(n)) & (1<<(1<<(n)))-1)) #define P_getbits_SS(v,i,n) ((int)((long)(v) << (SETBITS - (((i)+1) << (n))) >> (SETBITS-(1<<(n))))) #define P_putbits_US(v,i,x,n) ((v) |= (x) << ((i) << (n))) #define P_putbits_SS(v,i,x,n) ((v) |= ((x) & (1<<(1<<(n)))-1) << ((i)<<(n))) #define P_clrbits_S(v,i,n) ((v) &= ~( ((1<<(1<<(n)))-1) << ((i)<<(n)) )) #define P_max(a,b) ((a) > (b) ? (a) : (b)) #define P_min(a,b) ((a) < (b) ? (a) : (b)) /* Fix toupper/tolower on Suns and other stupid BSD systems */ #ifdef toupper # undef toupper # undef tolower # define toupper(c) my_toupper(c) # define tolower(c) my_tolower(c) #endif #ifndef _toupper # if 'A' == 65 && 'a' == 97 # define _toupper(c) ((c)-'a'+'A') # define _tolower(c) ((c)-'A'+'a') # else # define _toupper(c) toupper(c) # define _tolower(c) tolower(c) # endif #endif #endif /* P2C_H */ /* End. */ /* end of p2c.h */ /* the following is p2clib.c */ /* Run-time library for use with "p2c", the Pascal to C translator */ /* "p2c" Copyright (C) 1989 Dave Gillespie. * This file may be copied, modified, etc. in any way. It is not restricted * by the licence agreement accompanying p2c itself. */ /* #define LACK_LABS */ /* Define these if necessary */ /* #define LACK_MEMMOVE */ #ifndef NO_TIME # include #endif #define Isspace(c) isspace(c) /* or "((c) == ' ')" if preferred */ int P_argc; char **P_argv; short P_escapecode; int P_ioresult; long EXCP_LINE; /* Used by Pascal workstation system */ Anyptr __MallocTemp__; __p2c_jmp_buf *__top_jb; void PASCAL_MAIN(argc, argv) int argc; char **argv; { P_argc = argc; P_argv = argv; __top_jb = NULL; #ifdef LOCAL_INIT LOCAL_INIT(); #endif } /* In case your system lacks these... */ #ifdef LACK_LABS long labs(x) long x; { return((x > 0) ? x : -x); } #endif #ifdef LACK_MEMMOVE Anyptr memmove(d, s, n) Anyptr d, s; register long n; { if (d < s || d - s >= n) { memcpy(d, s, n); return d; } else if (n > 0) { register char *dd = d + n, *ss = s + n; while (--n >= 0) *--dd = *--ss; } return d; } #endif int my_toupper(c) int c; { if (islower(c)) return _toupper(c); else return c; } int my_tolower(c) int c; { if (isupper(c)) return _tolower(c); else return c; } long ipow(a, b) long a, b; { long v; if (a == 0 || a == 1) return a; if (a == -1) return (b & 1) ? -1 : 1; if (b < 0) return 0; if (a == 2) return 1 << b; v = (b & 1) ? a : 1; while ((b >>= 1) > 0) { a *= a; if (b & 1) v *= a; } return v; } /* Common string functions: */ /* Store in "ret" the substring of length "len" starting from "pos" (1-based). Store a shorter or null string if out-of-range. Return "ret". */ char *strsub(ret, s, pos, len) register char *ret, *s; register int pos, len; { register char *s2; if (--pos < 0 || len <= 0) { *ret = 0; return ret; } while (pos > 0) { if (!*s++) { *ret = 0; return ret; } pos--; } s2 = ret; while (--len >= 0) { if (!(*s2++ = *s++)) return ret; } *s2 = 0; return ret; } /* Return the index of the first occurrence of "pat" as a substring of "s", starting at index "pos" (1-based). Result is 1-based, 0 if not found. */ int strpos2(s, pat, pos) char *s; register char *pat; register int pos; { register char *cp, ch; register int slen; if (--pos < 0) return 0; slen = strlen(s) - pos; cp = s + pos; if (!(ch = *pat++)) return 0; pos = strlen(pat); slen -= pos; while (--slen >= 0) { if (*cp++ == ch && !strncmp(cp, pat, pos)) return cp - s; } return 0; } /* Case-insensitive version of strcmp. */ int strcicmp(s1, s2) register char *s1, *s2; { register unsigned char c1, c2; while (*s1) { if (*s1++ != *s2++) { if (!s2[-1]) return 1; c1 = toupper(s1[-1]); c2 = toupper(s2[-1]); if (c1 != c2) return c1 - c2; } } if (*s2) return -1; return 0; } /* HP and Turbo Pascal string functions: */ /* Trim blanks at left end of string. */ char *strltrim(s) register char *s; { while (Isspace(*s++)) ; return s - 1; } /* Trim blanks at right end of string. */ char *strrtrim(s) register char *s; { register char *s2 = s; if (!*s) return s; while (*++s2) ; while (s2 > s && Isspace(*--s2)) *s2 = 0; return s; } /* Store in "ret" "num" copies of string "s". Return "ret". */ char *strrpt(ret, s, num) char *ret; register char *s; register int num; { register char *s2 = ret; register char *s1; while (--num >= 0) { s1 = s; while ((*s2++ = *s1++)) ; s2--; } return ret; } /* Store in "ret" string "s" with enough pad chars added to reach "size". */ char *strpad(ret, s, padchar, num) char *ret; register char *s; register int padchar, num; { register char *d = ret; if (s == d) { while (*d++) ; } else { while ((*d++ = *s++)) ; } num -= (--d - ret); while (--num >= 0) *d++ = padchar; *d = 0; return ret; } /* Copy the substring of length "len" from index "spos" of "s" (1-based) to index "dpos" of "d", lengthening "d" if necessary. Length and indices must be in-range. */ void strmove(len, s, spos, d, dpos) register char *s, *d; register int len, spos, dpos; { s += spos - 1; d += dpos - 1; while (*d && --len >= 0) *d++ = *s++; if (len > 0) { while (--len >= 0) *d++ = *s++; *d = 0; } } /* Delete the substring of length "len" at index "pos" from "s". Delete less if out-of-range. */ void strdelete(s, pos, len) register char *s; register int pos, len; { register int slen; if (--pos < 0) return; slen = strlen(s) - pos; if (slen <= 0) return; s += pos; if (slen <= len) { *s = 0; return; } while ((*s = s[len])) s++; } /* Insert string "src" at index "pos" of "dst". */ void strinsert(src, dst, pos) register char *src, *dst; register int pos; { register int slen, dlen; if (--pos < 0) return; dlen = strlen(dst); dst += dlen; dlen -= pos; if (dlen <= 0) { strcpy(dst, src); return; } slen = strlen(src); do { dst[slen] = *dst; --dst; } while (--dlen >= 0); dst++; while (--slen >= 0) *dst++ = *src++; } /* File functions */ /* Peek at next character of input stream; return EOF at end-of-file. */ int P_peek(f) FILE *f; { int ch; ch = getc(f); if (ch == EOF) return EOF; ungetc(ch, f); return (ch == '\n') ? ' ' : ch; } /* Check if at end of file, using Pascal "eof" semantics. End-of-file for stdin is broken; remove the special case for it to be broken in a different way. */ int P_eof(f) FILE *f; { register int ch; if (feof(f)) return 1; if (f == stdin) return 0; /* not safe to look-ahead on the keyboard! */ ch = getc(f); if (ch == EOF) return 1; ungetc(ch, f); return 0; } /* Check if at end of line (or end of entire file). */ int P_eoln(f) FILE *f; { register int ch; ch = getc(f); if (ch == EOF) return 1; ungetc(ch, f); return (ch == '\n'); } /* Read a packed array of characters from a file. */ Void P_readpaoc(f, s, len) FILE *f; char *s; int len; { int ch; for (;;) { if (len <= 0) return; ch = getc(f); if (ch == EOF || ch == '\n') break; *s++ = ch; --len; } while (--len >= 0) *s++ = ' '; if (ch != EOF) ungetc(ch, f); } Void P_readlnpaoc(f, s, len) FILE *f; char *s; int len; { int ch; for (;;) { ch = getc(f); if (ch == EOF || ch == '\n') break; if (len > 0) { *s++ = ch; --len; } } while (--len >= 0) *s++ = ' '; } /* Compute maximum legal "seek" index in file (0-based). */ long P_maxpos(f) FILE *f; { long savepos = ftell(f); long val; if (fseek(f, 0L, SEEK_END)) return -1; val = ftell(f); if (fseek(f, savepos, SEEK_SET)) return -1; return val; } /* Use packed array of char for a file name. */ Char *P_trimname(fn, len) register Char *fn; register int len; { static Char fnbuf[256]; register Char *cp = fnbuf; while (--len >= 0 && *fn && !isspace(*fn)) *cp++ = *fn++; return fnbuf; } /* Pascal's "memavail" doesn't make much sense in Unix with virtual memory. We fix memory size as 10Meg as a reasonable compromise. */ long memavail() { return 10000000; /* worry about this later! */ } long maxavail() { return memavail(); } /* Sets are stored as an array of longs. S[0] is the size of the set; S[N] is the N'th 32-bit chunk of the set. S[0] equals the maximum I such that S[I] is nonzero. S[0] is zero for an empty set. Within each long, bits are packed from lsb to msb. The first bit of the set is the element with ordinal value 0. (Thus, for a "set of 5..99", the lowest five bits of the first long are unused and always zero.) */ /* (Sets with 32 or fewer elements are normally stored as plain longs.) */ long *P_setunion(d, s1, s2) /* d := s1 + s2 */ register long *d, *s1, *s2; { long *dbase = d++; register int sz1 = *s1++, sz2 = *s2++; while (sz1 > 0 && sz2 > 0) { *d++ = *s1++ | *s2++; sz1--, sz2--; } while (--sz1 >= 0) *d++ = *s1++; while (--sz2 >= 0) *d++ = *s2++; *dbase = d - dbase - 1; return dbase; } long *P_setint(d, s1, s2) /* d := s1 * s2 */ register long *d, *s1, *s2; { long *dbase = d++; register int sz1 = *s1++, sz2 = *s2++; while (--sz1 >= 0 && --sz2 >= 0) *d++ = *s1++ & *s2++; while (--d > dbase && !*d) ; *dbase = d - dbase; return dbase; } long *P_setdiff(d, s1, s2) /* d := s1 - s2 */ register long *d, *s1, *s2; { long *dbase = d++; register int sz1 = *s1++, sz2 = *s2++; while (--sz1 >= 0 && --sz2 >= 0) *d++ = *s1++ & ~*s2++; if (sz1 >= 0) { while (sz1-- >= 0) *d++ = *s1++; } while (--d > dbase && !*d) ; *dbase = d - dbase; return dbase; } long *P_setxor(d, s1, s2) /* d := s1 / s2 */ register long *d, *s1, *s2; { long *dbase = d++; register int sz1 = *s1++, sz2 = *s2++; while (sz1 > 0 && sz2 > 0) { *d++ = *s1++ ^ *s2++; sz1--, sz2--; } while (--sz1 >= 0) *d++ = *s1++; while (--sz2 >= 0) *d++ = *s2++; while (--d > dbase && !*d) ; *dbase = d - dbase; return dbase; } int P_inset(val, s) /* val IN s */ register unsigned val; register long *s; { register int bit; bit = val % SETBITS; val /= SETBITS; if (val < *s++ && ((1< size) { s += size; while (val > size) *++s = 0, size++; *sbase = size; } else s += val; *s |= 1< (int)v2) return sbase; b1 = v1 % SETBITS; v1 /= SETBITS; b2 = v2 % SETBITS; v2 /= SETBITS; size = *s; v1++; if (++v2 > size) { while (v2 > size) s[++size] = 0; s[v2] = 0; *s = v2; } s += v1; if (v1 == v2) { *s |= (~((-2)<<(b2-b1))) << b1; } else { *s++ |= (-1) << b1; while (++v1 < v2) *s++ = -1; *s |= ~((-2) << b2); } return sbase; } long *P_remset(s, val) /* s := s - [val] */ register long *s; register unsigned val; { register int bit; bit = val % SETBITS; val /= SETBITS; if (++val <= *s) { if (!(s[val] &= ~(1<= 0) { if (*s1++ != *s2++) return 0; } return 1; } int P_subset(s1, s2) /* s1 <= s2 */ register long *s1, *s2; { register int sz1 = *s1++, sz2 = *s2++; if (sz1 > sz2) return 0; while (--sz1 >= 0) { if (*s1++ & ~*s2++) return 0; } return 1; } long *P_setcpy(d, s) /* d := s */ register long *d, *s; { register long *save_d = d; #ifdef SETCPY_MEMCPY memcpy(d, s, (*s + 1) * sizeof(long)); #else register int i = *s + 1; while (--i >= 0) *d++ = *s++; #endif return save_d; } /* s is a "smallset", i.e., a 32-bit or less set stored directly in a long. */ long *P_expset(d, s) /* d := s */ register long *d; register long s; { if (s) { d[1] = s; *d = 1; } else *d = 0; return d; } long P_packset(s) /* convert s to a small-set */ register long *s; { if (*s++) return *s; else return 0; } /* Oregon Software Pascal extensions, courtesy of William Bader */ int P_getcmdline(l, h, line) int l, h; Char *line; { int i, len; char *s; h = h - l + 1; len = 0; for(i = 1; i < P_argc; i++) { s = P_argv[i]; while (*s) { if (len >= h) return len; line[len++] = *s++; } if (len >= h) return len; line[len++] = ' '; } return len; } Void TimeStamp(Day, Month, Year, Hour, Min, Sec) int *Day, *Month, *Year, *Hour, *Min, *Sec; { #ifndef NO_TIME struct tm *tm; long clock; time(&clock); tm = localtime(&clock); *Day = tm->tm_mday; *Month = tm->tm_mon + 1; /* Jan = 0 */ *Year = tm->tm_year; if (*Year < 1900) *Year += 1900; /* year since 1900 */ *Hour = tm->tm_hour; *Min = tm->tm_min; *Sec = tm->tm_sec; #endif } /* SUN Berkeley Pascal extensions */ Void P_sun_argv(s, len, n) register char *s; register int len, n; { register char *cp; if ((unsigned)n < P_argc) cp = P_argv[n]; else cp = ""; while (*cp && --len >= 0) *s++ = *cp++; while (--len >= 0) *s++ = ' '; } int _OutMem() { return _Escape(-2); } int _CaseCheck() { return _Escape(-9); } int _NilCheck() { return _Escape(-3); } /* The following is suitable for the HP Pascal operating system. It might want to be revised when emulating another system. */ char *_ShowEscape(buf, code, ior, prefix) char *buf, *prefix; int code, ior; { char *bufp; if (prefix && *prefix) { strcpy(buf, prefix); strcat(buf, ": "); bufp = buf + strlen(buf); } else { bufp = buf; } if (code == -10) { sprintf(bufp, "Pascal system I/O error %d", ior); switch (ior) { case 3: strcat(buf, " (illegal I/O request)"); break; case 7: strcat(buf, " (bad file name)"); break; case FileNotFound: /*10*/ strcat(buf, " (file not found)"); break; case FileNotOpen: /*13*/ strcat(buf, " (file not open)"); break; case BadInputFormat: /*14*/ strcat(buf, " (bad input format)"); break; case 24: strcat(buf, " (not open for reading)"); break; case 25: strcat(buf, " (not open for writing)"); break; case 26: strcat(buf, " (not open for direct access)"); break; case 28: strcat(buf, " (string subscript out of range)"); break; case EndOfFile: /*30*/ strcat(buf, " (end-of-file)"); break; case FileWriteError: /*38*/ strcat(buf, " (file write error)"); break; } } else { sprintf(bufp, "Pascal system error %d", code); switch (code) { case -2: strcat(buf, " (out of memory)"); break; case -3: strcat(buf, " (reference to NIL pointer)"); break; case -4: strcat(buf, " (integer overflow)"); break; case -5: strcat(buf, " (divide by zero)"); break; case -6: strcat(buf, " (real math overflow)"); break; case -8: strcat(buf, " (value range error)"); break; case -9: strcat(buf, " (CASE value range error)"); break; case -12: strcat(buf, " (bus error)"); break; case -20: strcat(buf, " (stopped by user)"); break; } } return buf; } int _Escape(code) int code; { char buf[100]; P_escapecode = code; if (__top_jb) { __p2c_jmp_buf *jb = __top_jb; __top_jb = jb->next; longjmp(jb->jbuf, 1); } if (code == 0) exit(0); if (code == -1) exit(1); fprintf(stderr, "%s\n", _ShowEscape(buf, P_escapecode, P_ioresult, "")); exit(1); } int _EscIO(code) int code; { P_ioresult = code; return _Escape(-10); } /* End. */ /* p2c: meteorJ.p, line 1: * Note: Unexpected name "magfile" in program header [262] */ /* p2c: meteorJ.p, line 1: * Note: Unexpected name "coefile" in program header [262] */ /* version I: Wed Jun 27 13:36:06 EDT 1990 */ /* copyright Prof. K. Steiglitz Dept. of Computer Science Princeton University Princeton, NJ 08544 */ /* Constraint-based design of linear-phase fir filters with upper and lower bounds, and convexity constraints. Finds minimum length, or optimizes fixed length, or pushes band-edges. If L is the filter length, the models are odd-length cosine: sum ( i from 0 to (L-1)/2 ) coeff[i]*cos(i*omega) sine: sum ( i from 0 to (L-3)/2 ) coeff[i]*sin((i+1)*omega) even-length cosine: sum ( i from 0 to L/2-1 ) coeff[i]*cos((i+.5)*omega) sine: sum ( i from 0 to L/2-1 ) coeff[i]*sin((i+.5)*omega) */ #define maxpivots 1000 /* maximum no. of pivots */ #define small 1.0e-8 /* small number used in defining band-edges */ #define large 1.0e+31 /* large number used in search for minimum cost column */ #define mmax 64 /* max. no. of coefficients */ #define Lmax 129 /* max. filter length, L = 2*m-1, 2*m, or 2*m+1 */ #define mmaxm 63 /* mmax - 1 */ #define nmax 640 /* max. size of n, where there are n+1 grid-points */ #define ncolmax 6000 /* max. no. of columns allowed in tableau */ #define nspecmax 20 /* max. no. of specifications */ #define eps 1.0e-8 /* for testing for zero */ Static FILE *magfile, *coefile; /* files for output */ Static boolean done, unbounded, optimal; /* flags for simplex */ Static enum { toomanycols, unbdual, infdual, toomanypivots, opt, infprimal } result; /* result of simplex */ Static long iteration; /* iteration count, index */ Static Char ch; /* character for input */ Static char m; /* no. of coefficients, left and right half m */ Static uchar L; /* filter length = 2*m-1, 2*m, 2*m+1 */ Static short n; /* there are n+1 grid-points from 0 to pi */ Static long numpivots; /* pivot count */ Static long pivotcol, pivotrow; /* pivot column and row */ Static double pivotel; /* pivot element */ Static double pi; /* 3.14159265358979... */ Static double cbar; /* price when searching for entering column */ Static double coeff[mmaxm + 1]; /* coefficients */ Static double carry[mmax + 2][mmax + 2]; /* inverse-basis matrix of the revised simplex method */ Static char phase; /* phase */ Static double price[mmax + 1]; /* shadow prices = row -1 of carry = -dual variables = -coefficients */ Static long basis[mmax + 1]; /* basis columns, negative integers artificial */ Static char nspec; /* no. of bands */ Static enum { con, lim } spectype[nspecmax]; /* type of band */ Static double left[nspecmax], right[nspecmax]; /* bandedges as read in */ Static double freq[ncolmax]; /* frequencies at grid points */ Static double bound1[nspecmax], bound2[nspecmax]; /* left and right bounds */ Static Char sense[nspecmax]; /* sense of constraint, + up, - down */ Static Char interpolate[nspecmax]; /* g=geometric, a=arithmetic */ Static short firstcol[nspecmax]; /* leftmost column of spec */ Static short lastcol[nspecmax]; /* rightmost column of spec */ Static boolean foundfeas; /* found feasible solution */ Static long mlargest, msmallest; /* range of m */ Static long Llargest, Lsmallest; /* range of L = 2*m-1, 2*m, or 2*m+1 */ Static short ncol; /* number of columns */ Static double tab[mmax + 1][ncolmax]; /* tableau */ Static double d[ncolmax]; /* current cost vector */ Static double c[ncolmax]; /* cost in original problem */ Static double curcol[mmax + 2]; /* current column */ Static double curcost; /* current cost */ Static long bestm; /* best order */ Static Char hug[nspecmax]; /* allow this constraint to be hugged? */ Static enum { findlen, maxdist, pushedge } whattodo; /* type of optimization */ Static long npushed; /* number of bandedges pushed */ Static enum { rr, ll } whichway; /* push which way? */ Static long bandpushed[nspecmax]; /* bandedges pushed */ Static double lowlim; /* lower limit for finding if primal is feasible */ Static boolean oddlength; /* odd-length filters? */ Static enum { cosine, sine } symtype; /* cosine or sine symmetry */ Static Void makebands(i) long i; { /* fills in frequencies to make grid - frequencies are kept as reals in radians, and each band has equally spaced grid points */ long k, kmax; if (i == 1) firstcol[i - 1] = 1; else firstcol[i - 1] = lastcol[i - 2] + 1; kmax = (long)((right[i - 1] - left[i - 1]) * n / 0.5 + small); /* kmax+1 cols. in this band */ if (kmax == 0) freq[firstcol[i - 1] - 1] = 2.0 * pi * left[i - 1]; else { for (k = 0; k <= kmax; k++) freq[firstcol[i - 1] + k - 1] = 2.0 * pi * (left[i - 1] + (right[i - 1] - left[i - 1]) * k / kmax); } lastcol[i - 1] = firstcol[i - 1] + kmax; } /* makebands */ Static double trig0(i, f) long i; double f; { /* trig function in filter transfer function */ double Result; if (oddlength) { if (symtype == cosine) Result = cos(i * f); else Result = sin((i + 1) * f); } if (oddlength) return Result; if (symtype == cosine) return cos((i + 0.5) * f); else { return sin((i + 0.5) * f); } /* return Result; removed to avoid warning */ } /* trig0 */ Static double trig2(i, f) long i; double f; { /* second derivative of trig function in filter transfer function */ double Result; if (oddlength) { if (symtype == cosine) Result = -i * i * cos(i * f); else Result = -(i + 1) * (i + 1) * sin((i + 1) * f); } if (oddlength) return Result; if (symtype == cosine) return (-(i + 0.5) * (i + 0.5) * cos((i + 0.5) * f)); else { return (-(i + 0.5) * (i + 0.5) * sin((i + 0.5) * f)); } /* return Result; removed to avoid warning */ } /* trig2 */ Static Void convex(i) long i; { /* sets up tableau columns for convexity constraints on magnitude */ long row, col; /* gridpoint, side of polygon, row, col */ long FORLIM, FORLIM1; makebands(i); FORLIM = lastcol[i - 1]; for (col = firstcol[i - 1] - 1; col < FORLIM; col++) { /* for all frequencies in band */ c[col] = 0.0; FORLIM1 = m; for (row = 0; row < FORLIM1; row++) { tab[row][col] = trig2(row, freq[col]); /* normal constraint is <= */ if (sense[i - 1] == '+') tab[row][col] = -tab[row][col]; } /* for row */ tab[m][col] = 0.0; } /*for col */ } /* convex */ Static Void limit(i) long i; { /* sets up tableau columns for upper or lower bounds on transfer function for specification i; the bound is linearly interpolated between the start and end of the band */ long row, col; /* gridpoint, side of polygon, row, col */ long FORLIM, FORLIM1; makebands(i); FORLIM = lastcol[i - 1]; for (col = firstcol[i - 1] - 1; col < FORLIM; col++) { /* for all frequencies in band */ if (firstcol[i - 1] == lastcol[i - 1]) c[col] = bound1[i - 1]; else { if (interpolate[i - 1] == 'g') c[col] = bound1[i - 1] * exp((double)(col - firstcol[i - 1] + 1) / (lastcol[i - 1] - firstcol[i - 1]) * log( fabs(bound2[i - 1] / bound1[i - 1]))); if (interpolate[i - 1] == 'a') c[col] = bound1[i - 1] + (double)(col - firstcol[i - 1] + 1) / (lastcol[i - 1] - firstcol[i - 1]) * (bound2[i - 1] - bound1[i - 1]); } if (sense[i - 1] == '-') c[col] = -c[col]; FORLIM1 = m; for (row = 0; row < FORLIM1; row++) { tab[row][col] = trig0(row, freq[col]); if (sense[i - 1] == '-') tab[row][col] = -tab[row][col]; } if (hug[i - 1] != 'h') tab[m][col] = 1.0; else tab[m][col] = 0.0; } /* for col */ } /* limit */ Static Void readdata() { /* reads in problem data */ /* not meant to be interactive; aborts on bad filter lengths*/ long i; boolean allhugged; /* to see if all constraints are hugged */ long FORLIM; scanf("%ld%ld%*[^\n]", &Lsmallest, &Llargest); getchar(); if (Lsmallest < 1 || Llargest > Lmax) { printf("Lsmallest or Llargest out of range: quitting\n"); _Escape(0); } if ((Lsmallest & 1) != (Llargest & 1)) { printf("parity of Lsmallest and Llargest unequal: quitting\n"); _Escape(0); } scanf("%c%*[^\n]", &ch); getchar(); if (ch == '\n') ch = ' '; if (ch == 'c') symtype = cosine; else symtype = sine; if (symtype == cosine) printf("cosine symmetry\n"); else printf("sine symmetry\n"); oddlength = Lsmallest & 1; if (oddlength) { if (symtype == cosine) { msmallest = (Lsmallest + 1) / 2; mlargest = (Llargest + 1) / 2; } else { msmallest = (Lsmallest - 1) / 2; mlargest = (Llargest - 1) / 2; } } if (!oddlength) { msmallest = Lsmallest / 2; mlargest = Llargest / 2; } if (Lsmallest != Llargest) { whattodo = findlen; printf("finding minimum length: range %4ld to %4ld\n", Lsmallest, Llargest); } if (Lsmallest == Llargest) { m = msmallest; L = Lsmallest; printf("fixed length of %4d\n", L); scanf("%c%*[^\n]", &ch); getchar(); /* right, left, or neither: edges to be pushed? */ if (ch == '\n') ch = ' '; if (ch == 'n') { whattodo = maxdist; printf("maximizing distance from constraints not marked hugged\n"); } else { whattodo = pushedge; if (ch == 'r') whichway = rr; else whichway = ll; scanf("%ld%*[^\n]", &npushed); getchar(); FORLIM = npushed; for (i = 0; i < FORLIM; i++) scanf("%ld", &bandpushed[i]); scanf("%*[^\n]"); getchar(); if (whichway == rr) printf("pushing bandedges right\n"); if (whichway == ll) printf("pushing bandedges left\n"); printf("constraint numbers: "); FORLIM = npushed; for (i = 0; i < FORLIM; i++) printf("%3ld ", bandpushed[i]); putchar('\n'); } } scanf("%hd%*[^\n]", &n); getchar(); /* there are n+1 grid points between 0 and pi */ nspec = 0; scanf("%c%*[^\n]", &ch); getchar(); if (ch == '\n') ch = ' '; while (ch != 'e') { /* 'e' for end */ nspec++; i = nspec; if (ch == 'c') { spectype[i - 1] = con; scanf("%c%*[^\n]", &sense[i - 1]); getchar(); if (sense[i - 1] == '\n') sense[i - 1] = ' '; scanf("%lg%lg%*[^\n]", &left[i - 1], &right[i - 1]); getchar(); printf("constraint %ld: convex, sense %c\n", i, sense[i - 1]); printf(" bandedges: %10.4f %10.4f\n", left[i - 1], right[i - 1]); } if (ch == 'l') { spectype[i - 1] = lim; scanf("%c%*[^\n]", &sense[i - 1]); getchar(); if (sense[i - 1] == '\n') sense[i - 1] = ' '; scanf("%c%*[^\n]", &interpolate[i - 1]); getchar(); if (interpolate[i - 1] == '\n') interpolate[i - 1] = ' '; scanf("%c%*[^\n]", &hug[i - 1]); getchar(); if (hug[i - 1] == '\n') hug[i - 1] = ' '; scanf("%lg%lg%*[^\n]", &left[i - 1], &right[i - 1]); getchar(); scanf("%lg%lg%*[^\n]", &bound1[i - 1], &bound2[i - 1]); getchar(); if (interpolate[i - 1] == 'g' && bound1[i - 1] * bound2[i - 1] == 0.0) { printf("geometrically interpolated bandedge in constraint %5ld is zero\n", i); _Escape(0); } if (sense[i - 1] == '+') printf("constraint %2ld: upper limit\n", i); if (sense[i - 1] == '-') printf("constraint %2ld: lower limit\n", i); if (interpolate[i - 1] == 'g') printf(" geometric interpolation\n"); if (interpolate[i - 1] == 'a') printf(" arithmetic interpolation\n"); if (hug[i - 1] == 'h') printf(" this constraint will be hugged\n"); else printf(" this constraint will be optimized\n"); printf(" bandedges: %10.4f %10.4f\n", left[i - 1], right[i - 1]); printf(" bounds: %10.4f %10.4f\n", bound1[i - 1], bound2[i - 1]); } makebands(i); printf(" initial columns: %10d %10d\n", firstcol[i - 1], lastcol[i - 1]); scanf("%c%*[^\n]", &ch); getchar(); /* next */ if (ch == '\n') ch = ' '; } /* while */ ncol = lastcol[nspec - 1]; printf("number of specs= %5d\n", nspec); printf("initial number of columns= %5d\n", ncol); allhugged = true; /* check to see if all limit constraints are hugged */ FORLIM = nspec; for (i = 0; i < FORLIM; i++) { if (spectype[i] == lim && hug[i] != 'h') allhugged = false; } if (allhugged) { printf("all constraints are hugged: ill-posed problem\n"); _Escape(0); } } /* readdata */ Static Void setup() { /* initializes constraints */ long i, FORLIM; pi = 4.0 * atan(1.0); FORLIM = nspec; for (i = 1; i <= FORLIM; i++) { switch (spectype[i - 1]) { case con: convex(i); break; case lim: limit(i); break; }/* case */ } /*for */ ncol = lastcol[nspec - 1]; } /* setup */ Static Void columnsearch() { /* looks for favorable column to enter basis. returns lowest cost and its column number, or turns on the flag optimal */ long i, col; double tempcost; /* minimum cost, temporary cost of column */ long FORLIM, FORLIM1; FORLIM = m; for (i = 0; i <= FORLIM; i++) /* set up price vector */ price[i] = -carry[0][i + 1]; optimal = false; cbar = large; pivotcol = 0; FORLIM = ncol; for (col = 1; col <= FORLIM; col++) { tempcost = d[col - 1]; FORLIM1 = m; for (i = 0; i <= FORLIM1; i++) tempcost -= price[i] * tab[i][col - 1]; if (cbar > tempcost) { cbar = tempcost; pivotcol = col; } } /* for col */ if (cbar > -eps) optimal = true; } /* columnsearch */ Static Void rowsearch() { /* looks for pivot row. returns pivot row number, or turns on the flag unbounded */ long i, j; double ratio, minratio; /* ratio and minimum ratio for ratio test */ long FORLIM, FORLIM1; FORLIM = m + 1; for (i = 1; i <= FORLIM; i++) { /* generate column */ curcol[i] = 0.0; /* current column = B inverse * original col. */ FORLIM1 = m; for (j = 0; j <= FORLIM1; j++) curcol[i] += carry[i][j + 1] * tab[j][pivotcol - 1]; } curcol[0] = cbar; /* first element in current column */ pivotrow = -1; minratio = large; FORLIM = m; for (i = 0; i <= FORLIM; i++) { /* ratio test */ if (curcol[i + 1] > eps) { ratio = carry[i + 1][0] / curcol[i + 1]; if (minratio > ratio) { /* favorable row */ minratio = ratio; pivotrow = i; pivotel = curcol[i + 1]; } else { /* break tie with max pivot */ if (minratio == ratio && pivotel < curcol[i + 1]) { pivotrow = i; pivotel = curcol[i + 1]; } } } /* curcol > eps */ } /* for i */ if (pivotrow == -1) unbounded = true; /* nothing found */ else unbounded = false; } /* rowsearch */ Static Void pivot() { /* pivots */ long i, j, FORLIM, FORLIM1; basis[pivotrow] = pivotcol; FORLIM = m + 1; for (j = 0; j <= FORLIM; j++) carry[pivotrow + 1][j] /= pivotel; FORLIM = m + 1; for (i = 0; i <= FORLIM; i++) { if (i - 1 != pivotrow) { FORLIM1 = m + 1; for (j = 0; j <= FORLIM1; j++) carry[i][j] -= carry[pivotrow + 1][j] * curcol[i]; } } curcost = -carry[0][0]; } /* pivot */ Static Void changephase() { /* changes phase from 1 to 2, by switching to original cost vector */ long i, j, b, FORLIM, FORLIM1; phase = 2; FORLIM = m; for (i = 0; i <= FORLIM; i++) { if (basis[i] <= 0) printf("...artificial basis element %5ld remains in basis after phase 1\n", basis[i]); } FORLIM = ncol; for (j = 0; j < FORLIM; j++) /* switch to original cost vector */ d[j] = c[j]; FORLIM = m + 1; for (j = 0; j <= FORLIM; j++) { carry[0][j] = 0.0; FORLIM1 = m; for (i = 0; i <= FORLIM1; i++) { b = basis[i]; /* ignore artificial basis elements that are */ if (b >= 1) /* still in basis */ carry[0][j] -= c[b - 1] * carry[i + 1][j]; } /* for i */ } /* for j */ curcost = -carry[0][0]; } /* changephase */ Static double mag(f) double f; { /* computes magnitude function, given radian frequency f */ long i; double temp; long FORLIM; temp = 0.0; FORLIM = m; for (i = 0; i < FORLIM; i++) temp += coeff[i] * trig0(i, f); return temp; } /* mag */ Static Void wrapup() { /* prints results, final frequency response */ long i, k, FORLIM; if (magfile != NULL) magfile = freopen("magfile", "wb", magfile); else magfile = fopen("magfile", "wb"); if (magfile == NULL) _EscIO(FileNotFound); if (coefile != NULL) { /* open files for writing */ coefile = freopen("coefile", "wb", coefile); } else coefile = fopen("coefile", "wb"); if (coefile == NULL) _EscIO(FileNotFound); if (oddlength && symtype == cosine) { /* write coeffs to coefile */ fprintf(coefile, "%4d cosine symmetry\n", m * 2 - 1); /* L = length, odd */ for (i = m - 1; i >= 1; i--) fprintf(coefile, "% .5E\n", coeff[i] / 2); fprintf(coefile, "% .5E\n", coeff[0]); FORLIM = m; for (i = 1; i < FORLIM; i++) fprintf(coefile, "% .5E\n", coeff[i] / 2); } /* odd, cosine */ if (!oddlength && symtype == cosine) { fprintf(coefile, "%4d cosine symmetry\n", m * 2); /* L = length, even */ for (i = m - 1; i >= 0; i--) fprintf(coefile, "% .5E\n", coeff[i] / 2); FORLIM = m; for (i = 0; i < FORLIM; i++) fprintf(coefile, "% .5E\n", coeff[i] / 2); } /* even, cosine */ if (oddlength && symtype == sine) { fprintf(coefile, "%4d sine symmetry\n", m * 2 + 1); /* L = length, odd */ for (i = m - 1; i >= 0; i--) /* negative of the first m coefs. */ fprintf(coefile, "% .5E\n", coeff[i] / -2); fprintf(coefile, " 0.0\n"); /* middle coefficient is always 0 */ FORLIM = m; for (i = 0; i < FORLIM; i++) fprintf(coefile, "% .5E\n", coeff[i] / 2); } /* odd, sine */ if (!oddlength && symtype == sine) { fprintf(coefile, "%4d sine symmetry\n", m * 2); /* L = length, even */ for (i = m - 1; i >= 0; i--) /* negative of the first m coefs. */ fprintf(coefile, "% .5E\n", coeff[i] / -2); FORLIM = m; for (i = 0; i < FORLIM; i++) fprintf(coefile, "% .5E\n", coeff[i] / 2); } /* even, sine */ FORLIM = n; for (k = 0; k <= FORLIM; k++) /* magnitude on regular grid */ fprintf(magfile, "%10.4f % .5E\n", 0.5 * k / n, mag(k * pi / n)); fprintf(magfile, "\n magnitude at bandedges\n\n"); FORLIM = nspec; for (i = 0; i < FORLIM; i++) { if (spectype[i] == lim) { fprintf(magfile, "%10.4f % .5E\n", freq[firstcol[i] - 1] * 0.5 / pi, mag(freq[firstcol[i] - 1])); fprintf(magfile, "%10.4f % .5E\n", freq[lastcol[i] - 1] * 0.5 / pi, mag(freq[lastcol[i] - 1])); putchar('\n'); } } } /* wrapup */ Static Void simplex() { /* simplex for linear programming */ long i, j, col, row, FORLIM, FORLIM1; done = false; phase = 1; FORLIM = m + 1; for (i = 0; i <= FORLIM; i++) { for (j = 0; j <= mmax + 1; j++) carry[i][j] = 0.0; } FORLIM = m + 1; for (i = 1; i <= FORLIM; i++) /* artificial basis */ carry[i][i] = 1.0; carry[0][0] = -1.0; /* - initial cost */ curcost = -carry[0][0]; carry[m + 1][0] = 1.0; /* variable minimized in primal */ FORLIM = m; for (i = 0; i <= FORLIM; i++) /* initial, artificial basis */ basis[i] = -i; if (ncol <= ncolmax) { /* check number of columns */ FORLIM = ncol; for (col = 0; col < FORLIM; col++) { /* initialize cost for phase 1 */ d[col] = 0.0; FORLIM1 = m; for (row = 0; row <= FORLIM1; row++) d[col] -= tab[row][col]; } } else { printf("...termination: too many columns for storage\n"); done = true; result = toomanycols; } numpivots = 0; while (numpivots < maxpivots && !done && (curcost > lowlim || phase == 1)) { columnsearch(); if (!optimal) { /* not optimal */ rowsearch(); if (unbounded) { done = true; result = unbdual; /* dual of problem is unbounded */ } else { pivot(); numpivots++; if (numpivots == 1 || numpivots % 10 == 0) printf("pivot %4ld cost= % .5E\n", numpivots, curcost); /* p2c: meteorJ.p, line 595: * Note: Using % for possibly-negative arguments [317] */ } continue; } /* not optimal */ if (phase == 1) { if (curcost > eps) { done = true; /* dual of problem is infeasible */ result = infdual; /* this happens if all specs are hugged */ } else { if (numpivots != 1 && numpivots % 10 != 0) printf("pivot %4ld cost= % .5E\n", numpivots, curcost); /* p2c: meteorJ.p, line 609: * Note: Using % for possibly-negative arguments [317] */ printf("phase 1 successfully completed\n"); changephase(); } continue; } /* if phase = 1 */ if (numpivots != 1 && numpivots % 10 != 0) printf("pivot %4ld cost= % .5E\n", numpivots, curcost); /* p2c: meteorJ.p, line 617: * Note: Using % for possibly-negative arguments [317] */ printf("phase 2 successfully completed\n"); done = true; result = opt; } /* while */ if (curcost <= lowlim && phase == 2) { if (numpivots != 1 && numpivots % 10 != 0) printf("pivot %4ld cost= % .5E\n", numpivots, curcost); /* p2c: meteorJ.p, line 626: * Note: Using % for possibly-negative arguments [317] */ result = infprimal; } if (numpivots >= maxpivots) { printf("...termination: maximum number of pivots exceeded\n"); result = toomanypivots; } /* optimal */ } /* simplex */ Static Void printresult() { /* prints enumerated type */ switch (result) { case toomanycols: printf("too many columns in specifications\n"); break; case unbdual: /* unbounded dual */ printf("infeasible (unbounded dual)\n"); break; case infdual: /* infeasible dual */ printf("infeasible or unbounded\n"); break; case toomanypivots: printf("too many pivots\n"); break; case opt: printf("optimum obtained\n"); break; case infprimal: printf("infeasible \n"); break; }/* case */ } /* printresult */ Static Void getm() { /* find best order (and hence length) */ char leftm, rightm; long i; boolean foundm, checkedleft, checkedright; long FORLIM; foundfeas = false; /* flag set when a feasible solution is found */ iteration = 0; leftm = msmallest; rightm = mlargest; foundm = false; checkedleft = false; checkedright = false; while (!foundm) { if (iteration == 0) /* first time through */ m = leftm + (rightm - leftm) / 2; iteration++; printf("\niteration %2ld\n", iteration); if (oddlength) { if (symtype == cosine) printf("L= %4d\n", m * 2 - 1); else printf("L= %4d\n", m * 2 + 1); } if (!oddlength) printf("L= %4d\n", m * 2); setup(); simplex(); printresult(); if (result == opt) { foundfeas = true; rightm = m; bestm = m; checkedright = true; /* right side of bracket has been checked */ if (oddlength) { if (symtype == cosine) printf("new best length L= %4ld\n", bestm * 2 - 1); else printf("new best length L= %4ld\n", bestm * 2 + 1); } if (!oddlength) printf("new best length L= %4ld\n", bestm * 2); FORLIM = m; for (i = 0; i < FORLIM; i++) coeff[i] = -carry[0][i + 1]; } if (result != opt) { leftm = m; checkedleft = true; /* left side of bracket has been checked */ } if (rightm > leftm + 1) m = leftm + (rightm - leftm) / 2; if (rightm == leftm + 1) { if (!checkedleft) { m = leftm; checkedleft = true; } else if (!checkedright) { m = rightm; checkedright = true; } else foundm = true; } if (rightm == leftm) { foundm = true; } } /* while */ if (!foundfeas) { printf("no feasible solution found\n"); _Escape(0); } m = bestm; putchar('\n'); if (oddlength) { if (symtype == cosine) printf("best length L= %4ld\n", bestm * 2 - 1); else printf("best length L= %4ld\n", bestm * 2 + 1); } if (!oddlength) { printf("best length L= %4ld\n", bestm * 2); } } /* getm */ Static Void getedge() { /* optimize bandedge */ double lefte, righte, newe, beste, onespace, stopspace; /* nominal grid spacing, stop criterion */ long i, FORLIM; onespace = 0.5 / n; /* space between grid points */ stopspace = onespace / 10.0; /* stop criterion is 1/10 nominal grid spacing */ if (whichway == rr) { lefte = left[bandpushed[0] - 1]; /* start with rightmost leftedge */ FORLIM = npushed; for (i = 1; i < FORLIM; i++) { if (left[bandpushed[i] - 1] > lefte) lefte = left[bandpushed[i] - 1]; } righte = 0.5; } else { lefte = 0.0; righte = right[bandpushed[0] - 1]; FORLIM = npushed; for (i = 1; i < FORLIM; i++) { /* start with leftmost rightedge */ if (right[bandpushed[i] - 1] < righte) righte = right[bandpushed[i] - 1]; } } iteration = 0; foundfeas = false; /* flag set when a feasible solution is found */ while (righte - lefte > stopspace) { iteration++; newe = (righte + lefte) / 2.0; printf("\niteration %4ld\n", iteration); printf("trying new edge = %10.4f\n", newe); FORLIM = npushed; for (i = 0; i < FORLIM; i++) { if (whichway == rr) right[bandpushed[i] - 1] = newe; else left[bandpushed[i] - 1] = newe; } setup(); simplex(); printresult(); if (result == opt) { if (whichway == rr) lefte = newe; else righte = newe; foundfeas = true; beste = newe; FORLIM = m; for (i = 0; i < FORLIM; i++) coeff[i] = -carry[0][i + 1]; } if (result != opt) { if (whichway == rr) righte = newe; else lefte = newe; } } /* while */ putchar('\n'); if (!foundfeas) { printf("no feasible bandedge found\n"); _Escape(0); } printf("found edge= %10.4f\n", beste); FORLIM = npushed; for (i = 0; i < FORLIM; i++) { if (whichway == rr) right[bandpushed[i] - 1] = beste; else left[bandpushed[i] - 1] = beste; } FORLIM = nspec; for (i = 1; i <= FORLIM; i++) makebands(i); } /* getedge */ Static Void getmaxdist() { /* maximizes distance from constraints */ long i, FORLIM; printf("optimization: maximize distance from constraints\n"); setup(); simplex(); printresult(); if (result != opt) _Escape(0); printf("final cost= distance from constraints= % .5E\n", curcost); FORLIM = m; for (i = 0; i < FORLIM; i++) /* record coefficients */ coeff[i] = -carry[0][i + 1]; /* don't go back if unsuccessful */ } /* getmaxdist */ main(argc, argv) int argc; Char *argv[]; { /* main */ PASCAL_MAIN(argc, argv); coefile = NULL; magfile = NULL; printf("welcome to meteor:\n"); printf("constraint-based, linear-phase FIR filter design\n"); readdata(); lowlim = -eps; /* dual cost negative => primal infeasible */ switch (whattodo) { /* case */ case findlen: getm(); break; case pushedge: getedge(); break; case maxdist: getmaxdist(); break; } /* no return here if unsuccessful */ wrapup(); if (magfile != NULL) fclose(magfile); if (coefile != NULL) fclose(coefile); exit(0); } /* main */ /* End. */