home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 5 Edit / 05-Edit.zip / e20313sr.zip / emacs / 20.3.1 / src / floatfns.c < prev    next >
C/C++ Source or Header  |  1999-07-31  |  28KB  |  1,076 lines
  1. /* Primitive operations on floating point for GNU Emacs Lisp interpreter.
  2.    Copyright (C) 1988, 1993, 1994 Free Software Foundation, Inc.
  3.  
  4. This file is part of GNU Emacs.
  5.  
  6. GNU Emacs is free software; you can redistribute it and/or modify
  7. it under the terms of the GNU General Public License as published by
  8. the Free Software Foundation; either version 2, or (at your option)
  9. any later version.
  10.  
  11. GNU Emacs is distributed in the hope that it will be useful,
  12. but WITHOUT ANY WARRANTY; without even the implied warranty of
  13. MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14. GNU General Public License for more details.
  15.  
  16. You should have received a copy of the GNU General Public License
  17. along with GNU Emacs; see the file COPYING.  If not, write to
  18. the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  19. Boston, MA 02111-1307, USA.  */
  20.  
  21.  
  22. /* ANSI C requires only these float functions:
  23.    acos, asin, atan, atan2, ceil, cos, cosh, exp, fabs, floor, fmod,
  24.    frexp, ldexp, log, log10, modf, pow, sin, sinh, sqrt, tan, tanh.
  25.  
  26.    Define HAVE_INVERSE_HYPERBOLIC if you have acosh, asinh, and atanh.
  27.    Define HAVE_CBRT if you have cbrt.
  28.    Define HAVE_RINT if you have a working rint.
  29.    If you don't define these, then the appropriate routines will be simulated.
  30.  
  31.    Define HAVE_MATHERR if on a system supporting the SysV matherr callback.
  32.    (This should happen automatically.)
  33.  
  34.    Define FLOAT_CHECK_ERRNO if the float library routines set errno.
  35.    This has no effect if HAVE_MATHERR is defined.
  36.  
  37.    Define FLOAT_CATCH_SIGILL if the float library routines signal SIGILL.
  38.    (What systems actually do this?  Please let us know.)
  39.  
  40.    Define FLOAT_CHECK_DOMAIN if the float library doesn't handle errors by
  41.    either setting errno, or signaling SIGFPE/SIGILL.  Otherwise, domain and
  42.    range checking will happen before calling the float routines.  This has
  43.    no effect if HAVE_MATHERR is defined (since matherr will be called when
  44.    a domain error occurs.)
  45.  */
  46.  
  47. #include <signal.h>
  48.  
  49. #include <config.h>
  50. #include "lisp.h"
  51. #include "syssignal.h"
  52.  
  53. #ifdef LISP_FLOAT_TYPE
  54.  
  55. #if STDC_HEADERS
  56. #include <float.h>
  57. #endif
  58.  
  59. /* If IEEE_FLOATING_POINT isn't defined, default it from FLT_*. */
  60. #ifndef IEEE_FLOATING_POINT
  61. #if (FLT_RADIX == 2 && FLT_MANT_DIG == 24 \
  62.      && FLT_MIN_EXP == -125 && FLT_MAX_EXP == 128)
  63. #define IEEE_FLOATING_POINT 1
  64. #else
  65. #define IEEE_FLOATING_POINT 0
  66. #endif
  67. #endif
  68.  
  69. /* Work around a problem that happens because math.h on hpux 7
  70.    defines two static variables--which, in Emacs, are not really static,
  71.    because `static' is defined as nothing.  The problem is that they are
  72.    defined both here and in lread.c.
  73.    These macros prevent the name conflict.  */
  74. #if defined (HPUX) && !defined (HPUX8)
  75. #define _MAXLDBL floatfns_maxldbl
  76. #define _NMAXLDBL floatfns_nmaxldbl
  77. #endif
  78.  
  79. #include <math.h>
  80.  
  81. /* This declaration is omitted on some systems, like Ultrix.  */
  82. #if !defined (HPUX) && defined (HAVE_LOGB) && !defined (logb)
  83. extern double logb ();
  84. #endif /* not HPUX and HAVE_LOGB and no logb macro */
  85.  
  86. #if defined(DOMAIN) && defined(SING) && defined(OVERFLOW)
  87.     /* If those are defined, then this is probably a `matherr' machine. */
  88. # ifndef HAVE_MATHERR
  89. #  define HAVE_MATHERR
  90. # endif
  91. #endif
  92.  
  93. #ifdef NO_MATHERR
  94. #undef HAVE_MATHERR
  95. #endif
  96.  
  97. #ifdef HAVE_MATHERR
  98. # ifdef FLOAT_CHECK_ERRNO
  99. #  undef FLOAT_CHECK_ERRNO
  100. # endif
  101. # ifdef FLOAT_CHECK_DOMAIN
  102. #  undef FLOAT_CHECK_DOMAIN
  103. # endif
  104. #endif
  105.  
  106. #ifndef NO_FLOAT_CHECK_ERRNO
  107. #define FLOAT_CHECK_ERRNO
  108. #endif
  109.  
  110. #ifdef FLOAT_CHECK_ERRNO
  111. # include <errno.h>
  112.  
  113. extern int errno;
  114. #endif
  115.  
  116. /* Avoid traps on VMS from sinh and cosh.
  117.    All the other functions set errno instead.  */
  118.  
  119. #ifdef VMS
  120. #undef cosh
  121. #undef sinh
  122. #define cosh(x) ((exp(x)+exp(-x))*0.5)
  123. #define sinh(x) ((exp(x)-exp(-x))*0.5)
  124. #endif /* VMS */
  125.  
  126. static SIGTYPE float_error ();
  127.  
  128. /* Nonzero while executing in floating point.
  129.    This tells float_error what to do.  */
  130.  
  131. static int in_float;
  132.  
  133. /* If an argument is out of range for a mathematical function,
  134.    here is the actual argument value to use in the error message.
  135.    These variables are used only across the floating point library call
  136.    so there is no need to staticpro them.  */
  137.  
  138. static Lisp_Object float_error_arg, float_error_arg2;
  139.  
  140. static char *float_error_fn_name;
  141.  
  142. /* Evaluate the floating point expression D, recording NUM
  143.    as the original argument for error messages.
  144.    D is normally an assignment expression.
  145.    Handle errors which may result in signals or may set errno.
  146.  
  147.    Note that float_error may be declared to return void, so you can't
  148.    just cast the zero after the colon to (SIGTYPE) to make the types
  149.    check properly.  */
  150.  
  151. #ifdef FLOAT_CHECK_ERRNO
  152. #define IN_FLOAT(d, name, num)                \
  153.   do {                            \
  154.     float_error_arg = num;                \
  155.     float_error_fn_name = name;                \
  156.     in_float = 1; errno = 0; (d); in_float = 0;        \
  157.     switch (errno) {                    \
  158.     case 0: break;                    \
  159.     case EDOM:     domain_error (float_error_fn_name, float_error_arg);    \
  160.     case ERANGE: range_error (float_error_fn_name, float_error_arg);    \
  161.     default:     arith_error (float_error_fn_name, float_error_arg);    \
  162.     }                            \
  163.   } while (0)
  164. #define IN_FLOAT2(d, name, num, num2)            \
  165.   do {                            \
  166.     float_error_arg = num;                \
  167.     float_error_arg2 = num2;                \
  168.     float_error_fn_name = name;                \
  169.     in_float = 1; errno = 0; (d); in_float = 0;        \
  170.     switch (errno) {                    \
  171.     case 0: break;                    \
  172.     case EDOM:     domain_error (float_error_fn_name, float_error_arg);    \
  173.     case ERANGE: range_error (float_error_fn_name, float_error_arg);    \
  174.     default:     arith_error (float_error_fn_name, float_error_arg);    \
  175.     }                            \
  176.   } while (0)
  177. #else
  178. #define IN_FLOAT(d, name, num) (in_float = 1, (d), in_float = 0)
  179. #define IN_FLOAT2(d, name, num, num2) (in_float = 1, (d), in_float = 0)
  180. #endif
  181.  
  182. /* Convert float to Lisp_Int if it fits, else signal a range error
  183.    using the given arguments.  */
  184. #define FLOAT_TO_INT(x, i, name, num)                    \
  185.   do                                    \
  186.     {                                    \
  187.       if ((x) >= (((EMACS_INT) 1) << (VALBITS-1)) ||            \
  188.       (x) <= - (((EMACS_INT) 1) << (VALBITS-1)) - 1)        \
  189.     range_error (name, num);                    \
  190.       XSETINT (i,  (EMACS_INT)(x));                    \
  191.     }                                    \
  192.   while (0)
  193. #define FLOAT_TO_INT2(x, i, name, num1, num2)                \
  194.   do                                    \
  195.     {                                    \
  196.       if ((x) >= (((EMACS_INT) 1) << (VALBITS-1)) ||            \
  197.       (x) <= - (((EMACS_INT) 1) << (VALBITS-1)) - 1)        \
  198.     range_error2 (name, num1, num2);                \
  199.       XSETINT (i,  (EMACS_INT)(x));                    \
  200.     }                                    \
  201.   while (0)
  202.  
  203. #define arith_error(op,arg) \
  204.   Fsignal (Qarith_error, Fcons (build_string ((op)), Fcons ((arg), Qnil)))
  205. #define range_error(op,arg) \
  206.   Fsignal (Qrange_error, Fcons (build_string ((op)), Fcons ((arg), Qnil)))
  207. #define range_error2(op,a1,a2) \
  208.   Fsignal (Qrange_error, Fcons (build_string ((op)), \
  209.                 Fcons ((a1), Fcons ((a2), Qnil))))
  210. #define domain_error(op,arg) \
  211.   Fsignal (Qdomain_error, Fcons (build_string ((op)), Fcons ((arg), Qnil)))
  212. #define domain_error2(op,a1,a2) \
  213.   Fsignal (Qdomain_error, Fcons (build_string ((op)), \
  214.                  Fcons ((a1), Fcons ((a2), Qnil))))
  215.  
  216. /* Extract a Lisp number as a `double', or signal an error.  */
  217.  
  218. double
  219. extract_float (num)
  220.      Lisp_Object num;
  221. {
  222.   CHECK_NUMBER_OR_FLOAT (num, 0);
  223.  
  224.   if (FLOATP (num))
  225.     return XFLOAT (num)->data;
  226.   return (double) XINT (num);
  227. }
  228.  
  229. /* Trig functions.  */
  230.  
  231. DEFUN ("acos", Facos, Sacos, 1, 1, 0,
  232.   "Return the inverse cosine of ARG.")
  233.   (arg)
  234.      register Lisp_Object arg;
  235. {
  236.   double d = extract_float (arg);
  237. #ifdef FLOAT_CHECK_DOMAIN
  238.   if (d > 1.0 || d < -1.0)
  239.     domain_error ("acos", arg);
  240. #endif
  241.   IN_FLOAT (d = acos (d), "acos", arg);
  242.   return make_float (d);
  243. }
  244.  
  245. DEFUN ("asin", Fasin, Sasin, 1, 1, 0,
  246.   "Return the inverse sine of ARG.")
  247.   (arg)
  248.      register Lisp_Object arg;
  249. {
  250.   double d = extract_float (arg);
  251. #ifdef FLOAT_CHECK_DOMAIN
  252.   if (d > 1.0 || d < -1.0)
  253.     domain_error ("asin", arg);
  254. #endif
  255.   IN_FLOAT (d = asin (d), "asin", arg);
  256.   return make_float (d);
  257. }
  258.  
  259. DEFUN ("atan", Fatan, Satan, 1, 1, 0,
  260.   "Return the inverse tangent of ARG.")
  261.   (arg)
  262.      register Lisp_Object arg;
  263. {
  264.   double d = extract_float (arg);
  265.   IN_FLOAT (d = atan (d), "atan", arg);
  266.   return make_float (d);
  267. }
  268.  
  269. DEFUN ("cos", Fcos, Scos, 1, 1, 0,
  270.   "Return the cosine of ARG.")
  271.   (arg)
  272.      register Lisp_Object arg;
  273. {
  274.   double d = extract_float (arg);
  275.   IN_FLOAT (d = cos (d), "cos", arg);
  276.   return make_float (d);
  277. }
  278.  
  279. DEFUN ("sin", Fsin, Ssin, 1, 1, 0,
  280.   "Return the sine of ARG.")
  281.   (arg)
  282.      register Lisp_Object arg;
  283. {
  284.   double d = extract_float (arg);
  285.   IN_FLOAT (d = sin (d), "sin", arg);
  286.   return make_float (d);
  287. }
  288.  
  289. DEFUN ("tan", Ftan, Stan, 1, 1, 0,
  290.   "Return the tangent of ARG.")
  291.   (arg)
  292.      register Lisp_Object arg;
  293. {
  294.   double d = extract_float (arg);
  295.   double c = cos (d);
  296. #ifdef FLOAT_CHECK_DOMAIN
  297.   if (c == 0.0)
  298.     domain_error ("tan", arg);
  299. #endif
  300.   IN_FLOAT (d = sin (d) / c, "tan", arg);
  301.   return make_float (d);
  302. }
  303.  
  304. #if 0 /* Leave these out unless we find there's a reason for them.  */
  305.  
  306. DEFUN ("bessel-j0", Fbessel_j0, Sbessel_j0, 1, 1, 0,
  307.   "Return the bessel function j0 of ARG.")
  308.   (arg)
  309.      register Lisp_Object arg;
  310. {
  311.   double d = extract_float (arg);
  312.   IN_FLOAT (d = j0 (d), "bessel-j0", arg);
  313.   return make_float (d);
  314. }
  315.  
  316. DEFUN ("bessel-j1", Fbessel_j1, Sbessel_j1, 1, 1, 0,
  317.   "Return the bessel function j1 of ARG.")
  318.   (arg)
  319.      register Lisp_Object arg;
  320. {
  321.   double d = extract_float (arg);
  322.   IN_FLOAT (d = j1 (d), "bessel-j1", arg);
  323.   return make_float (d);
  324. }
  325.  
  326. DEFUN ("bessel-jn", Fbessel_jn, Sbessel_jn, 2, 2, 0,
  327.   "Return the order N bessel function output jn of ARG.\n\
  328. The first arg (the order) is truncated to an integer.")
  329.   (n, arg)
  330.      register Lisp_Object n, arg;
  331. {
  332.   int i1 = extract_float (n);
  333.   double f2 = extract_float (arg);
  334.  
  335.   IN_FLOAT (f2 = jn (i1, f2), "bessel-jn", n);
  336.   return make_float (f2);
  337. }
  338.  
  339. DEFUN ("bessel-y0", Fbessel_y0, Sbessel_y0, 1, 1, 0,
  340.   "Return the bessel function y0 of ARG.")
  341.   (arg)
  342.      register Lisp_Object arg;
  343. {
  344.   double d = extract_float (arg);
  345.   IN_FLOAT (d = y0 (d), "bessel-y0", arg);
  346.   return make_float (d);
  347. }
  348.  
  349. DEFUN ("bessel-y1", Fbessel_y1, Sbessel_y1, 1, 1, 0,
  350.   "Return the bessel function y1 of ARG.")
  351.   (arg)
  352.      register Lisp_Object arg;
  353. {
  354.   double d = extract_float (arg);
  355.   IN_FLOAT (d = y1 (d), "bessel-y0", arg);
  356.   return make_float (d);
  357. }
  358.  
  359. DEFUN ("bessel-yn", Fbessel_yn, Sbessel_yn, 2, 2, 0,
  360.   "Return the order N bessel function output yn of ARG.\n\
  361. The first arg (the order) is truncated to an integer.")
  362.   (n, arg)
  363.      register Lisp_Object n, arg;
  364. {
  365.   int i1 = extract_float (n);
  366.   double f2 = extract_float (arg);
  367.  
  368.   IN_FLOAT (f2 = yn (i1, f2), "bessel-yn", n);
  369.   return make_float (f2);
  370. }
  371.  
  372. #endif
  373.  
  374. #if 0 /* Leave these out unless we see they are worth having.  */
  375.  
  376. DEFUN ("erf", Ferf, Serf, 1, 1, 0,
  377.   "Return the mathematical error function of ARG.")
  378.   (arg)
  379.      register Lisp_Object arg;
  380. {
  381.   double d = extract_float (arg);
  382.   IN_FLOAT (d = erf (d), "erf", arg);
  383.   return make_float (d);
  384. }
  385.  
  386. DEFUN ("erfc", Ferfc, Serfc, 1, 1, 0,
  387.   "Return the complementary error function of ARG.")
  388.   (arg)
  389.      register Lisp_Object arg;
  390. {
  391.   double d = extract_float (arg);
  392.   IN_FLOAT (d = erfc (d), "erfc", arg);
  393.   return make_float (d);
  394. }
  395.  
  396. DEFUN ("log-gamma", Flog_gamma, Slog_gamma, 1, 1, 0,
  397.   "Return the log gamma of ARG.")
  398.   (arg)
  399.      register Lisp_Object arg;
  400. {
  401.   double d = extract_float (arg);
  402.   IN_FLOAT (d = lgamma (d), "log-gamma", arg);
  403.   return make_float (d);
  404. }
  405.  
  406. DEFUN ("cube-root", Fcube_root, Scube_root, 1, 1, 0,
  407.   "Return the cube root of ARG.")
  408.   (arg)
  409.      register Lisp_Object arg;
  410. {
  411.   double d = extract_float (arg);
  412. #ifdef HAVE_CBRT
  413.   IN_FLOAT (d = cbrt (d), "cube-root", arg);
  414. #else
  415.   if (d >= 0.0)
  416.     IN_FLOAT (d = pow (d, 1.0/3.0), "cube-root", arg);
  417.   else
  418.     IN_FLOAT (d = -pow (-d, 1.0/3.0), "cube-root", arg);
  419. #endif
  420.   return make_float (d);
  421. }
  422.  
  423. #endif
  424.  
  425. DEFUN ("exp", Fexp, Sexp, 1, 1, 0,
  426.   "Return the exponential base e of ARG.")
  427.   (arg)
  428.      register Lisp_Object arg;
  429. {
  430.   double d = extract_float (arg);
  431. #ifdef FLOAT_CHECK_DOMAIN
  432.   if (d > 709.7827)   /* Assume IEEE doubles here */
  433.     range_error ("exp", arg);
  434.   else if (d < -709.0)
  435.     return make_float (0.0);
  436.   else
  437. #endif
  438.     IN_FLOAT (d = exp (d), "exp", arg);
  439.   return make_float (d);
  440. }
  441.  
  442. DEFUN ("expt", Fexpt, Sexpt, 2, 2, 0,
  443.   "Return the exponential ARG1 ** ARG2.")
  444.   (arg1, arg2)
  445.      register Lisp_Object arg1, arg2;
  446. {
  447.   double f1, f2;
  448.  
  449.   CHECK_NUMBER_OR_FLOAT (arg1, 0);
  450.   CHECK_NUMBER_OR_FLOAT (arg2, 0);
  451.   if (INTEGERP (arg1)     /* common lisp spec */
  452.       && INTEGERP (arg2)) /* don't promote, if both are ints */
  453.     {                /* this can be improved by pre-calculating */
  454.       EMACS_INT acc, x, y;    /* some binary powers of x then accumulating */
  455.       Lisp_Object val;
  456.  
  457.       x = XINT (arg1);
  458.       y = XINT (arg2);
  459.       acc = 1;
  460.       
  461.       if (y < 0)
  462.     {
  463.       if (x == 1)
  464.         acc = 1;
  465.       else if (x == -1)
  466.         acc = (y & 1) ? -1 : 1;
  467.       else
  468.         acc = 0;
  469.     }
  470.       else
  471.     {
  472.       while (y > 0)
  473.         {
  474.           if (y & 1)
  475.         acc *= x;
  476.           x *= x;
  477.           y = (unsigned)y >> 1;
  478.         }
  479.     }
  480.       XSETINT (val, acc);
  481.       return val;
  482.     }
  483.   f1 = FLOATP (arg1) ? XFLOAT (arg1)->data : XINT (arg1);
  484.   f2 = FLOATP (arg2) ? XFLOAT (arg2)->data : XINT (arg2);
  485.   /* Really should check for overflow, too */
  486.   if (f1 == 0.0 && f2 == 0.0)
  487.     f1 = 1.0;
  488. #ifdef FLOAT_CHECK_DOMAIN
  489.   else if ((f1 == 0.0 && f2 < 0.0) || (f1 < 0 && f2 != floor(f2)))
  490.     domain_error2 ("expt", arg1, arg2);
  491. #endif
  492.   IN_FLOAT2 (f1 = pow (f1, f2), "expt", arg1, arg2);
  493.   return make_float (f1);
  494. }
  495.  
  496. DEFUN ("log", Flog, Slog, 1, 2, 0,
  497.   "Return the natural logarithm of ARG.\n\
  498. If second optional argument BASE is given, return log ARG using that base.")
  499.   (arg, base)
  500.      register Lisp_Object arg, base;
  501. {
  502.   double d = extract_float (arg);
  503.  
  504. #ifdef FLOAT_CHECK_DOMAIN
  505.   if (d <= 0.0)
  506.     domain_error2 ("log", arg, base);
  507. #endif
  508.   if (NILP (base))
  509.     IN_FLOAT (d = log (d), "log", arg);
  510.   else
  511.     {
  512.       double b = extract_float (base);
  513.  
  514. #ifdef FLOAT_CHECK_DOMAIN
  515.       if (b <= 0.0 || b == 1.0)
  516.     domain_error2 ("log", arg, base);
  517. #endif
  518.       if (b == 10.0)
  519.     IN_FLOAT2 (d = log10 (d), "log", arg, base);
  520.       else
  521.     IN_FLOAT2 (d = log (d) / log (b), "log", arg, base);
  522.     }
  523.   return make_float (d);
  524. }
  525.  
  526. DEFUN ("log10", Flog10, Slog10, 1, 1, 0,
  527.   "Return the logarithm base 10 of ARG.")
  528.   (arg)
  529.      register Lisp_Object arg;
  530. {
  531.   double d = extract_float (arg);
  532. #ifdef FLOAT_CHECK_DOMAIN
  533.   if (d <= 0.0)
  534.     domain_error ("log10", arg);
  535. #endif
  536.   IN_FLOAT (d = log10 (d), "log10", arg);
  537.   return make_float (d);
  538. }
  539.  
  540. DEFUN ("sqrt", Fsqrt, Ssqrt, 1, 1, 0,
  541.   "Return the square root of ARG.")
  542.   (arg)
  543.      register Lisp_Object arg;
  544. {
  545.   double d = extract_float (arg);
  546. #ifdef FLOAT_CHECK_DOMAIN
  547.   if (d < 0.0)
  548.     domain_error ("sqrt", arg);
  549. #endif
  550.   IN_FLOAT (d = sqrt (d), "sqrt", arg);
  551.   return make_float (d);
  552. }
  553.  
  554. #if 0 /* Not clearly worth adding.  */
  555.  
  556. DEFUN ("acosh", Facosh, Sacosh, 1, 1, 0,
  557.   "Return the inverse hyperbolic cosine of ARG.")
  558.   (arg)
  559.      register Lisp_Object arg;
  560. {
  561.   double d = extract_float (arg);
  562. #ifdef FLOAT_CHECK_DOMAIN
  563.   if (d < 1.0)
  564.     domain_error ("acosh", arg);
  565. #endif
  566. #ifdef HAVE_INVERSE_HYPERBOLIC
  567.   IN_FLOAT (d = acosh (d), "acosh", arg);
  568. #else
  569.   IN_FLOAT (d = log (d + sqrt (d*d - 1.0)), "acosh", arg);
  570. #endif
  571.   return make_float (d);
  572. }
  573.  
  574. DEFUN ("asinh", Fasinh, Sasinh, 1, 1, 0,
  575.   "Return the inverse hyperbolic sine of ARG.")
  576.   (arg)
  577.      register Lisp_Object arg;
  578. {
  579.   double d = extract_float (arg);
  580. #ifdef HAVE_INVERSE_HYPERBOLIC
  581.   IN_FLOAT (d = asinh (d), "asinh", arg);
  582. #else
  583.   IN_FLOAT (d = log (d + sqrt (d*d + 1.0)), "asinh", arg);
  584. #endif
  585.   return make_float (d);
  586. }
  587.  
  588. DEFUN ("atanh", Fatanh, Satanh, 1, 1, 0,
  589.   "Return the inverse hyperbolic tangent of ARG.")
  590.   (arg)
  591.      register Lisp_Object arg;
  592. {
  593.   double d = extract_float (arg);
  594. #ifdef FLOAT_CHECK_DOMAIN
  595.   if (d >= 1.0 || d <= -1.0)
  596.     domain_error ("atanh", arg);
  597. #endif
  598. #ifdef HAVE_INVERSE_HYPERBOLIC
  599.   IN_FLOAT (d = atanh (d), "atanh", arg);
  600. #else
  601.   IN_FLOAT (d = 0.5 * log ((1.0 + d) / (1.0 - d)), "atanh", arg);
  602. #endif
  603.   return make_float (d);
  604. }
  605.  
  606. DEFUN ("cosh", Fcosh, Scosh, 1, 1, 0,
  607.   "Return the hyperbolic cosine of ARG.")
  608.   (arg)
  609.      register Lisp_Object arg;
  610. {
  611.   double d = extract_float (arg);
  612. #ifdef FLOAT_CHECK_DOMAIN
  613.   if (d > 710.0 || d < -710.0)
  614.     range_error ("cosh", arg);
  615. #endif
  616.   IN_FLOAT (d = cosh (d), "cosh", arg);
  617.   return make_float (d);
  618. }
  619.  
  620. DEFUN ("sinh", Fsinh, Ssinh, 1, 1, 0,
  621.   "Return the hyperbolic sine of ARG.")
  622.   (arg)
  623.      register Lisp_Object arg;
  624. {
  625.   double d = extract_float (arg);
  626. #ifdef FLOAT_CHECK_DOMAIN
  627.   if (d > 710.0 || d < -710.0)
  628.     range_error ("sinh", arg);
  629. #endif
  630.   IN_FLOAT (d = sinh (d), "sinh", arg);
  631.   return make_float (d);
  632. }
  633.  
  634. DEFUN ("tanh", Ftanh, Stanh, 1, 1, 0,
  635.   "Return the hyperbolic tangent of ARG.")
  636.   (arg)
  637.      register Lisp_Object arg;
  638. {
  639.   double d = extract_float (arg);
  640.   IN_FLOAT (d = tanh (d), "tanh", arg);
  641.   return make_float (d);
  642. }
  643. #endif
  644.  
  645. DEFUN ("abs", Fabs, Sabs, 1, 1, 0,
  646.   "Return the absolute value of ARG.")
  647.   (arg)
  648.      register Lisp_Object arg;
  649. {
  650.   CHECK_NUMBER_OR_FLOAT (arg, 0);
  651.  
  652.   if (FLOATP (arg))
  653.     IN_FLOAT (arg = make_float (fabs (XFLOAT (arg)->data)), "abs", arg);
  654.   else if (XINT (arg) < 0)
  655.     XSETINT (arg, - XINT (arg));
  656.  
  657.   return arg;
  658. }
  659.  
  660. DEFUN ("float", Ffloat, Sfloat, 1, 1, 0,
  661.   "Return the floating point number equal to ARG.")
  662.   (arg)
  663.      register Lisp_Object arg;
  664. {
  665.   CHECK_NUMBER_OR_FLOAT (arg, 0);
  666.  
  667.   if (INTEGERP (arg))
  668.     return make_float ((double) XINT (arg));
  669.   else                /* give 'em the same float back */
  670.     return arg;
  671. }
  672.  
  673. DEFUN ("logb", Flogb, Slogb, 1, 1, 0,
  674.   "Returns largest integer <= the base 2 log of the magnitude of ARG.\n\
  675. This is the same as the exponent of a float.")
  676.      (arg)
  677.      Lisp_Object arg;
  678. {
  679.   Lisp_Object val;
  680.   EMACS_INT value;
  681.   double f = extract_float (arg);
  682.  
  683.   if (f == 0.0)
  684.     value = -(VALMASK >> 1);
  685.   else
  686.     {
  687. #ifdef HAVE_LOGB
  688.       IN_FLOAT (value = logb (f), "logb", arg);
  689. #else
  690. #ifdef HAVE_FREXP
  691.       int ivalue;
  692.       IN_FLOAT (frexp (f, &ivalue), "logb", arg);
  693.       value = ivalue - 1;
  694. #else
  695.       int i;
  696.       double d;
  697.       if (f < 0.0)
  698.     f = -f;
  699.       value = -1;
  700.       while (f < 0.5)
  701.     {
  702.       for (i = 1, d = 0.5; d * d >= f; i += i)
  703.         d *= d;
  704.       f /= d;
  705.       value -= i;
  706.     }
  707.       while (f >= 1.0)
  708.     {
  709.       for (i = 1, d = 2.0; d * d <= f; i += i)
  710.         d *= d;
  711.       f /= d;
  712.       value += i;
  713.     }
  714. #endif
  715. #endif
  716.     }
  717.   XSETINT (val, value);
  718.   return val;
  719. }
  720.  
  721. #endif /* LISP_FLOAT_TYPE */
  722.  
  723.  
  724. /* the rounding functions  */
  725.  
  726. static Lisp_Object
  727. rounding_driver (arg, divisor, double_round, int_round2, name)
  728.      register Lisp_Object arg, divisor;
  729.      double (*double_round) ();
  730.      EMACS_INT (*int_round2) ();
  731.      char *name;
  732. {
  733.   CHECK_NUMBER_OR_FLOAT (arg, 0);
  734.  
  735.   if (! NILP (divisor))
  736.     {
  737.       EMACS_INT i1, i2;
  738.  
  739.       CHECK_NUMBER_OR_FLOAT (divisor, 1);
  740.  
  741. #ifdef LISP_FLOAT_TYPE
  742.       if (FLOATP (arg) || FLOATP (divisor))
  743.     {
  744.       double f1, f2;
  745.  
  746.       f1 = FLOATP (arg) ? XFLOAT (arg)->data : XINT (arg);
  747.       f2 = (FLOATP (divisor) ? XFLOAT (divisor)->data : XINT (divisor));
  748.       if (! IEEE_FLOATING_POINT && f2 == 0)
  749.         Fsignal (Qarith_error, Qnil);
  750.  
  751.       IN_FLOAT2 (f1 = (*double_round) (f1 / f2), name, arg, divisor);
  752.       FLOAT_TO_INT2 (f1, arg, name, arg, divisor);
  753.       return arg;
  754.     }
  755. #endif
  756.  
  757.       i1 = XINT (arg);
  758.       i2 = XINT (divisor);
  759.  
  760.       if (i2 == 0)
  761.     Fsignal (Qarith_error, Qnil);
  762.  
  763.       XSETINT (arg, (*int_round2) (i1, i2));
  764.       return arg;
  765.     }
  766.  
  767. #ifdef LISP_FLOAT_TYPE
  768.   if (FLOATP (arg))
  769.     {
  770.       double d;
  771.  
  772.       IN_FLOAT (d = (*double_round) (XFLOAT (arg)->data), name, arg);
  773.       FLOAT_TO_INT (d, arg, name, arg);
  774.     }
  775. #endif
  776.  
  777.   return arg;
  778. }
  779.  
  780. /* With C's /, the result is implementation-defined if either operand
  781.    is negative, so take care with negative operands in the following
  782.    integer functions.  */
  783.  
  784. static EMACS_INT
  785. ceiling2 (i1, i2)
  786.      EMACS_INT i1, i2;
  787. {
  788.   return (i2 < 0
  789.       ? (i1 < 0  ?  ((-1 - i1) / -i2) + 1  :  - (i1 / -i2))
  790.       : (i1 <= 0  ?  - (-i1 / i2)  :  ((i1 - 1) / i2) + 1));
  791. }
  792.  
  793. static EMACS_INT
  794. floor2 (i1, i2)
  795.      EMACS_INT i1, i2;
  796. {
  797.   return (i2 < 0
  798.       ? (i1 <= 0  ?  -i1 / -i2  :  -1 - ((i1 - 1) / -i2))
  799.       : (i1 < 0  ?  -1 - ((-1 - i1) / i2)  :  i1 / i2));
  800. }
  801.  
  802. static EMACS_INT
  803. truncate2 (i1, i2)
  804.      EMACS_INT i1, i2;
  805. {
  806.   return (i2 < 0
  807.       ? (i1 < 0  ?  -i1 / -i2  :  - (i1 / -i2))
  808.       : (i1 < 0  ?  - (-i1 / i2)  :  i1 / i2));
  809. }
  810.  
  811. static EMACS_INT
  812. round2 (i1, i2)
  813.      EMACS_INT i1, i2;
  814. {
  815.   /* The C language's division operator gives us one remainder R, but
  816.      we want the remainder R1 on the other side of 0 if R1 is closer
  817.      to 0 than R is; because we want to round to even, we also want R1
  818.      if R and R1 are the same distance from 0 and if C's quotient is
  819.      odd.  */
  820.   EMACS_INT q = i1 / i2;
  821.   EMACS_INT r = i1 % i2;
  822.   EMACS_INT abs_r = r < 0 ? -r : r;
  823.   EMACS_INT abs_r1 = (i2 < 0 ? -i2 : i2) - abs_r;
  824.   return q + (abs_r + (q & 1) <= abs_r1 ? 0 : (i2 ^ r) < 0 ? -1 : 1);
  825. }
  826.  
  827. /* The code uses emacs_rint, so that it works to undefine HAVE_RINT
  828.    if `rint' exists but does not work right.  */
  829. #ifdef HAVE_RINT
  830. #define emacs_rint rint
  831. #else
  832. static double
  833. emacs_rint (d)
  834.      double d;
  835. {
  836.   return floor (d + 0.5);
  837. }
  838. #endif
  839.  
  840. static double
  841. double_identity (d)
  842.      double d;
  843. {
  844.   return d;
  845. }
  846.  
  847. DEFUN ("ceiling", Fceiling, Sceiling, 1, 2, 0,
  848.   "Return the smallest integer no less than ARG.  (Round toward +inf.)\n\
  849. With optional DIVISOR, return the smallest integer no less than ARG/DIVISOR.")
  850.   (arg, divisor)
  851.      Lisp_Object arg, divisor;
  852. {
  853.   return rounding_driver (arg, divisor, ceil, ceiling2, "ceiling");
  854. }
  855.  
  856. DEFUN ("floor", Ffloor, Sfloor, 1, 2, 0,
  857.   "Return the largest integer no greater than ARG.  (Round towards -inf.)\n\
  858. With optional DIVISOR, return the largest integer no greater than ARG/DIVISOR.")
  859.   (arg, divisor)
  860.      Lisp_Object arg, divisor;
  861. {
  862.   return rounding_driver (arg, divisor, floor, floor2, "floor");
  863. }
  864.  
  865. DEFUN ("round", Fround, Sround, 1, 2, 0,
  866.   "Return the nearest integer to ARG.\n\
  867. With optional DIVISOR, return the nearest integer to ARG/DIVISOR.")
  868.   (arg, divisor)
  869.      Lisp_Object arg, divisor;
  870. {
  871.   return rounding_driver (arg, divisor, emacs_rint, round2, "round");
  872. }
  873.  
  874. DEFUN ("truncate", Ftruncate, Struncate, 1, 2, 0,
  875.        "Truncate a floating point number to an int.\n\
  876. Rounds ARG toward zero.\n\
  877. With optional DIVISOR, truncate ARG/DIVISOR.")
  878.   (arg, divisor)
  879.      Lisp_Object arg, divisor;
  880. {
  881.   return rounding_driver (arg, divisor, double_identity, truncate2,
  882.               "truncate");
  883. }
  884.  
  885. #ifdef LISP_FLOAT_TYPE
  886.  
  887. Lisp_Object
  888. fmod_float (x, y)
  889.      register Lisp_Object x, y;
  890. {
  891.   double f1, f2;
  892.  
  893.   f1 = FLOATP (x) ? XFLOAT (x)->data : XINT (x);
  894.   f2 = FLOATP (y) ? XFLOAT (y)->data : XINT (y);
  895.  
  896.   if (! IEEE_FLOATING_POINT && f2 == 0)
  897.     Fsignal (Qarith_error, Qnil);
  898.  
  899.   /* If the "remainder" comes out with the wrong sign, fix it.  */
  900.   IN_FLOAT2 ((f1 = fmod (f1, f2),
  901.           f1 = (f2 < 0 ? f1 > 0 : f1 < 0) ? f1 + f2 : f1),
  902.          "mod", x, y);
  903.   return make_float (f1);
  904. }
  905.  
  906. /* It's not clear these are worth adding.  */
  907.  
  908. DEFUN ("fceiling", Ffceiling, Sfceiling, 1, 1, 0,
  909.   "Return the smallest integer no less than ARG, as a float.\n\
  910. \(Round toward +inf.\)")
  911.   (arg)
  912.      register Lisp_Object arg;
  913. {
  914.   double d = extract_float (arg);
  915.   IN_FLOAT (d = ceil (d), "fceiling", arg);
  916.   return make_float (d);
  917. }
  918.  
  919. DEFUN ("ffloor", Fffloor, Sffloor, 1, 1, 0,
  920.   "Return the largest integer no greater than ARG, as a float.\n\
  921. \(Round towards -inf.\)")
  922.   (arg)
  923.      register Lisp_Object arg;
  924. {
  925.   double d = extract_float (arg);
  926.   IN_FLOAT (d = floor (d), "ffloor", arg);
  927.   return make_float (d);
  928. }
  929.  
  930. DEFUN ("fround", Ffround, Sfround, 1, 1, 0,
  931.   "Return the nearest integer to ARG, as a float.")
  932.   (arg)
  933.      register Lisp_Object arg;
  934. {
  935.   double d = extract_float (arg);
  936.   IN_FLOAT (d = emacs_rint (d), "fround", arg);
  937.   return make_float (d);
  938. }
  939.  
  940. DEFUN ("ftruncate", Fftruncate, Sftruncate, 1, 1, 0,
  941.   "Truncate a floating point number to an integral float value.\n\
  942. Rounds the value toward zero.")
  943.   (arg)
  944.      register Lisp_Object arg;
  945. {
  946.   double d = extract_float (arg);
  947.   if (d >= 0.0)
  948.     IN_FLOAT (d = floor (d), "ftruncate", arg);
  949.   else
  950.     IN_FLOAT (d = ceil (d), "ftruncate", arg);
  951.   return make_float (d);
  952. }
  953.  
  954. #ifdef FLOAT_CATCH_SIGILL
  955. static SIGTYPE
  956. float_error (signo)
  957.      int signo;
  958. {
  959.   if (! in_float)
  960.     fatal_error_signal (signo);
  961.  
  962. #ifdef BSD_SYSTEM
  963. #ifdef BSD4_1
  964.   sigrelse (SIGILL);
  965. #else /* not BSD4_1 */
  966.   sigsetmask (SIGEMPTYMASK);
  967. #endif /* not BSD4_1 */
  968. #else
  969.   /* Must reestablish handler each time it is called.  */
  970.   signal (SIGILL, float_error);
  971. #endif /* BSD_SYSTEM */
  972.  
  973.   in_float = 0;
  974.  
  975.   Fsignal (Qarith_error, Fcons (float_error_arg, Qnil));
  976. }
  977.  
  978. /* Another idea was to replace the library function `infnan'
  979.    where SIGILL is signaled.  */
  980.  
  981. #endif /* FLOAT_CATCH_SIGILL */
  982.  
  983. #ifdef HAVE_MATHERR
  984. int 
  985. matherr (x)
  986.      struct exception *x;
  987. {
  988.   Lisp_Object args;
  989.   if (! in_float)
  990.     /* Not called from emacs-lisp float routines; do the default thing. */
  991.     return 0;
  992.   if (!strcmp (x->name, "pow"))
  993.     x->name = "expt";
  994.  
  995.   args
  996.     = Fcons (build_string (x->name),
  997.          Fcons (make_float (x->arg1),
  998.             ((!strcmp (x->name, "log") || !strcmp (x->name, "pow"))
  999.              ? Fcons (make_float (x->arg2), Qnil)
  1000.              : Qnil)));
  1001.   switch (x->type)
  1002.     {
  1003.     case DOMAIN:    Fsignal (Qdomain_error, args);        break;
  1004.     case SING:        Fsignal (Qsingularity_error, args);    break;
  1005.     case OVERFLOW:    Fsignal (Qoverflow_error, args);    break;
  1006.     case UNDERFLOW:    Fsignal (Qunderflow_error, args);    break;
  1007.     default:        Fsignal (Qarith_error, args);        break;
  1008.     }
  1009.   return (1);    /* don't set errno or print a message */
  1010. }
  1011. #endif /* HAVE_MATHERR */
  1012.  
  1013. void
  1014. init_floatfns ()
  1015. {
  1016. #ifdef FLOAT_CATCH_SIGILL
  1017.   signal (SIGILL, float_error);
  1018. #endif 
  1019.   in_float = 0;
  1020. }
  1021.  
  1022. #else /* not LISP_FLOAT_TYPE */
  1023.  
  1024. init_floatfns ()
  1025. {}
  1026.  
  1027. #endif /* not LISP_FLOAT_TYPE */
  1028.  
  1029. void
  1030. syms_of_floatfns ()
  1031. {
  1032. #ifdef LISP_FLOAT_TYPE
  1033.   defsubr (&Sacos);
  1034.   defsubr (&Sasin);
  1035.   defsubr (&Satan);
  1036.   defsubr (&Scos);
  1037.   defsubr (&Ssin);
  1038.   defsubr (&Stan);
  1039. #if 0
  1040.   defsubr (&Sacosh);
  1041.   defsubr (&Sasinh);
  1042.   defsubr (&Satanh);
  1043.   defsubr (&Scosh);
  1044.   defsubr (&Ssinh);
  1045.   defsubr (&Stanh);
  1046.   defsubr (&Sbessel_y0);
  1047.   defsubr (&Sbessel_y1);
  1048.   defsubr (&Sbessel_yn);
  1049.   defsubr (&Sbessel_j0);
  1050.   defsubr (&Sbessel_j1);
  1051.   defsubr (&Sbessel_jn);
  1052.   defsubr (&Serf);
  1053.   defsubr (&Serfc);
  1054.   defsubr (&Slog_gamma);
  1055.   defsubr (&Scube_root);
  1056. #endif
  1057.   defsubr (&Sfceiling);
  1058.   defsubr (&Sffloor);
  1059.   defsubr (&Sfround);
  1060.   defsubr (&Sftruncate);
  1061.   defsubr (&Sexp);
  1062.   defsubr (&Sexpt);
  1063.   defsubr (&Slog);
  1064.   defsubr (&Slog10);
  1065.   defsubr (&Ssqrt);
  1066.  
  1067.   defsubr (&Sabs);
  1068.   defsubr (&Sfloat);
  1069.   defsubr (&Slogb);
  1070. #endif /* LISP_FLOAT_TYPE */
  1071.   defsubr (&Sceiling);
  1072.   defsubr (&Sfloor);
  1073.   defsubr (&Sround);
  1074.   defsubr (&Struncate);
  1075. }
  1076.