/*--------------------------------------------------------------------
 *	$Id: gmtmath_main.c,v 1.35 2004/09/30 18:22:06 pwessel Exp $
 *
 *	Copyright (c) 1991-2004 by P. Wessel and W. H. F. Smith
 *	See COPYING file for copying and redistribution conditions.
 *
 *	This program is free software; you can redistribute it and/or modify
 *	it under the terms of the GNU General Public License as published by
 *	the Free Software Foundation; version 2 of the License.
 *
 *	This program is distributed in the hope that it will be useful,
 *	but WITHOUT ANY WARRANTY; without even the implied warranty of
 *	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *	GNU General Public License for more details.
 *
 *	Contact info: gmt.soest.hawaii.edu
 *--------------------------------------------------------------------*/
/*
 * gmtmath.c is a reverse polish calculator that operates on table files 
 * (and constants) and perform basic mathematical operations
 * on them like add, multiply, etc.
 * Some operators only work on one operand (e.g., log, exp)
 *
 * Author:	Paul Wessel
 * Date:	10-NOV-1998
 * Version:	1.0 based on 3.1 grdmath and sample1d
 *		3.1.2 PW 03/26/99 Added -H capability
 *		3.1.2 PW 04/07/99 Added -Q for quick scalar calculator
 *		3.3.2 PW 09/10/99 Added erfinv
 *		3.3.3 PW 12/10/99 Added RAND and NRAND
 *		3.3.4 PW 03/14/00 Fix problems with EXCH and DUP for constants
 *		3.3.5 PW 07/14/00 Changed STEP to STEPT and added new STEP
 *		3.3.5 PW 07/31/00 Empty -T means there are no time columns.
 *		3.3.6 PW 08/16/00 Add LT, LE, EQ, GE, GT, NAN, CHIDIST, FDIST, TDIST
 *		3.3.6 PW 08/18/00 Add INT, ISNAN, XOR, MODE, MAD, LMSSCL, SUM
 *		3.3.6 PW 08/18/00 Added -S to just return first row
 *		3.3.6 PW 08/23/00 Added LOWER and UPPER
 *		3.4   PW 03/01/01
 *		4.0   PW 11/28/01 Added Critical values for Chi2, F, T, and Z distributions
 *		      PW 11/30/01 Added LSQFIT to solve a general least squares system
 *		      PW 12/08/01 Allow for -T<filename> with irregular time coordinates
 *		      PW 01/27/04 Added SINC and NEQ
 *		      PW 03/24/04 Added ROOTS
 *		      PW 03/28/04 Added FLIPUD, ROTT
 *		      PW 07/01/04 Added LRAND
 *		      PW 07/17/04 Added LOG2
 *
 */
 
#include "gmt.h"

#define N_OPERATORS	106
#define STACK_SIZE	50

#define ARG_IS_FILE	0
#define ARG_IS_NUMBER	1
#define ARG_IS_PI	2
#define ARG_IS_E	3
#define ARG_IS_T_MATRIX	4
#define ARG_IS_OPERATOR	100
#define N_SPECIAL_ARGS	9

#include "gmtmath_def.h"

PFV call_operator[N_OPERATORS];

int nm = 0, consumed_operands[N_OPERATORS], produced_operands[N_OPERATORS];

struct GMT_HASH localhashnode[HASH_SIZE];

double *t_coordinates;

struct TABLE_HEADER {
	int n_row;	/* Number of time-nodes (rows) */
	int n_col;	/* Number of columns */
	double t_min;	/* Minimum t value */
	double t_max;	/* Maximum t value */
	double t_inc;	/* t increment */
} header;

int decode_argument(char *txt), get_operator(char *choice);
void gmtmath_init(PFV ops[], int n_args[], int n_out[]);
void GMT_read_table (char *file, struct TABLE_HEADER *h, double ***p, int t_col);
void GMT_write_table (char *file, struct TABLE_HEADER *h, double **p, int first_last_all);
void new_table (double ***s, int n_col, int n);
void decode_columns (char *txt, int *skip, int n_col, int t_col);
void solve_LSQFIT (double **stack[], int last, int n_col, int n_row, int skip[], char *file);

BOOLEAN very_first = TRUE, irregular = FALSE, roots_found = FALSE;
char head_record[BUFSIZ];
int n_roots, r_col;

char *special_arg[N_SPECIAL_ARGS] = {
	"PI",
	"pi",
	"E",
	"e",
	"x",
	"T",
	"+",
	"-",
	"^"
};

int GMT_mode_selection = 0, GMT_n_multiples = 0;

main (int argc, char **argv)
{
	int i, j, arg, op = 0, nstack = 0, new_stack = -1, last_arg, ok = 1, type, *skip;
	int n_col = 2, t_col = 0, use_t_col, first_last_all = 0;
	
	BOOLEAN constant[STACK_SIZE], error = FALSE, set_t = FALSE, set_n = FALSE, got_t_from_file = FALSE;
	BOOLEAN set_q = FALSE, set_a = FALSE, no_time = FALSE, set_t_table = FALSE;
	
	double **stack[STACK_SIZE], **rhs, **start;

	double factor[STACK_SIZE], t_noise;
	
	char *outfile = CNULL, file[BUFSIZ];
	
	struct TABLE_HEADER tbl[STACK_SIZE], rhs_header;
	
	FILE *fp;
	
	argc = GMT_begin (argc, argv);
	
	if (argc == 2 && !strcmp (argv[1], "-")) error = GMT_quick = TRUE;
	
	if (argc == 1 || GMT_quick) {
		fprintf (stderr, "gmtmath %s - Reverse Polish Notation (RPN) calculator for table data\n\n", GMT_VERSION);
		fprintf (stderr, "usage: gmtmath [-A<t_f(t).d>] [-C<cols>] [-H[<nrec>]] [-N<n_col>/<t_col>] [-Q]\n");
		fprintf (stderr, "\t[-S[f|l]] [-T[<tmin/tmax/t_inc>]] [-V] [-bi[s][<n>]]\n\t[-bo[s][<n>]] [-f[i|o]<colinfo>] A B op C op ... = [outfile]\n\n");
		
		if (GMT_quick) exit (EXIT_FAILURE);
		
		fprintf (stderr, "\tA, B, etc are table files, constants, or symbols (see below)\n");
		fprintf (stderr, "\tThe stack can hold up to %d entries (given enough memory)\n", STACK_SIZE);
		fprintf (stderr, "\tTrigonometric operators expect radians.  The operators are:\n\n");
		fprintf (stderr, "\tName	#args	Returns:\n");
		fprintf (stderr, "\t-----------------------\n");
#include "gmtmath_explain.h"
		fprintf (stderr, "\n\tThe special symbols are:\n\n");
		fprintf (stderr, "\t  PI	= 3.1415926...\n");
		fprintf (stderr, "\t  E	= 2.7182818...\n");
		fprintf (stderr, "\t  T	= table with t-coordinates\n");
		fprintf (stderr, "\n\tOPTIONS:\n\n");
		fprintf (stderr, "\t-A Requires -N and will initialize table with file containing t and f(t) only.\n");
		fprintf (stderr, "\t   t goes into column <t_col> while f(t) goes into column <n_col> - 1.\n");
		fprintf (stderr, "\t-C change which columns to operate on [Default is all except time]\n");
		fprintf (stderr, "\t   -C reverts to the default, -Cr toggles current settings, and -Ca selects all columns.\n");
		GMT_explain_option ('H');
		fprintf (stderr, "\t-N sets the number of columns and the id of the time column (0 is first) [2/0]\n");
		fprintf (stderr, "\t-Q quick scalar calculator. Shorthand for -Ca -N1/0 -T0/0/1\n");
		fprintf (stderr, "\t-S Only write first row upon completion of calculations [write all rows]\n");
		fprintf (stderr, "\t   Optionally, append l for last row or f for first row [Default]\n");
		fprintf (stderr, "\t-T Set domain from t_min to t_max in steps of t_inc\n");
		fprintf (stderr, "\t   If a filename is given instead we read t coordinates from first column.\n");
		fprintf (stderr, "\t   If no domain is given we assume no time, i.e. only data columns are present.\n");
		fprintf (stderr, "\t   This choice also implies -Ca.\n");
		GMT_explain_option ('V');
		GMT_explain_option ('i');
		GMT_explain_option ('n');
		GMT_explain_option ('o');
		GMT_explain_option ('n');
		GMT_explain_option ('f');
		exit (EXIT_FAILURE);
	}
	
	if (! ((argc > 1 && argv[argc-2][0] == '=') || argv[argc-1][0] == '=')) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Usage is <operations> = outfile\n", GMT_program);
		exit (EXIT_FAILURE);
	}

	GMT_put_history (argc, argv);	/* Update .gmtcommands4 */

	/* Determine if the user have files with stupid names that will cause conflict */

	for (op = 0; op < N_OPERATORS; op++) {
		if (!access (operator[op], R_OK))
			fprintf (stderr, "%s Warning: Your file %s may be confused with a %s operator!\n", GMT_program, operator[op], GMT_program);
	}
	for (i = 0; i < N_SPECIAL_ARGS; i++) {
		if (!access (special_arg[i], R_OK))
			fprintf (stderr, "%s Warning: Your file %s may be confused with a %s operator!\n", GMT_program, special_arg[i], GMT_program);
	}

	GMT_hash_init (localhashnode, operator, HASH_SIZE, N_OPERATORS);

	for (i = 0; i < STACK_SIZE; i++) {
		constant[i] = FALSE;
		factor[i] = 0.0;
		tbl[i].n_col = tbl[i].n_row = 0;
		stack[i] = (double **)NULL;
	}
		
	if (!strcmp (argv[argc-2], "=")) {	/* Have  = file */
		outfile = argv[argc-1];
		last_arg = argc - 2;
	}
	else {
		outfile = NULL;
		last_arg = argc - 1;
	}
	
	/* Must first scan command line for -b, -f, -H switches before reading any file */

	for (arg = 1; arg < last_arg; arg++) {
		if (!strncmp (argv[arg], "-H", 2)) error += GMT_get_common_args (argv[arg], NULL, NULL, NULL, NULL);
		if (!strncmp (argv[arg], "-b", 2)) error += GMT_get_common_args (argv[arg], NULL, NULL, NULL, NULL);
		if (!strncmp (argv[arg], "-f", 2)) error += GMT_get_common_args (argv[arg], NULL, NULL, NULL, NULL);
		if (!strcmp (argv[arg], "-T")) GMT_io.skip_if_NaN[0] = GMT_io.skip_if_NaN[1] = FALSE;	/* Turn off default GMT NaN-handling in x/y columns */
	}
	/* Get header from one file so we can allocate space */
	
	for (arg = 1; nm == 0 && arg < last_arg; arg++) {
	
		if (argv[arg][0] == '-') continue;	/* Command line option */
		if (decode_argument (argv[arg]) != ARG_IS_FILE) continue;
		
		strcpy (file, argv[arg]);
		GMT_read_table (argv[arg], &header, &start, t_col);
		
		nm = header.n_row * header.n_col;
		n_col = header.n_col;
		got_t_from_file = TRUE;
		use_t_col = t_col;
	}
	
	/* Scan command line for -A, -T, -N, -Q, -S, -V */

	for (arg = 1; arg < last_arg; arg++) {
		if (argv[arg][0] == '-') {

			switch (argv[arg][1]) {

				case 'V':
					error += GMT_get_common_args (argv[arg], NULL, NULL, NULL, NULL);
					break;

				case 'A':	/* y(x) table for LSQFIT operations */
					GMT_read_table (&argv[arg][2], &rhs_header, &rhs, 0);	/* Always store as t and f(t) in cols 0 and 1 */
					if (rhs_header.n_col != 2) {
						fprintf (stderr, "%s: GMT SYNTAX ERROR:  -A must take a file with 2 (t,f(t)) colums\n", GMT_program);
						exit (EXIT_FAILURE);
					}
					set_a = TRUE;
					break;
				case 'N':
					sscanf (&argv[arg][2], "%d/%d", &n_col, &t_col);
					GMT_io.skip_if_NaN[0] = GMT_io.skip_if_NaN[1] = FALSE;
					GMT_io.skip_if_NaN[t_col] = TRUE;
					set_n = TRUE;
					break;
				case 'Q':	/* Quick for -Ca -N1/0 -T0/0/1 */
					n_col = 1;
					t_col = 0;
					set_n = set_t = set_q = TRUE;
					header.t_min = header.t_max = 0;
					header.t_inc = 1.0;
					break;
				case 'S':	/* Only want one row (first or last) */
					first_last_all = (argv[arg][2] == 'l' || argv[arg][2] == 'L') ? +1 : -1;
					break;
				case 'T':	/* Either get a file with time coordinate or a min/max/dt setting */
					if (!(access (&argv[arg][2], R_OK))) {	/* File can be opened */
						if (got_t_from_file) {
							fprintf (stderr, "%s: GMT SYNTAX ERROR:  Cannot use -T when data files are specified\n", GMT_program);
							exit (EXIT_FAILURE);
						}
						GMT_read_table (&argv[arg][2], &header, &start, 0);
						use_t_col = 0;
						got_t_from_file = set_t_table = TRUE;
					}
					else {
						i = sscanf (&argv[arg][2], "%lf/%lf/%lf", &header.t_min, &header.t_max, &header.t_inc);
						if (i == 3)
							set_t = TRUE;
						else
							no_time = TRUE;
					}
					break;
			}
		}
	}
	if (GMT_io.binary[GMT_IN] && gmtdefs.io_header[GMT_IN]) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR.  Binary input data cannot have header -H\n", GMT_program);
		error++;
	}
	if (set_n && (n_col <= 0 || t_col < 0 || t_col >= n_col)) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  -N must have positive n_cols and 0 <= t_col < n_col\n", GMT_program);
		exit (EXIT_FAILURE);
	}
	if (nm && set_t) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Cannot use -T when data files are specified\n", GMT_program);
		exit (EXIT_FAILURE);
	}
	if (set_t && !set_q) {
		switch (GMT_minmaxinc_verify (header.t_min, header.t_max, header.t_inc, SMALL)) {
			case 1:
				fprintf (stderr, "%s: GMT SYNTAX ERROR -T:  (max - min) is not a whole multiple of inc\n", GMT_program);
				exit (EXIT_FAILURE);
				break;
			case 2:
				if (header.t_inc != 1.0) {	/* Allow for somebody explicitly saying -T0/0/1 */
					fprintf (stderr, "%s: GMT SYNTAX ERROR -T:  (max - min) is <= 0\n", GMT_program);
					exit (EXIT_FAILURE);
				}
				break;
			case 3:
				fprintf (stderr, "%s: GMT SYNTAX ERROR -T:  inc is <= 0\n", GMT_program);
				exit (EXIT_FAILURE);
				break;
			default:	/* OK */
				break;
		}
		
		header.n_row = irint ((header.t_max - header.t_min) / header.t_inc) + 1;
		header.n_col = n_col;
		nm = header.n_row * header.n_col;
	}
	if (set_a) {	/* Get number of rows and time from the file, but not n_cols (that takes -N, which defaults to 2) */
		header.n_row = rhs_header.n_row;
		header.n_col = n_col;
		header.t_min = rhs_header.t_min;
		header.t_max = rhs_header.t_max;
		header.t_inc = rhs_header.t_inc;
		nm = header.n_row * header.n_col;
	}
	if (set_q) header.n_row = header.n_col = nm = 1;
	if (set_t_table) {
		header.n_col = n_col;
		nm = header.n_row * header.n_col;
	}
	if (nm == 0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Expression must contain at least one table file or -T [and -N]\n", GMT_program);
		exit (EXIT_FAILURE);
	}

	new_table (&stack[0], header.n_col, header.n_row);

	skip = (int *) GMT_memory (VNULL, (size_t)header.n_col, sizeof (int), GMT_program);
	if (!no_time) skip[t_col] = (set_q) ? FALSE : TRUE;

	/* Get t vector */

	t_coordinates = (double *) GMT_memory (VNULL, (size_t)header.n_row, sizeof (double), GMT_program);
	if (got_t_from_file) {
		memcpy ((void *)t_coordinates, (void *)start[use_t_col], (size_t)(header.n_row * sizeof (double)));
		header.t_inc = (start[use_t_col][header.n_row-1] - start[use_t_col][0]) / (header.n_row - 1);
		t_noise = fabs (SMALL * header.t_inc);
		for (i = 1; i < header.n_row && !irregular; i++) if (fabs (fabs (t_coordinates[i] - t_coordinates[i-1]) - fabs (header.t_inc)) > t_noise) irregular = TRUE;
		for (i = 0; i < header.n_col; i++) GMT_free ((void *)start[i]);
		GMT_free ((void *)start);
	}
	else {
		for (i = 0; i < header.n_row; i++) t_coordinates[i] = (i == (header.n_row-1)) ? header.t_max: header.t_min + i * header.t_inc;
		t_noise = fabs (SMALL * header.t_inc);
	}
	memcpy ((void *)stack[0][t_col], (void *)t_coordinates, (size_t)(header.n_row * sizeof (double)));
	if (set_a) {
		memcpy ((void *)stack[0][n_col-1], (void *)rhs[1], (size_t)(header.n_row * sizeof (double)));
		GMT_free ((void *)rhs[0]);
		GMT_free ((void *)rhs[1]);
		nstack = 1;
	}
	else
		nstack = 0;

	gmtmath_init (call_operator, consumed_operands, produced_operands);

	
	for (arg = 1; !error && arg < last_arg; arg++) {
	
		/* First check if we should skip optional arguments */

		if (!(strncmp (argv[arg], "-T", 2) && strncmp (argv[arg], "-b", 2) && strncmp (argv[arg], "-f", 2) && strncmp (argv[arg], "-N", 2))) continue;
		if (!(strncmp (argv[arg], "-H", 2)  && strncmp (argv[arg], "-Q", 2) && strncmp (argv[arg], "-S", 2) && strncmp (argv[arg], "-V", 2))) continue;
		if (!strncmp (argv[arg], "-A", 2)) continue;

		if (!strncmp (argv[arg], "-C", 2)) {	/* Change affected columns */
			decode_columns (&argv[arg][2], skip, n_col, t_col);
			continue;
		}
		
		if ((type = decode_argument (argv[arg])) != ARG_IS_OPERATOR) {	/* File name or factor */
		
			if (nstack == STACK_SIZE) {	/* Stack overflow */
				error = TRUE;
				continue;
			}

			if (type == ARG_IS_NUMBER) {
				constant[nstack] = TRUE;
				ok = sscanf (argv[arg], "%lf", &factor[nstack]);
				error = !ok;
				if (gmtdefs.verbose) fprintf (stderr, "%lg ", factor[nstack]);
				nstack++;
				continue;
			}
			else if (type == ARG_IS_PI) {
				constant[nstack] = TRUE;
				factor[nstack] = M_PI;
				if (gmtdefs.verbose) fprintf (stderr, "%lg ", factor[nstack]);
				nstack++;
				continue;
			}
			else if (type == ARG_IS_E) {
				constant[nstack] = TRUE;
				factor[nstack] = M_E;
				if (gmtdefs.verbose) fprintf (stderr, "%lg ", factor[nstack]);
				nstack++;
				continue;
			}

			/* Here we need a matrix */

			if (!stack[nstack]) new_table (&stack[nstack], header.n_col, header.n_row);

			constant[nstack] = FALSE;

			if (type == ARG_IS_T_MATRIX) {	/* Need to set up matrix of t-values */
				if (no_time) {
					fprintf (stderr, "%s: T is not defined for plain data files!\n", GMT_program);
					exit (EXIT_FAILURE);
				}
				if (gmtdefs.verbose) fprintf (stderr, "T ");
				for (j = 0; j < header.n_col; j++) memcpy ((void *)stack[nstack][j], (void *)t_coordinates, (size_t)(header.n_row * sizeof (double)));
			}
			else if (type == ARG_IS_FILE) {		/* Filename given */
				if (gmtdefs.verbose) fprintf (stderr, "%s ", argv[arg]);
				GMT_read_table (argv[arg], &tbl[nstack], &stack[nstack], t_col);
				if (tbl[nstack].n_row != header.n_row || tbl[nstack].n_col != header.n_col) {
					fprintf (stderr, "%s: tables not of same size!\n", GMT_program);
					exit (EXIT_FAILURE);
				}
				else if (!no_time && (fabs (tbl[nstack].t_min - header.t_min) > t_noise || fabs (tbl[nstack].t_max - header.t_max) > t_noise)) {
					fprintf (stderr, "%s: tables do not cover the same domain!\n", GMT_program);
					exit (EXIT_FAILURE);
				}
			}
			nstack++;
			continue;
		}
		
		/* Here we have an operator */
		
		if ((op = get_operator (argv[arg])) < 0) {
			fprintf (stderr, "%s: GMT SYNTAX ERROR:  Unrecognized operator %s\n", GMT_program, argv[arg]);
			exit (EXIT_FAILURE);
		}

		if (!strncmp (argv[arg], "ROOTS", 5) && !((arg+1) == last_arg && argv[arg+1][0] == '=')) {
			fprintf (stderr, "%s: GMT SYNTAX ERROR:  Only = may follow operator ROOTS\n", GMT_program);
			exit (EXIT_FAILURE);
		}
		
		if ((new_stack = nstack - consumed_operands[op] + produced_operands[op]) >= STACK_SIZE) {
			error = TRUE;
			continue;
		}

		if (nstack < consumed_operands[op]) {
			fprintf (stderr, "%s: GMT SYNTAX ERROR:  Operation \"%s\" requires %d operands\n", GMT_program, operator[op], consumed_operands[op]);
			exit (EXIT_FAILURE);
		}
		
		if (gmtdefs.verbose) fprintf (stderr, "%s ", operator[op]);

		for (i = produced_operands[op] - consumed_operands[op]; i > 0; i--) {
			 if (stack[nstack+i-1])	continue;

			/* Must make space for more */

			new_table (&stack[nstack+i-1], header.n_col, header.n_row);
		}

		/* If operators operates on constants only we may have to make space as well */

		for (j = 0, i = nstack - consumed_operands[op]; j < produced_operands[op]; j++, i++) {
			if (constant[i] && !stack[i]) new_table (&stack[i], header.n_col, header.n_row);
		}

		if (!strcmp (operator[op], "LSQFIT")) {	/* Special case, solve LSQ system and exit */
			solve_LSQFIT (stack, nstack - 1, n_col, header.n_row, skip, outfile);
			exit (EXIT_SUCCESS);
		}

		for (j = 0; j < n_col; j++) {
			if (skip[j]) continue;
			(*call_operator[op]) (stack, constant, factor, nstack - 1, j, header.n_row);	/* Do it */
		}

		nstack = new_stack;
		
		for (i = 1; i <= produced_operands[op]; i++)
			constant[nstack-i] = FALSE;	/* Now filled with table */
	}
	
	if (error && !ok) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Unable to decode constant %s (File not found?)\n", GMT_program, argv[i-1]);
		exit (EXIT_FAILURE);
	}
	
	if (error) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Stack overflow (%s)\n", GMT_program, argv[i-1]);
		exit (EXIT_FAILURE);
	}
	
	if (gmtdefs.verbose) {
		(outfile) ? fprintf (stderr, "= %s", outfile) : fprintf (stderr, "= <stdout>");
	}

	if (new_stack < 0 && constant[0]) {	/* Only a constant provided, set table accordingly */
		for (j = 0; j < header.n_col; j++) {
			if (skip[j]) continue;
			for (i = 0; i < header.n_row; i++) stack[0][j][i] = factor[0];
		}
	}

	if (roots_found) {	/* Special treatment of root finding */
		if (outfile) {
			if ((fp = GMT_fopen (outfile, "w")) == NULL) {
				fprintf (stderr, "%s: GMT ERROR:  Could not create file (%s)\n", GMT_program, outfile);
				exit (EXIT_FAILURE);
			}
		}
		else
			fp = GMT_stdout;
		for (i = 0; i < n_roots; i++) GMT_output (fp, 1, &stack[0][r_col][i]);
		if (fp != GMT_stdout) GMT_fclose (fp);
	}
	else
		GMT_write_table (outfile, &header, stack[0], first_last_all);
	
	for (i = 0; i < STACK_SIZE; i++) if (stack[i]) {
		for (j = 0; j < header.n_col; j++) GMT_free ((void *)stack[i][j]);
		GMT_free ((void *)stack[i]);
	}
	GMT_free ((void *)t_coordinates);
	GMT_free ((void *)skip);
	
	if (gmtdefs.verbose) fprintf (stderr, "\n");

	if (nstack > 1) fprintf (stderr, "%s: Warning: %d more operands left on the stack!\n", GMT_program, nstack-1);

	GMT_end (argc, argv);
}

void GMT_read_table (char *file, struct TABLE_HEADER *h, double ***p, int t_col)
{
	int j, n = 0, n_alloc, n_expected_fields, n_fields;
	double *in, **table;
	char buffer[BUFSIZ];
	FILE *fp;
	
	n_expected_fields = (GMT_io.binary[GMT_IN]) ? GMT_io.ncol[GMT_IN] : BUFSIZ;

#ifdef SET_IO_MODE
	GMT_setmode (0);
#endif

	if (file && (fp = GMT_fopen (file, GMT_io.r_mode)) == NULL) {
		fprintf (stderr, "%s: Error opening file %s\n", GMT_program, file);
		exit (EXIT_FAILURE);
	}

	for (j = 0; gmtdefs.io_header[GMT_IN] && j < gmtdefs.n_header_recs; j++) {
		GMT_fgets (buffer, BUFSIZ, fp);
		if (very_first && j == 0) strcpy (head_record, buffer);
		very_first = FALSE;
	}

	GMT_input (fp, &n_expected_fields, &in);
	if (GMT_io.status & GMT_IO_EOF) {
		fprintf (stderr, "%s: Error reading 1st record of file %s\n", GMT_program, file);
		exit (EXIT_FAILURE);
	}

	h->n_col = n_expected_fields;
	n_alloc = GMT_CHUNK;
	new_table (&table, n_expected_fields, n_alloc);

	do {
		for (j = 0; j < h->n_col; j++) table[j][n] = in[j];
		n++;
		if (n == n_alloc) {
			n_alloc += GMT_CHUNK;
			for (j = 0; j < h->n_col; j++) table[j] = (double *) GMT_memory ((void *)table[j], (size_t)n_alloc, sizeof (double), GMT_program);
		}
	} while ((n_fields = GMT_input (fp, &n_expected_fields, &in)) >= 0 && !(GMT_io.status & GMT_IO_EOF));

	h->t_min = table[t_col][0];
	h->t_max = table[t_col][n-1];
	h->t_inc = (h->t_max - h->t_min) / (n-1);

	if (p)
		*p = table;
	else {	/* Discard */
		for (j = 0; j < h->n_col; j++) GMT_free ((void *)table[j]);
		GMT_free ((void *)table);
	}
	h->n_row = n;
}

		
void GMT_write_table (char *file, struct TABLE_HEADER *h, double **p, int first_last_all)
{	/* first_last_all will write first [-1], last [+1], or all [0] records */
	int i, j, start, stop;
	double *out;
	FILE *fp;

	if (file && (fp = GMT_fopen (file, GMT_io.w_mode)) == NULL) {
		fprintf (stderr, "%s: Error creating file %s\n", GMT_program, file);
		exit (EXIT_FAILURE);
	}
	else if (!file)
		fp = GMT_stdout;

#ifdef SET_IO_MODE
	GMT_setmode (1);
#endif

	out = (double *) GMT_memory (VNULL, (size_t)h->n_col, sizeof (double), GMT_program);

	if (gmtdefs.io_header[GMT_OUT] && !GMT_io.binary[GMT_OUT]) {
		fprintf (fp, "%s", head_record);
		for (i = 1; i < gmtdefs.n_header_recs; i++) fprintf (fp, "# gmtmath header record\n");
	}
	
	start = (first_last_all == +1) ? h->n_row - 1 : 0;
	stop = (first_last_all == -1) ? 1 : h->n_row;
	for (i = start; i < stop; i++) {
		for (j = 0; j < h->n_col; j++) out[j] = p[j][i];
		GMT_output (fp, h->n_col, out);
	}

	if (fp != GMT_stdout) GMT_fclose (fp);
	GMT_free ((void *)out);
}

void new_table (double ***s, int n_col, int n)
{
	int j;
	double **p;
	p = (double **) GMT_memory (VNULL, (size_t)n_col, sizeof (double *), GMT_program);
	for (j = 0; j < n_col; j++) p[j] = (double *) GMT_memory (VNULL, (size_t)n, sizeof (double), GMT_program);
	*s = p;
}

void decode_columns (char *txt, int *skip, int n_col, int t_col)
{
	int i, start, stop;
	char line[BUFSIZ], *p;

	if (!txt) {	/* Reset to default */
		for (i = 0; i < n_col; i++) skip[i] = FALSE;
		skip[t_col] = TRUE;
	}
	if (txt[0] == 'r' && txt[1] == '\0') {	/* Reverse all settings */
		for (i = 0; i < n_col; i++) skip[i] = !skip[i];
	}
	else if (txt[0] == 'a') {	/* Select all columns */
		for (i = 0; i < n_col; i++) skip[i] = FALSE;
	}
	else {	/* Set the selected columns */
		for (i = 0; i < n_col; i++) skip[i] = TRUE;
		strcpy (line, txt);
		p = strtok (line, ",");
		while (p) {
			if (strchr (p, '-'))
				sscanf (p, "%d-%d", &start, &stop);
			else {
				sscanf (p, "%d", &start);
				stop = start;
			}
			stop = MIN (stop, n_col-1);
			for (i = start; i <= stop; i++) skip[i] = FALSE;
			p = strtok (NULL, ",");
		}
	}
}


/* -----------------------------------------------------------------
 *              Definitions of all operator functions
 * -----------------------------------------------------------------*/

void table_ABS (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, operand == 0!\n", GMT_program);
	if (constant[last]) a = fabs (factor[last]);
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : fabs (stack[last][col][i]);
}

void table_ACOS (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (gmtdefs.verbose && constant[last] && fabs (factor[last]) > 1.0) fprintf (stderr, "%s: Warning, |operand| > 1 for ACOS!\n", GMT_program);
	if (constant[last]) a = d_acos (factor[last]);
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : d_acos (stack[last][col][i]);
}

void table_ACOSH (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (gmtdefs.verbose && constant[last] && fabs (factor[last]) > 1.0) fprintf (stderr, "%s: Warning, operand < 1 for ACOSH!\n", GMT_program);
	if (constant[last]) a = acosh (factor[last]);
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : acosh (stack[last][col][i]);
}

void table_ADD (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	if (gmtdefs.verbose && constant[prev] && factor[prev] == 0.0) fprintf (stderr, "%s: Warning, operand one == 0!\n", GMT_program);
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, operand two == 0!\n", GMT_program);
	for (i = 0; i < n_row; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][col][i];
		b = (constant[last]) ? factor[last] : stack[last][col][i];
		stack[prev][col][i] = a + b;
	}
}

void table_AND (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	for (i = 0; i < n_row; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][col][i];
		b = (constant[last]) ? factor[last] : stack[last][col][i];
		stack[prev][col][i] = (GMT_is_dnan (a)) ? b : a;
	}
}

void table_ASIN (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (gmtdefs.verbose && constant[last] && fabs (factor[last]) > 1.0) fprintf (stderr, "%s: Warning, |operand| > 1 for ASIN!\n", GMT_program);
	if (constant[last]) a = d_asin (factor[last]);
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : d_asin (stack[last][col][i]);
}

void table_ASINH (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (constant[last]) a = asinh (factor[last]);
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : asinh (stack[last][col][i]);
}

void table_ATAN (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (constant[last]) a = atan (factor[last]);
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : atan (stack[last][col][i]);
}

void table_ATAN2 (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	if (gmtdefs.verbose && constant[prev] && factor[prev] == 0.0) fprintf (stderr, "%s: Warning, operand one == 0 for ATAN2!\n", GMT_program);
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, operand two == 0 for ATAN2!\n", GMT_program);
	for (i = 0; i < n_row; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][col][i];
		b = (constant[last]) ? factor[last] : stack[last][col][i];
		stack[prev][col][i] = d_atan2 (a, b);
	}
}

void table_ATANH (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (gmtdefs.verbose && constant[last] && fabs (factor[last]) >= 1.0) fprintf (stderr, "%s: Warning, |operand| >= 1 for ATANH!\n", GMT_program);
	if (constant[last]) a = atanh (factor[last]);
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : atanh (stack[last][col][i]);
}

void table_BEI (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (constant[last]) a = GMT_bei (fabs (factor[last]));
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : GMT_bei (fabs (stack[last][col][i]));
}

void table_BER (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (constant[last]) a = GMT_ber (fabs (factor[last]));
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : GMT_ber (fabs (stack[last][col][i]));
}

void table_CEIL (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (constant[last]) a = ceil (factor[last]);
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : ceil (stack[last][col][i]);
}

void table_CHICRIT (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	if (gmtdefs.verbose && constant[prev] && factor[prev] == 0.0) fprintf (stderr, "%s: Warning, operand one == 0 for CHICRIT!\n", GMT_program);
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, operand two == 0 for CHICRIT!\n", GMT_program);
	for (i = 0; i < n_row; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][col][i];
		b = (constant[last]) ? factor[last] : stack[last][col][i];
		stack[prev][col][i] = GMT_chi2crit (a, b);
	}
}

void table_CHIDIST (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	if (gmtdefs.verbose && constant[prev] && factor[prev] == 0.0) fprintf (stderr, "%s: Warning, operand one == 0 for CHIDIST!\n", GMT_program);
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, operand two == 0 for CHIDIST!\n", GMT_program);
	for (i = 0; i < n_row; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][col][i];
		b = (constant[last]) ? factor[last] : stack[last][col][i];
		GMT_chi2 (a, b, &stack[prev][col][i]);
	}
}

void table_COL (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i, k, prev;
	
	if (!constant[last]) {
		fprintf (stderr, "%s: Error, argument to COL must be a constant column  number (0 <= k < n_col)!\n", GMT_program);
		exit (EXIT_FAILURE);
	}
	prev = last - 1;
	k = irint (factor[last]);
	for (i = 0; i < n_row; i++) {
		stack[last][col][i] = stack[prev][k][i];
	}
}


void table_COS (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (constant[last]) a = cos (factor[last]);
	for (i = 0; i < n_row; i++) {
		stack[last][col][i] = (constant[last]) ? a : cos (stack[last][col][i]);
	}
}

void table_COSD (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (constant[last]) a = cosd (factor[last]);
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : cosd (stack[last][col][i]);
}

void table_COSH (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (constant[last]) a = cosh (factor[last]);
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : cosh (stack[last][col][i]);
}

void table_DDT (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double c, left, next_left;

	/* Central 1st difference in t */

	if (irregular) fprintf (stderr, "%s: Warning, DDT called on irregularly spaced data (not supported)!\n", GMT_program);
	if (gmtdefs.verbose && constant[last]) fprintf (stderr, "%s: Warning, operand to DDT is constant!\n", GMT_program);

	c = 0.5 / header.t_inc;
	next_left = 2.0 * stack[last][col][0] - stack[last][col][1];
	for (i = 0; i < n_row - 1; i++) {
		left = next_left;
		next_left = stack[last][col][i];
		stack[last][col][i] = (constant[last]) ? 0.0 : c * (stack[last][col][i+1] - left);
	}
	stack[last][col][i] = (constant[last]) ? 0.0 : 2.0 * c * (stack[last][col][i] - next_left);
}

void table_D2DT2 (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double c, left, next_left;

	/* Central 2nd difference in t */

	if (irregular) fprintf (stderr, "%s: Warning, D2DT2 called on irregularly spaced data (not supported)!\n", GMT_program);
	if (gmtdefs.verbose && constant[last]) fprintf (stderr, "%s: Warning, operand to D2DT2 is constant!\n", GMT_program);

	c = 1.0 / (header.t_inc * header.t_inc);
	next_left = stack[last][col][0];
	stack[last][col][0] = 0.0;
	for (i = 1; i < n_row - 1; i++) {
		left = next_left;
		next_left = stack[last][col][i];
		stack[last][col][i] = (constant[last]) ? 0.0 : c * (stack[last][col][i+1] - 2 * stack[last][col][i] + left);
	}
	stack[last][col][i] = 0.0;
}

void table_D2R (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (constant[last]) a = factor[last] * D2R;
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : stack[last][col][i] * D2R;
}

void table_DILOG (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (constant[last]) a = GMT_dilog (factor[last]);
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : GMT_dilog (stack[last][col][i]);
}

void table_DIV (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	
	if (constant[last] && factor[last] == 0.0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Cannot divide by zero\n", GMT_program);
		exit (EXIT_FAILURE);
	}
	if (constant[last]) {	/* Turn divide into multiply */
		a = factor[last];	/* Save old factor */
		factor[last] = 1.0 / factor[last];
		table_MUL (stack, constant, factor, last, col, n_row);
		factor[last] = a;	/* Restore factor to original value */
		return;
	}
	
	for (i = 0; i < n_row; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][col][i];
		b = (constant[last]) ? factor[last] : stack[last][col][i];
		stack[prev][col][i] = a / b;
	}
}

void table_DUP (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int next, i;
	
	next = last + 1;
	factor[next] = factor[last];
	constant[next] = constant[last];
	if (constant[last])
		for (i = 0; i < n_row; i++) stack[next][col][i] = stack[last][col][i] = factor[next];
	else
		memcpy ((void *)stack[next][col], (void *)stack[last][col], (size_t)(n_row * sizeof (double)));
}

void table_ERF (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (constant[last]) a = erf (factor[last]);
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : erf (stack[last][col][i]);
}

void table_ERFC (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (constant[last]) a = erfc (factor[last]);
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : erfc (stack[last][col][i]);
}

void table_ERFINV (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (constant[last]) a = GMT_erfinv (factor[last]);
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : GMT_erfinv (stack[last][col][i]);
}

void table_EQ (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	for (i = 0; i < n_row; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][col][i];
		b = (constant[last]) ? factor[last] : stack[last][col][i];
		stack[prev][col][i] = (double)(a == b);
	}
}

void table_EXCH (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i, prev;
	
	prev = last - 1;
	
	for (i = 0; i < n_row; i++) {
		if (constant[last]) stack[last][col][i] = factor[last];
		if (constant[prev]) stack[prev][col][i] = factor[prev];
		d_swap (stack[last][col][i], stack[prev][col][i]);
	}
	i_swap (constant[last], constant[prev]);
	d_swap (factor[last], factor[prev]);
}

void table_EXP (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, operand == 0!\n", GMT_program);
	if (constant[last]) a = exp (factor[last]);
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : exp (stack[last][col][i]);
}

void table_FCRIT (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i, nu1, nu2, prev1, prev2;
	double alpha;
	
	prev1 = last - 1;
	prev2 = last - 2;
	if (gmtdefs.verbose && constant[prev2] && factor[prev2] == 0.0) fprintf (stderr, "%s: Warning, operand one == 0 for FCRIT!\n", GMT_program);
	if (gmtdefs.verbose && constant[prev1] && factor[prev1] == 0.0) fprintf (stderr, "%s: Warning, operand two == 0 for FCRIT!\n", GMT_program);
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, operand three == 0 for FCRIT!\n", GMT_program);
	for (i = 0; i < n_row; i++) {
		alpha = (constant[prev2]) ? factor[prev2] : stack[prev2][col][i];
		nu1 = irint ((double)((constant[prev1]) ? factor[prev1] : stack[prev1][col][i]));
		nu2 = irint ((double)((constant[last]) ? factor[last] : stack[last][col][i]));
		stack[prev2][col][i] = GMT_Fcrit (alpha, (double)nu1, (double)nu2);
	}
}

void table_FDIST (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i, nu1, nu2, prev1, prev2;
	double F, chisq1, chisq2 = 1.0;
	
	prev1 = last - 1;
	prev2 = last - 2;
	if (gmtdefs.verbose && constant[prev1] && factor[prev1] == 0.0) fprintf (stderr, "%s: Warning, operand two == 0 for FDIST!\n", GMT_program);
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, operand three == 0 for FDIST!\n", GMT_program);
	for (i = 0; i < n_row; i++) {
		F = (constant[prev2]) ? factor[prev2] : stack[prev2][col][i];
		nu1 = irint ((double)((constant[prev1]) ? factor[prev1] : stack[prev1][col][i]));
		nu2 = irint ((double)((constant[last]) ? factor[last] : stack[last][col][i]));
		/* Since GMT_f_q needs chisq1 and chisq2, we set chisq2 = 1 and solve for chisq1 */
		chisq1 = F * nu1 / nu2;
		(void) GMT_f_q (chisq1, nu1, chisq2, nu2, &stack[prev2][col][i]);
	}
}

void table_FLIPUD (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i, k;
	/* Reverse the order of points in a column */	
	if (constant[last]) return;
	for (i = 0, k = n_row-1; i < n_row/2; i++, k--) d_swap (stack[last][col][i], stack[last][col][k]);
}

void table_FLOOR (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (constant[last]) a = floor (factor[last]);
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : floor (stack[last][col][i]);
}

void table_FMOD (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	if (gmtdefs.verbose && constant[prev] && factor[prev] == 0.0) fprintf (stderr, "%s: Warning, operand one == 0!\n", GMT_program);
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, operand two == 0!\n", GMT_program);
	for (i = 0; i < n_row; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][col][i];
		b = (constant[last]) ? factor[last] : stack[last][col][i];
		stack[prev][col][i] = fmod (a, b);
	}
}

void table_GE (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	for (i = 0; i < n_row; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][col][i];
		b = (constant[last]) ? factor[last] : stack[last][col][i];
		stack[prev][col][i] = (double)(a >= b);
	}
}

void table_GT (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	for (i = 0; i < n_row; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][col][i];
		b = (constant[last]) ? factor[last] : stack[last][col][i];
		stack[prev][col][i] = (double)(a > b);
	}
}

void table_HYPOT (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	if (gmtdefs.verbose && constant[prev] && factor[prev] == 0.0) fprintf (stderr, "%s: Warning, operand one == 0!\n", GMT_program);
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, operand two == 0!\n", GMT_program);
	for (i = 0; i < n_row; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][col][i];
		b = (constant[last]) ? factor[last] : stack[last][col][i];
		stack[prev][col][i] = hypot (a, b);
	}
}

void table_I0 (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (constant[last]) a = GMT_i0 (factor[last]);
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : GMT_i0 (stack[last][col][i]);
}

void table_I1 (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (constant[last]) a = GMT_i1 (factor[last]);
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : GMT_i1 (stack[last][col][i]);
}

void table_IN (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i, prev, order;
	BOOLEAN simple = FALSE;
	double b;
	
	prev = last - 1;
	if (gmtdefs.verbose && constant[last] && factor[last] < 0.0) fprintf (stderr, "%s: Warning, order < 0 for IN!\n", GMT_program);
	if (gmtdefs.verbose && constant[last] && fabs (rint(factor[last]) - factor[last]) > SMALL) fprintf (stderr, "%s: Warning, order not an integer for IN!\n", GMT_program);
	if (constant[last]) order = irint (fabs (factor[last]));
	if (constant[last] && constant[prev]) {
		b = GMT_in (order, fabs (factor[prev]));
		simple = TRUE;
	}
	for (i = 0; i < n_row; i++) {
		if (simple)
			stack[prev][col][i] = b;
		else {
			if (!constant[last]) order = irint (fabs (stack[last][col][i]));
			stack[prev][col][i] = GMT_in (order, fabs (stack[prev][col][i]));
		}
	}
}

void table_INT (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double f, left, right, sum;
	
	if (constant[last]) {	/* Trivial case */
		sum = factor[last] * header.t_inc;
		for (i = 0; i < n_row; i++) stack[last][col][i] = i * sum;
		return;
	}
	
	/* We use dumb trapezoidal rule - one day we will replace with more sophisticated rules */

	left = stack[last][col][0];
	stack[last][col][0] = sum = 0.0;
	if (!irregular) f = 0.5 * header.t_inc;
	for (i = 1; i < n_row; i++) {	/* Dumb trapezoidal rule */
		if (irregular) f = 0.5 * (t_coordinates[i] - t_coordinates[i-1]);
		right = stack[last][col][i];
		sum += f * (left + right);
		stack[last][col][i] = sum;
		left = right;
	}
}

void table_INV (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (constant[last] && factor[last] == 0.0) {
		fprintf (stderr, "%s: Error, Cannot take inverse of zero!\n", GMT_program);
		exit (EXIT_FAILURE);
	}
	if (constant[last]) a = 1.0 / factor[last];
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : 1.0 / stack[last][col][i];
}

void table_ISNAN (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (constant[last]) a = (double)GMT_is_dnan (factor[last]);
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : (double)GMT_is_dnan (stack[last][col][i]);
}

void table_J0 (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (constant[last]) a = j0 (factor[last]);
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : j0 (stack[last][col][i]);
}

void table_J1 (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (constant[last]) a = j1 (fabs (factor[last]));
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : j1 (fabs (stack[last][col][i]));
}

void table_JN (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i, prev, order;
	BOOLEAN simple = FALSE;
	double b;
	
	prev = last - 1;
	if (gmtdefs.verbose && constant[last] && factor[last] < 0.0) fprintf (stderr, "%s: Warning, order < 0 for JN!\n", GMT_program);
	if (gmtdefs.verbose && constant[last] && fabs (rint(factor[last]) - factor[last]) > SMALL) fprintf (stderr, "%s: Warning, order not an integer for JN!\n", GMT_program);
	if (constant[last]) order = irint (fabs (factor[last]));
	if (constant[last] && constant[prev]) {
		b = jn (order, fabs (factor[prev]));
		simple = TRUE;
	}
	for (i = 0; i < n_row; i++) {
		if (simple)
			stack[prev][col][i] = b;
		else {
			if (!constant[last]) order = irint (fabs (stack[last][col][i]));
			stack[prev][col][i] = jn (order, fabs (stack[prev][col][i]));
		}
	}
}

void table_K0 (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (constant[last]) a = GMT_k0 (factor[last]);
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : GMT_k0 (stack[last][col][i]);
}

void table_K1 (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (constant[last]) a = GMT_k1 (factor[last]);
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : GMT_k1 (stack[last][col][i]);
}

void table_KN (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i, prev, order;
	BOOLEAN simple = FALSE;
	double b;
	
	prev = last - 1;
	if (gmtdefs.verbose && constant[last] && factor[last] < 0.0) fprintf (stderr, "%s: Warning, order < 0 for KN!\n", GMT_program);
	if (gmtdefs.verbose && constant[last] && fabs (rint(factor[last]) - factor[last]) > SMALL) fprintf (stderr, "%s: Warning, order not an integer for KN!\n", GMT_program);
	if (constant[last]) order = irint (fabs (factor[last]));
	if (constant[last] && constant[prev]) {
		b = GMT_kn (order, fabs (factor[prev]));
		simple = TRUE;
	}
	for (i = 0; i < n_row; i++) {
		if (simple)
			stack[prev][col][i] = b;
		else {
			if (!constant[last]) order = irint (fabs (stack[last][col][i]));
			stack[prev][col][i] = GMT_kn (order, fabs (stack[prev][col][i]));
		}
	}
}

void table_KEI (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (constant[last]) a = GMT_kei (fabs (factor[last]));
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : GMT_kei (fabs (stack[last][col][i]));
}

void table_KER (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (constant[last]) a = GMT_ker (fabs (factor[last]));
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : GMT_ker (fabs (stack[last][col][i]));
}

void table_LE (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	for (i = 0; i < n_row; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][col][i];
		b = (constant[last]) ? factor[last] : stack[last][col][i];
		stack[prev][col][i] = (double)(a <= b);
	}
}

void table_LMSSCL (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double lmsscl, mode;
	
	if (constant[last]) {	/* Trivial case */
		memset ((void *)stack[last][col], 0, (size_t)(n_row * sizeof (double)));
		return;
	}

	qsort ((void *)stack[last][col], (size_t)n_row, sizeof (double), GMT_comp_double_asc);
	for (i = n_row; GMT_is_fnan (stack[last][col][i-1]) && i > 1; i--);
	if (i) {
		GMT_mode (stack[last][col], i, i/2, 0, GMT_mode_selection, &GMT_n_multiples, &mode);
		GMT_getmad (stack[last][col], i, mode, &lmsscl);
	}
	else
		lmsscl = GMT_d_NaN;

	for (i = 0; i < n_row; i++) stack[last][col][i] = lmsscl;
	if (GMT_n_multiples > 0) {
		fprintf (stderr, "%s: WARNING: %d Multiple modes found\n", GMT_program, GMT_n_multiples);
		GMT_n_multiples = 0;
	}
}

void table_LOG (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, argument to log = 0!\n", GMT_program);

	if (constant[last]) a = d_log (fabs (factor[last]));
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : d_log (fabs (stack[last][col][i]));
}

void table_LOG10 (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, argument to log10 = 0!\n", GMT_program);

	if (constant[last]) a = d_log10 (fabs (factor[last]));
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : d_log10 (fabs (stack[last][col][i]));
}

void table_LOG1P (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (gmtdefs.verbose && constant[last] && factor[last] < 0.0) fprintf (stderr, "%s: Warning, argument to log1p < 0!\n", GMT_program);

	if (constant[last]) a = d_log1p (fabs (factor[last]));
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : d_log1p (fabs (stack[last][col][i]));
}

void table_LOG2 (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, argument to log2 = 0!\n", GMT_program);

	if (constant[last]) a = d_log (fabs (factor[last])) * M_LN2_INV;
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : d_log (fabs (stack[last][col][i])) * M_LN2_INV;
}

void table_LOWER (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double low;

	if (constant[last]) {	/* Trivial case */
		for (i = 0; i < n_row; i++) stack[last][col][i] = factor[last];
		return;
	}
	
	for (i = 0, low = DBL_MAX; i < n_row; i++) {
		if (GMT_is_fnan (stack[last][col][i])) continue;
		if (stack[last][col][i] < low) low = stack[last][col][i];
	}
	for (i = 0; i < n_row; i++) if (!GMT_is_fnan (stack[last][col][i])) stack[last][col][i] = low;
}

void table_LRAND (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	if (constant[prev]) a = factor[prev];
	if (constant[last]) b = factor[last];
	for (i = 0; i < n_row; i++) {
		if (!constant[prev]) a = stack[prev][col][i];
		if (!constant[last]) b = stack[last][col][i];
		stack[prev][col][i] = a + b * GMT_lrand ();
	}
}

void table_LSQFIT (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	/* Dummy routine needed since the automatically generated include file will have table_LSQFIT
	 * with these parameters just like any other function.  However, when we find LSQFIT we will
	 * instead call solve_LSQFIT which can be found at the end of these functions */
}

void table_LT (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	for (i = 0; i < n_row; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][col][i];
		b = (constant[last]) ? factor[last] : stack[last][col][i];
		stack[prev][col][i] = (double)(a < b);
	}
}

void table_MAD (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double mad, med;
	
	if (constant[last]) {	/* Trivial case */
		memset ((void *)stack[last][col], 0, (size_t)(n_row * sizeof (double)));
		return;
	}

	qsort ((void *)stack[last][col], (size_t)n_row, sizeof (double), GMT_comp_double_asc);
	for (i = n_row; GMT_is_fnan (stack[last][col][i-1]) && i > 1; i--);
	if (i) {
		med = (i%2) ? stack[last][col][i/2] : 0.5 * (stack[last][col][(i-1)/2] + stack[last][col][i/2]);
		GMT_getmad (stack[last][col], i, med, &mad);
	}
	else
		mad = GMT_d_NaN;

	for (i = 0; i < n_row; i++) stack[last][col][i] = mad;
}

void table_MAX (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	for (i = 0; i < n_row; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][col][i];
		b = (constant[last]) ? factor[last] : stack[last][col][i];
		stack[prev][col][i] = MAX (a, b);
	}
}

void table_MEAN (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i, n_a = 0;
	double sum_a = 0.0;

	if (constant[last]) {	/* Trivial case */
		for (i = 0; i < n_row; i++) stack[last][col][i] = factor[last];
		return;
	}
	
	for (i = 0; i < n_row; i++) {
		if (GMT_is_fnan (stack[last][col][i])) continue;
		sum_a += stack[last][col][i];
		n_a++;
	}
	sum_a = (n_a) ? sum_a / n_a : 0.0;
	for (i = 0; i < n_row; i++) stack[last][col][i] = sum_a;
}

void table_MED (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double med;
	
	if (constant[last]) {	/* Trivial case */
		for (i = 0; i < n_row; i++) stack[last][col][i] = factor[last];
		return;
	}

	qsort ((void *)stack[last][col], (size_t)n_row, sizeof (double), GMT_comp_double_asc);
	for (i = n_row; GMT_is_fnan (stack[last][col][i-1]) && i > 1; i--);
	if (i)
		med = (i%2) ? stack[last][col][i/2] : 0.5 * (stack[last][col][(i-1)/2] + stack[last][col][i/2]);
	else
		med = GMT_d_NaN;

	for (i = 0; i < n_row; i++) stack[last][col][i] = med;
}

void table_MIN (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	for (i = 0; i < n_row; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][col][i];
		b = (constant[last]) ? factor[last] : stack[last][col][i];
		stack[prev][col][i] = MIN (a, b);
	}
}

void table_MODE (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double mode;
	
	if (constant[last]) {	/* Trivial case */
		for (i = 0; i < n_row; i++) stack[last][col][i] = factor[last];
		return;
	}

	qsort ((void *)stack[last][col], (size_t)n_row, sizeof (double), GMT_comp_double_asc);
	for (i = n_row; GMT_is_fnan (stack[last][col][i-1]) && i > 1; i--);
	if (i)
		GMT_mode (stack[last][col], i, i/2, 0, GMT_mode_selection, &GMT_n_multiples, &mode);
	else
		mode = GMT_d_NaN;

	for (i = 0; i < n_row; i++) stack[last][col][i] = mode;
	if (GMT_n_multiples > 0) {
		fprintf (stderr, "%s: WARNING: %d Multiple modes found\n", GMT_program, GMT_n_multiples);
		GMT_n_multiples = 0;
	}
}

void table_MUL (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	if (gmtdefs.verbose && constant[prev] && factor[prev] == 0.0) fprintf (stderr, "%s: Warning, operand one == 0!\n", GMT_program);
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, operand two == 0!\n", GMT_program);
	for (i = 0; i < n_row; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][col][i];
		b = (constant[last]) ? factor[last] : stack[last][col][i];
		stack[prev][col][i] = a * b;
	}
}

void table_NAN (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	if (constant[prev]) a = factor[prev];
	if (constant[last]) b = factor[last];
	for (i = 0; i < n_row; i++) {
		if (!constant[prev]) a = stack[prev][col][i];
		if (!constant[last]) b = stack[last][col][i];
		stack[prev][col][i] = ((a == b) ? GMT_d_NaN : a);
	}
}

void table_NEG (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, operand == 0!\n", GMT_program);
	if (constant[last]) a = -factor[last];
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : -stack[last][col][i];
}

void table_NEQ (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	for (i = 0; i < n_row; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][col][i];
		b = (constant[last]) ? factor[last] : stack[last][col][i];
		stack[prev][col][i] = (double)(a != b);
	}
}

void table_NRAND (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	if (constant[prev]) a = factor[prev];
	if (constant[last]) b = factor[last];
	for (i = 0; i < n_row; i++) {
		if (!constant[prev]) a = stack[prev][col][i];
		if (!constant[last]) b = stack[last][col][i];
		stack[prev][col][i] = a + b * GMT_nrand ();
	}
}

void table_OR (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	for (i = 0; i < n_row; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][col][i];
		b = (constant[last]) ? factor[last] : stack[last][col][i];
		stack[prev][col][i] = (GMT_is_dnan (a) || GMT_is_dnan (b)) ? GMT_d_NaN : a;
	}
}

void table_PLM (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i, prev, first, L, M;
	double a;
				/* last holds the order M */
	prev  = last - 1;	/* prev holds the degree L */
	first = prev - 1;	/* first holds the argument x = cos(colat) */

	if (!(constant[prev] && constant[last])) {
		fprintf (stderr, "%s: L and M must be constants in PLM!\n", GMT_program);
		exit (EXIT_FAILURE);
	}
	L = irint (factor[prev]);
	M = irint (factor[last]);
	if (gmtdefs.verbose && constant[first] && (factor[first] < -1.0 || factor[first] > 1.0)) fprintf (stderr, "%s: Warning, argument to PLM outside domain!\n", GMT_program);

	if (constant[first]) a = GMT_plm (L, M, factor[first]);
	for (i = 0; i < n_row; i++) stack[first][col][i] = (constant[first]) ? a : GMT_plm (L, M, stack[first][col][i]);
}

void table_POP (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	/* Dummy routine that does nothing but consume the top element of stack */
}

void table_POW (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;

	if (gmtdefs.verbose && constant[prev] && factor[prev] == 0.0) fprintf (stderr, "%s: Warning, operand one == 0!\n", GMT_program);
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, operand two == 0!\n", GMT_program);
	for (i = 0; i < n_row; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][col][i];
		b = (constant[last]) ? factor[last] : stack[last][col][i];
		stack[prev][col][i] = pow (a, b);
	}
}

void table_R2 (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	if (gmtdefs.verbose && constant[prev] && factor[prev] == 0.0) fprintf (stderr, "%s: Warning, operand one == 0!\n", GMT_program);
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, operand two == 0!\n", GMT_program);
	if (constant[prev]) factor[prev] *= factor[prev];
	if (constant[last]) factor[last] *= factor[last];
	for (i = 0; i < n_row; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][col][i] * stack[prev][col][i];
		b = (constant[last]) ? factor[last] : stack[last][col][i] * stack[last][col][i];
		stack[prev][col][i] = a + b;
	}
}

void table_R2D (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (constant[last]) a = factor[last] * R2D;
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : stack[last][col][i] * R2D;
}

void table_RAND (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	if (constant[prev]) a = factor[prev];
	if (constant[last]) b = factor[last];
	for (i = 0; i < n_row; i++) {
		if (!constant[prev]) a = stack[prev][col][i];
		if (!constant[last]) b = stack[last][col][i];
		stack[prev][col][i] = a + GMT_rand () * (b - a);
	}
}

void table_RINT (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (constant[last]) a = rint (factor[last]);
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : rint (stack[last][col][i]);
}

void table_ROTT (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i, shift, prev;
	double *z;
	
	if (!constant[last]) {
		fprintf (stderr, "%s: T-shift must be a constant in ROTT\n", GMT_program);
		exit (EXIT_FAILURE);
	}
	prev = last - 1;
	shift = irint (factor[last] / header.t_inc);
	if (constant[prev] || !shift) return;	/* Easy, constant or no shift */
	if (shift < 0) shift += n_row;		/* Same thing */
	
	z = (double *) GMT_memory (VNULL, (size_t)(n_row), sizeof (double), GMT_program);
	
	for (i = 0; i < n_row; i++) z[(i+shift)%n_row] = stack[prev][col][i];
	memcpy ((void *)stack[prev][col], (void *)z, (size_t)(header.n_row * sizeof (double)));
	GMT_free ((void *)z);
}

void table_SIGN (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, operand == 0!\n", GMT_program);
	if (constant[last]) a = copysign (1.0, factor[last]);
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : copysign (1.0, stack[last][col][i]);
}

void table_SIN (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (constant[last]) a = sin (factor[last]);
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : sin (stack[last][col][i]);
}

void table_SINC (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (constant[last]) a = GMT_sinc (factor[last]);
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : GMT_sinc (stack[last][col][i]);
}

void table_SIND (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (constant[last]) a = sind (factor[last]);
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : sind (stack[last][col][i]);
}

void table_SINH (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (constant[last]) a = sinh (factor[last]);
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : sinh (stack[last][col][i]);
}

void table_SQRT (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, operand one == 0!\n", GMT_program);
	if (constant[last]) a = sqrt (factor[last]);
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : sqrt (stack[last][col][i]);
}

void table_STD (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i, n_a = 0;
	double sum_a = 0.0, sum_a2 = 0.0, std;
	
	if (constant[last]) {	/* Trivial case */
		memset ((void *)stack[last][col], 0, (size_t)(n_row * sizeof (double)));
		return;
	}

	for (i = 0; i < n_row; i++) {
		if (GMT_is_fnan (stack[last][col][i])) continue;
		sum_a += stack[last][col][i];
		sum_a2 += (stack[last][col][i] * stack[last][col][i]);
		n_a++;
	}
	if (n_a > 1)
		std = sqrt ((n_a * sum_a2 - sum_a * sum_a) / (n_a * (n_a - 1.0)));
	else
		std = 0.0;
	for (i = 0; i < n_row; i++) stack[last][col][i] = std;
}

void table_STEP (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	for (i = 0; i < n_row; i++) {
		a = (constant[last]) ? factor[last] : stack[last][col][i];
		if (a == 0.0)
			stack[last][col][i] = 0.5;
		else
			stack[last][col][i] = (a < 0.0) ? 0.0 : 1.0;
	}
}

void table_STEPT (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	for (i = 0; i < n_row; i++) {
		a = t_coordinates[i] - ((constant[last]) ? factor[last] : stack[last][col][i]);
		if (a == 0.0)
			stack[last][col][i] = 0.5;
		else
			stack[last][col][i] = (a < 0.0) ? 0.0 : 1.0;
	}
}

void table_SUB (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	if (gmtdefs.verbose && constant[prev] && factor[prev] == 0.0) fprintf (stderr, "%s: Warning, operand one == 0!\n", GMT_program);
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, operand two == 0!\n", GMT_program);
	for (i = 0; i < n_row; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][col][i];
		b = (constant[last]) ? factor[last] : stack[last][col][i];
		stack[prev][col][i] = a - b;
	}
}

void table_SUM (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a, sum = 0.0;
	
	if (constant[last]) a = factor[last];
	for (i = 0; i < n_row; i++) {
		if (!constant[last]) a = stack[last][col][i];
		if (!GMT_is_dnan (a)) sum += a;
		stack[last][col][i] = sum;
	}
}

void table_TAN (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (constant[last]) a = tan (factor[last]);
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : tan (stack[last][col][i]);
}

void table_TAND (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (constant[last]) a = tand (factor[last]);
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : tand (stack[last][col][i]);
}

void table_TANH (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (constant[last]) a = tanh (factor[last]);
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : tanh (stack[last][col][i]);
}

void table_TCRIT (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	if (gmtdefs.verbose && constant[prev] && factor[prev] == 0.0) fprintf (stderr, "%s: Warning, operand one == 0 for TCRIT!\n", GMT_program);
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, operand two == 0 for TCRIT!\n", GMT_program);
	for (i = 0; i < n_row; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][col][i];
		b = (constant[last]) ? factor[last] : stack[last][col][i];
		stack[prev][col][i] = GMT_tcrit (a, b);
	}
}

void table_TDIST (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i, b, prev;
	double a;
	
	prev = last - 1;
	if (gmtdefs.verbose && constant[prev] && factor[prev] == 0.0) fprintf (stderr, "%s: Warning, operand one == 0 for TDIST!\n", GMT_program);
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, operand two == 0 for TDIST!\n", GMT_program);
	for (i = 0; i < n_row; i++) {
		a = (constant[prev]) ? factor[prev] : stack[prev][col][i];
		b = irint ((constant[last]) ? factor[last] : stack[last][col][i]);
		(void) GMT_student_t_a (a, b, &stack[prev][col][i]);
	}
}

void table_UPPER (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double high;

	if (constant[last]) {	/* Trivial case */
		for (i = 0; i < n_row; i++) stack[last][col][i] = factor[last];
		return;
	}
	
	for (i = 0, high = -DBL_MAX; i < n_row; i++) {
		if (GMT_is_fnan (stack[last][col][i])) continue;
		if (stack[last][col][i] > high) high = stack[last][col][i];
	}
	for (i = 0; i < n_row; i++) if (!GMT_is_fnan (stack[last][col][i])) stack[last][col][i] = high;
}

void table_XOR (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i, prev;
	double a, b;
	
	prev = last - 1;
	if (constant[prev]) a = factor[prev];
	if (constant[last]) b = factor[last];
	for (i = 0; i < n_row; i++) {
		if (!constant[prev]) a = stack[prev][col][i];
		if (!constant[last]) b = stack[last][col][i];
		stack[prev][col][i] = (GMT_is_dnan (a)) ? b : a;
	}
}

void table_Y0 (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, operand = 0 for Y0!\n", GMT_program);
	if (constant[last]) a = y0 (fabs (factor[last]));
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : y0 (fabs (stack[last][col][i]));
}

void table_Y1 (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, operand = 0 for Y1!\n", GMT_program);
	if (constant[last]) a = y1 (fabs (factor[last]));
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : y1 (fabs (stack[last][col][i]));
}

void table_YN (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i, prev, order;
	BOOLEAN simple = FALSE;
	double b;
	
	prev = last - 1;
	if (gmtdefs.verbose && constant[last] && factor[last] < 0.0) fprintf (stderr, "%s: Warning, order < 0 for YN!\n", GMT_program);
	if (gmtdefs.verbose && constant[last] && fabs (rint(factor[last]) - factor[last]) > SMALL) fprintf (stderr, "%s: Warning, order not an integer for YN!\n", GMT_program);
	if (gmtdefs.verbose && constant[prev] && factor[prev] == 0.0) fprintf (stderr, "%s: Warning, argument = 0 for YN!\n", GMT_program);
	if (constant[last]) order = irint (fabs (factor[last]));
	if (constant[last] && constant[prev]) {
		b = yn (order, fabs (factor[prev]));
		simple = TRUE;
	}
	for (i = 0; i < n_row; i++) {
		if (simple)
			stack[prev][col][i] = b;
		else {
			if (!constant[last]) order = irint (fabs (stack[last][col][i]));
			stack[prev][col][i] = yn (order, fabs (stack[prev][col][i]));
		}
	}
}

void table_ZCRIT (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (constant[last]) a = GMT_zcrit (factor[last]);
	for (i = 0; i < n_row; i++) stack[last][col][i] = (constant[last]) ? a : GMT_zcrit (stack[last][col][i]);
}

void table_ROOTS (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i, prev;
	double *roots;
	
	/* Treats the chosen column (at there is only one) as f(t) and solves for t that makes f(t) == 0.
	 * For now we only solve using a linear spline but in the future this should depend on the users
	 * choice of INTERPOLANT. */

	if (roots_found) return;	/* Already been here */
	if (!constant[last]) {
		fprintf (stderr, "%s: Argument to operator ROOTS must be a constant: the column number. Reset to 0\n", GMT_program);
		r_col = 0;
	}
	else
		r_col = irint (factor[last]);
	if (r_col < 0 || r_col >= header.n_col) {
		fprintf (stderr, "%s: Argument to operator ROOTS must be a column number 0 < col < %d. Reset to 0\n", GMT_program, header.n_col);
		r_col = 0;
	}
	roots = (double *) GMT_memory (VNULL, (size_t)(n_row), sizeof (double), GMT_program);
	n_roots = 0;
	prev = last - 1;
	if (stack[prev][r_col][0] == 0.0) roots[n_roots++] = t_coordinates[0]; 
	for (i = 1; i < n_row; i++) {
		if (stack[prev][r_col][i] == 0.0) {
			roots[n_roots++] = t_coordinates[i];
			continue;
		}
		
		if ((stack[prev][r_col][i-1] * stack[prev][r_col][i]) < 0.0) {	/* Crossing 0 */
			roots[n_roots] = t_coordinates[i-1] - stack[prev][r_col][i-1] * (t_coordinates[i] - t_coordinates[i-1]) / (stack[prev][r_col][i] - stack[prev][r_col][i-1]);
			n_roots++;
		}
	}
	for (i = 0; i < n_roots; i++) stack[prev][r_col][i] = roots[i];
	GMT_free ((void *)roots);
	roots_found = TRUE;
}

/* ---------------------- end operator functions --------------------- */

void solve_LSQFIT (double **stack[], int last, int n_col, int n_row, int skip[], char *file)
{
	/* Consider the current table the augmented matrix [A | b], making up the linear system Ax = b.
	 * We will set up the normal equations, solve for x, and output the solution before quitting.
	 * This function is special since it operates across columns and returns n_col scalars.
	 * We try to solve this positive definite & symmetric matrix with Cholsky methods; if that fails
	 * we do a full SVD decomposition and set small eigenvalues to zero, yielding an approximate solution.
	 */
	 
	int i, j, k, k0, i2, j2, rhs, n, ier;
	double cond, *N, *B, *d, *x, *b, *z, *v, *lambda;
	FILE *fp;
	
	for (i = n = 0; i < n_col; i++) if (!skip[i]) n++;	/* Need to find how many active columns we have */
	if (n < 2) {
		fprintf (stderr, "%s: Error, LSQSOL requires at least 2 active columns!\n", GMT_program);
		exit (EXIT_FAILURE);
	}
	rhs = n_col - 1;
	while (skip[rhs] && rhs > 0) rhs--; 	/* Get last active col number as the rhs vector b */
	n--;					/* Account for b, the rhs vector, to get row & col dimensions of normal matrix N */
	
	N = (double *) GMT_memory (VNULL, (size_t)(n*n), sizeof (double), GMT_program);
	B = (double *) GMT_memory (VNULL, (size_t)n_row, sizeof (double), GMT_program);
	
	/* Do the row & col dot products, skipping inactive columns as we go along */
	for (j = j2 = 0; j < n; j2++) {	/* j2 is table column, j is row in N matrix */
		if (skip[j2]) continue;
		for (i = i2 = 0; i < n; i2++) {	/* i2 is table column, i is column in N matrix */
			if (skip[i2]) continue;
			k0 = j * n + i;
			N[k0] = 0.0;
			for (k = 0; k < n_row; k++) N[k0] += stack[last][j][k] * stack[last][i][k];
			i++;
		}
		B[j] = 0.0;
		for (k = 0; k < n_row; k++) B[j] += stack[last][j][k] * stack[last][rhs][k];
		j++;
	}
	
	d = (double *) GMT_memory (VNULL, (size_t)n, sizeof(double), GMT_program);
	x = (double *) GMT_memory (VNULL, (size_t)n, sizeof(double), GMT_program);
	if ( (ier = GMT_chol_dcmp (N, d, &cond, n, n) ) != 0) {	/* Decomposition failed, use SVD method */
		int nrots;
		GMT_chol_recover (N, d, n, n, ier, TRUE);		/* Restore to former matrix N */
		/* Solve instead using GMT_jacobi */
		lambda = (double *) GMT_memory (VNULL, (size_t)n, sizeof(double), GMT_program);
		b = (double *) GMT_memory (VNULL, (size_t)n, sizeof(double), GMT_program);
		z = (double *) GMT_memory (VNULL, (size_t)n, sizeof(double), GMT_program);
		v = (double *) GMT_memory (VNULL, (size_t)n*n, sizeof(double), GMT_program);

		if (GMT_jacobi (N, &n, &n, lambda, v, b, z, &nrots)) {
			fprintf (stderr, "%s: Eigenvalue routine failed to converge in 50 sweeps.\n", GMT_program);
			fprintf (stderr, "%s: The reported L2 positions might be garbage.\n", GMT_program);
		}
		/* Solution x = v * lambda^-1 * v' * B */
		
		/* First do d = V' * B, so x = v * lambda^-1 * d */
		for (j = 0; j < n; j++) for (k = 0, d[j] = 0.0; k < n; k++) d[j] += v[k*n+j] * B[k];
		/* Then do d = lambda^-1 * d by setting small lambda's to zero */
		for (j = k = 0; j < n; j++) {
			if (lambda[j] < 1.0e7) {
				d[j] = 0.0;
				k++;
			}
			else
				d[j] /= lambda[j];
		}
		if (k) fprintf (stderr,"%s: %d eigenvalues < 1.0e-7 set to zero to yield a stable solution\n", GMT_program, k);

		/* Finally do x = v * d */
		for (j = 0; j < n; j++) for (k = 0; k < n; k++) x[j] += v[j*n+k] * d[k];

		GMT_free ((void *)b);
		GMT_free ((void *)z);
		GMT_free ((void *)v);
	}
	else {	/* Decomposition worked, now solve system */
		GMT_chol_solv (N, x, B, n, n);
	}
	
	if (file && (fp = GMT_fopen (file, GMT_io.w_mode)) == NULL) {
		fprintf (stderr, "%s: Error creating file %s\n", GMT_program, file);
		exit (EXIT_FAILURE);
	}
	else if (!file)
		fp = GMT_stdout;

#ifdef SET_IO_MODE
	GMT_setmode (1);
#endif
	GMT_output (fp, n, x);
	if (fp != GMT_stdout) GMT_fclose (fp);
	
	GMT_free ((void *)x);
	GMT_free ((void *)d);
	GMT_free ((void *)N);
	GMT_free ((void *)B);
}

/* ---------------------- start convenience functions --------------------- */

int decode_argument (txt)
char *txt; {
	int number, minus, plus, exp, dec, n_digits = 0;
	char *s;
	
	if (!strcmp (txt, "/")) return ARG_IS_OPERATOR;	/* Special test since / is also a directory */

	/* First see if argument it can be opened as a file */
	
	if (!access (txt, R_OK)) return ARG_IS_FILE; 	/* returns FALSE if file exists and can be read */

	/* Next look for symbols with special meaning */

	if (!(strcmp (txt, "PI") && strcmp (txt, "pi"))) return ARG_IS_PI;
	if (!(strcmp (txt, "E") && strcmp (txt, "e"))) return ARG_IS_E;
	if (!(strcmp (txt, "T") && strcmp (txt, "t"))) return ARG_IS_T_MATRIX;

	/* Here we must check if argument is a numerical value */

	s = txt;
	if (*s == '-' || *s == '+') s++;	/* Skip leading sign */
	
	minus = plus = exp = dec = 0;
	number = TRUE;

	while (number && *s) {
		if (isdigit ((int)*s))
			n_digits++;
		else {
			switch (*s) {
				case '-':
					minus++;
					break;
				case '+':
					plus++;
					break;
				case 'E':
				case 'e':
					exp++;
					break;
				case '.':
					dec++;
					break;
				default:
					number = FALSE;
					break;
			}
		}
		if (minus > 1 || exp > 1 || dec > 1) number = FALSE;
		s++;
	}
	if (txt[0] == '-' && !number) {	/* Probably a bad commandline option */
		fprintf (stderr, "%s: ERROR: Option %s not recognized\n", GMT_program, txt);
		exit (EXIT_FAILURE);
	}
	
	return ((number && n_digits > 0) ? ARG_IS_NUMBER : ARG_IS_OPERATOR);
	
}

int get_operator (choice)
char *choice; {
	int op;
	
	/* Returns -1 if not a registered operator */

	op = GMT_hash_lookup (choice, localhashnode, HASH_SIZE, HASH_SIZE);

	if (op < 0 && strlen (choice) == 1) {	/* Check for old-style operators */
	
		switch (choice[0]) {
			case '+':
				op = ADD;
				break;
			case '-':
				op = SUB;
				break;
			case 'x':
				op = MUL;
				break;
			case '/':
				op = DIV;
				break;
			case '^':
				op = RAISE;
				break;
		}
	}

	return (op);
}

void gmtmath_init (PFV ops[], int n_args[], int n_out[])
{

	/* Operator function		# of operands  		# of outputs */

	ops[0] = table_ABS;		n_args[0] = 1;		n_out[0] = 1;
	ops[1] = table_ACOS;		n_args[1] = 1;		n_out[1] = 1;
	ops[2] = table_ACOSH;		n_args[2] = 1;		n_out[2] = 1;
	ops[3] = table_ADD;		n_args[3] = 2;		n_out[3] = 1;
	ops[4] = table_AND;		n_args[4] = 2;		n_out[4] = 1;
	ops[5] = table_ASIN;		n_args[5] = 1;		n_out[5] = 1;
	ops[6] = table_ASINH;		n_args[6] = 1;		n_out[6] = 1;
	ops[7] = table_ATAN;		n_args[7] = 1;		n_out[7] = 1;
	ops[8] = table_ATAN2;		n_args[8] = 2;		n_out[8] = 1;
	ops[9] = table_ATANH;		n_args[9] = 1;		n_out[9] = 1;
	ops[10] = table_BEI;		n_args[10] = 1;		n_out[10] = 1;
	ops[11] = table_BER;		n_args[11] = 1;		n_out[11] = 1;
	ops[12] = table_CEIL;		n_args[12] = 1;		n_out[12] = 1;
	ops[13] = table_CHICRIT;		n_args[13] = 2;		n_out[13] = 1;
	ops[14] = table_CHIDIST;		n_args[14] = 2;		n_out[14] = 1;
	ops[15] = table_COL;		n_args[15] = 1;		n_out[15] = 1;
	ops[16] = table_COS;		n_args[16] = 1;		n_out[16] = 1;
	ops[17] = table_COSD;		n_args[17] = 1;		n_out[17] = 1;
	ops[18] = table_COSH;		n_args[18] = 1;		n_out[18] = 1;
	ops[19] = table_D2DT2;		n_args[19] = 1;		n_out[19] = 1;
	ops[20] = table_D2R;		n_args[20] = 1;		n_out[20] = 1;
	ops[21] = table_DILOG;		n_args[21] = 1;		n_out[21] = 1;
	ops[22] = table_DIV;		n_args[22] = 2;		n_out[22] = 1;
	ops[23] = table_DDT;		n_args[23] = 1;		n_out[23] = 1;
	ops[24] = table_DUP;		n_args[24] = 1;		n_out[24] = 2;
	ops[25] = table_ERF;		n_args[25] = 1;		n_out[25] = 1;
	ops[26] = table_ERFC;		n_args[26] = 1;		n_out[26] = 1;
	ops[27] = table_ERFINV;		n_args[27] = 1;		n_out[27] = 1;
	ops[28] = table_EQ;		n_args[28] = 2;		n_out[28] = 1;
	ops[29] = table_EXCH;		n_args[29] = 2;		n_out[29] = 2;
	ops[30] = table_EXP;		n_args[30] = 1;		n_out[30] = 1;
	ops[31] = table_FCRIT;		n_args[31] = 3;		n_out[31] = 1;
	ops[32] = table_FDIST;		n_args[32] = 3;		n_out[32] = 1;
	ops[33] = table_FLIPUD;		n_args[33] = 1;		n_out[33] = 1;
	ops[34] = table_FLOOR;		n_args[34] = 1;		n_out[34] = 1;
	ops[35] = table_FMOD;		n_args[35] = 2;		n_out[35] = 1;
	ops[36] = table_GE;		n_args[36] = 2;		n_out[36] = 1;
	ops[37] = table_GT;		n_args[37] = 2;		n_out[37] = 1;
	ops[38] = table_HYPOT;		n_args[38] = 2;		n_out[38] = 1;
	ops[39] = table_I0;		n_args[39] = 1;		n_out[39] = 1;
	ops[40] = table_I1;		n_args[40] = 1;		n_out[40] = 1;
	ops[41] = table_IN;		n_args[41] = 2;		n_out[41] = 1;
	ops[42] = table_INT;		n_args[42] = 1;		n_out[42] = 1;
	ops[43] = table_INV;		n_args[43] = 1;		n_out[43] = 1;
	ops[44] = table_ISNAN;		n_args[44] = 1;		n_out[44] = 1;
	ops[45] = table_J0;		n_args[45] = 1;		n_out[45] = 1;
	ops[46] = table_J1;		n_args[46] = 1;		n_out[46] = 1;
	ops[47] = table_JN;		n_args[47] = 2;		n_out[47] = 1;
	ops[48] = table_K0;		n_args[48] = 1;		n_out[48] = 1;
	ops[49] = table_K1;		n_args[49] = 1;		n_out[49] = 1;
	ops[50] = table_KN;		n_args[50] = 2;		n_out[50] = 1;
	ops[51] = table_KEI;		n_args[51] = 1;		n_out[51] = 1;
	ops[52] = table_KER;		n_args[52] = 1;		n_out[52] = 1;
	ops[53] = table_LE;		n_args[53] = 2;		n_out[53] = 1;
	ops[54] = table_LMSSCL;		n_args[54] = 1;		n_out[54] = 1;
	ops[55] = table_LOG;		n_args[55] = 1;		n_out[55] = 1;
	ops[56] = table_LOG10;		n_args[56] = 1;		n_out[56] = 1;
	ops[57] = table_LOG1P;		n_args[57] = 1;		n_out[57] = 1;
	ops[58] = table_LOG2;		n_args[58] = 1;		n_out[58] = 1;
	ops[59] = table_LOWER;		n_args[59] = 1;		n_out[59] = 1;
	ops[60] = table_LRAND;		n_args[60] = 2;		n_out[60] = 1;
	ops[61] = table_LSQFIT;		n_args[61] = 1;		n_out[61] = 0;
	ops[62] = table_LT;		n_args[62] = 2;		n_out[62] = 1;
	ops[63] = table_MAD;		n_args[63] = 1;		n_out[63] = 1;
	ops[64] = table_MAX;		n_args[64] = 2;		n_out[64] = 1;
	ops[65] = table_MEAN;		n_args[65] = 1;		n_out[65] = 1;
	ops[66] = table_MED;		n_args[66] = 1;		n_out[66] = 1;
	ops[67] = table_MIN;		n_args[67] = 2;		n_out[67] = 1;
	ops[68] = table_MODE;		n_args[68] = 1;		n_out[68] = 1;
	ops[69] = table_MUL;		n_args[69] = 2;		n_out[69] = 1;
	ops[70] = table_NAN;		n_args[70] = 2;		n_out[70] = 1;
	ops[71] = table_NEG;		n_args[71] = 1;		n_out[71] = 1;
	ops[72] = table_NEQ;		n_args[72] = 2;		n_out[72] = 1;
	ops[73] = table_NRAND;		n_args[73] = 2;		n_out[73] = 1;
	ops[74] = table_OR;		n_args[74] = 2;		n_out[74] = 1;
	ops[75] = table_PLM;		n_args[75] = 3;		n_out[75] = 1;
	ops[76] = table_POP;		n_args[76] = 1;		n_out[76] = 0;
	ops[77] = table_POW;		n_args[77] = 2;		n_out[77] = 1;
	ops[78] = table_R2;		n_args[78] = 2;		n_out[78] = 1;
	ops[79] = table_R2D;		n_args[79] = 1;		n_out[79] = 1;
	ops[80] = table_RAND;		n_args[80] = 2;		n_out[80] = 1;
	ops[81] = table_RINT;		n_args[81] = 1;		n_out[81] = 1;
	ops[82] = table_ROOTS;		n_args[82] = 2;		n_out[82] = 1;
	ops[83] = table_ROTT;		n_args[83] = 2;		n_out[83] = 1;
	ops[84] = table_SIGN;		n_args[84] = 1;		n_out[84] = 1;
	ops[85] = table_SIN;		n_args[85] = 1;		n_out[85] = 1;
	ops[86] = table_SINC;		n_args[86] = 1;		n_out[86] = 1;
	ops[87] = table_SIND;		n_args[87] = 1;		n_out[87] = 1;
	ops[88] = table_SINH;		n_args[88] = 1;		n_out[88] = 1;
	ops[89] = table_SQRT;		n_args[89] = 1;		n_out[89] = 1;
	ops[90] = table_STD;		n_args[90] = 1;		n_out[90] = 1;
	ops[91] = table_STEP;		n_args[91] = 1;		n_out[91] = 1;
	ops[92] = table_STEPT;		n_args[92] = 1;		n_out[92] = 1;
	ops[93] = table_SUB;		n_args[93] = 2;		n_out[93] = 1;
	ops[94] = table_SUM;		n_args[94] = 1;		n_out[94] = 1;
	ops[95] = table_TAN;		n_args[95] = 1;		n_out[95] = 1;
	ops[96] = table_TAND;		n_args[96] = 1;		n_out[96] = 1;
	ops[97] = table_TANH;		n_args[97] = 1;		n_out[97] = 1;
	ops[98] = table_TCRIT;		n_args[98] = 2;		n_out[98] = 1;
	ops[99] = table_TDIST;		n_args[99] = 2;		n_out[99] = 1;
	ops[100] = table_UPPER;		n_args[100] = 1;		n_out[100] = 1;
	ops[101] = table_XOR;		n_args[101] = 2;		n_out[101] = 1;
	ops[102] = table_Y0;		n_args[102] = 1;		n_out[102] = 1;
	ops[103] = table_Y1;		n_args[103] = 1;		n_out[103] = 1;
	ops[104] = table_YN;		n_args[104] = 2;		n_out[104] = 1;
	ops[105] = table_ZCRIT;		n_args[105] = 1;		n_out[105] = 1;
}
