#include	<math.h>
#include	"misc.h"
#define		RMAX	100	/* maximum number of registers */
#define		CMDS	500	/* maximum number of commands */
#define		PI	3.141592654/* pi */


/* OPSTRM -- Operate on an input stream. */

opstrm (input, argc, argv)
char	input[];
int	*argc;
char	*argv[];
{
    float   focasrand(), func3(), func4(), pois();
    float   x[RMAX], z, store[10], value;
    int     steps, flag, cflag, lseed;
    int     i, j, k, m, n;
    char    cmds[CMDS][20], tmp[CMDS][20];
    char    op[20], y[100], frmt[20], prgm[100];
    FILE    *fd, *fin, *fout;

    if (strcmp (input, "STDIN") == 0)
	fd = getstdin;
    else if ((fd = fopen (input, "r")) == NULL) {
	printf ("OPSTRM: Can't read file %s\n", input);
	return;
    }

    flag = 0;
    steps = 0;
    cflag = 0;
    lseed = 1;
    if (*argc == 0) {
	cflag = 1;
	steps = 2;
    } else {
	for (i = 0; i < *argc; i++, steps++) {
	    switch (argv[i][0]) {
	    case 'P': 
		sprintf (prgm, "%s.opstrm", &argv[i][1]);
		if ((fout = fopen (prgm, "w")) == NULL)
		    focaserr (1, "opstrm: Can't create ", prgm);
		for (;;) {
		    fscanf (fd, "%s", op);
		    if (feof (fd))
			goto done;
		    strncpy (cmds[0], op, 20);
		    fprintf (fout, "%s\n", op);
		}
	    case 'R': 
		if (argv[i][0] == 'R')
		    flag++;
	    default: 
		strncpy (cmds[steps], argv[i], 20);
	    }
	}
    }

    j = -1;
    for (; (flag != 0) && (cflag == 0);) {
	for (i = 0; i < steps; i++)
	    strcpy (tmp[i], cmds[i]);
	for (i = 0, j = steps, steps = 0, flag = 0; i < j; i++, steps++) {
	    switch (tmp[i][0]) {
		case 'R': 
		    sprintf (prgm, "%s.opstrm", &tmp[i][1]);
		    if ((fin = fopen (prgm, "r")) == NULL)
			focaserr (2, "opstrm: Can't find ", prgm);
		    for (;; steps++) {
			fscanf (fin, "%s", cmds[steps]);
			if (feof (fin))
			    break;
			if (cmds[steps][0] == 'R')
			    flag++;
		    }
		    fclose (fin);
		    steps--;
		    break;
		default: 
		    strcpy (cmds[steps], tmp[i]);
	    }
	}
    }

    strcpy (frmt, "%g ");
    for (m = 1;; m++) {
	for (i = 0, j = -1; i < steps; i++) {
	    if (j > (RMAX - 1))
		goto reg;
	    if (cflag) {
		i = 0;
		fscanf (fd, "%s", cmds[i]);
		if (feof (fd))
		    goto done;
	    }
	    switch (cmds[i][0]) {
		case 'I': 
		    if (m != 1)
			++i;
		    break;
		case 'q': 
		    printf ("\n");
		    goto done;
		case 'r': 
		    fscanf (fd, "%hf", &x[++j]);
		    if (feof (fd))
			goto done;
		    break;
		case 'c': 
		    fscanf (fd, "%s", y);
		    if (feof (fd))
			goto done;
		    printf ("%s ", y);
		    if (cflag)
			fflush (stdout);
		    break;
		case 'p': 
		    if (cmds[i][1] == 'i')
			x[++j] = PI;
		    else {
			sscanf (cmds[i], "p%d", &n);
			if ((n >= 0) && (n <= 9))
			    store[n] = x[j];
			else
			    goto stre;
		    }
		    break;
		case 'g': 
		    switch (cmds[i][1]) {
			case 'c': 
			    x[++j] = m;
			    break;
			case 'l': 
			    x[++j] = log10 ((double) m);
			    break;
			default: 
			    sscanf (cmds[i], "g%d", &n);
			    if ((n >= 0) && (n <= 9))
				x[++j] = store[n];
			    else
				goto stre;
		    }
		    break;
		case 'd': 
		    if (j > 0)
			j--;
		    else
			goto stack;
		    break;
		case '+': 
		    if (j >= 1) {
			x[j - 1] = x[j - 1] + x[j];
			j--;
		    }
		    else
			goto stack;
		    break;
		case '-': 
		    if (cmds[i][1] == 0) {
			if (j >= 1) {
			    x[j - 1] = x[j - 1] - x[j];
			    j--;
			}
			else
			    goto stack;
		    }
		    else
			sscanf (cmds[i], "%hf", &x[++j]);
		    break;
		case 'x': 
		    if (j >= 1) {
			x[j - 1] = x[j - 1] * x[j];
			j--;
		    }
		    else
			goto stack;
		    break;
		case '/': 
		    if (x[j] == 0.)
			goto divide;
		    if (j >= 1) {
			x[j - 1] = x[j - 1] / x[j];
			j--;
		    }
		    else
			goto stack;
		    break;
		case 'i': 
		    if (j >= 1) {
			value = x[j - 1];
			x[j - 1] = x[j];
			x[j] = value;
		    }
		    else
			goto stack;
		    break;
		case 's': 
		    fscanf (fd, "%*s");
		    break;
		case 'e': 
		    x[j+1] = x[j];
		    j++;
		    break;
		case '=': 
		    if (j >= 0) {
			printf (frmt, x[j]);
			if (cflag)
			    fflush (stdout);
		    } else
			goto stack;
		    break;
		case 'S': 
		    printf ("%s ", cmds[i] + 1);
		    if (cflag)
			fflush (stdout);
		    break;
		case 'f': 
		    sprintf (frmt, "%%%s ", cmds[i] + 1);
		    break;
		case 'F': 
		    if (!strcmp (cmds[i], "Fmax")) {
			if (j >= 1) {
			    if (x[j] < x[j - 1]) {
				value = x[j];
				x[j] = x[j - 1];
				x[j - 1] = value;
			    }
			}
			else
			    goto stack;
		    }
		    if (!strcmp (cmds[i], "Fabs"))
			x[j] = fabs (x[j]);
		    if (!strcmp (cmds[i], "Fran"))
			x[++j] = focasrand (&lseed);
		    if (!strcmp (cmds[i], "Fsq"))
			x[j] = x[j] * x[j];
		    if (!strcmp (cmds[i], "Fpow")) {
			x[j-1] = pow (x[j-1], x[j]);
			j--;
		    }
		    if (!strcmp (cmds[i], "Fsqrt"))
			x[j] = sqrt (x[j]);
		    if (!strcmp (cmds[i], "Fln"))
			x[j] = log (x[j]);
		    if (!strcmp (cmds[i], "Fexp"))
			x[j] = exp (x[j]);
		    if (!strcmp (cmds[i], "Flog"))
			x[j] = log10 (x[j]);
		    if (!strcmp (cmds[i], "Fdex"))
			x[j] = pow (10., x[j]);
		    if (!strcmp (cmds[i], "Fmag"))
			x[j] = -2.5 * log10 (x[j]);
		    if (!strcmp (cmds[i], "Flum"))
			x[j] = pow (10., -.4 * x[j]);
		    if (!strcmp (cmds[i], "Fcos"))
			x[j] = cos (x[j]);
		    if (!strcmp (cmds[i], "Fsin"))
			x[j] = sin (x[j]);
		    if (!strcmp (cmds[i], "Frec")) {
			if (j >= 1) {
			    z = x[j] * sin (x[j - 1]);
			    x[j] = x[j] * cos (x[j - 1]);
			    x[j - 1] = z;
			}
			else
			    goto stack;
		    }
		    if (!strcmp (cmds[i], "Fpolar")) {
			if (j >= 1) {
			    z = sqrt (x[j] * x[j] + x[j - 1] * x[j - 1]);
			    x[j - 1] = atan2 (x[j - 1], x[j]);
			    x[j] = z;
			}
			else
			    goto stack;
		    }
		    if (!strcmp (cmds[i], "Fpolar1")) {
			if (j >= 1) {
			    z = sqrt (x[j] * x[j] + x[j - 1] * x[j - 1]);
			    x[j - 1] = atan2 (x[j - 1], x[j]);
			    if (x[j - 1] < -PI / 2.)
				x[j - 1] += PI;
			    if (x[j - 1] > PI / 2.)
				x[j - 1] -= PI;
			    x[j] = z;
			}
			else
			    goto stack;
		    }
		    if (!strcmp (cmds[i], "Fpois"))
			x[j] = pois (x + j - 1);
		    break;
		default: 
		    sscanf (cmds[i], "%hf", &x[++j]);
	    }
	}
	printf ("\n");
	fflush (stdout);
    }
done:
    fclose (fd);
    return;
reg: 
    focaserr (1, "opstrm: register overflow", "");
stack: 
    focaserr (2, "opstrm: stack empty", "");
divide: 
    focaserr (3, "opstrm: divide by zero", "");
stre: 
    focaserr (4, "opstrm: invalid store - ", cmds[i]);
}

float   pois (x)
float  *x;
{
    int     i, N;
    double  y, lambda;

    N = x[1];
    lambda = x[0];
    for (y = 1, i = 2; i <= N; i++)
	y *= i;
    return (pow (lambda, x[1]) * exp (-lambda) / y);
}