#include <stdio.h>
#include <math.h>
#include <errno.h>
#if defined (sgi) || defined (sun) || defined (ultrix)
#include <values.h>
#include <nan.h>
#endif
#ifdef AIX
#include <fp.h>
#endif
#ifndef NOSTDLIB_H
#include <stdlib.h>
#endif

#include "fudgit.h"
#include "head.h"

extern int errno;
#ifndef __HAVE_68881__
extern double cbrt(double), atanh(double), asinh(double), acosh(double);
#ifndef AIX   /* AIX define the following as macros whenever possible */
extern double log(double), log10(double), acos(double);
extern double asin(double), atan(double);
extern double exp(double), sqrt(double);
#endif
#endif
#if defined(NeXT) || defined(OSF)
extern double lgamma(double);
#define gamma(x) lgamma(x)
#else
extern double gamma(double);
#endif
#ifndef __HAVE_68881__
extern double pow(double, double);
extern double cosh(double), sinh(double), tanh(double);
#endif
extern void Ft_matherror (char *s1, char *s2, int lino);
double Ft_errcheck(double x, char *str);

double Ft_Rand(void)
{
#ifdef NODRAND48
/* For systems which do not have a better random number generator  */
#include <stdlib.h>
    double x = RAND_MAX;

    return((double)rand()/x);
#else
	extern double drand48(void);
    return(drand48());
#endif
}

double Ft_Srand(double x)
{
    long xx;
#ifdef NODRAND48
#define srand48(x) srand(x)
#else
extern void srand48(long);
#endif
    xx = (long) x;
    srand48(xx);
    return((double)x);
}

double Ft_Log(double x)
{
    errno = 0;
    return(Ft_errcheck(log(x), "ln"));
}

double Ft_Log10(double x)
{
    errno = 0;
    return(Ft_errcheck(log10(x), "log"));
}

double Ft_Lgamma(double x)
{
    errno = 0;
    return(Ft_errcheck(gamma(x), "lgamma"));
}

double Ft_Exp(double x)
{
    errno = 0;
    return(Ft_errcheck(exp(x), "exp"));
}

double Ft_Sqrt(double x)
{
    errno = 0;
    return(Ft_errcheck(sqrt(x), "sqrt"));
}

double Ft_Pow(double x, double y)
{
    errno = 0;
    return(Ft_errcheck(pow(x, y), "exponentiation"));
}

double Ft_Cbrt(double x)
{
    errno = 0;
    return(Ft_errcheck(cbrt(x), "cbrt"));
}

double Ft_integer(double x)
{
    return((double)(long)x);
}

double Ft_Cosh(double x)
{
    errno = 0;
    return(Ft_errcheck(cosh(x), "cosh"));
}

double Ft_Sinh(double x)
{
    errno = 0;
    return(Ft_errcheck(sinh(x), "sinh"));
}

double Ft_Tanh(double x)
{
    errno = 0;
    return(Ft_errcheck(tanh(x), "tanh"));
}

double Ft_Acosh(double x)
{
    errno = 0;
    return(Ft_errcheck(acosh(x), "acosh"));
}

double Ft_Acos(double x)
{
    errno = 0;
    return(Ft_errcheck(acos(x), "acos"));
}

double Ft_Asin(double x)
{
    errno = 0;
    return(Ft_errcheck(asin(x), "asin"));
}

double Ft_Asinh(double x)
{
    errno = 0;
    return(Ft_errcheck(asinh(x), "asinh"));
}

double Ft_Atanh(double x)
{
    errno = 0;
    return(Ft_errcheck(atanh(x), "atanh"));
}

double Ft_Coth(double x)
{
    errno = 0;
    return(Ft_errcheck(1.0/tanh(x), "coth"));
}

double Ft_Csch(double x)
{
    errno = 0;
    return(Ft_errcheck(1.0/sinh(x), "csch"));
}

double Ft_Sech(double x)
{
    errno = 0;
    return(Ft_errcheck(1.0/cosh(x), "sech"));
}

double Ft_Cot(double x)
{
    errno = 0;
    return(Ft_errcheck(1.0/tan(x), "cot"));
}

double Ft_Hypot(double x, double y)
{
    errno = 0;
    return(Ft_errcheck(hypot(x, y), "hypot"));
}

double Ft_Atan2(double x, double y)
{
    errno = 0;
    return(Ft_errcheck(atan2(x, y), "atan2"));
}

double Ft_Atan(double x)
{
    errno = 0;
    return(Ft_errcheck(atan(x), "atan"));
}

double Ft_Tan(double x)
{
    errno = 0;
    return(Ft_errcheck(tan(x), "tan"));
}

double Ft_Csc(double x)
{
    errno = 0;
    return(Ft_errcheck(1.0/sin(x), "csc"));
}

double Ft_Sec(double x)
{
    errno = 0;
    return(Ft_errcheck(1.0/cos(x), "sec"));
}

double Ft_Y0(double d)
{
#ifdef AMIGA
    Ft_matherror("%s: Function y0() unavailable", NULL, 0);
    errno = EDOM;
    return 0;
#else
    errno = 0;
    return(Ft_errcheck(y0(d), "besy0"));
#endif
}

double Ft_Y1(double d)
{
#ifdef AMIGA
    Ft_matherror("%s: Function y1() unavailable", NULL, 0);
    errno = EDOM;
    return 0;
#else
    errno = 0;
    return(Ft_errcheck(y1(d), "besy1"));
#endif
}

double Ft_Yn(double i, double d)
{
#ifdef AMIGA
    Ft_matherror("%s: Function yn() unavailable", NULL, 0);
    errno = EDOM;
    return 0;
#else
    errno = 0;
    return(Ft_errcheck(yn((int)i, d), "besyn"));
#endif
}

double Ft_Jn(double i, double d)
{
#ifdef AMIGA
    Ft_matherror("%s: Function jn() unavailable", NULL, 0);
    errno = EDOM;
    return 0;
#else
    errno = 0;
    return(Ft_errcheck(jn((int)i, d), "besjn"));
#endif
}

double Ft_errcheck(double x, char *str)
{
	extern int Ft_Check;

	if (!Ft_Check)
		return(x);
    if (errno == EDOM && Ft_Check & EDOM_CHK) {
        errno = 0;
        Ft_matherror("%s: Argument out of domain.", str, 0);
    }
    else if (errno == ERANGE && Ft_Check & ERANGE_CHK) {
        errno = 0;
        Ft_matherror("%s: Result out of range.", str, 0);
    }
#if defined (IsNaNorINF)
	else if (IsNaNorINF(x) && (Ft_Check & NAN_CHK || Ft_Check & INF_CHK)) {
# if defined (IsINF)
		if (IsINF(x) && Ft_Check & INF_CHK) {
        	Ft_matherror("%s: Result infinite.", str, 0);
		}
# endif  /* IsINF */
        Ft_matherror("%s: Result not a number.", str, 0);
	}
#else  /* IsNaNorINF */
# if defined (NaN)
	else if (NaN(x) && Ft_Check & NAN_CHK) {
        Ft_matherror("%s: Result not a number.", str, 0);
	}
# else  /* NaN */
#  if defined (IS_NAN)
	else if (IS_NAN(x) && Ft_Check & NAN_CHK) {
        Ft_matherror("%s: Result not a number.", str, 0);
	}
#  endif  /* IS_NAN */
/* nothing */ /* SUNOS 3.5 does not support #elif */
# endif  /* NaN */
# if defined (IS_INF)
	else if (IS_INF(x) && Ft_Check & INF_CHK) {
        Ft_matherror("%s: Result infinite.", str, 0);
	}
# endif  /* IS_INF */
/* nothing */ /* SUNOS 3.5 does not support #elif */
#endif  /* IsNaNorINF */
	else {
    	return(x);
	}
	Ft_matherror("Internal error: Impossible case in Ft_errcheck().", NULL, 0);
	return(ERRR); /* DUMMY */
}

double Ft_dbscan(char *s1, char *s2)
{
	double dd;

	if (sscanf(s1, s2, &dd) != 1) {
		Ft_matherror("scan: Wrong assignment \"%s\".", s2, 0);
	}
	return(dd);
}

double Ft_octal(double x)
{
	long i = (long) rint(x);
	char str[128];

	sprintf(str, "%d", i);
	i = strtol(str, (char **)NULL, 8);
	return((double)i);
}

double Ft_minimum(double x, double y)
{
	if (x < y)
		return(x);
	return(y);
}

double Ft_maximum(double x, double y)
{
	if (x > y)
		return(x);
	return(y);
}

#include "symbol.h"

double Ft_sum(double *vec)
{
	int i, ndata;
	double val=0.0;
	extern double *Ft_Data;

	ndata = (int) *Ft_Data;
	for (i=1; i<=ndata; i++)
		val += vec[i];	

	return(val);
}

#include <string.h>

double Ft_vread(void)
{
	double value;
	int ret = 0;
	extern char Ft_Inname[];
	extern FILE *Ft_Inread;

    if (Ft_Inread == stdin) {
		while (ret != 1) {
            fputs("vread? ", stderr);
        	ret = fscanf(Ft_Inread, "%lf", &value);
			if (ret < 0)
				Ft_matherror("vread: stdin EOF encountered.", NULL, 0);
			fprintf(stderr, "vread: Bad entry: Flushing...\n");
			fflush(stdin);
        }
    }
	else {
		ret = fscanf(Ft_Inread, "%lf", &value);
		if (ret < 0) {
			fprintf(stderr,"vread: Reached end of file \"%s\".", Ft_Inname);
			Ft_Inread = stdin;
			strcpy(Ft_Inname, "stdin");
			Ft_matherror("Resetting input to stdin...", NULL, 0);
		}
		else if (ret == 0) {
			Ft_matherror("vread: Could not read value from file \"%s\"\n",
			Ft_Inname, 0);
		}
	}
	return(value);
}