comparison src/floatfns.c @ 2094:c77607f8e32d

(FLOAT_CHECK_ERRNO): Define unless NO_FLOAT_CHECK_ERRNO. Changes from Lucid: (HAVE_MATHERR, FLOAT_CHECK_ERRNO, FLOAT_CATCH_SIGILL): New parm macros. (FLOAT_CHECK_DOMAIN, HAVE_RINT): New parm macros. (HAVE_INVERSE_HYPERBOLIC, HAVE_CBRT): New parm macros. [!HAVE_RINT]: Define rint as macro. (IN_FLOAT): Major rewrite; several alternate versions. (IN_FLOAT2): New macro. (arith_error, range_error, domain_error, domain_error2): New macros. (Facos, Fasin, Fatan, Fcos, Fsin, Ftan, Fexp, Fexpt, Flog): Changed. (Flog10, Fsqrt, Fabs, Ffloat, Flogb): Changed. (Ffloor, Fceiling, Fround, Ftruncate): Changed. (Fcube_root): Renamed from Fcbrt. (matherr): New function. (float_error): Only if FLOAT_CATCH_SIGILL.

2093:ce8bad247b1a

You should have received a copy of the GNU General Public License
along with GNU Emacs; see the file COPYING.  If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.  */


#include <signal.h>

#include "config.h"
#include "lisp.h"
#include "syssignal.h"

Lisp_Object Qarith_error;

#ifdef LISP_FLOAT_TYPE

#include <math.h>
#include <errno.h>

extern int errno;

/* Avoid traps on VMS from sinh and cosh.
   All the other functions set errno instead.  */

#ifdef VMS

#undef sinh
#define cosh(x) ((exp(x)+exp(-x))*0.5)
#define sinh(x) ((exp(x)-exp(-x))*0.5)
#endif /* VMS */

static SIGTYPE float_error ();

/* Nonzero while executing in floating point.
   This tells float_error what to do.  */

static int in_float;

/* If an argument is out of range for a mathematical function,
   here is the actual argument value to use in the error message.  */

static Lisp_Object float_error_arg;

/* Evaluate the floating point expression D, recording NUM
   as the original argument for error messages.
   D is normally an assignment expression.
   Handle errors which may result in signals or may set errno.

   Note that float_error may be declared to return void, so you can't
   just cast the zero after the colon to (SIGTYPE) to make the types
   check properly.  */

#define IN_FLOAT(D, NUM) \
(in_float = 1, errno = 0, float_error_arg = NUM, (D),			\
 (errno == ERANGE || errno == EDOM ? (float_error (),0) : 0),		\
 in_float = 0)

/* Extract a Lisp number as a `double', or signal an error.  */

double
extract_float (num)


/* Trig functions.  */

DEFUN ("acos", Facos, Sacos, 1, 1, 0,
  "Return the inverse cosine of ARG.")
  (num)
     register Lisp_Object num;
{
  double d = extract_float (num);
  IN_FLOAT (d = acos (d), num);
  return make_float (d);
}

DEFUN ("asin", Fasin, Sasin, 1, 1, 0,
  "Return the inverse sine of ARG.")
  (num)
     register Lisp_Object num;
{
  double d = extract_float (num);
  IN_FLOAT (d = asin (d), num);
  return make_float (d);
}

DEFUN ("atan", Fatan, Satan, 1, 1, 0,
  "Return the inverse tangent of ARG.")
  (num)
     register Lisp_Object num;
{
  double d = extract_float (num);
  IN_FLOAT (d = atan (d), num);
  return make_float (d);
}

DEFUN ("cos", Fcos, Scos, 1, 1, 0,
  "Return the cosine of ARG.")
  (num)
     register Lisp_Object num;
{
  double d = extract_float (num);
  IN_FLOAT (d = cos (d), num);
  return make_float (d);
}

DEFUN ("sin", Fsin, Ssin, 1, 1, 0,
  "Return the sine of ARG.")
  (num)
     register Lisp_Object num;
{
  double d = extract_float (num);
  IN_FLOAT (d = sin (d), num);
  return make_float (d);
}

DEFUN ("tan", Ftan, Stan, 1, 1, 0,
  "Return the tangent of ARG.")
  (num)
     register Lisp_Object num;
{
  double d = extract_float (num);
  IN_FLOAT (d = tan (d), num);
  return make_float (d);
}

#if 0 /* Leave these out unless we find there's a reason for them.  */

DEFUN ("bessel-j0", Fbessel_j0, Sbessel_j0, 1, 1, 0,
  "Return the bessel function j0 of ARG.")
  (num)
     register Lisp_Object num;
{
  double d = extract_float (num);
  IN_FLOAT (d = j0 (d), num);
  return make_float (d);
}

DEFUN ("bessel-j1", Fbessel_j1, Sbessel_j1, 1, 1, 0,
  "Return the bessel function j1 of ARG.")
  (num)
     register Lisp_Object num;
{
  double d = extract_float (num);
  IN_FLOAT (d = j1 (d), num);
  return make_float (d);
}

DEFUN ("bessel-jn", Fbessel_jn, Sbessel_jn, 2, 2, 0,
  "Return the order N bessel function output jn of ARG.\n\
The first arg (the order) is truncated to an integer.")
  (num1, num2)
     register Lisp_Object num1, num2;
{
  int i1 = extract_float (num1);
  double f2 = extract_float (num2);

  IN_FLOAT (f2 = jn (i1, f2), num1);
  return make_float (f2);
}

DEFUN ("bessel-y0", Fbessel_y0, Sbessel_y0, 1, 1, 0,
  "Return the bessel function y0 of ARG.")
  (num)
     register Lisp_Object num;
{
  double d = extract_float (num);
  IN_FLOAT (d = y0 (d), num);
  return make_float (d);
}

DEFUN ("bessel-y1", Fbessel_y1, Sbessel_y1, 1, 1, 0,
  "Return the bessel function y1 of ARG.")
  (num)
     register Lisp_Object num;
{
  double d = extract_float (num);
  IN_FLOAT (d = y1 (d), num);
  return make_float (d);
}

DEFUN ("bessel-yn", Fbessel_yn, Sbessel_yn, 2, 2, 0,
  "Return the order N bessel function output yn of ARG.\n\
The first arg (the order) is truncated to an integer.")
  (num1, num2)
     register Lisp_Object num1, num2;
{
  int i1 = extract_float (num1);
  double f2 = extract_float (num2);

  IN_FLOAT (f2 = yn (i1, f2), num1);
  return make_float (f2);
}

#endif

#if 0 /* Leave these out unless we see they are worth having.  */

DEFUN ("erf", Ferf, Serf, 1, 1, 0,
  "Return the mathematical error function of ARG.")
  (num)
     register Lisp_Object num;
{
  double d = extract_float (num);
  IN_FLOAT (d = erf (d), num);
  return make_float (d);
}

DEFUN ("erfc", Ferfc, Serfc, 1, 1, 0,
  "Return the complementary error function of ARG.")
  (num)
     register Lisp_Object num;
{
  double d = extract_float (num);
  IN_FLOAT (d = erfc (d), num);
  return make_float (d);
}

DEFUN ("log-gamma", Flog_gamma, Slog_gamma, 1, 1, 0,
  "Return the log gamma of ARG.")
  (num)
     register Lisp_Object num;
{
  double d = extract_float (num);
  IN_FLOAT (d = lgamma (d), num);
  return make_float (d);
}

DEFUN ("cbrt", Fcbrt, Scbrt, 1, 1, 0,
  "Return the cube root of ARG.")
  (num)
     register Lisp_Object num;
{
  double d = extract_float (num);
  IN_FLOAT (d = cbrt (d), num);
  return make_float (d);
}

#endif

DEFUN ("exp", Fexp, Sexp, 1, 1, 0,
  "Return the exponential base e of ARG.")
  (num)
     register Lisp_Object num;
{
  double d = extract_float (num);
  IN_FLOAT (d = exp (d), num);
  return make_float (d);
}

DEFUN ("expt", Fexpt, Sexpt, 2, 2, 0,
  "Return the exponential X ** Y.")
  (num1, num2)
     register Lisp_Object num1, num2;
{
  double f1, f2;

  CHECK_NUMBER_OR_FLOAT (num1, 0);
  CHECK_NUMBER_OR_FLOAT (num2, 0);
  if ((XTYPE (num1) == Lisp_Int) && /* common lisp spec */
      (XTYPE (num2) == Lisp_Int)) /* don't promote, if both are ints */
    {				/* this can be improved by pre-calculating */
      int acc, x, y;		/* some binary powers of x then acumulating */
      /* these, therby saving some time. -wsr */
      x = XINT (num1);
      y = XINT (num2);
      acc = 1;
      
      if (y < 0)
	{
	  for (; y < 0; y++)
	    acc /= x;
	}
      else
	{
	  for (; y > 0; y--)
	    acc *= x;
	}
      XFASTINT (x) = acc;
      return x;
    }
  f1 = (XTYPE (num1) == Lisp_Float) ? XFLOAT (num1)->data : XINT (num1);
  f2 = (XTYPE (num2) == Lisp_Float) ? XFLOAT (num2)->data : XINT (num2);
  IN_FLOAT (f1 = pow (f1, f2), num1);
  return make_float (f1);
}

DEFUN ("log", Flog, Slog, 1, 2, 0,
  "Return the natural logarithm of NUM.\n\
If second optional argument BASE is given, return log NUM using that base.")
  (num, base)
     register Lisp_Object num, base;
{
  double d = extract_float (num);

  if (NILP (base))
    IN_FLOAT (d = log (d), num);
  else
    {
      double b = extract_float (base);

      IN_FLOAT (d = log (num) / log (b), num);
    }
  return make_float (d);
}

DEFUN ("log10", Flog10, Slog10, 1, 1, 0,
  "Return the logarithm base 10 of ARG.")
  (num)
     register Lisp_Object num;
{
  double d = extract_float (num);
  IN_FLOAT (d = log10 (d), num);
  return make_float (d);
}

DEFUN ("sqrt", Fsqrt, Ssqrt, 1, 1, 0,
  "Return the square root of ARG.")
  (num)
     register Lisp_Object num;
{
  double d = extract_float (num);
  IN_FLOAT (d = sqrt (d), num);
  return make_float (d);
}

#if 0 /* Not clearly worth adding.  */

DEFUN ("acosh", Facosh, Sacosh, 1, 1, 0,
  "Return the inverse hyperbolic cosine of ARG.")
  (num)
     register Lisp_Object num;
{
  double d = extract_float (num);
  IN_FLOAT (d = acosh (d), num);
  return make_float (d);
}

DEFUN ("asinh", Fasinh, Sasinh, 1, 1, 0,
  "Return the inverse hyperbolic sine of ARG.")
  (num)
     register Lisp_Object num;
{
  double d = extract_float (num);
  IN_FLOAT (d = asinh (d), num);
  return make_float (d);
}

DEFUN ("atanh", Fatanh, Satanh, 1, 1, 0,
  "Return the inverse hyperbolic tangent of ARG.")
  (num)
     register Lisp_Object num;
{
  double d = extract_float (num);
  IN_FLOAT (d = atanh (d), num);
  return make_float (d);
}

DEFUN ("cosh", Fcosh, Scosh, 1, 1, 0,
  "Return the hyperbolic cosine of ARG.")
  (num)
     register Lisp_Object num;
{
  double d = extract_float (num);
  IN_FLOAT (d = cosh (d), num);
  return make_float (d);
}

DEFUN ("sinh", Fsinh, Ssinh, 1, 1, 0,
  "Return the hyperbolic sine of ARG.")
  (num)
     register Lisp_Object num;
{
  double d = extract_float (num);
  IN_FLOAT (d = sinh (d), num);
  return make_float (d);
}

DEFUN ("tanh", Ftanh, Stanh, 1, 1, 0,
  "Return the hyperbolic tangent of ARG.")
  (num)
     register Lisp_Object num;
{
  double d = extract_float (num);
  IN_FLOAT (d = tanh (d), num);
  return make_float (d);
}
#endif

DEFUN ("abs", Fabs, Sabs, 1, 1, 0,
  "Return the absolute value of ARG.")
  (num)
     register Lisp_Object num;
{
  CHECK_NUMBER_OR_FLOAT (num, 0);

  if (XTYPE (num) == Lisp_Float)
    IN_FLOAT (num = make_float (fabs (XFLOAT (num)->data)), num);
  else if (XINT (num) < 0)
    XSETINT (num, - XFASTINT (num));

  return num;
}

DEFUN ("float", Ffloat, Sfloat, 1, 1, 0,
  "Return the floating point number equal to ARG.")
  (num)
     register Lisp_Object num;
{
  CHECK_NUMBER_OR_FLOAT (num, 0);

  if (XTYPE (num) == Lisp_Int)
    return make_float ((double) XINT (num));
  else				/* give 'em the same float back */
    return num;
}

DEFUN ("logb", Flogb, Slogb, 1, 1, 0,
  "Returns the integer that is the base 2 log of ARG.\n\
This is the same as the exponent of a float.")
     (num)
Lisp_Object num;
{
  /* System V apparently doesn't have a `logb' function.  It might be
     better to use it on systems that have it, but Ultrix (at least)
     doesn't declare it properly in <math.h>; does anyone really care? */
  return Flog (num, make_number (2));
}

/* the rounding functions  */

DEFUN ("ceiling", Fceiling, Sceiling, 1, 1, 0,
  "Return the smallest integer no less than ARG.  (Round toward +inf.)")
  (num)
     register Lisp_Object num;
{
  CHECK_NUMBER_OR_FLOAT (num, 0);

  if (XTYPE (num) == Lisp_Float)
    IN_FLOAT (XSET (num, Lisp_Int, ceil (XFLOAT (num)->data)), num);

  return num;
}

DEFUN ("floor", Ffloor, Sfloor, 1, 1, 0,
  "Return the largest integer no greater than ARG.  (Round towards -inf.)")
  (num)
     register Lisp_Object num;
{
  CHECK_NUMBER_OR_FLOAT (num, 0);

  if (XTYPE (num) == Lisp_Float)
    IN_FLOAT (XSET (num, Lisp_Int, floor (XFLOAT (num)->data)), num);

  return num;
}

DEFUN ("round", Fround, Sround, 1, 1, 0,
  "Return the nearest integer to ARG.")
  (num)
     register Lisp_Object num;
{
  CHECK_NUMBER_OR_FLOAT (num, 0);

  if (XTYPE (num) == Lisp_Float)
    {
      /* Screw the prevailing rounding mode.  */
      IN_FLOAT (XSET (num, Lisp_Int, floor (XFLOAT (num)->data + 0.5)), num);

      /* It used to be that on non-USG systems we'd use the `rint'
	 function.  But that seems not to be declared properly in
	 <math.h> on Ultrix, I don't want to declare it myself because
	 that might conflict with <math.h> on other systems, and I
	 don't see what's wrong with the code above anyway.  */
    }

  return num;
}

DEFUN ("truncate", Ftruncate, Struncate, 1, 1, 0,
       "Truncate a floating point number to an int.\n\
Rounds the value toward zero.")
  (num)
     register Lisp_Object num;
{
  CHECK_NUMBER_OR_FLOAT (num, 0);

  if (XTYPE (num) == Lisp_Float)
    XSET (num, Lisp_Int, (int) XFLOAT (num)->data);

  return num;
}

static SIGTYPE
float_error (signo)
     int signo;
{
  if (! in_float)

  in_float = 0;

  Fsignal (Qarith_error, Fcons (float_error_arg, Qnil));
}

init_floatfns ()
{
  signal (SIGILL, float_error);
  in_float = 0;
}

syms_of_floatfns ()
{

  defsubr (&Sbessel_j1);
  defsubr (&Sbessel_jn);
  defsubr (&Serf);
  defsubr (&Serfc);
  defsubr (&Slog_gamma);
  defsubr (&Scbrt);
#endif
  defsubr (&Sexp);
  defsubr (&Sexpt);
  defsubr (&Slog);
  defsubr (&Slog10);

2094:c77607f8e32d

You should have received a copy of the GNU General Public License
along with GNU Emacs; see the file COPYING.  If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.  */


/* ANSI C requires only these float functions:
   acos, asin, atan, atan2, ceil, cos, cosh, exp, fabs, floor, fmod,
   frexp, ldexp, log, log10, modf, pow, sin, sinh, sqrt, tan, tanh.

   Define HAVE_INVERSE_HYPERBOLIC if you have acosh, asinh, and atanh.
   Define HAVE_CBRT if you have cbrt.
   Define HAVE_RINT if you have rint.
   If you don't define these, then the appropriate routines will be simulated.

   Define HAVE_MATHERR if on a system supporting the SysV matherr callback.
   (This should happen automatically.)

   Define FLOAT_CHECK_ERRNO if the float library routines set errno.
   This has no effect if HAVE_MATHERR is defined.

   Define FLOAT_CATCH_SIGILL if the float library routines signal SIGILL.
   (What systems actually do this?  Please let us know.)

   Define FLOAT_CHECK_DOMAIN if the float library doesn't handle errors by
   either setting errno, or signalling SIGFPE/SIGILL.  Otherwise, domain and
   range checking will happen before calling the float routines.  This has
   no effect if HAVE_MATHERR is defined (since matherr will be called when
   a domain error occurs.)
 */

#include <signal.h>

#include "config.h"
#include "lisp.h"
#include "syssignal.h"

Lisp_Object Qarith_error;

#ifdef LISP_FLOAT_TYPE

#include <math.h>

#if defined(DOMAIN) && defined(SING) && defined(OVERFLOW)
    /* If those are defined, then this is probably a `matherr' machine. */
# ifndef HAVE_MATHERR
#  define HAVE_MATHERR
# endif
#endif

#ifdef HAVE_MATHERR
# ifdef FLOAT_CHECK_ERRNO
#  undef FLOAT_CHECK_ERRNO
# endif
# ifdef FLOAT_CHECK_DOMAIN
#  undef FLOAT_CHECK_DOMAIN
# endif
#endif

#ifndef NO_FLOAT_CHECK_ERRNO
#define FLOAT_CHECK_ERRNO
#endif

#ifdef FLOAT_CHECK_ERRNO
# include <errno.h>

extern int errno;
#endif

/* Avoid traps on VMS from sinh and cosh.
   All the other functions set errno instead.  */

#ifdef VMS

#undef sinh
#define cosh(x) ((exp(x)+exp(-x))*0.5)
#define sinh(x) ((exp(x)-exp(-x))*0.5)
#endif /* VMS */

#ifndef HAVE_RINT
#define rint(x) (floor((x)+0.5))
#endif

static SIGTYPE float_error ();

/* Nonzero while executing in floating point.
   This tells float_error what to do.  */

static int in_float;

/* If an argument is out of range for a mathematical function,
   here is the actual argument value to use in the error message.  */

static Lisp_Object float_error_arg, float_error_arg2;

static char *float_error_fn_name;

/* Evaluate the floating point expression D, recording NUM
   as the original argument for error messages.
   D is normally an assignment expression.
   Handle errors which may result in signals or may set errno.

   Note that float_error may be declared to return void, so you can't
   just cast the zero after the colon to (SIGTYPE) to make the types
   check properly.  */

#ifdef FLOAT_CHECK_ERRNO
#define IN_FLOAT(d, name, num)				\
  do {							\
    float_error_arg = num;				\
    float_error_fn_name = name;				\
    in_float = 1; errno = 0; (d); in_float = 0;		\
    switch (errno) {					\
    case 0: break;					\
    case EDOM:	 domain_error (float_error_fn_name, float_error_arg);	\
    case ERANGE: range_error (float_error_fn_name, float_error_arg);	\
    default:	 arith_error (float_error_fn_name, float_error_arg);	\
    }							\
  } while (0)
#define IN_FLOAT2(d, name, num, num2)			\
  do {							\
    float_error_arg = num;				\
    float_error_arg2 = num2;				\
    float_error_fn_name = name;				\
    in_float = 1; errno = 0; (d); in_float = 0;		\
    switch (errno) {					\
    case 0: break;					\
    case EDOM:	 domain_error (float_error_fn_name, float_error_arg);	\
    case ERANGE: range_error (float_error_fn_name, float_error_arg);	\
    default:	 arith_error (float_error_fn_name, float_error_arg);	\
    }							\
  } while (0)
#else
#define IN_FLOAT2(d, name, num, num2) (in_float = 1, (d), in_float = 0)
#endif

#define arith_error(op,arg) \
  Fsignal (Qarith_error, Fcons (build_string ((op)), Fcons ((arg), Qnil)))
#define range_error(op,arg) \
  Fsignal (Qrange_error, Fcons (build_string ((op)), Fcons ((arg), Qnil)))
#define domain_error(op,arg) \
  Fsignal (Qdomain_error, Fcons (build_string ((op)), Fcons ((arg), Qnil)))
#define domain_error2(op,a1,a2) \
  Fsignal (Qdomain_error, Fcons (build_string ((op)), Fcons ((a1), Fcons ((a2), Qnil))))

/* Extract a Lisp number as a `double', or signal an error.  */

double
extract_float (num)


/* Trig functions.  */

DEFUN ("acos", Facos, Sacos, 1, 1, 0,
  "Return the inverse cosine of ARG.")
  (arg)
     register Lisp_Object arg;
{
  double d = extract_float (arg);
#ifdef FLOAT_CHECK_DOMAIN
  if (d > 1.0 || d < -1.0)
    domain_error ("acos", arg);
#endif
  IN_FLOAT (d = acos (d), "acos", arg);
  return make_float (d);
}

DEFUN ("asin", Fasin, Sasin, 1, 1, 0,
  "Return the inverse sine of ARG.")
  (arg)
     register Lisp_Object arg;
{
  double d = extract_float (arg);
#ifdef FLOAT_CHECK_DOMAIN
  if (d > 1.0 || d < -1.0)
    domain_error ("asin", arg);
#endif
  IN_FLOAT (d = asin (d), "asin", arg);
  return make_float (d);
}

DEFUN ("atan", Fatan, Satan, 1, 1, 0,
  "Return the inverse tangent of ARG.")
  (arg)
     register Lisp_Object arg;
{
  double d = extract_float (arg);
  IN_FLOAT (d = atan (d), "atan", arg);
  return make_float (d);
}

DEFUN ("cos", Fcos, Scos, 1, 1, 0,
  "Return the cosine of ARG.")
  (arg)
     register Lisp_Object arg;
{
  double d = extract_float (arg);
  IN_FLOAT (d = cos (d), "cos", arg);
  return make_float (d);
}

DEFUN ("sin", Fsin, Ssin, 1, 1, 0,
  "Return the sine of ARG.")
  (arg)
     register Lisp_Object arg;
{
  double d = extract_float (arg);
  IN_FLOAT (d = sin (d), "sin", arg);
  return make_float (d);
}

DEFUN ("tan", Ftan, Stan, 1, 1, 0,
  "Return the tangent of ARG.")
  (arg)
     register Lisp_Object arg;
{
  double d = extract_float (arg);
  double c = cos (d);
#ifdef FLOAT_CHECK_DOMAIN
  if (c == 0.0)
    domain_error ("tan", arg);
#endif
  IN_FLOAT (d = sin (d) / c, "tan", arg);
  return make_float (d);
}

#if 0 /* Leave these out unless we find there's a reason for them.  */

DEFUN ("bessel-j0", Fbessel_j0, Sbessel_j0, 1, 1, 0,
  "Return the bessel function j0 of ARG.")
  (arg)
     register Lisp_Object arg;
{
  double d = extract_float (arg);
  IN_FLOAT (d = j0 (d), "bessel-j0", arg);
  return make_float (d);
}

DEFUN ("bessel-j1", Fbessel_j1, Sbessel_j1, 1, 1, 0,
  "Return the bessel function j1 of ARG.")
  (arg)
     register Lisp_Object arg;
{
  double d = extract_float (arg);
  IN_FLOAT (d = j1 (d), "bessel-j1", arg);
  return make_float (d);
}

DEFUN ("bessel-jn", Fbessel_jn, Sbessel_jn, 2, 2, 0,
  "Return the order N bessel function output jn of ARG.\n\
The first arg (the order) is truncated to an integer.")
  (arg1, arg2)
     register Lisp_Object arg1, arg2;
{
  int i1 = extract_float (arg1);
  double f2 = extract_float (arg2);

  IN_FLOAT (f2 = jn (i1, f2), "bessel-jn", arg1);
  return make_float (f2);
}

DEFUN ("bessel-y0", Fbessel_y0, Sbessel_y0, 1, 1, 0,
  "Return the bessel function y0 of ARG.")
  (arg)
     register Lisp_Object arg;
{
  double d = extract_float (arg);
  IN_FLOAT (d = y0 (d), "bessel-y0", arg);
  return make_float (d);
}

DEFUN ("bessel-y1", Fbessel_y1, Sbessel_y1, 1, 1, 0,
  "Return the bessel function y1 of ARG.")
  (arg)
     register Lisp_Object arg;
{
  double d = extract_float (arg);
  IN_FLOAT (d = y1 (d), "bessel-y0", arg);
  return make_float (d);
}

DEFUN ("bessel-yn", Fbessel_yn, Sbessel_yn, 2, 2, 0,
  "Return the order N bessel function output yn of ARG.\n\
The first arg (the order) is truncated to an integer.")
  (arg1, arg2)
     register Lisp_Object arg1, arg2;
{
  int i1 = extract_float (arg1);
  double f2 = extract_float (arg2);

  IN_FLOAT (f2 = yn (i1, f2), "bessel-yn", arg1);
  return make_float (f2);
}

#endif

#if 0 /* Leave these out unless we see they are worth having.  */

DEFUN ("erf", Ferf, Serf, 1, 1, 0,
  "Return the mathematical error function of ARG.")
  (arg)
     register Lisp_Object arg;
{
  double d = extract_float (arg);
  IN_FLOAT (d = erf (d), "erf", arg);
  return make_float (d);
}

DEFUN ("erfc", Ferfc, Serfc, 1, 1, 0,
  "Return the complementary error function of ARG.")
  (arg)
     register Lisp_Object arg;
{
  double d = extract_float (arg);
  IN_FLOAT (d = erfc (d), "erfc", arg);
  return make_float (d);
}

DEFUN ("log-gamma", Flog_gamma, Slog_gamma, 1, 1, 0,
  "Return the log gamma of ARG.")
  (arg)
     register Lisp_Object arg;
{
  double d = extract_float (arg);
  IN_FLOAT (d = lgamma (d), "log-gamma", arg);
  return make_float (d);
}

DEFUN ("cube-root", Fcube_root, Scube_root, 1, 1, 0,
  "Return the cube root of ARG.")
  (arg)
     register Lisp_Object arg;
{
  double d = extract_float (arg);
#ifdef HAVE_CBRT
  IN_FLOAT (d = cbrt (d), "cube-root", arg);
#else
  if (d >= 0.0)
    IN_FLOAT (d = pow (d, 1.0/3.0), "cube-root", arg);
  else
    IN_FLOAT (d = -pow (-d, 1.0/3.0), "cube-root", arg);
#endif
  return make_float (d);
}

#endif

DEFUN ("exp", Fexp, Sexp, 1, 1, 0,
  "Return the exponential base e of ARG.")
  (arg)
     register Lisp_Object arg;
{
  double d = extract_float (arg);
#ifdef FLOAT_CHECK_DOMAIN
  if (d > 709.7827)   /* Assume IEEE doubles here */
    range_error ("exp", arg);
  else if (d < -709.0)
    return make_float (0.0);
  else
#endif
    IN_FLOAT (d = exp (d), "exp", arg);
  return make_float (d);
}

DEFUN ("expt", Fexpt, Sexpt, 2, 2, 0,
  "Return the exponential X ** Y.")
  (arg1, arg2)
     register Lisp_Object arg1, arg2;
{
  double f1, f2;

  CHECK_NUMBER_OR_FLOAT (arg1, 0);
  CHECK_NUMBER_OR_FLOAT (arg2, 0);
  if ((XTYPE (arg1) == Lisp_Int) && /* common lisp spec */
      (XTYPE (arg2) == Lisp_Int)) /* don't promote, if both are ints */
    {				/* this can be improved by pre-calculating */
      int acc, x, y;		/* some binary powers of x then acumulating */
      /* these, therby saving some time. -wsr */
      x = XINT (arg1);
      y = XINT (arg2);
      acc = 1;
      
      if (y < 0)
	{
	  if (x == 1)
	    acc = 1;
	  else if (x == -1)
	    acc = (y & 1) ? -1 : 1;
	  else
	    acc = 0;
	}
      else
	{
	  for (; y > 0; y--)
	  while (y > 0)
	    {
	      if (y & 1)
		acc *= x;
	      x *= x;
	      y = (unsigned)y >> 1;
	    }
	}
      XSET (x, Lisp_Int, acc);
      return x;
    }
  f1 = (XTYPE (arg1) == Lisp_Float) ? XFLOAT (arg1)->data : XINT (arg1);
  f2 = (XTYPE (arg2) == Lisp_Float) ? XFLOAT (arg2)->data : XINT (arg2);
  /* Really should check for overflow, too */
  if (f1 == 0.0 && f2 == 0.0)
    f1 = 1.0;
#ifdef FLOAT_CHECK_DOMAIN
  else if ((f1 == 0.0 && f2 < 0.0) || (f1 < 0 && f2 != floor(f2)))
    domain_error2 ("expt", arg1, arg2);
#endif
  IN_FLOAT (f1 = pow (f1, f2), "expt", arg1);
  return make_float (f1);
}

DEFUN ("log", Flog, Slog, 1, 2, 0,
  "Return the natural logarithm of ARG.\n\
If second optional argument BASE is given, return log ARG using that base.")
  (arg, base)
     register Lisp_Object arg, base;
{
  double d = extract_float (arg);

#ifdef FLOAT_CHECK_DOMAIN
  if (d <= 0.0)
    domain_error2 ("log", arg, base);
#endif
  if (NILP (base))
    IN_FLOAT (d = log (d), "log", arg);
  else
    {
      double b = extract_float (base);

#ifdef FLOAT_CHECK_DOMAIN
      if (b <= 0.0 || b == 1.0)
	domain_error2 ("log", arg, base);
#endif
      if (b == 10.0)
	IN_FLOAT2 (d = log10 (d), "log", arg, base);
      else
	IN_FLOAT2 (d = log (arg) / log (b), "log", arg, base);
    }
  return make_float (d);
}

DEFUN ("log10", Flog10, Slog10, 1, 1, 0,
  "Return the logarithm base 10 of ARG.")
  (arg)
     register Lisp_Object arg;
{
  double d = extract_float (arg);
#ifdef FLOAT_CHECK_DOMAIN
  if (d <= 0.0)
    domain_error ("log10", arg);
#endif
  IN_FLOAT (d = log10 (d), "log10", arg);
  return make_float (d);
}

DEFUN ("sqrt", Fsqrt, Ssqrt, 1, 1, 0,
  "Return the square root of ARG.")
  (arg)
     register Lisp_Object arg;
{
  double d = extract_float (arg);
#ifdef FLOAT_CHECK_DOMAIN
  if (d < 0.0)
    domain_error ("sqrt", arg);
#endif
  IN_FLOAT (d = sqrt (d), "sqrt", arg);
  return make_float (d);
}

#if 0 /* Not clearly worth adding.  */

DEFUN ("acosh", Facosh, Sacosh, 1, 1, 0,
  "Return the inverse hyperbolic cosine of ARG.")
  (arg)
     register Lisp_Object arg;
{
  double d = extract_float (arg);
#ifdef FLOAT_CHECK_DOMAIN
  if (d < 1.0)
    domain_error ("acosh", arg);
#endif
#ifdef HAVE_INVERSE_HYPERBOLIC
  IN_FLOAT (d = acosh (d), "acosh", arg);
#else
  IN_FLOAT (d = log (d + sqrt (d*d - 1.0)), "acosh", arg);
#endif
  return make_float (d);
}

DEFUN ("asinh", Fasinh, Sasinh, 1, 1, 0,
  "Return the inverse hyperbolic sine of ARG.")
  (arg)
     register Lisp_Object arg;
{
  double d = extract_float (arg);
#ifdef HAVE_INVERSE_HYPERBOLIC
  IN_FLOAT (d = asinh (d), "asinh", arg);
#else
  IN_FLOAT (d = log (d + sqrt (d*d + 1.0)), "asinh", arg);
#endif
  return make_float (d);
}

DEFUN ("atanh", Fatanh, Satanh, 1, 1, 0,
  "Return the inverse hyperbolic tangent of ARG.")
  (arg)
     register Lisp_Object arg;
{
  double d = extract_float (arg);
#ifdef FLOAT_CHECK_DOMAIN
  if (d >= 1.0 || d <= -1.0)
    domain_error ("atanh", arg);
#endif
#ifdef HAVE_INVERSE_HYPERBOLIC
  IN_FLOAT (d = atanh (d), "atanh", arg);
#else
  IN_FLOAT (d = 0.5 * log ((1.0 + d) / (1.0 - d)), "atanh", arg);
#endif
  return make_float (d);
}

DEFUN ("cosh", Fcosh, Scosh, 1, 1, 0,
  "Return the hyperbolic cosine of ARG.")
  (arg)
     register Lisp_Object arg;
{
  double d = extract_float (arg);
#ifdef FLOAT_CHECK_DOMAIN
  if (d > 710.0 || d < -710.0)
    range_error ("cosh", arg);
#endif
  IN_FLOAT (d = cosh (d), "cosh", arg);
  return make_float (d);
}

DEFUN ("sinh", Fsinh, Ssinh, 1, 1, 0,
  "Return the hyperbolic sine of ARG.")
  (arg)
     register Lisp_Object arg;
{
  double d = extract_float (arg);
#ifdef FLOAT_CHECK_DOMAIN
  if (d > 710.0 || d < -710.0)
    range_error ("sinh", arg);
#endif
  IN_FLOAT (d = sinh (d), "sinh", arg);
  return make_float (d);
}

DEFUN ("tanh", Ftanh, Stanh, 1, 1, 0,
  "Return the hyperbolic tangent of ARG.")
  (arg)
     register Lisp_Object arg;
{
  double d = extract_float (arg);
  IN_FLOAT (d = tanh (d), "tanh", arg);
  return make_float (d);
}
#endif

DEFUN ("abs", Fabs, Sabs, 1, 1, 0,
  "Return the absolute value of ARG.")
  (arg)
     register Lisp_Object arg;
{
  CHECK_NUMBER_OR_FLOAT (arg, 0);

  if (XTYPE (arg) == Lisp_Float)
    IN_FLOAT (arg = make_float (fabs (XFLOAT (arg)->data)), "abs", arg);
  else if (XINT (arg) < 0)
    XSETINT (arg, - XFASTINT (arg));

  return arg;
}

DEFUN ("float", Ffloat, Sfloat, 1, 1, 0,
  "Return the floating point number equal to ARG.")
  (arg)
     register Lisp_Object arg;
{
  CHECK_NUMBER_OR_FLOAT (arg, 0);

  if (XTYPE (arg) == Lisp_Int)
    return make_float ((double) XINT (arg));
  else				/* give 'em the same float back */
    return arg;
}

DEFUN ("logb", Flogb, Slogb, 1, 1, 0,
  "Returns the integer that is the base 2 log of ARG.\n\
This is the same as the exponent of a float.")
     (arg)
     Lisp_Object arg;
{
  /* System V apparently doesn't have a `logb' function.  It might be
     better to use it on systems that have it, but Ultrix (at least)
     doesn't declare it properly in <math.h>; does anyone really care? */
  return Flog (arg, make_number (2));
}

/* the rounding functions  */

DEFUN ("ceiling", Fceiling, Sceiling, 1, 1, 0,
  "Return the smallest integer no less than ARG.  (Round toward +inf.)")
  (arg)
     register Lisp_Object arg;
{
  CHECK_NUMBER_OR_FLOAT (arg, 0);

  if (XTYPE (arg) == Lisp_Float)
    IN_FLOAT (XSET (arg, Lisp_Int, ceil (XFLOAT (arg)->data)), "celing", arg);

  return arg;
}

DEFUN ("floor", Ffloor, Sfloor, 1, 1, 0,
  "Return the largest integer no greater than ARG.  (Round towards -inf.)")
  (arg)
     register Lisp_Object arg;
{
  CHECK_NUMBER_OR_FLOAT (arg, 0);

  if (XTYPE (arg) == Lisp_Float)
    IN_FLOAT (XSET (arg, Lisp_Int, floor (XFLOAT (arg)->data)), "floor", arg);

  return arg;
}

DEFUN ("round", Fround, Sround, 1, 1, 0,
  "Return the nearest integer to ARG.")
  (arg)
     register Lisp_Object arg;
{
  CHECK_NUMBER_OR_FLOAT (arg, 0);

  if (XTYPE (arg) == Lisp_Float)
    /* Screw the prevailing rounding mode.  */
    IN_FLOAT (XSET (arg, Lisp_Int, rint (XFLOAT (arg)->data)), "round", arg);

  return arg;
}

DEFUN ("truncate", Ftruncate, Struncate, 1, 1, 0,
       "Truncate a floating point number to an int.\n\
Rounds the value toward zero.")
  (arg)
     register Lisp_Object arg;
{
  CHECK_NUMBER_OR_FLOAT (arg, 0);

  if (XTYPE (arg) == Lisp_Float)
    XSET (arg, Lisp_Int, (int) XFLOAT (arg)->data);

  return arg;
}

#if 0
/* It's not clear these are worth adding.  */

DEFUN ("fceiling", Ffceiling, Sfceiling, 1, 1, 0,
  "Return the smallest integer no less than ARG, as a float.\n\
\(Round toward +inf.\)")
  (arg)
     register Lisp_Object arg;
{
  double d = extract_float (arg);
  IN_FLOAT (d = ceil (d), "fceiling", arg);
  return make_float (d);
}

DEFUN ("ffloor", Fffloor, Sffloor, 1, 1, 0,
  "Return the largest integer no greater than ARG, as a float.\n\
\(Round towards -inf.\)")
  (arg)
     register Lisp_Object arg;
{
  double d = extract_float (arg);
  IN_FLOAT (d = floor (d), "ffloor", arg);
  return make_float (d);
}

DEFUN ("fround", Ffround, Sfround, 1, 1, 0,
  "Return the nearest integer to ARG, as a float.")
  (arg)
     register Lisp_Object arg;
{
  double d = extract_float (arg);
  IN_FLOAT (d = rint (XFLOAT (arg)->data), "fround", arg);
  return make_float (d);
}

DEFUN ("ftruncate", Fftruncate, Sftruncate, 1, 1, 0,
       "Truncate a floating point number to an integral float value.\n\
Rounds the value toward zero.")
  (arg)
     register Lisp_Object arg;
{
  double d = extract_float (arg);
  if (d >= 0.0)
    IN_FLOAT (d = floor (d), "ftruncate", arg);
  else
    IN_FLOAT (d = ceil (d), arg);
  return make_float (d);
}
#endif

#ifdef FLOAT_CATCH_SIGILL
static SIGTYPE
float_error (signo)
     int signo;
{
  if (! in_float)

  in_float = 0;

  Fsignal (Qarith_error, Fcons (float_error_arg, Qnil));
}

/* Another idea was to replace the library function `infnan'
   where SIGILL is signaled.  */

#endif /* FLOAT_CATCH_SIGILL */

#ifdef HAVE_MATHERR
int 
matherr (x)
     struct exception *x;
{
  Lisp_Object args;
  if (! in_float)
    /* Not called from emacs-lisp float routines; do the default thing. */
    return 0;
  if (!strcmp (x->name, "pow"))
    x->name = "expt";

  args
    = Fcons (build_string (x->name),
	     Fcons (make_float (x->arg1),
		    ((!strcmp (x->name, "log") || !strcmp (x->name, "pow"))
		     ? Fcons (make_float (x->arg2), Qnil)
		     : Qnil)));
  switch (x->type)
    {
    case DOMAIN:	Fsignal (Qdomain_error, args);		break;
    case SING:		Fsignal (Qsingularity_error, args);	break;
    case OVERFLOW:	Fsignal (Qoverflow_error, args);	break;
    case UNDERFLOW:	Fsignal (Qunderflow_error, args);	break;
    default:		Fsignal (Qarith_error, args);		break;
    }
  return (1);	/* don't set errno or print a message */
}
#endif /* HAVE_MATHERR */

init_floatfns ()
{
#ifdef FLOAT_CATCH_SIGILL
  signal (SIGILL, float_error);
#endif 
  in_float = 0;
}

syms_of_floatfns ()
{

  defsubr (&Sbessel_j1);
  defsubr (&Sbessel_jn);
  defsubr (&Serf);
  defsubr (&Serfc);
  defsubr (&Slog_gamma);
  defsubr (&Scube_root);
  defsubr (&Sfceiling);
  defsubr (&Sffloor);
  defsubr (&Sfround);
  defsubr (&Sftruncate);
#endif
  defsubr (&Sexp);
  defsubr (&Sexpt);
  defsubr (&Slog);
  defsubr (&Slog10);