/*
 * PROJECT C Library, X68000 PROGRAMMING INTERFACE DEFINITION
 * --------------------------------------------------------------------
 * This file is written by the Project C Library Group,  and completely
 * in public domain. You can freely use, copy, modify, and redistribute
 * the whole contents, without this notice.
 * --------------------------------------------------------------------
 * $Id: $
 */

#ifndef __sys_math_f_h__
#define __sys_math_f_h__

#ifndef __GNUC__
#error You lose. This file can be compiled only by GNU-C compiler.
#endif

#ifndef __cdecl_h__
#include <cdecl.h>
#endif


#define __INLINE_SQRT(x)				\
    (__extension__ ({					\
	double f;					\
	__asm ("fsqrt.x %1,%0" : "=f" (f) : "f" (x));	\
	f;						\
    }))

#define __INLINE_LOG(x)					\
    (__extension__ ({					\
	double f;					\
	__asm ("flogn.x %1,%0" : "=f" (f) : "f" (x));	\
	f;						\
    }))


static __inline __const double __inline_acos (double x)
{
    double value;

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

    return value;
}

static __inline __const double __inline_acosh (double x)
{
    return __INLINE_LOG (x + __INLINE_SQRT ((x * x) - 1.0));
}

static __inline __const double __inline_asin (double x)
{
    double value;

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

    return value;
}

static __inline __const double __inline_asinh (double x)
{
    return __INLINE_LOG (x + __INLINE_SQRT ((x * x) + 1.0));
}

static __inline __const double __inline_atan (double x)
{
    double value;

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

    return value;
}

static __inline __const double __inline_atan2 (double y, double x)
{
    double pi;

    if (x > 0.0) {
	__asm ("fatan.x %1,%0" : "=f" (x) : "f" (y / x));
	return x;
    }

    __asm ("fmovecr.x #$0,%0" : "=f" (pi));

    if (x < 0.0) {
	if (y > 0.0) {
	    __asm ("fatan.x %1,%0" : "=f" (x):"f" (y / x));
	    return x + pi;
	}
	else if (y < 0.0) {
	    __asm ("fatan.x %1,%0" : "=f" (x):"f" (y / x));
	    return x - pi;
	}
	else
	    return pi;
    }
    else {
	double pi_2;

	__asm ("fscale.w #-1,%0" : "=f" (pi_2) : "0" (pi));

	if (y < 0.0)
	    return -pi_2;
	else if (y > 0.0)
	    return +pi_2;
	else
	    return _NAN;
    }
}

static __inline __const double __inline_atanh (double x)
{
    return 0.50 * __INLINE_LOG ((1.0 + x) / (1.0 - x));
}

static __inline __const double __inline_ceil (double x)
{
    double value;

    __asm ("fmove.l fpcr,d0\n\t"
	   "fmove.l %2,fpcr\n\t"
	   "fint.x %1,%0\n\t"
	   "fmove.l d0,fpcr"
	   : "=f" (value)
	   : "f" (x), "dmi" (0x30)
	   : "d0");

    return value;
}

static __inline __const double __inline_cos (double x)
{
    double value;

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

    return value;
}

static __inline __const double __inline_cosh (double x)
{
    double value;

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

    return value;
}

static __inline __const double __inline_exp (double x)
{
    double value;

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

    return value;
}

static __inline __const double __inline_fabs (double x)
{
    double value;

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

    return value;
}

static __inline __const double __inline_floor (double x)
{
    double value;

    __asm ("fmove.l fpcr,d0\n\t"
	   "fmove.l %2,fpcr\n\t"
	   "fint.x %1,%0\n\t"
	   "fmove.l d0,fpcr"
	   : "=f" (value)
	   : "f" (x), "dmi" (0x20)
	   : "d0");

    return value;
}

static __inline __const double __inline_fmod (double x, double y)
{
    double value;

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

    return value;
}

static __inline __const double __inline_ldexp (double x, int_ exp)
{
    double value;

    __asm ("fscale.l %2,%0" : "=f" (value) : "0" (x), "g" (exp));

    return value;
}

static __inline __const double __inline_log (double x)
{
    double value;

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

    return value;
}

static __inline __const double __inline_log10 (double x)
{
    double value;

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

    return value;
}

static __inline __const double __inline_pow (double x, double y)
{
    if (x > 0.0)
	return exp (y * __INLINE_LOG (x));

    else if (x == 0.0)
	return (y > 0.0) ? x : _NAN;

    else {
	double yi;

	__asm ("fintrz.x %1,%0" : "=f" (yi) : "f" (y));

	if (yi == y)
	    return ((int) y) & 1 ? - exp (y * __INLINE_LOG (-x))
				 :   exp (y * __INLINE_LOG (-x));
	else
	    return _NAN;
    }
}

static __inline __const double __inline_sin (double x)
{
    double value;

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

    return value;
}

static __inline __const double __inline_sinh (double x)
{
    double value;

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

    return value;
}

static __inline __const double __inline_sqrt (double x)
{
    double value;

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

    return value;
}

static __inline __const double __inline_tan (double x)
{
    double value;

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

    return value;
}

static __inline __const double __inline_tanh (double x)
{
    double value;

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

    return value;
}

static __inline double __inline_frexp (double x, int_ *exp)
{
    xdouble_t value;

    if ((value.d = x) == 0.0)
	*exp = 0;
    else {
	*exp = value.ieee.exp - 0x03fe;
	value.ieee.exp = 0x3fe;
    }

    return value.d;
}

static __inline double __inline_modf (double x, double *iptr)
{
    double value;

    __asm ("fintrz.x %1,%0" : "=f" (value) : "f" (x));

    *iptr = value;
    return x - value;
}

static __inline __const int_ __inline_isinf (double x)
{
    int_ result;

    __asm ("ftst.x %1\n\t"
	   "fmove fpsr,%0\n\t"
	   "bfextu %0{#6:#1},%0"
		: "=d" (result)
		: "f" (x));

    return result;
}

static __inline __const int_ __inline_isnan (double x)
{
    int_ result;

    __asm ("ftst.x %1\n\t"
	   "fmove fpsr,%0\n\t"
	   "bfextu %0{#7:#1},%0"
		: "=d" (result)
		: "f" (x));

    return result;
}

#endif
