/* sin, cos, etc, for S-Lang */
/* Copyright (c) 1992, 1995 John E. Davis
 * All rights reserved.
 *
 * You may distribute under the terms of either the GNU General Public
 * License or the Perl Artistic License.
 */

#include "config.h"
#include "sl-feat.h"

#include <stdio.h>
#include <math.h>

#include "slang.h"
#include "_slang.h"

#ifdef PI
# undef PI
#endif
#define PI 3.14159265358979323846


#ifndef HAVE_STDLIB_H
  extern double atof ();
#endif

#if defined(__unix__)
#include <signal.h>
#include <errno.h>

#define SIGNAL  SLsignal

static void math_floating_point_exception (int sig)
{
   sig = errno;
   SLang_Error = SL_INTRINSIC_ERROR;
   (void) SIGNAL (SIGFPE, math_floating_point_exception);
   errno = sig;
}
#endif

double SLmath_hypot (double x, double y)
{
   double fr, fi, ratio;

   fr = fabs(x);
   fi = fabs(y);

   if (fr > fi)
     {
	ratio = y / x;
	x = fr * sqrt (1.0 + ratio * ratio);
     }
   else if (fi == 0.0) x = 0.0;
   else
     {
	ratio = x / y;
	x = fi * sqrt (1.0 + ratio * ratio);
     }

   return x;
}

/* usage here is a1 a2 ... an n x ==> a1x^n + a2 x ^(n - 1) + ... + an */
static double math_poly (void)
{
   int n;
   int dum1, dum2;
   double xn = 1.0, sum = 0.0;
   double an, x;

   if ((SLang_pop_double(&x, &dum1, &dum2))
       || (SLang_pop_integer(&n))) return(0.0);

   while (n-- > 0)
     {
	if (SLang_pop_double(&an, &dum1, &dum2)) break;
	(void) dum1; (void) dum2;
	sum += an * xn;
	xn = xn * x;
     }
   return (double) sum;
}

static int double_math_op_result (int op, unsigned char a, unsigned char *b)
{
   (void) a;
   (void) op;
   *b = SLANG_DOUBLE_TYPE;
   return 1;
}

static int double_math_op (int op,
			   unsigned char type, VOID_STAR ap, unsigned int na,
			   VOID_STAR bp)
{
   double *a, *b;
   unsigned int i;

   (void) type;
   a = (double *) ap;
   b = (double *) bp;

   switch (op)
     {
      default:
	return 0;

      case SLMATH_SINH:
	for (i = 0; i < na; i++)
	  b[i] = sinh (a[i]);
	break;

#ifdef HAVE_ASINH
      case SLMATH_ASINH:
	for (i = 0; i < na; i++)
	  b[i] = asinh (a[i]);
	break;
#endif

      case SLMATH_COSH:
	for (i = 0; i < na; i++)
	  b[i] = cosh (a[i]);
	break;
	
#ifdef HAVE_ACOSH
      case SLMATH_ACOSH:
	for (i = 0; i < na; i++)
	  b[i] = acosh (a[i]);
	break;
#endif
      case SLMATH_TANH:
	for (i = 0; i < na; i++)
	  b[i] = tanh (a[i]);
	break;
#ifdef HAVE_ATANH
      case SLMATH_ATANH:
	for (i = 0; i < na; i++)
	  b[i] = atanh (a[i]);
	break;
#endif
      case SLMATH_TAN:
	for (i = 0; i < na; i++)
	  b[i] = tan (a[i]);
	break;

      case SLMATH_ASIN:
	for (i = 0; i < na; i++)
	  b[i] = asin (a[i]);
	break;

      case SLMATH_ACOS:
	for (i = 0; i < na; i++)
	  b[i] = acos (a[i]);
	break;

      case SLMATH_ATAN:
	for (i = 0; i < na; i++)
	  b[i] = atan (a[i]);
	break;

      case SLMATH_EXP:
	for (i = 0; i < na; i++)
	  b[i] = exp (a[i]);
	break;

      case SLMATH_LOG:
	for (i = 0; i < na; i++)
	  b[i] = log (a[i]);
	break;

      case SLMATH_LOG10:
	for (i = 0; i < na; i++)
	  b[i] = log10 (a[i]);
	break;

      case SLMATH_SQRT:
	for (i = 0; i < na; i++)
	  b[i] = sqrt (a[i]);
	break;

      case SLMATH_SIN:
	for (i = 0; i < na; i++)
	  b[i] = sin (a[i]);
	break;

      case SLMATH_COS:
	for (i = 0; i < na; i++)
	  b[i] = cos (a[i]);
	break;

      case SLMATH_REAL:
	for (i = 0; i < na; i++)
	  b[i] = a[i];
	break;

      case SLMATH_IMAG:
	for (i = 0; i < na; i++)
	  b[i] = 0.0;
	break;
     }

   return 1;
}

static int int_math_op (int op,
			unsigned char type, VOID_STAR ap, unsigned int na,
			VOID_STAR bp)
{
   double *b;
   int *a;
   unsigned int i;

   (void) type;
   a = (int *) ap;
   b = (double *) bp;

   switch (op)
     {
      default:
	return 0;

      case SLMATH_SINH:
	for (i = 0; i < na; i++)
	  b[i] = sinh ((double) a[i]);
	break;
#ifdef HAVE_ASINH
      case SLMATH_ASINH:
	for (i = 0; i < na; i++)
	  b[i] = asinh ((double) a[i]);
	break;
#endif
      case SLMATH_COSH:
	for (i = 0; i < na; i++)
	  b[i] = cosh ((double) a[i]);
	break;
#ifdef HAVE_ACOSH
      case SLMATH_ACOSH:
	for (i = 0; i < na; i++)
	  b[i] = acosh ((double) a[i]);
	break;
#endif
      case SLMATH_TANH:
	for (i = 0; i < na; i++)
	  b[i] = tanh ((double) a[i]);
	break;
#ifdef HAVE_ATANH
      case SLMATH_ATANH:
	for (i = 0; i < na; i++)
	  b[i] = atanh ((double) a[i]);
	break;
#endif
      case SLMATH_TAN:
	for (i = 0; i < na; i++)
	  b[i] = tan ((double) a[i]);
	break;

      case SLMATH_ASIN:
	for (i = 0; i < na; i++)
	  b[i] = asin ((double) a[i]);
	break;

      case SLMATH_ACOS:
	for (i = 0; i < na; i++)
	  b[i] = acos ((double) a[i]);
	break;

      case SLMATH_ATAN:
	for (i = 0; i < na; i++)
	  b[i] = atan ((double) a[i]);
	break;

      case SLMATH_EXP:
	for (i = 0; i < na; i++)
	  b[i] = exp ((double) a[i]);
	break;

      case SLMATH_LOG:
	for (i = 0; i < na; i++)
	  b[i] = log ((double) a[i]);
	break;

      case SLMATH_LOG10:
	for (i = 0; i < na; i++)
	  b[i] = log10 ((double) a[i]);
	break;

      case SLMATH_SQRT:
	for (i = 0; i < na; i++)
	  b[i] = sqrt ((double) a[i]);
	break;

      case SLMATH_SIN:
	for (i = 0; i < na; i++)
	  b[i] = sin ((double) a[i]);
	break;

      case SLMATH_COS:
	for (i = 0; i < na; i++)
	  b[i] = cos ((double) a[i]);
	break;

      case SLMATH_REAL:
	for (i = 0; i < na; i++)
	  b[i] = a[i];
	break;

      case SLMATH_IMAG:
	for (i = 0; i < na; i++)
	  b[i] = 0.0;
	break;
     }

   return 1;
}

#if SLANG_HAS_COMPLEX
static int complex_math_op_result (int op, unsigned char a, unsigned char *b)
{
   (void) a;
   switch (op)
     {
      default:
	*b = SLANG_COMPLEX_TYPE;
	break;

      case SLMATH_REAL:
      case SLMATH_IMAG:
	*b = SLANG_DOUBLE_TYPE;
	break;
     }
   return 1;
}

static int complex_math_op (int op,
			    unsigned char type, VOID_STAR ap, unsigned int na,
			    VOID_STAR bp)
{
   double *a, *b;
   unsigned int i;

   (void) type;
   a = (double *) ap;
   b = (double *) bp;

   switch (op)
     {
      default:
      case SLMATH_ATANH:
      case SLMATH_ACOSH:
      case SLMATH_ASINH:
	return 0;

      case SLMATH_REAL:
	for (i = 0; i < na; i++)
	  b[i] = a[2 * i];
	break;

      case SLMATH_IMAG:
	for (i = 0; i < na; i++)
	  b[i] = a[2 * i + 1];
	break;

      case SLMATH_EXP:
	for (i = 0; i < na; i += 2)
	  SLcomplex_exp (b + i, a + i);
	break;

      case SLMATH_LOG:
	for (i = 0; i < na; i += 2)
	  SLcomplex_log (b + i, a + i);
	break;

      case SLMATH_LOG10:
	for (i = 0; i < na; i += 2)
	  SLcomplex_log10 (b + i, a + i);
	break;

      case SLMATH_SQRT:
	for (i = 0; i < na; i += 2)
	  SLcomplex_sqrt (b + i, a + i);
	break;

      case SLMATH_SIN:
	for (i = 0; i < na; i += 2)
	  SLcomplex_sin (b + i, a + i);
	break;

      case SLMATH_COS:
	for (i = 0; i < na; i += 2)
	  SLcomplex_cos (b + i, a + i);
	break;

      case SLMATH_SINH:
	for (i = 0; i < na; i += 2)
	  SLcomplex_sinh (b + i, a + i);
	break;

      case SLMATH_COSH:
	for (i = 0; i < na; i += 2)
	  SLcomplex_cosh (b + i, a + i);
	break;

      case SLMATH_TANH:
	for (i = 0; i < na; i += 2)
	  SLcomplex_tanh (b + i, a + i);
	break;

      case SLMATH_TAN:
	for (i = 0; i < na; i += 2)
	  SLcomplex_tan (b + i, a + i);
	break;

      case SLMATH_ASIN:
	for (i = 0; i < na; i += 2)
	  SLcomplex_asin (b + i, a + i);
	break;

      case SLMATH_ACOS:
	for (i = 0; i < na; i += 2)
	  SLcomplex_acos (b + i, a + i);
	break;

      case SLMATH_ATAN:
	for (i = 0; i < na; i += 2)
	  SLcomplex_atan (b + i, a + i);
	break;
     }

   return 1;
}
#endif

static double Const_E =  2.718281828459045;
static double Const_Pi = 3.141592653589793;

static SLang_Math_Unary_Type SLmath_Table [] =
{
   MAKE_MATH_UNARY("sinh", SLMATH_SINH),
   MAKE_MATH_UNARY("asinh", SLMATH_ASINH),
   MAKE_MATH_UNARY("cosh", SLMATH_COSH),
   MAKE_MATH_UNARY("acosh", SLMATH_ACOSH),
   MAKE_MATH_UNARY("tanh", SLMATH_TANH),
   MAKE_MATH_UNARY("atanh", SLMATH_ATANH),
   MAKE_MATH_UNARY("sin", SLMATH_SIN),
   MAKE_MATH_UNARY("cos", SLMATH_COS),
   MAKE_MATH_UNARY("tan", SLMATH_TAN),
   MAKE_MATH_UNARY("atan", SLMATH_ATAN),
   MAKE_MATH_UNARY("acos", SLMATH_ACOS),
   MAKE_MATH_UNARY("asin", SLMATH_ASIN),
   MAKE_MATH_UNARY("exp", SLMATH_EXP),
   MAKE_MATH_UNARY("log", SLMATH_LOG),
   MAKE_MATH_UNARY("sqrt", SLMATH_SQRT),
   MAKE_MATH_UNARY("log10", SLMATH_LOG10),
#if SLANG_HAS_COMPLEX
   MAKE_MATH_UNARY("Real", SLMATH_REAL),
   MAKE_MATH_UNARY("Imag", SLMATH_IMAG),
#endif
   SLANG_END_TABLE
};

int SLang_init_slmath (void)
{
#if defined(__unix__)
   (void) SIGNAL (SIGFPE, math_floating_point_exception);
#endif
   if ((-1 == SLclass_add_math_op (SLANG_INT_TYPE, int_math_op, double_math_op_result))
       || (-1 == SLclass_add_math_op (SLANG_DOUBLE_TYPE, double_math_op, double_math_op_result))
#if SLANG_HAS_COMPLEX
       || (-1 == SLclass_add_math_op (SLANG_COMPLEX_TYPE, complex_math_op, complex_math_op_result))
#endif
       )
     return -1;

   if ((-1 == SLadd_math_unary_table (SLmath_Table, "__SLMATH__"))
       || (-1 == SLadd_intrinsic_variable ("E", (VOID_STAR)&Const_E, SLANG_DOUBLE_TYPE, 1))
       || (-1 == SLadd_intrinsic_variable ("PI", (VOID_STAR)&Const_Pi, SLANG_DOUBLE_TYPE, 1))
       || (-1 == SLadd_intrinsic_function ("polynom", (FVOID_STAR) math_poly, SLANG_DOUBLE_TYPE, 0)))
     return -1;
   
   return 0;
}