1: /* Primitive operations on floating point for GNU Emacs Lisp interpreter.
    2:    Copyright (C) 1988, 1993, 1994, 1999, 2001, 2002, 2003, 2004,
    3:                  2005, 2006, 2007  Free Software Foundation, Inc.
    4: 
    5: This file is part of GNU Emacs.
    6: 
    7: GNU Emacs is free software; you can redistribute it and/or modify
    8: it under the terms of the GNU General Public License as published by
    9: the Free Software Foundation; either version 2, or (at your option)
   10: any later version.
   11: 
   12: GNU Emacs is distributed in the hope that it will be useful,
   13: but WITHOUT ANY WARRANTY; without even the implied warranty of
   14: MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   15: GNU General Public License for more details.
   16: 
   17: You should have received a copy of the GNU General Public License
   18: along with GNU Emacs; see the file COPYING.  If not, write to
   19: the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
   20: Boston, MA 02110-1301, USA.  */
   21: 
   22: 
   23: /* ANSI C requires only these float functions:
   24:    acos, asin, atan, atan2, ceil, cos, cosh, exp, fabs, floor, fmod,
   25:    frexp, ldexp, log, log10, modf, pow, sin, sinh, sqrt, tan, tanh.
   26: 
   27:    Define HAVE_INVERSE_HYPERBOLIC if you have acosh, asinh, and atanh.
   28:    Define HAVE_CBRT if you have cbrt.
   29:    Define HAVE_RINT if you have a working rint.
   30:    If you don't define these, then the appropriate routines will be simulated.
   31: 
   32:    Define HAVE_MATHERR if on a system supporting the SysV matherr callback.
   33:    (This should happen automatically.)
   34: 
   35:    Define FLOAT_CHECK_ERRNO if the float library routines set errno.
   36:    This has no effect if HAVE_MATHERR is defined.
   37: 
   38:    Define FLOAT_CATCH_SIGILL if the float library routines signal SIGILL.
   39:    (What systems actually do this?  Please let us know.)
   40: 
   41:    Define FLOAT_CHECK_DOMAIN if the float library doesn't handle errors by
   42:    either setting errno, or signaling SIGFPE/SIGILL.  Otherwise, domain and
   43:    range checking will happen before calling the float routines.  This has
   44:    no effect if HAVE_MATHERR is defined (since matherr will be called when
   45:    a domain error occurs.)
   46:  */
   47: 
   48: #include <config.h>
   49: #include <signal.h>
   50: #include "lisp.h"
   51: #include "syssignal.h"
   52: 
   53: #if STDC_HEADERS
   54: #include <float.h>
   55: #endif
   56: 
   57: /* If IEEE_FLOATING_POINT isn't defined, default it from FLT_*. */
   58: #ifndef IEEE_FLOATING_POINT
   59: #if (FLT_RADIX == 2 && FLT_MANT_DIG == 24 \
   60:      && FLT_MIN_EXP == -125 && FLT_MAX_EXP == 128)
   61: #define IEEE_FLOATING_POINT 1
   62: #else
   63: #define IEEE_FLOATING_POINT 0
   64: #endif
   65: #endif
   66: 
   67: /* Work around a problem that happens because math.h on hpux 7
   68:    defines two static variables--which, in Emacs, are not really static,
   69:    because `static' is defined as nothing.  The problem is that they are
   70:    defined both here and in lread.c.
   71:    These macros prevent the name conflict.  */
   72: #if defined (HPUX) && !defined (HPUX8)
   73: #define _MAXLDBL floatfns_maxldbl
   74: #define _NMAXLDBL floatfns_nmaxldbl
   75: #endif
   76: 
   77: #include <math.h>
   78: 
   79: /* This declaration is omitted on some systems, like Ultrix.  */
   80: #if !defined (HPUX) && defined (HAVE_LOGB) && !defined (logb)
   81: extern double logb ();
   82: #endif /* not HPUX and HAVE_LOGB and no logb macro */
   83: 
   84: #if defined(DOMAIN) && defined(SING) && defined(OVERFLOW)
   85:     /* If those are defined, then this is probably a `matherr' machine. */
   86: # ifndef HAVE_MATHERR
   87: #  define HAVE_MATHERR
   88: # endif
   89: #endif
   90: 
   91: #ifdef NO_MATHERR
   92: #undef HAVE_MATHERR
   93: #endif
   94: 
   95: #ifdef HAVE_MATHERR
   96: # ifdef FLOAT_CHECK_ERRNO
   97: #  undef FLOAT_CHECK_ERRNO
   98: # endif
   99: # ifdef FLOAT_CHECK_DOMAIN
  100: #  undef FLOAT_CHECK_DOMAIN
  101: # endif
  102: #endif
  103: 
  104: #ifndef NO_FLOAT_CHECK_ERRNO
  105: #define FLOAT_CHECK_ERRNO
  106: #endif
  107: 
  108: #ifdef FLOAT_CHECK_ERRNO
  109: # include <errno.h>
  110: 
  111: #ifndef USE_CRT_DLL
  112: extern int errno;
  113: #endif
  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: #ifdef FLOAT_CATCH_SIGILL
  127: static SIGTYPE float_error ();
  128: #endif
  129: 
  130: /* Nonzero while executing in floating point.
  131:    This tells float_error what to do.  */
  132: 
  133: static int in_float;
  134: 
  135: /* If an argument is out of range for a mathematical function,
  136:    here is the actual argument value to use in the error message.
  137:    These variables are used only across the floating point library call
  138:    so there is no need to staticpro them.  */
  139: 
  140: static Lisp_Object float_error_arg, float_error_arg2;
  141: 
  142: static char *float_error_fn_name;
  143: 
  144: /* Evaluate the floating point expression D, recording NUM
  145:    as the original argument for error messages.
  146:    D is normally an assignment expression.
  147:    Handle errors which may result in signals or may set errno.
  148: 
  149:    Note that float_error may be declared to return void, so you can't
  150:    just cast the zero after the colon to (SIGTYPE) to make the types
  151:    check properly.  */
  152: 
  153: #ifdef FLOAT_CHECK_ERRNO
  154: #define IN_FLOAT(d, name, num)                          \
  155:   do {                                                  \
  156:     float_error_arg = num;                              \
  157:     float_error_fn_name = name;                         \
  158:     in_float = 1; errno = 0; (d); in_float = 0;         \
  159:     switch (errno) {                                    \
  160:     case 0: break;                                      \
  161:     case EDOM:   domain_error (float_error_fn_name, float_error_arg);    \
  162:     case ERANGE: range_error (float_error_fn_name, float_error_arg);    \
  163:     default:     arith_error (float_error_fn_name, float_error_arg);       \
  164:     }                                                   \
  165:   } while (0)
  166: #define IN_FLOAT2(d, name, num, num2)                   \
  167:   do {                                                  \
  168:     float_error_arg = num;                              \
  169:     float_error_arg2 = num2;                            \
  170:     float_error_fn_name = name;                         \
  171:     in_float = 1; errno = 0; (d); in_float = 0;         \
  172:     switch (errno) {                                    \
  173:     case 0: break;                                      \
  174:     case EDOM:   domain_error (float_error_fn_name, float_error_arg);    \
  175:     case ERANGE: range_error (float_error_fn_name, float_error_arg);    \
  176:     default:     arith_error (float_error_fn_name, float_error_arg);       \
  177:     }                                                   \
  178:   } while (0)
  179: #else
  180: #define IN_FLOAT(d, name, num) (in_float = 1, (d), in_float = 0)
  181: #define IN_FLOAT2(d, name, num, num2) (in_float = 1, (d), in_float = 0)
  182: #endif
  183: 
  184: /* Convert float to Lisp_Int if it fits, else signal a range error
  185:    using the given arguments.  */
  186: #define FLOAT_TO_INT(x, i, name, num)                                   \
  187:   do                                                                    \
  188:     {                                                                   \
  189:       if (FIXNUM_OVERFLOW_P (x))                                        \
  190:         range_error (name, num);                                       \
  191:       XSETINT (i,  (EMACS_INT)(x));                                     \
  192:     }                                                                   \
  193:   while (0)
  194: #define FLOAT_TO_INT2(x, i, name, num1, num2)                           \
  195:   do                                                                    \
  196:     {                                                                   \
  197:       if (FIXNUM_OVERFLOW_P (x))                                        \
  198:         range_error2 (name, num1, num2);                               \
  199:       XSETINT (i,  (EMACS_INT)(x));                                     \
  200:     }                                                                   \
  201:   while (0)
  202: 
  203: #define arith_error(op,arg) \
  204:   xsignal2 (Qarith_error, build_string ((op)), (arg))
  205: #define range_error(op,arg) \
  206:   xsignal2 (Qrange_error, build_string ((op)), (arg))
  207: #define range_error2(op,a1,a2) \
  208:   xsignal3 (Qrange_error, build_string ((op)), (a1), (a2))
  209: #define domain_error(op,arg) \
  210:   xsignal2 (Qdomain_error, build_string ((op)), (arg))
  211: #define domain_error2(op,a1,a2) \
  212:   xsignal3 (Qdomain_error, build_string ((op)), (a1), (a2))
  213: 
  214: /* Extract a Lisp number as a `double', or signal an error.  */
  215: 
  216: double
  217: extract_float (num)
  218:      Lisp_Object num;
  219: {
  220:   CHECK_NUMBER_OR_FLOAT (num);
  221: 
  222:   if (FLOATP (num))
  223:     return XFLOAT_DATA (num);
  224:   return (double) XINT (num);
  225: }
  226: ^L
  227: /* Trig functions.  */
  228: 
  229: DEFUN ("acos", Facos, Sacos, 1, 1, 0,
  230:        doc: /* Return the inverse cosine of ARG.  */)
  231:      (arg)
  232:      register Lisp_Object arg;
  233: {
  234:   double d = extract_float (arg);
  235: #ifdef FLOAT_CHECK_DOMAIN
  236:   if (d > 1.0 || d < -1.0)
  237:     domain_error ("acos", arg);
  238: #endif
  239:   IN_FLOAT (d = acos (d), "acos", arg);
  240:   return make_float (d);
  241: }
  242: 
  243: DEFUN ("asin", Fasin, Sasin, 1, 1, 0,
  244:        doc: /* Return the inverse sine of ARG.  */)
  245:      (arg)
  246:      register Lisp_Object arg;
  247: {
  248:   double d = extract_float (arg);
  249: #ifdef FLOAT_CHECK_DOMAIN
  250:   if (d > 1.0 || d < -1.0)
  251:     domain_error ("asin", arg);
  252: #endif
  253:   IN_FLOAT (d = asin (d), "asin", arg);
  254:   return make_float (d);
  255: }
  256: 
  257: DEFUN ("atan", Fatan, Satan, 1, 2, 0,
  258:        doc: /* Return the inverse tangent of the arguments.
  259: If only one argument Y is given, return the inverse tangent of Y.
  260: If two arguments Y and X are given, return the inverse tangent of Y
  261: divided by X, i.e. the angle in radians between the vector (X, Y)
  262: and the x-axis.  */)
  263:      (y, x)
  264:      register Lisp_Object y, x;
  265: {
  266:   double d = extract_float (y);
  267: 
  268:   if (NILP (x))
  269:     IN_FLOAT (d = atan (d), "atan", y);
  270:   else
  271:     {
  272:       double d2 = extract_float (x);
  273: 
  274:       IN_FLOAT2 (d = atan2 (d, d2), "atan", y, x);
  275:     }
  276:   return make_float (d);
  277: }
  278: 
  279: DEFUN ("cos", Fcos, Scos, 1, 1, 0,
  280:        doc: /* Return the cosine of ARG.  */)
  281:      (arg)
  282:      register Lisp_Object arg;
  283: {
  284:   double d = extract_float (arg);
  285:   IN_FLOAT (d = cos (d), "cos", arg);
  286:   return make_float (d);
  287: }
  288: 
  289: DEFUN ("sin", Fsin, Ssin, 1, 1, 0,
  290:        doc: /* Return the sine of ARG.  */)
  291:      (arg)
  292:      register Lisp_Object arg;
  293: {
  294:   double d = extract_float (arg);
  295:   IN_FLOAT (d = sin (d), "sin", arg);
  296:   return make_float (d);
  297: }
  298: 
  299: DEFUN ("tan", Ftan, Stan, 1, 1, 0,
  300:        doc: /* Return the tangent of ARG.  */)
  301:      (arg)
  302:      register Lisp_Object arg;
  303: {
  304:   double d = extract_float (arg);
  305:   double c = cos (d);
  306: #ifdef FLOAT_CHECK_DOMAIN
  307:   if (c == 0.0)
  308:     domain_error ("tan", arg);
  309: #endif
  310:   IN_FLOAT (d = sin (d) / c, "tan", arg);
  311:   return make_float (d);
  312: }
  313: ^L
  314: #if 0 /* Leave these out unless we find there's a reason for them.  */
  315: 
  316: DEFUN ("bessel-j0", Fbessel_j0, Sbessel_j0, 1, 1, 0,
  317:        doc: /* Return the bessel function j0 of ARG.  */)
  318:      (arg)
  319:      register Lisp_Object arg;
  320: {
  321:   double d = extract_float (arg);
  322:   IN_FLOAT (d = j0 (d), "bessel-j0", arg);
  323:   return make_float (d);
  324: }
  325: 
  326: DEFUN ("bessel-j1", Fbessel_j1, Sbessel_j1, 1, 1, 0,
  327:        doc: /* Return the bessel function j1 of ARG.  */)
  328:      (arg)
  329:      register Lisp_Object arg;
  330: {
  331:   double d = extract_float (arg);
  332:   IN_FLOAT (d = j1 (d), "bessel-j1", arg);
  333:   return make_float (d);
  334: }
  335: 
  336: DEFUN ("bessel-jn", Fbessel_jn, Sbessel_jn, 2, 2, 0,
  337:        doc: /* Return the order N bessel function output jn of ARG.
  338: The first arg (the order) is truncated to an integer.  */)
  339:      (n, arg)
  340:      register Lisp_Object n, arg;
  341: {
  342:   int i1 = extract_float (n);
  343:   double f2 = extract_float (arg);
  344: 
  345:   IN_FLOAT (f2 = jn (i1, f2), "bessel-jn", n);
  346:   return make_float (f2);
  347: }
  348: 
  349: DEFUN ("bessel-y0", Fbessel_y0, Sbessel_y0, 1, 1, 0,
  350:        doc: /* Return the bessel function y0 of ARG.  */)
  351:      (arg)
  352:      register Lisp_Object arg;
  353: {
  354:   double d = extract_float (arg);
  355:   IN_FLOAT (d = y0 (d), "bessel-y0", arg);
  356:   return make_float (d);
  357: }
  358: 
  359: DEFUN ("bessel-y1", Fbessel_y1, Sbessel_y1, 1, 1, 0,
  360:        doc: /* Return the bessel function y1 of ARG.  */)
  361:      (arg)
  362:      register Lisp_Object arg;
  363: {
  364:   double d = extract_float (arg);
  365:   IN_FLOAT (d = y1 (d), "bessel-y0", arg);
  366:   return make_float (d);
  367: }
  368: 
  369: DEFUN ("bessel-yn", Fbessel_yn, Sbessel_yn, 2, 2, 0,
  370:        doc: /* Return the order N bessel function output yn of ARG.
  371: The first arg (the order) is truncated to an integer.  */)
  372:      (n, arg)
  373:      register Lisp_Object n, arg;
  374: {
  375:   int i1 = extract_float (n);
  376:   double f2 = extract_float (arg);
  377: 
  378:   IN_FLOAT (f2 = yn (i1, f2), "bessel-yn", n);
  379:   return make_float (f2);
  380: }
  381: 
  382: #endif
  383: ^L
  384: #if 0 /* Leave these out unless we see they are worth having.  */
  385: 
  386: DEFUN ("erf", Ferf, Serf, 1, 1, 0,
  387:        doc: /* Return the mathematical error function of ARG.  */)
  388:      (arg)
  389:      register Lisp_Object arg;
  390: {
  391:   double d = extract_float (arg);
  392:   IN_FLOAT (d = erf (d), "erf", arg);
  393:   return make_float (d);
  394: }
  395: 
  396: DEFUN ("erfc", Ferfc, Serfc, 1, 1, 0,
  397:        doc: /* Return the complementary error function of ARG.  */)
  398:      (arg)
  399:      register Lisp_Object arg;
  400: {
  401:   double d = extract_float (arg);
  402:   IN_FLOAT (d = erfc (d), "erfc", arg);
  403:   return make_float (d);
  404: }
  405: 
  406: DEFUN ("log-gamma", Flog_gamma, Slog_gamma, 1, 1, 0,
  407:        doc: /* Return the log gamma of ARG.  */)
  408:      (arg)
  409:      register Lisp_Object arg;
  410: {
  411:   double d = extract_float (arg);
  412:   IN_FLOAT (d = lgamma (d), "log-gamma", arg);
  413:   return make_float (d);
  414: }
  415: 
  416: DEFUN ("cube-root", Fcube_root, Scube_root, 1, 1, 0,
  417:        doc: /* Return the cube root of ARG.  */)
  418:      (arg)
  419:      register Lisp_Object arg;
  420: {
  421:   double d = extract_float (arg);
  422: #ifdef HAVE_CBRT
  423:   IN_FLOAT (d = cbrt (d), "cube-root", arg);
  424: #else
  425:   if (d >= 0.0)
  426:     IN_FLOAT (d = pow (d, 1.0/3.0), "cube-root", arg);
  427:   else
  428:     IN_FLOAT (d = -pow (-d, 1.0/3.0), "cube-root", arg);
  429: #endif
  430:   return make_float (d);
  431: }
  432: 
  433: #endif
  434: ^L
  435: DEFUN ("exp", Fexp, Sexp, 1, 1, 0,
  436:        doc: /* Return the exponential base e of ARG.  */)
  437:      (arg)
  438:      register Lisp_Object arg;
  439: {
  440:   double d = extract_float (arg);
  441: #ifdef FLOAT_CHECK_DOMAIN
  442:   if (d > 709.7827)   /* Assume IEEE doubles here */
  443:     range_error ("exp", arg);
  444:   else if (d < -709.0)
  445:     return make_float (0.0);
  446:   else
  447: #endif
  448:     IN_FLOAT (d = exp (d), "exp", arg);
  449:   return make_float (d);
  450: }
  451: 
  452: DEFUN ("expt", Fexpt, Sexpt, 2, 2, 0,
  453:        doc: /* Return the exponential ARG1 ** ARG2.  */)
  454:      (arg1, arg2)
  455:      register Lisp_Object arg1, arg2;
  456: {
  457:   double f1, f2;
  458: 
  459:   CHECK_NUMBER_OR_FLOAT (arg1);
  460:   CHECK_NUMBER_OR_FLOAT (arg2);
  461:   if (INTEGERP (arg1)     /* common lisp spec */
  462:       && INTEGERP (arg2)   /* don't promote, if both are ints, and */
  463:       && 0 <= XINT (arg2)) /* we are sure the result is not fractional */
  464:     {                           /* this can be improved by pre-calculating */
  465:       EMACS_INT acc, x, y;      /* some binary powers of x then accumulating */
  466:       Lisp_Object val;
  467: 
  468:       x = XINT (arg1);
  469:       y = XINT (arg2);
  470:       acc = 1;
  471: 
  472:       if (y < 0)
  473:         {
  474:           if (x == 1)
  475:             acc = 1;
  476:           else if (x == -1)
  477:             acc = (y & 1) ? -1 : 1;
  478:           else
  479:             acc = 0;
  480:         }
  481:       else
  482:         {
  483:           while (y > 0)
  484:             {
  485:               if (y & 1)
  486:                 acc *= x;
  487:               x *= x;
  488:               y = (unsigned)y >> 1;
  489:             }
  490:         }
  491:       XSETINT (val, acc);
  492:       return val;
  493:     }
  494:   f1 = FLOATP (arg1) ? XFLOAT_DATA (arg1) : XINT (arg1);
  495:   f2 = FLOATP (arg2) ? XFLOAT_DATA (arg2) : XINT (arg2);
  496:   /* Really should check for overflow, too */
  497:   if (f1 == 0.0 && f2 == 0.0)
  498:     f1 = 1.0;
  499: #ifdef FLOAT_CHECK_DOMAIN
  500:   else if ((f1 == 0.0 && f2 < 0.0) || (f1 < 0 && f2 != floor(f2)))
  501:     domain_error2 ("expt", arg1, arg2);
  502: #endif
  503:   IN_FLOAT2 (f1 = pow (f1, f2), "expt", arg1, arg2);
  504:   return make_float (f1);
  505: }
  506: 
  507: DEFUN ("log", Flog, Slog, 1, 2, 0,
  508:        doc: /* Return the natural logarithm of ARG.
  509: If the optional argument BASE is given, return log ARG using that base.  */)
  510:      (arg, base)
  511:      register Lisp_Object arg, base;
  512: {
  513:   double d = extract_float (arg);
  514: 
  515: #ifdef FLOAT_CHECK_DOMAIN
  516:   if (d <= 0.0)
  517:     domain_error2 ("log", arg, base);
  518: #endif
  519:   if (NILP (base))
  520:     IN_FLOAT (d = log (d), "log", arg);
  521:   else
  522:     {
  523:       double b = extract_float (base);
  524: 
  525: #ifdef FLOAT_CHECK_DOMAIN
  526:       if (b <= 0.0 || b == 1.0)
  527:         domain_error2 ("log", arg, base);
  528: #endif
  529:       if (b == 10.0)
  530:         IN_FLOAT2 (d = log10 (d), "log", arg, base);
  531:       else
  532:         IN_FLOAT2 (d = log (d) / log (b), "log", arg, base);
  533:     }
  534:   return make_float (d);
  535: }
  536: 
  537: DEFUN ("log10", Flog10, Slog10, 1, 1, 0,
  538:        doc: /* Return the logarithm base 10 of ARG.  */)
  539:      (arg)
  540:      register Lisp_Object arg;
  541: {
  542:   double d = extract_float (arg);
  543: #ifdef FLOAT_CHECK_DOMAIN
  544:   if (d <= 0.0)
  545:     domain_error ("log10", arg);
  546: #endif
  547:   IN_FLOAT (d = log10 (d), "log10", arg);
  548:   return make_float (d);
  549: }
  550: 
  551: DEFUN ("sqrt",