/**
 * -- IEEE 754/854 floating point extensions
 *
 * Copyright (C) 2005 David N. Williams
 *
 * @see LGPL 
 * @author David N. Williams              @(#) %derived_by: dnw %
 * @version %version: 0.5.2%
 *   (%date_modified: Mon Nov 07 08:39:00 2005 %)
 *
 * If you take advantage of the option in the LGPL to put a
 * particular version of this library under the GPL, the author
 * would regard it as polite if you would put any direct
 * modifications under the LGPL as well, and include a copy of
 * this request near the beginning of the modified library
 * source.  A "direct modification" is one that enhances or
 * extends the library in line with its original concept, as
 * opposed to developing a distinct application or library which
 * might use it.
 * 
 * @description
 *   This extension word set is an interface to some of the C99,
 *   IEEE 754 floating-point math extensions, including some
 *   floating-point environment, exception, and rounding-mode
 *   functions.  The system we use (darwin) does not comply with
 *   the C99 standard return values for all of the environment
 *   and exception functions, so we do not use those returns.
 *
 *   Remember that in pfe, the default floating-point size is
 *   that of a double.  So where relevant, that's the C99
 *   function that is used.
 *
 *   This implementation assumes that TOS and FTOS are in
 *   memory, not in registers.  It also requires a version of pfe
 *   recent enough to use the ROUND_MODE session variable to
 *   prevent longjmp() from clobbering the C99 rounding mode.
 */

#define _P4_SOURCE 1
#include <pfe/pfe-base.h>
#include <pfe/incl-sup.h>  /* for BITSOF in pfe-0.33.58 */

#include <float.h>
#include <math.h>
#include <fenv.h>
#if __STDC_VERSION__+0 > 199900L
#pragma STDC FENV_ACCESS ON  /* still ignored by gcc? */
#endif

/* Align BSD? with C99 */
#if !defined FE_DEFL_ENV &&  defined FE_DFL_ENV
#define FE_DEFL_ENV   FE_DFL_ENV
#endif

#include <pfe/double-sub.h>  /* p4_d_negate() */
#include "fpieee-ext.h"

#define CELLBITS   BITSOF (p4cell)

/****************************************************************/
/* elementary functions new in C99				*/
/****************************************************************/

/** F*+  (f: x y z -- y*z+x )
 */
FCode (p4_f_star_plus)
{
  FP[2] = fma ( FP[1], FP[0], FP[2]);
  FP += 2;
}

/** FCBRT  (f: x -- x^[1/3] )
 */
FCode (p4_f_cbrt)  { *FP = cbrt (*FP); }

/** FDIM  (f: x y -- x-y | +0 )
 */
FCode (p4_f_dim)
{
  FP[1] = fdim ( FP[1], FP[0]);
  FP++;
}

/** FEXP2  (f: x -- 2^x )
 */
FCode (p4_f_exp2)
{ *FP = exp2 (*FP); }

/** FLOG2  (f: x -- log.base.2[x] )
 */
FCode (p4_f_log2)
{ *FP = log2 (*FP); }

/** FERF  (f: x -- erf[x] )
 */
FCode (p4_f_erf)
{ *FP = erf (*FP); }

/** FERFC  (f: x -- 1-erf[x] )
 */
FCode (p4_f_erfc)
{ *FP = erfc (*FP); }

/** FGAMMA  (f: x -- gamma[x] )
 */
FCode (p4_f_gamma)
{ *FP = tgamma (*FP); }

/** FLNGAMMA  (f: x -- ln[gamma[x]] )
 */
FCode (p4_f_ln_gamma)
{ *FP = lgamma (*FP); }

/****************************************************************/
/* floating-point auxiliary functions				*/
/****************************************************************/

/** FPCLASSIFY (f: x -- s: type )
 * Return the appropriate value FP_INFINITE, FP_NAN, FP_NORMAL,
 * FP_SUBNORMAL, or FP_ZERO provided by the host C99 compiler. 
 * The corresponding Forth constants are FP-INFINITE, etc.
 */
FCode (p4_fpclassify)
{ *--SP = fpclassify (*FP++); }

/** ISFINITE (f: x -- s: flag )
 */
FCode (p4_isfinite)
{ *--SP = isfinite (*FP++) ? ~0 : 0 ; }

/** ISINF (f: x -- s: flag )
 */
FCode (p4_isinf)
{ *--SP = isinf (*FP++) ? ~0 : 0 ; }

/** ISNAN (f: x -- s: flag )
 */
FCode (p4_isnan)
{ *--SP = isnan (*FP++) ? ~0 : 0 ; }

/** ISNORMAL (f: x -- s: flag )
 */
FCode (p4_isnormal)
{ *--SP = isnormal (*FP++) ? ~0 : 0 ; }

/** COPYSIGN  (f: x y -- |x|*sgn[y] )
 */
FCode (p4_copysign)
{
  FP[1] = copysign ( FP[1], FP[0]);
  FP++ ;
}

/** NEXTAFTER (f: x y -- x.next )
 * Return the next machine representable number after x in the y
 * direction.
 */
FCode (p4_nextafter)
{
  FP[1] = nextafter ( FP[1], FP[0] );
  FP++;
}

/** SIGNBIT (f: x -- s: minus? )
 */
FCode (p4_signbit)
{ signbit (*FP++) ? (*--SP = ~0) : (*--SP = 0); }

/** FREXP  (f: x -- frac s: n )
 * Return the fractional part of x as a floating-point number,
 * frac = 0.0 or 0.5 <= frac < 1.0, and the power n such that
 * frac*2^n = x.
 */
FCode (p4_frexp)
{ *FP = frexp ( *FP, --SP ); }

/** LDEXP  (f: x s: n -- f: x*2^n )
 */
FCode (p4_ldexp)
{ *FP = ldexp ( *FP, *SP++ ); }

/** MODF  (f: x -- frac integ )
 * Return the integer part of x as a double, and the fractional
 * part f of x with |f| < 1, both having the same sign as x, so
 * that x = frac + integ.
 */
FCode (p4_modf)
{
  double integ;

  *FP = modf (*FP, &integ);
  *--FP = integ;
}

/** MODFD  (f: x -- frac s: d )
 * A version of MODF that returns the integer part as a
 * double-cell signed integer.  Not in C99.
 * TENTATIVE!
 */
FCode (p4_modfd)
{
  /* based on F>D in floating-ext.c */
  double integ, hi, lo;
  int sign;

  *FP = modf (*FP, &integ);

  if (integ < 0)
    sign = 1, integ = -integ;
  else
    sign = 0;
  lo = modf (ldexp (integ, -CELLBITS), &hi);
  SP -= 2;
  SP[0] = (p4ucell) hi;
  SP[1] = (p4ucell) ldexp (lo, CELLBITS);
  if (sign)
    p4_d_negate ((p4dcell *) &SP[0]);
}

/** SCALBN  (f: x s: n -- f: x*b^n )
 * The output is efficiently scaled by b^n, with b = FLT_RADIX,
 * which is also the value of the Forth constant FLT-RADIX.
 */
FCode (p4_scalbn)
{ *FP = scalbn ( *FP, *SP++ ); }

/** LOGB  (f: x -- f: b.exponent )
 * Leave the exponent (radix b = FLT_RADIX) of the
 * floating-point representation.
 */
FCode (p4_logb)
{ *FP = logb (*FP); }

/** ILOGB  (f: x -- s: b.exponent )
 * Leave the exponent (radix b = FLT_RADIX) of the
 * floating-point representation, with special values:
 *    x      exp
 *    0.0    FP_ILOGB0
 *    Inf    INT_MAX
 *   -Inf    INT_MAX
 *    Nan    FP_ILOGBNAN
 * The special output constants are defined elsewhere with Forth
 * names FP-ILOGB0, FP-ILOGBNAN, INT-MAX.
 */
FCode (p4_ilogb)
{ *--SP = logb (*FP++); }

/****************************************************************/
/* floating-point environment					*/
/****************************************************************/

/* In the following stack effects, "fenv" indicates
 * implementation-defined control mode and exception bit data in
 * a C99 fenv_t type for the floating-point state.  An ambiguous
 * condition exists if fenv_t does not fit into a cell.  The
 * four C99 functions fegetenv(), fesetenv(), feholdexcept(),
 * and feupdateenv() return zero if successful and a nonzero
 * value if not, which we ignore, because some older systems
 * (BSD?) return void instead.
 */

/** FE-DEFL-ENV  ( -- default.fenv)
 * Leave the default floating-point environment. 
 */
FCode (p4_fe_defl_env)
{ *--SP = (fenv_t) *FE_DEFL_ENV; }

/** FEGETENV  ( -- fenv )
 * Leave the current floating-point environment.
 */
FCode (p4_fegetenv)
{ fegetenv ( (fenv_t*) --SP ); }

/** FESETENV  ( fenv -- )
 * Set the floating-point environment to fenv.
 */
FCode (p4_fesetenv)
{ fesetenv ( (fenv_t*) SP++ ); }

/** FEHOLDEXCEPT  ( -- fenv )
 * Leave the current floating-point environment for later
 * restoration, and install an environment that turns off
 * floating-point exceptions.
 */
FCode (p4_feholdexcept)
{ feholdexcept ( (fenv_t*) --SP ); }

/** FEUPDATEENV  ( fenv -- )
 * Set the floating-point environment to fenv, and raise the
 * exceptions that were up just before the call.
 */
FCode (p4_feupdateenv)
{ feupdateenv ( (fenv_t*) SP++ ); }

/****************************************************************/
/* floating-point exceptions					*/
/****************************************************************/

/* In the following stack effects, "fexcept" indicates an
 * implementation-defined floating-point exception bit mask in a
 * C99 fexcept_t type, while "excepts" indicates exception flag
 * data corresponding to a combination of the masks in a C int
 * type.  An ambiguous condition exists if either fexcept_t or
 * int does not fit into a cell.  The four C99 functions
 * fegetexceptflag(), fesetexceptflag(), feraiseexcept(), and
 * feclearexcept() return zero if successful and a nonzero value
 * if not, which we ignore, because some older systems (BSD?)
 * return void instead.
 *
 * Forth constants FE-DIVBYZERO, FE-INEXACT, FE-INVALID,
 * FE-OVERFLOW, FE-UNDERFLOW, FE-ALL-EXCEPT, corresponding to
 * the C99 exception masks named with "_" instead of "-". are
 * defined under P4_LISTWORDS.
 */

/** FEGETEXCEPTFLAG  ( excepts -- fenv )
 */
FCode (p4_fegetexceptflag)
{ fegetexceptflag ( (fexcept_t*) SP, *SP ); }

/** FESETEXCEPTFLAG  ( fenv excepts -- )
 */
FCode (p4_fesetexceptflag)
{
  fesetexceptflag ( (fexcept_t*) &SP[1], SP[0] );
  SP += 2;
}

/** FETESTEXCEPT  ( excepts -- flags.excepts )
 */
FCode (p4_fetestexcept)
{ *SP = fetestexcept (*SP); }

/** FERAISEEXCEPT  ( excepts -- )
 */
FCode (p4_feraiseexcept)
{ feraiseexcept (*SP++); }

/** FECLEAREXCEPT ( excepts -- )
 */
FCode (p4_feclearexcept)
{ feclearexcept (*SP++); }

/****************************************************************/
/* floating-point rounding modes				*/
/****************************************************************/

/* The Forth constants FE-DOWNWARD, FE-UPWARD, FE-TONEAREST,
 * FE-TOWARDZERO, corresponding to C99 rounding mode values
 * with "_" instead of "-" in their names, are defined under
 * P4_LISTWORDS.  They represent the legal rounding modes.
 */

/** FEGETROUND  ( -- round )
 * Leave the current rounding mode value, or a negative number
 * if unsuccessful.
 */
FCode (p4_fegetround)
{
  *--SP = fegetround ();
}

/** FESETROUND  ( round -- failure.flag )
 * Set the rounding mode to round.  In C99, a zero return
 * means success.  If round is illegal, the rounding mode is
 * unchanged, and a nonzero value is returned.
 *
 * When successful, the new rounding mode is cached in the
 * session variable ROUND_MODE for anticlobbering in case gcc's
 * longjmp() misbehaves in the pfe interpreter.  The user is not
 * supposed to access this variable.
 */
FCode (p4_fesetround)
{
  if ( !(*SP = fesetround (*SP)) )
    ROUND_MODE = fegetround ();
}


P4_LISTWORDS (fpieee) =
{
  P4_NEED ("floating-ext"),
  P4_INTO ("EXTENSIONS", 0),
  /* elementary functions new in C99 */
  P4_FXco ("F*+",		p4_f_star_plus),
  P4_FXco ("FCBRT",		p4_f_cbrt),
  P4_FXco ("FDIM",		p4_f_dim),
  P4_FXco ("FEXP2",		p4_f_exp2),
  P4_FXco ("FLOG2",		p4_f_log2),
  P4_FXco ("FERF",		p4_f_erf),
  P4_FXco ("FERFC",		p4_f_erfc),
  P4_FXco ("FLNGAMMA",		p4_f_ln_gamma),
  P4_FXco ("FGAMMA",		p4_f_gamma),
  /* floating-point auxiliary functions */
  P4_FXco ("FPCLASSIFY",	p4_fpclassify),
  P4_OCoN ("FP-INFINITE",	FP_INFINITE),
  P4_OCoN ("FP-NAN",		FP_NAN),
  P4_OCoN ("FP-NORMAL",		FP_NORMAL),
  P4_OCoN ("FP-SUBNORMAL",	FP_SUBNORMAL),
  P4_OCoN ("FP-ZERO",		FP_ZERO),
  P4_OCoN ("FLT-RADIX",		FLT_RADIX),
  P4_OCoN ("FP-ILOGB0",		FP_ILOGB0),
  P4_OCoN ("INT-MAX",		INT_MAX),
  P4_OCoN ("FP-ILOGBNAN",	FP_ILOGBNAN),
  P4_FXco ("ISFINITE",		p4_isfinite),
  P4_FXco ("ISINF",		p4_isinf),
  P4_FXco ("ISNAN",		p4_isnan),
  P4_FXco ("ISNORMAL",		p4_isnormal),
  P4_FXco ("COPYSIGN",		p4_copysign),
  P4_FXco ("NEXTAFTER",		p4_nextafter),
  P4_FXco ("SIGNBIT",		p4_signbit),
  P4_FXco ("FREXP",		p4_frexp),
  P4_FXco ("LDEXP",		p4_ldexp),
  P4_FXco ("MODF",		p4_modf),
  P4_FXco ("MODFD",		p4_modfd),
  P4_FXco ("SCALBN",		p4_scalbn),
  P4_FXco ("LOGB",		p4_logb),
  P4_FXco ("ILOGB",		p4_ilogb),
  /* floating-point environment */
  P4_FXco ("FE-DEFL-ENV",	p4_fe_defl_env),
  P4_FXco ("FEGETENV",		p4_fegetenv),
  P4_FXco ("FESETENV",		p4_fesetenv),
  P4_FXco ("FEHOLDEXCEPT",	p4_feholdexcept),
  P4_FXco ("FEUPDATEENV",	p4_feupdateenv),
  /* floating-point exceptions */
  P4_OCoN ("FE-DIVBYZERO",	FE_DIVBYZERO),
  P4_OCoN ("FE-INEXACT",	FE_INEXACT),
  P4_OCoN ("FE-INVALID",	FE_INVALID),
  P4_OCoN ("FE-OVERFLOW",	FE_OVERFLOW),
  P4_OCoN ("FE-UNDERFLOW",	FE_UNDERFLOW),
  P4_OCoN ("FE-ALL-EXCEPT",	FE_ALL_EXCEPT),
  P4_FXco ("FEGETEXCEPTFLAG",	p4_fegetexceptflag),
  P4_FXco ("FESETEXCEPTFLAG",	p4_fesetexceptflag),
  P4_FXco ("FETESTEXCEPT",	p4_fetestexcept),
  P4_FXco ("FERAISEEXCEPT",	p4_feraiseexcept),
  P4_FXco ("FECLEAREXCEPT",	p4_feclearexcept),
  /* floating-point rounding modes */
  P4_OCoN ("FE-DOWNWARD",	FE_DOWNWARD),
  P4_OCoN ("FE-UPWARD",		FE_UPWARD),
  P4_OCoN ("FE-TONEAREST",	FE_TONEAREST),
  P4_OCoN ("FE-TOWARDZERO",	FE_TOWARDZERO),
  P4_FXco ("FEGETROUND",	p4_fegetround),
  P4_FXco ("FESETROUND",	p4_fesetround),

  P4_INTO ("ENVIRONMENT",	0),
  P4_OCoN ("FPIEEE-EXT",	2005),
};

P4_COUNTWORDS (fpieee, "C99 IEEE 754/854 interface");