/* I included parts of math_68881.h from gcc 1.37.1 in here. HP 24-Jun-90 */
/* This file is missing the following functions, which i didn't want to
   include:
   atanh, expm1, log1p, rint, drem, scalb,
   logb, ldexp, frexp
*/
/* Please note that using gcc -traditional will disable inline expansion
   of these functions, since they are defined with ansi-style delarations. HP */
/* Please note also that using -ansi -pedantic will give warnings because
   I do try to use all gcc extensions giving more speed. HP */

#ifndef __MATH_H__
#define __MATH_H__

extern int errno;

#ifdef __GNUC__

#define fabs(x) __builtin_fabs(x)
#define _FABS

#if defined(__mc68881__) && defined(__mc68020__) && defined(__STDC__)

extern double atof(const char *);

#ifndef _FABS
__inline static const double fabs (double x)
{
  double value;

  __asm ("fabs%.x %1,%0"
	 : "=f" (value)
	 : "f" (x));
  return value;
}
#endif

__inline static const double sin (double x)
{
  double value;

  __asm ("fsin%.x %1,%0"
	 : "=f" (value)
	 : "f" (x));
  return value;
}

__inline static const double sinh (double x)
{
  double value;

  __asm ("fsinh%.x %1,%0"
	 : "=f" (value)
	 : "f" (x));
  return value;
}

__inline static const double cos (double x)
{
  double value;

  __asm ("fcos%.x %1,%0"
	 : "=f" (value)
	 : "f" (x));
  return value;
}

__inline static const double cosh (double x)
{
  double value;

  __asm ("fcosh%.x %1,%0"
	 : "=f" (value)
	 : "f" (x));
  return value;
}

__inline static const double tan (double x)
{
  double value;

  __asm ("ftan%.x %1,%0"
	 : "=f" (value)
	 : "f" (x));
  return value;
}

__inline static const double tanh (double x)
{
  double value;

  __asm ("ftanh%.x %1,%0"
	 : "=f" (value)
	 : "f" (x));
  return value;
}

__inline static const double asin (double x)
{
  double value;

  __asm ("fasin%.x %1,%0"
	 : "=f" (value)
	 : "f" (x));
  return value;
}

__inline static const double acos (double x)
{
  double value;

  __asm ("facos%.x %1,%0"
	 : "=f" (value)
	 : "f" (x));
  return value;
}

__inline static const double atan (double x)
{
  double value;

  __asm ("fatan%.x %1,%0"
	 : "=f" (value)
	 : "f" (x));
  return value;
}

__inline static const double atan2 (double y, double x)
{
  double pi, pi_over_2;

  __asm ("fmovecr%.x %#0,%0"		/* extended precision pi */
	 : "=f" (pi)
	 : /* no inputs */ );
  __asm ("fscale%.b %#-1,%0"		/* no loss of accuracy */
	 : "=f" (pi_over_2)
	 : "0" (pi));
  if (x > 0)
    {
      if (y > 0)
	{
	  if (x > y)
	    return atan (y / x);
	  else
	    return pi_over_2 - atan (x / y);
	}
      else
	{
	  if (x > -y)
	    return atan (y / x);
	  else
	    return - pi_over_2 - atan (x / y);
	}
    }
  else
    {
      if (y > 0)
	{
	  if (-x > y)
	    return pi + atan (y / x);
	  else
	    return pi_over_2 - atan (x / y);
	}
      else
	{
	  if (-x > -y)
	    return - pi + atan (y / x);
	  else if (y < 0)
	    return - pi_over_2 - atan (x / y);
	  else
	    {
	      double value;

	      errno = 0x43;   /* EDOM */
	      __asm ("fmove%.d %#0x7fffffffffffffff,%0" 	/* quiet NaN */
		     : "=f" (value)
		     : /* no inputs */);
	      return value;
	    }
	}
    }
}

__inline static const double exp (double x)
{
  double value;

  __asm ("fetox%.x %1,%0"
	 : "=f" (value)
	 : "f" (x));
  return value;
}


__inline static const double log (double x)
{
  double value;

  __asm ("flogn%.x %1,%0"
	 : "=f" (value)
	 : "f" (x));
  return value;
}

__inline static const double log10 (double x)
{
  double value;

  __asm ("flog10%.x %1,%0"
	 : "=f" (value)
	 : "f" (x));
  return value;
}

__inline static const double fmod (double x, double y)
{
  double value;

  __asm ("fmod%.x %2,%0"
	 : "=f" (value)
	 : "0" (x),
	   "f" (y));
  return value;
}

__inline static const double sqrt (double x)
{
  double value;

  __asm ("fsqrt%.x %1,%0"
	 : "=f" (value)
	 : "f" (x));
  return value;
}

__inline static const double pow (const double x, const double y)
{
  if (x > 0)
    return exp (y * log (x));
  else if (x == 0)
    {
      if (y > 0)
	return 0.0;
      else
#ifdef __OSK__
	/* We currently don't handle this error to be compatible with
	   the library routine. I may change this in the future.
	   HP 24-Jun-90 */
	return 1.0;
#else
	{
	  double value;

	  errno = EDOM;
	  __asm ("fmove%.d %#0rnan,%0"		/* quiet NaN */
		 : "=f" (value)
		 : /* no inputs */);
	  return value;
	}
#endif
    }
  else
    {
      double temp;

      __asm ("fintrz%.x %1,%0"
	     : "=f" (temp)			/* integer-valued float */
	     : "f" (y));
      if (y == temp)
        {
	  int i = (int) y;

	  temp = exp (y * log (-x));
	  if ((i & 1) == 0)			/* even */
	    return temp;
	  else
	    return -temp;
        }
      else
#ifdef __OSK__
	/* We currently don't handle this error to be compatible with
	   the library routine. I may change this in the future.
	   HP 24-Jun-90 */
	return 0.0;
#else
        {
	  double value;

	  errno = EDOM;
	  __asm ("fmove%.d %#0rnan,%0"		/* quiet NaN */
		 : "=f" (value)
		 : /* no inputs */);
	  return value;
        }
#endif
    }
}

__inline static const double ceil (double x)
{
  int rounding_mode, round_up;
  double value;

  __asm volatile ("fmove%.l fpcr,%0"
		  : "=d" (rounding_mode)
		  : /* no inputs */ );
  round_up = rounding_mode | 0x30;
  __asm volatile ("fmove%.l %0,fpcr"
		  : /* no outputs */
		  : "di" (round_up));
  __asm volatile ("fint%.x %1,%0"
		  : "=f" (value)
		  : "f" (x));
  __asm volatile ("fmove%.l %0,fpcr"
		  : /* no outputs */
		  : "di" (rounding_mode));
  return value;
}

__inline static const double floor (double x)
{
  int rounding_mode, round_down;
  double value;

  __asm volatile ("fmove%.l fpcr,%0"
		  : "=d" (rounding_mode)
		  : /* no inputs */ );
  round_down = (rounding_mode & ~0x10)
		| 0x20;
  __asm volatile ("fmove%.l %0,fpcr"
		  : /* no outputs */
		  : "di" (round_down));
  __asm volatile ("fint%.x %1,%0"
		  : "=f" (value)
		  : "f" (x));
  __asm volatile ("fmove%.l %0,fpcr"
		  : /* no outputs */
		  : "di" (rounding_mode));
  return value;
}

__inline static double modf (double x, double *ip)
{
  double temp;

  __asm ("fintrz%.x %1,%0"
	 : "=f" (temp)			/* integer-valued float */
	 : "f" (x));
  *ip = temp;
  return x - temp;
}

#if 1
/* Standard definition of hypot. I don't change the formula to watch for
   unwarranted overflows since the calculation will be done in extended
   precision registers. This may fail if you use -ffloat-store ... HP */

__inline static const double hypot (double x, double y)
{
  register double temp1,temp2;

  __asm ("fmul%.x %0,%0"
	 : "=f" (temp1)
	 : "0" (x));
  __asm ("fmul%.x %0,%0"
	 : "=f" (temp2)
	 : "0" (y));
  __asm ("fadd%.x %2,%0"
	 : "=f" (temp1)
	 : "0" (temp1), "f" (temp2));
  __asm ("fsqrt%.x %1,%0"
	 : "=f" (temp1)
	 : "0" (temp1));

  return temp1;
}
#else
/* Modified version of hypot (untested). This should do the trick
   for -ffloat-store. Also its slower and longer. HP */
/* Oh well, I just tried to test it and found -ffloat-store to be
   buggy. So i suggest for the time being to not use -ffloat-store.
   So the standard version will be just perfect. HP */
__inline static const double hypot (double x,double y)
{
  double temp1=fabs(x),temp2=fabs(y);

  if( temp2 < temp1 )
  {
    double temp3;
    temp3 = temp2;
    temp2 = temp1;
    temp1 = temp3;
  }
  if( temp2 == 0.0 )
    return 0.0;
  return( sqrt( (temp1/temp2)*(temp1/temp2) + 1.0 ) * temp2 );
}
#endif
/* This is different from the other definitions since it generates
   infinity instead of a large real number. HP */
#define HUGE								\
({									\
  double huge_val;							\
									\
  __asm ("fmove%.d %#0x7ff0000000000000,%0"	/* Infinity */		\
	 : "=f" (huge_val)						\
	 : /* no inputs */);						\
  huge_val;								\
})

/* This might be faster or slower than a constant depending on the
   context. I define it anyway, since it may give higher precision. HP */
#define PI								\
({									\
  double pi;								\
									\
  __asm ("fmovecr%.x %#0,%0"	/* extended precision pi */		\
	 : "=f" (pi)							\
	 : /* no inputs */);						\
  pi;									\
})

#else /* NOT mc68881 */

#ifdef __STDC__
extern double atof(const char *);
extern const double sin(double),cos(double),tan(double),asin(double);
extern const double acos(double),atan(double),exp(double),log(double);
extern const double log10(double),pow(double,double),sqrt(double);
extern const double floor(double),ceil(double),hypot(double,double);
extern double modf(double,double *);
#else /* NOT __STDC__ i.e. -traditional */
extern double atof();
extern double sin(),cos(),tan(),asin(),acos(),atan();
extern double exp(),log(),log10(),pow(),sqrt();
extern double floor(),ceil(),modf(),hypot();
#endif /* __STDC__ */

#define HUGE	1.701411733192644270e38
#define PI	3.141592653589796

#endif

#else /* NOT GNU */

extern double atof();
extern double sin(),cos(),tan(),asin(),acos(),atan();
extern double exp(),log(),log10(),pow(),sqrt();
extern double floor(),ceil(),modf(),hypot();
extern double fabs();

#define HUGE	1.701411733192644270e38
#define PI	3.141592653589796

#endif /* NOT GNU */
#endif /* __MATH_H__ */