/*--------------------------------------------------------------------
 *	$Id: gmtmath_main.c,v 1.4.4.3 2002/02/27 17:58:55 pwessel Exp $
 *
 *	Copyright (c) 1991-2002 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
 *		3.4.1
 *
 */
 
#include "gmt.h"

#define N_OPERATORS	93
#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 hashnode[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);
void new_table (double ***s, int n_col, int n);
void decode_columns (char *txt, int *skip, int n_col, int t_col);

BOOLEAN very_first = TRUE;
char head_record[BUFSIZ];

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

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;
	
	BOOLEAN constant[STACK_SIZE], error = FALSE, set_t = FALSE, set_n = FALSE;
	BOOLEAN one_row_out = FALSE, set_q = FALSE, no_time = FALSE;
	
	double **stack[STACK_SIZE];

	double factor[STACK_SIZE], t_noise;
	
	char *outfile = CNULL, file[BUFSIZ];
	
	struct TABLE_HEADER tbl[STACK_SIZE];
	
	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 [-C<cols>] [-H[<nrec>]] [-N<n_col>/<t_col>] [-Q] [-S] [-T[<tmin/tmax/t_inc>]]\n");
		fprintf (stderr, "\t[-bi[s][<n>]] [-bo[s][<n>]] [-V] 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-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-T Set domain from t_min to t_max in steps of t_inc\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');
		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 .gmtcommands */

	/* 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 (hashnode, 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, -H switches before reading any file */

	for (arg = 1; arg < last_arg; arg++) {
		if (!strncmp (argv[arg], "-b", 2)) error += GMT_io_selection (&argv[arg][2]);
		if (!strncmp (argv[arg], "-H", 2)) error += GMT_get_common_args (argv[arg], NULL, NULL, NULL, NULL);
	}
	/* 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, (double ***)NULL, t_col);
		
		nm = header.n_row * header.n_col;
		n_col = header.n_col;
	}
	
	/* Scan command line for -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 'N':
					sscanf (&argv[arg][2], "%d/%d", &n_col, &t_col);
					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':
					one_row_out = TRUE;
					break;
				case 'T':
					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[0] && gmtdefs.io_header) {
		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:
				fprintf (stderr, "%s: GMT SYNTAX ERROR -T:  inc is <= 0\n", GMT_program);
				exit (EXIT_FAILURE);
				break;
			case 3:
				fprintf (stderr, "%s: GMT SYNTAX ERROR -T:  (max - min) 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_q) header.n_row = header.n_col = nm = 1;
	
	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 x and y vectors */

	t_coordinates = (double *) GMT_memory (VNULL, (size_t)header.n_row, sizeof (double), GMT_program);
	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)));

	gmtmath_init (call_operator, consumed_operands, produced_operands);

	nstack = 0;
	
	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], "-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], "-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!", 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 ((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);
		}

		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) fprintf (stderr, "= %s", outfile);

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

	if (one_row_out) header.n_row = 1;	/* Trick to get only the first row on output */
	GMT_write_table (outfile, &header, stack[0]);
	
	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[0]) ? GMT_io.ncol[0] : 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 && 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 i, j;
	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_io.binary[1]) {
		fprintf (fp, "%s", head_record);
		for (i = 1; i < gmtdefs.n_header_recs; i++) fprintf (fp, "# gmtmath header record\n");
	}
		
	for (i = 0; i < h->n_row; 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);
	for (i = 0; i < n_row; i++) {
		a = (constant[last]) ? factor[last] : stack[last][col][i];
		stack[last][col][i] = fabs (a);
	}
}

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);
	for (i = 0; i < n_row; i++) {
		a = (constant[last]) ? factor[last] : stack[last][col][i];
		stack[last][col][i] = d_acos (a);
	}
}

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);
	for (i = 0; i < n_row; i++) {
		a = (constant[last]) ? factor[last] : stack[last][col][i];
		stack[last][col][i] = acosh (a);
	}
}

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);
	for (i = 0; i < n_row; i++) {
		a = (constant[last]) ? factor[last] : stack[last][col][i];
		stack[last][col][i] = d_asin (a);
	}
}

void table_ASINH (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];
		stack[last][col][i] = asinh (a);
	}
}

void table_ATAN (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];
		stack[last][col][i] = atan (a);
	}
}

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);
	for (i = 0; i < n_row; i++) {
		a = (constant[last]) ? factor[last] : stack[last][col][i];
		stack[last][col][i] = atanh (a);
	}
}

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;
	
	for (i = 0; i < n_row; i++) {
		a = (constant[last]) ? factor[last] : stack[last][col][i];
		stack[last][col][i] = ceil (a);
	}
}

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_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;
	
	for (i = 0; i < n_row; i++) {
		a = (constant[last]) ? factor[last] : stack[last][col][i];
		stack[last][col][i] = cosd (a);
	}
}

void table_COSH (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];
		stack[last][col][i] = cosh (a);
	}
}

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 (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 (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;
	
	for (i = 0; i < n_row; i++) {
		a = (constant[last]) ? factor[last] : stack[last][col][i];
		stack[last][col][i] = a * 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 */
		factor[last] = 1.0 / factor[last];
		table_MUL (stack, constant, factor, last, col, n_row);
		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);
	for (i = 0; i < n_row; i++) {
		a = (constant[last]) ? factor[last] : stack[last][col][i];
		stack[last][col][i] = exp (a);
	}
}

void table_FDIST (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i, c, d, prev1, prev2, prev3;
	double a, b;
	
	prev1 = last - 1;
	prev2 = last - 2;
	prev3 = last - 3;
	if (gmtdefs.verbose && constant[prev3] && factor[prev3] == 0.0) fprintf (stderr, "%s: Warning, operand one == 0 for FDIST!\n", GMT_program);
	if (gmtdefs.verbose && constant[prev2] && factor[prev2] == 0.0) fprintf (stderr, "%s: Warning, operand two == 0 for FDIST!\n", GMT_program);
	if (gmtdefs.verbose && constant[prev1] && factor[prev1] == 0.0) fprintf (stderr, "%s: Warning, operand three == 0 for FDIST!\n", GMT_program);
	if (gmtdefs.verbose && constant[last] && factor[last] == 0.0) fprintf (stderr, "%s: Warning, operand four == 0 for FDIST!\n", GMT_program);
	for (i = 0; i < n_row; i++) {
		a = (constant[prev3]) ? factor[prev3] : stack[prev3][col][i];
		b = (constant[prev2]) ? factor[prev2] : stack[prev2][col][i];
		c = irint ((double)((constant[prev1]) ? factor[prev1] : stack[prev1][col][i]));
		d = irint ((double)((constant[last]) ? factor[last] : stack[last][col][i]));
		(void) GMT_f_q (a * c, c, b * d, d, &stack[prev3][col][i]);
	}
}

void table_FLOOR (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];
		stack[last][col][i] = floor (a);
	}
}

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;
	}
	
	left = stack[last][col][0];
	stack[last][col][0] = sum = 0.0;
	f = 0.5 * header.t_inc;
	for (i = 1; i < n_row; i++) {	/* Dumb trapezoidal rule */
		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]) factor[last] = 1.0 / factor[last];
	for (i = 0; i < n_row; i++) {
		a = (constant[last]) ? factor[last] : 1.0 / stack[last][col][i];
		stack[last][col][i] = a;
	}
}

void table_ISNAN (double **stack[], BOOLEAN *constant, double *factor, int last, int col, int n_row)
{
	int i;
	double a;
	
	if (constant[last]) a = factor[last];
	for (i = 0; i < n_row; i++) {
		if (!constant[last]) a = stack[last][col][i];
		stack[last][col][i] = (double)GMT_is_dnan (a);
	}
}
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, &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;
}

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_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_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, &mode);
	else
		mode = GMT_d_NaN;

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

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);
	for (i = 0; i < n_row; i++) {
		a = (constant[last]) ? factor[last] : stack[last][col][i];
		stack[last][col][i] = -a;
	}
}

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 rutine 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;
	
	for (i = 0; i < n_row; i++) {
		a = (constant[last]) ? factor[last] : stack[last][col][i];
		stack[last][col][i] = a * 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;
	
	for (i = 0; i < n_row; i++) {
		a = (constant[last]) ? factor[last] : stack[last][col][i];
		stack[last][col][i] = rint (a);
	}
}

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);
	for (i = 0; i < n_row; i++) {
		a = (constant[last]) ? factor[last] : stack[last][col][i];
		stack[last][col][i] = copysign (1.0, a);
	}
}

void table_SIN (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];
		stack[last][col][i] = sin (a);
	}
}

void table_SIND (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];
		stack[last][col][i] = sind (a);
	}
}

void table_SINH (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];
		stack[last][col][i] = sinh (a);
	}
}

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);
	for (i = 0; i < n_row; i++) {
		a = (constant[last]) ? factor[last] : stack[last][col][i];
		stack[last][col][i] = sqrt (a);
	}
}

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;
	
	for (i = 0; i < n_row; i++) {
		a = (constant[last]) ? factor[last] : stack[last][col][i];
		stack[last][col][i] = tan (a);
	}
}

void table_TAND (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];
		stack[last][col][i] = tand (a);
	}
}

void table_TANH (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];
		stack[last][col][i] = tanh (a);
	}
}

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]));
		}
	}
}

/* ---------------------- end operator 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, hashnode, 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_CHIDIST;		n_args[13] = 2;		n_out[13] = 1;
	ops[14] = table_COS;		n_args[14] = 1;		n_out[14] = 1;
	ops[15] = table_COSD;		n_args[15] = 1;		n_out[15] = 1;
	ops[16] = table_COSH;		n_args[16] = 1;		n_out[16] = 1;
	ops[17] = table_D2DT2;		n_args[17] = 1;		n_out[17] = 1;
	ops[18] = table_D2R;		n_args[18] = 1;		n_out[18] = 1;
	ops[19] = table_DILOG;		n_args[19] = 1;		n_out[19] = 1;
	ops[20] = table_DIV;		n_args[20] = 2;		n_out[20] = 1;
	ops[21] = table_DDT;		n_args[21] = 1;		n_out[21] = 1;
	ops[22] = table_DUP;		n_args[22] = 1;		n_out[22] = 2;
	ops[23] = table_ERF;		n_args[23] = 1;		n_out[23] = 1;
	ops[24] = table_ERFC;		n_args[24] = 1;		n_out[24] = 1;
	ops[25] = table_ERFINV;		n_args[25] = 1;		n_out[25] = 1;
	ops[26] = table_EQ;		n_args[26] = 2;		n_out[26] = 1;
	ops[27] = table_EXCH;		n_args[27] = 2;		n_out[27] = 2;
	ops[28] = table_EXP;		n_args[28] = 1;		n_out[28] = 1;
	ops[29] = table_FDIST;		n_args[29] = 4;		n_out[29] = 1;
	ops[30] = table_FLOOR;		n_args[30] = 1;		n_out[30] = 1;
	ops[31] = table_FMOD;		n_args[31] = 2;		n_out[31] = 1;
	ops[32] = table_GE;		n_args[32] = 2;		n_out[32] = 1;
	ops[33] = table_GT;		n_args[33] = 2;		n_out[33] = 1;
	ops[34] = table_HYPOT;		n_args[34] = 2;		n_out[34] = 1;
	ops[35] = table_I0;		n_args[35] = 1;		n_out[35] = 1;
	ops[36] = table_I1;		n_args[36] = 1;		n_out[36] = 1;
	ops[37] = table_IN;		n_args[37] = 2;		n_out[37] = 1;
	ops[38] = table_INT;		n_args[38] = 1;		n_out[38] = 1;
	ops[39] = table_INV;		n_args[39] = 1;		n_out[39] = 1;
	ops[40] = table_ISNAN;		n_args[40] = 1;		n_out[40] = 1;
	ops[41] = table_J0;		n_args[41] = 1;		n_out[41] = 1;
	ops[42] = table_J1;		n_args[42] = 1;		n_out[42] = 1;
	ops[43] = table_JN;		n_args[43] = 2;		n_out[43] = 1;
	ops[44] = table_K0;		n_args[44] = 1;		n_out[44] = 1;
	ops[45] = table_K1;		n_args[45] = 1;		n_out[45] = 1;
	ops[46] = table_KN;		n_args[46] = 2;		n_out[46] = 1;
	ops[47] = table_KEI;		n_args[47] = 1;		n_out[47] = 1;
	ops[48] = table_KER;		n_args[48] = 1;		n_out[48] = 1;
	ops[49] = table_LE;		n_args[49] = 2;		n_out[49] = 1;
	ops[50] = table_LMSSCL;		n_args[50] = 1;		n_out[50] = 1;
	ops[51] = table_LOG;		n_args[51] = 1;		n_out[51] = 1;
	ops[52] = table_LOG10;		n_args[52] = 1;		n_out[52] = 1;
	ops[53] = table_LOG1P;		n_args[53] = 1;		n_out[53] = 1;
	ops[54] = table_LOWER;		n_args[54] = 1;		n_out[54] = 1;
	ops[55] = table_LT;		n_args[55] = 2;		n_out[55] = 1;
	ops[56] = table_MAD;		n_args[56] = 1;		n_out[56] = 1;
	ops[57] = table_MAX;		n_args[57] = 2;		n_out[57] = 1;
	ops[58] = table_MEAN;		n_args[58] = 1;		n_out[58] = 1;
	ops[59] = table_MED;		n_args[59] = 1;		n_out[59] = 1;
	ops[60] = table_MIN;		n_args[60] = 2;		n_out[60] = 1;
	ops[61] = table_MODE;		n_args[61] = 1;		n_out[61] = 1;
	ops[62] = table_MUL;		n_args[62] = 2;		n_out[62] = 1;
	ops[63] = table_NAN;		n_args[63] = 2;		n_out[63] = 1;
	ops[64] = table_NEG;		n_args[64] = 1;		n_out[64] = 1;
	ops[65] = table_NRAND;		n_args[65] = 2;		n_out[65] = 1;
	ops[66] = table_OR;		n_args[66] = 2;		n_out[66] = 1;
	ops[67] = table_PLM;		n_args[67] = 3;		n_out[67] = 1;
	ops[68] = table_POP;		n_args[68] = 1;		n_out[68] = 0;
	ops[69] = table_POW;		n_args[69] = 2;		n_out[69] = 1;
	ops[70] = table_R2;		n_args[70] = 2;		n_out[70] = 1;
	ops[71] = table_R2D;		n_args[71] = 1;		n_out[71] = 1;
	ops[72] = table_RAND;		n_args[72] = 2;		n_out[72] = 1;
	ops[73] = table_RINT;		n_args[73] = 1;		n_out[73] = 1;
	ops[74] = table_SIGN;		n_args[74] = 1;		n_out[74] = 1;
	ops[75] = table_SIN;		n_args[75] = 1;		n_out[75] = 1;
	ops[76] = table_SIND;		n_args[76] = 1;		n_out[76] = 1;
	ops[77] = table_SINH;		n_args[77] = 1;		n_out[77] = 1;
	ops[78] = table_SQRT;		n_args[78] = 1;		n_out[78] = 1;
	ops[79] = table_STD;		n_args[79] = 1;		n_out[79] = 1;
	ops[80] = table_STEP;		n_args[80] = 1;		n_out[80] = 1;
	ops[81] = table_STEPT;		n_args[81] = 1;		n_out[81] = 1;
	ops[82] = table_SUB;		n_args[82] = 2;		n_out[82] = 1;
	ops[83] = table_SUM;		n_args[83] = 1;		n_out[83] = 1;
	ops[84] = table_TAN;		n_args[84] = 1;		n_out[84] = 1;
	ops[85] = table_TAND;		n_args[85] = 1;		n_out[85] = 1;
	ops[86] = table_TANH;		n_args[86] = 1;		n_out[86] = 1;
	ops[87] = table_TDIST;		n_args[87] = 2;		n_out[87] = 1;
	ops[88] = table_UPPER;		n_args[88] = 1;		n_out[88] = 1;
	ops[89] = table_XOR;		n_args[89] = 2;		n_out[89] = 1;
	ops[90] = table_Y0;		n_args[90] = 1;		n_out[90] = 1;
	ops[91] = table_Y1;		n_args[91] = 1;		n_out[91] = 1;
	ops[92] = table_YN;		n_args[92] = 2;		n_out[92] = 1;
}