floatfns.c 28.2 KB
Newer Older
Mike Rowan's avatar
Mike Rowan committed
1
/* Primitive operations on floating point for GNU Emacs Lisp interpreter.
2
   Copyright (C) 1988, 1993, 1994, 1999, 2001, 2002, 2003, 2004,
Glenn Morris's avatar
Glenn Morris committed
3
                 2005, 2006, 2007, 2008, 2009, 2010  Free Software Foundation, Inc.
Mike Rowan's avatar
Mike Rowan committed
4

5 6 7
Author: Wolfgang Rupprecht
(according to ack.texi)

Mike Rowan's avatar
Mike Rowan committed
8 9
This file is part of GNU Emacs.

10
GNU Emacs is free software: you can redistribute it and/or modify
Mike Rowan's avatar
Mike Rowan committed
11
it under the terms of the GNU General Public License as published by
12 13
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
Mike Rowan's avatar
Mike Rowan committed
14 15 16 17 18 19 20

GNU Emacs is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
21
along with GNU Emacs.  If not, see <http://www.gnu.org/licenses/>.  */
Mike Rowan's avatar
Mike Rowan committed
22 23


24 25 26 27 28 29
/* 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.
30
   Define HAVE_RINT if you have a working rint.
31 32 33 34 35 36 37 38 39 40 41 42
   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
Karl Heuer's avatar
Karl Heuer committed
43
   either setting errno, or signaling SIGFPE/SIGILL.  Otherwise, domain and
44 45 46 47 48
   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.)
 */

49
#include <config.h>
50
#include <signal.h>
51
#include <setjmp.h>
52 53 54
#include "lisp.h"
#include "syssignal.h"

55 56 57 58
#if STDC_HEADERS
#include <float.h>
#endif

59 60 61 62 63 64 65 66 67 68
/* If IEEE_FLOATING_POINT isn't defined, default it from FLT_*. */
#ifndef IEEE_FLOATING_POINT
#if (FLT_RADIX == 2 && FLT_MANT_DIG == 24 \
     && FLT_MIN_EXP == -125 && FLT_MAX_EXP == 128)
#define IEEE_FLOATING_POINT 1
#else
#define IEEE_FLOATING_POINT 0
#endif
#endif

Mike Rowan's avatar
Mike Rowan committed
69
#include <math.h>
70

71
/* This declaration is omitted on some systems, like Ultrix.  */
72
#if !defined (HPUX) && defined (HAVE_LOGB) && !defined (logb)
73
extern double logb ();
74
#endif /* not HPUX and HAVE_LOGB and no logb macro */
75

76 77 78 79 80 81 82
#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

Richard M. Stallman's avatar
Richard M. Stallman committed
83
#ifdef NO_MATHERR
84 85 86
#undef HAVE_MATHERR
#endif

87 88 89 90 91 92 93 94 95 96 97 98 99 100 101
#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>
102
#endif
Jim Blandy's avatar
Jim Blandy committed
103

104
#ifdef FLOAT_CATCH_SIGILL
Jim Blandy's avatar
Jim Blandy committed
105
static SIGTYPE float_error ();
106
#endif
Mike Rowan's avatar
Mike Rowan committed
107 108 109 110 111 112 113

/* 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,
Richard M. Stallman's avatar
Richard M. Stallman committed
114 115 116
   here is the actual argument value to use in the error message.
   These variables are used only across the floating point library call
   so there is no need to staticpro them.  */
Mike Rowan's avatar
Mike Rowan committed
117

118 119 120
static Lisp_Object float_error_arg, float_error_arg2;

static char *float_error_fn_name;
Mike Rowan's avatar
Mike Rowan committed
121

Jim Blandy's avatar
Jim Blandy committed
122 123 124
/* Evaluate the floating point expression D, recording NUM
   as the original argument for error messages.
   D is normally an assignment expression.
125 126 127 128 129
   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.  */
Jim Blandy's avatar
Jim Blandy committed
130

131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157
#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
158
#define IN_FLOAT(d, name, num) (in_float = 1, (d), in_float = 0)
159 160 161
#define IN_FLOAT2(d, name, num, num2) (in_float = 1, (d), in_float = 0)
#endif

162 163 164 165 166
/* Convert float to Lisp_Int if it fits, else signal a range error
   using the given arguments.  */
#define FLOAT_TO_INT(x, i, name, num)					\
  do									\
    {									\
167
      if (FIXNUM_OVERFLOW_P (x))					\
168
	range_error (name, num);					\
169
      XSETINT (i,  (EMACS_INT)(x));					\
170 171 172 173 174
    }									\
  while (0)
#define FLOAT_TO_INT2(x, i, name, num1, num2)				\
  do									\
    {									\
175
      if (FIXNUM_OVERFLOW_P (x))					\
176
	range_error2 (name, num1, num2);				\
177
      XSETINT (i,  (EMACS_INT)(x));					\
178 179 180
    }									\
  while (0)

181
#define arith_error(op,arg) \
182
  xsignal2 (Qarith_error, build_string ((op)), (arg))
183
#define range_error(op,arg) \
184
  xsignal2 (Qrange_error, build_string ((op)), (arg))
185
#define range_error2(op,a1,a2) \
186
  xsignal3 (Qrange_error, build_string ((op)), (a1), (a2))
187
#define domain_error(op,arg) \
188
  xsignal2 (Qdomain_error, build_string ((op)), (arg))
189
#define domain_error2(op,a1,a2) \
190
  xsignal3 (Qdomain_error, build_string ((op)), (a1), (a2))
Mike Rowan's avatar
Mike Rowan committed
191 192 193 194 195 196 197

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

double
extract_float (num)
     Lisp_Object num;
{
198
  CHECK_NUMBER_OR_FLOAT (num);
Mike Rowan's avatar
Mike Rowan committed
199

200
  if (FLOATP (num))
201
    return XFLOAT_DATA (num);
Mike Rowan's avatar
Mike Rowan committed
202 203
  return (double) XINT (num);
}
204 205

/* Trig functions.  */
Mike Rowan's avatar
Mike Rowan committed
206 207

DEFUN ("acos", Facos, Sacos, 1, 1, 0,
208 209
       doc: /* Return the inverse cosine of ARG.  */)
     (arg)
210
     register Lisp_Object arg;
Mike Rowan's avatar
Mike Rowan committed
211
{
212 213 214 215 216 217
  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);
Mike Rowan's avatar
Mike Rowan committed
218 219 220
  return make_float (d);
}

221
DEFUN ("asin", Fasin, Sasin, 1, 1, 0,
222 223
       doc: /* Return the inverse sine of ARG.  */)
     (arg)
224
     register Lisp_Object arg;
Mike Rowan's avatar
Mike Rowan committed
225
{
226 227 228 229 230 231
  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);
Mike Rowan's avatar
Mike Rowan committed
232 233 234
  return make_float (d);
}

235 236 237 238 239 240 241 242
DEFUN ("atan", Fatan, Satan, 1, 2, 0,
       doc: /* Return the inverse tangent of the arguments.
If only one argument Y is given, return the inverse tangent of Y.
If two arguments Y and X are given, return the inverse tangent of Y
divided by X, i.e. the angle in radians between the vector (X, Y)
and the x-axis.  */)
     (y, x)
     register Lisp_Object y, x;
Mike Rowan's avatar
Mike Rowan committed
243
{
244 245 246 247 248 249 250 251 252 253
  double d = extract_float (y);

  if (NILP (x))
    IN_FLOAT (d = atan (d), "atan", y);
  else
    {
      double d2 = extract_float (x);

      IN_FLOAT2 (d = atan2 (d, d2), "atan", y, x);
    }
Mike Rowan's avatar
Mike Rowan committed
254 255 256
  return make_float (d);
}

257
DEFUN ("cos", Fcos, Scos, 1, 1, 0,
258 259
       doc: /* Return the cosine of ARG.  */)
     (arg)
260
     register Lisp_Object arg;
Mike Rowan's avatar
Mike Rowan committed
261
{
262 263
  double d = extract_float (arg);
  IN_FLOAT (d = cos (d), "cos", arg);
Mike Rowan's avatar
Mike Rowan committed
264 265 266
  return make_float (d);
}

267
DEFUN ("sin", Fsin, Ssin, 1, 1, 0,
268 269
       doc: /* Return the sine of ARG.  */)
     (arg)
270
     register Lisp_Object arg;
Mike Rowan's avatar
Mike Rowan committed
271
{
272 273
  double d = extract_float (arg);
  IN_FLOAT (d = sin (d), "sin", arg);
Mike Rowan's avatar
Mike Rowan committed
274 275 276
  return make_float (d);
}

277
DEFUN ("tan", Ftan, Stan, 1, 1, 0,
278 279
       doc: /* Return the tangent of ARG.  */)
     (arg)
280 281 282 283 284 285 286 287 288
     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);
Mike Rowan's avatar
Mike Rowan committed
289 290
  return make_float (d);
}
291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354

#if defined HAVE_ISNAN && defined HAVE_COPYSIGN
DEFUN ("isnan", Fisnan, Sisnan, 1, 1, 0,
       doc: /* Return non nil iff argument X is a NaN.  */)
     (x)
     Lisp_Object x;
{
  CHECK_FLOAT (x);
  return isnan (XFLOAT_DATA (x)) ? Qt : Qnil;
}

DEFUN ("copysign", Fcopysign, Scopysign, 1, 2, 0,
       doc: /* Copy sign of X2 to value of X1, and return the result.
Cause an error if X1 or X2 is not a float.  */)
     (x1, x2)
     Lisp_Object x1, x2;
{
  double f1, f2;

  CHECK_FLOAT (x1);
  CHECK_FLOAT (x2);

  f1 = XFLOAT_DATA (x1);
  f2 = XFLOAT_DATA (x2);

  return make_float (copysign (f1, f2));
}

DEFUN ("frexp", Ffrexp, Sfrexp, 1, 1, 0,
       doc: /* Get significand and exponent of a floating point number.
Breaks the floating point number X into its binary significand SGNFCAND
\(a floating point value between 0.5 (included) and 1.0 (excluded))
and an integral exponent EXP for 2, such that:

  X = SGNFCAND * 2^EXP

The function returns the cons cell (SGNFCAND . EXP).
If X is zero, both parts (SGNFCAND and EXP) are zero.  */)
     (x)
     Lisp_Object x;
{
  double f = XFLOATINT (x);

  if (f == 0.0)
    return Fcons (make_float (0.0), make_number (0));
  else
    {
      int    exp;
      double sgnfcand = frexp (f, &exp);
      return Fcons (make_float (sgnfcand), make_number (exp));
    }
}

DEFUN ("ldexp", Fldexp, Sldexp, 1, 2, 0,
       doc: /* Construct number X from significand SGNFCAND and exponent EXP.
Returns the floating point value resulting from multiplying SGNFCAND
(the significand) by 2 raised to the power of EXP (the exponent).   */)
     (sgnfcand, exp)
     Lisp_Object sgnfcand, exp;
{
  CHECK_NUMBER (exp);
  return make_float (ldexp (XFLOATINT (sgnfcand), XINT (exp)));
}
#endif
Mike Rowan's avatar
Mike Rowan committed
355

356 357
#if 0 /* Leave these out unless we find there's a reason for them.  */

Mike Rowan's avatar
Mike Rowan committed
358
DEFUN ("bessel-j0", Fbessel_j0, Sbessel_j0, 1, 1, 0,
359 360
       doc: /* Return the bessel function j0 of ARG.  */)
     (arg)
361
     register Lisp_Object arg;
Mike Rowan's avatar
Mike Rowan committed
362
{
363 364
  double d = extract_float (arg);
  IN_FLOAT (d = j0 (d), "bessel-j0", arg);
Mike Rowan's avatar
Mike Rowan committed
365 366 367 368
  return make_float (d);
}

DEFUN ("bessel-j1", Fbessel_j1, Sbessel_j1, 1, 1, 0,
369 370
       doc: /* Return the bessel function j1 of ARG.  */)
     (arg)
371
     register Lisp_Object arg;
Mike Rowan's avatar
Mike Rowan committed
372
{
373 374
  double d = extract_float (arg);
  IN_FLOAT (d = j1 (d), "bessel-j1", arg);
Mike Rowan's avatar
Mike Rowan committed
375 376 377 378
  return make_float (d);
}

DEFUN ("bessel-jn", Fbessel_jn, Sbessel_jn, 2, 2, 0,
379 380 381
       doc: /* Return the order N bessel function output jn of ARG.
The first arg (the order) is truncated to an integer.  */)
     (n, arg)
382
     register Lisp_Object n, arg;
Mike Rowan's avatar
Mike Rowan committed
383
{
384 385
  int i1 = extract_float (n);
  double f2 = extract_float (arg);
Mike Rowan's avatar
Mike Rowan committed
386

387
  IN_FLOAT (f2 = jn (i1, f2), "bessel-jn", n);
Mike Rowan's avatar
Mike Rowan committed
388 389 390 391
  return make_float (f2);
}

DEFUN ("bessel-y0", Fbessel_y0, Sbessel_y0, 1, 1, 0,
392 393
       doc: /* Return the bessel function y0 of ARG.  */)
     (arg)
394
     register Lisp_Object arg;
Mike Rowan's avatar
Mike Rowan committed
395
{
396 397
  double d = extract_float (arg);
  IN_FLOAT (d = y0 (d), "bessel-y0", arg);
Mike Rowan's avatar
Mike Rowan committed
398 399 400 401
  return make_float (d);
}

DEFUN ("bessel-y1", Fbessel_y1, Sbessel_y1, 1, 1, 0,
402 403
       doc: /* Return the bessel function y1 of ARG.  */)
     (arg)
404
     register Lisp_Object arg;
Mike Rowan's avatar
Mike Rowan committed
405
{
406 407
  double d = extract_float (arg);
  IN_FLOAT (d = y1 (d), "bessel-y0", arg);
Mike Rowan's avatar
Mike Rowan committed
408 409 410 411
  return make_float (d);
}

DEFUN ("bessel-yn", Fbessel_yn, Sbessel_yn, 2, 2, 0,
412 413 414
       doc: /* Return the order N bessel function output yn of ARG.
The first arg (the order) is truncated to an integer.  */)
     (n, arg)
415
     register Lisp_Object n, arg;
Mike Rowan's avatar
Mike Rowan committed
416
{
417 418
  int i1 = extract_float (n);
  double f2 = extract_float (arg);
Mike Rowan's avatar
Mike Rowan committed
419

420
  IN_FLOAT (f2 = yn (i1, f2), "bessel-yn", n);
Mike Rowan's avatar
Mike Rowan committed
421 422 423
  return make_float (f2);
}

424 425 426
#endif

#if 0 /* Leave these out unless we see they are worth having.  */
Mike Rowan's avatar
Mike Rowan committed
427 428

DEFUN ("erf", Ferf, Serf, 1, 1, 0,
429 430
       doc: /* Return the mathematical error function of ARG.  */)
     (arg)
431
     register Lisp_Object arg;
Mike Rowan's avatar
Mike Rowan committed
432
{
433 434
  double d = extract_float (arg);
  IN_FLOAT (d = erf (d), "erf", arg);
Mike Rowan's avatar
Mike Rowan committed
435 436 437 438
  return make_float (d);
}

DEFUN ("erfc", Ferfc, Serfc, 1, 1, 0,
439 440
       doc: /* Return the complementary error function of ARG.  */)
     (arg)
441
     register Lisp_Object arg;
Mike Rowan's avatar
Mike Rowan committed
442
{
443 444
  double d = extract_float (arg);
  IN_FLOAT (d = erfc (d), "erfc", arg);
Mike Rowan's avatar
Mike Rowan committed
445 446 447 448
  return make_float (d);
}

DEFUN ("log-gamma", Flog_gamma, Slog_gamma, 1, 1, 0,
449 450
       doc: /* Return the log gamma of ARG.  */)
     (arg)
451
     register Lisp_Object arg;
Mike Rowan's avatar
Mike Rowan committed
452
{
453 454
  double d = extract_float (arg);
  IN_FLOAT (d = lgamma (d), "log-gamma", arg);
Mike Rowan's avatar
Mike Rowan committed
455 456 457
  return make_float (d);
}

458
DEFUN ("cube-root", Fcube_root, Scube_root, 1, 1, 0,
459 460
       doc: /* Return the cube root of ARG.  */)
     (arg)
461
     register Lisp_Object arg;
Mike Rowan's avatar
Mike Rowan committed
462
{
463 464 465 466 467 468 469 470 471
  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
Mike Rowan's avatar
Mike Rowan committed
472 473 474
  return make_float (d);
}

Richard M. Stallman's avatar
Richard M. Stallman committed
475 476
#endif

477
DEFUN ("exp", Fexp, Sexp, 1, 1, 0,
478 479
       doc: /* Return the exponential base e of ARG.  */)
     (arg)
480 481 482 483 484 485 486 487 488 489 490
     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);
Mike Rowan's avatar
Mike Rowan committed
491 492 493 494
  return make_float (d);
}

DEFUN ("expt", Fexpt, Sexpt, 2, 2, 0,
495 496
       doc: /* Return the exponential ARG1 ** ARG2.  */)
     (arg1, arg2)
497
     register Lisp_Object arg1, arg2;
Mike Rowan's avatar
Mike Rowan committed
498
{
499
  double f1, f2, f3;
Mike Rowan's avatar
Mike Rowan committed
500

501 502
  CHECK_NUMBER_OR_FLOAT (arg1);
  CHECK_NUMBER_OR_FLOAT (arg2);
503
  if (INTEGERP (arg1)     /* common lisp spec */
504 505
      && INTEGERP (arg2)   /* don't promote, if both are ints, and */
      && 0 <= XINT (arg2)) /* we are sure the result is not fractional */
Mike Rowan's avatar
Mike Rowan committed
506
    {				/* this can be improved by pre-calculating */
507
      EMACS_INT acc, x, y;	/* some binary powers of x then accumulating */
508 509
      Lisp_Object val;

510 511
      x = XINT (arg1);
      y = XINT (arg2);
Mike Rowan's avatar
Mike Rowan committed
512
      acc = 1;
513

Mike Rowan's avatar
Mike Rowan committed
514 515
      if (y < 0)
	{
516 517 518 519 520 521
	  if (x == 1)
	    acc = 1;
	  else if (x == -1)
	    acc = (y & 1) ? -1 : 1;
	  else
	    acc = 0;
Mike Rowan's avatar
Mike Rowan committed
522 523 524
	}
      else
	{
525 526 527 528 529 530 531
	  while (y > 0)
	    {
	      if (y & 1)
		acc *= x;
	      x *= x;
	      y = (unsigned)y >> 1;
	    }
Mike Rowan's avatar
Mike Rowan committed
532
	}
533
      XSETINT (val, acc);
534
      return val;
Mike Rowan's avatar
Mike Rowan committed
535
    }
536 537
  f1 = FLOATP (arg1) ? XFLOAT_DATA (arg1) : XINT (arg1);
  f2 = FLOATP (arg2) ? XFLOAT_DATA (arg2) : XINT (arg2);
538 539 540 541 542 543 544
  /* 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
545 546 547 548 549
  IN_FLOAT2 (f3 = pow (f1, f2), "expt", arg1, arg2);
  /* Check for overflow in the result.  */
  if (f1 != 0.0 && f3 == 0.0)
    range_error ("expt", arg1);
  return make_float (f3);
Mike Rowan's avatar
Mike Rowan committed
550
}
551

552
DEFUN ("log", Flog, Slog, 1, 2, 0,
553
       doc: /* Return the natural logarithm of ARG.
Richard M. Stallman's avatar
Richard M. Stallman committed
554
If the optional argument BASE is given, return log ARG using that base.  */)
555
     (arg, base)
556
     register Lisp_Object arg, base;
Mike Rowan's avatar
Mike Rowan committed
557
{
558
  double d = extract_float (arg);
559

560 561 562 563
#ifdef FLOAT_CHECK_DOMAIN
  if (d <= 0.0)
    domain_error2 ("log", arg, base);
#endif
564
  if (NILP (base))
565
    IN_FLOAT (d = log (d), "log", arg);
566 567 568 569
  else
    {
      double b = extract_float (base);

570 571 572 573 574 575 576
#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
577
	IN_FLOAT2 (d = log (d) / log (b), "log", arg, base);
578
    }
Mike Rowan's avatar
Mike Rowan committed
579 580 581
  return make_float (d);
}

582
DEFUN ("log10", Flog10, Slog10, 1, 1, 0,
583 584
       doc: /* Return the logarithm base 10 of ARG.  */)
     (arg)
585
     register Lisp_Object arg;
Mike Rowan's avatar
Mike Rowan committed
586
{
587 588 589 590 591 592
  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);
593 594 595
  return make_float (d);
}

Mike Rowan's avatar
Mike Rowan committed
596
DEFUN ("sqrt", Fsqrt, Ssqrt, 1, 1, 0,
597 598
       doc: /* Return the square root of ARG.  */)
     (arg)
599
     register Lisp_Object arg;
Mike Rowan's avatar
Mike Rowan committed
600
{
601 602 603 604 605 606
  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);
Mike Rowan's avatar
Mike Rowan committed
607 608
  return make_float (d);
}
609

Richard M. Stallman's avatar
Richard M. Stallman committed
610
#if 0 /* Not clearly worth adding.  */
Mike Rowan's avatar
Mike Rowan committed
611

612
DEFUN ("acosh", Facosh, Sacosh, 1, 1, 0,
613 614
       doc: /* Return the inverse hyperbolic cosine of ARG.  */)
     (arg)
615
     register Lisp_Object arg;
Mike Rowan's avatar
Mike Rowan committed
616
{
617 618 619 620 621 622 623 624 625 626
  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
627 628 629 630
  return make_float (d);
}

DEFUN ("asinh", Fasinh, Sasinh, 1, 1, 0,
631 632
       doc: /* Return the inverse hyperbolic sine of ARG.  */)
     (arg)
633
     register Lisp_Object arg;
634
{
635 636 637 638 639 640
  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
641 642 643 644
  return make_float (d);
}

DEFUN ("atanh", Fatanh, Satanh, 1, 1, 0,
645 646
       doc: /* Return the inverse hyperbolic tangent of ARG.  */)
     (arg)
647
     register Lisp_Object arg;
648
{
649 650 651 652 653 654 655 656 657 658
  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
659 660 661 662
  return make_float (d);
}

DEFUN ("cosh", Fcosh, Scosh, 1, 1, 0,
663 664
       doc: /* Return the hyperbolic cosine of ARG.  */)
     (arg)
665
     register Lisp_Object arg;
666
{
667 668 669 670 671 672
  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);
673 674 675 676
  return make_float (d);
}

DEFUN ("sinh", Fsinh, Ssinh, 1, 1, 0,
677 678
       doc: /* Return the hyperbolic sine of ARG.  */)
     (arg)
679
     register Lisp_Object arg;
680
{
681 682 683 684 685 686
  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);
Mike Rowan's avatar
Mike Rowan committed
687 688 689 690
  return make_float (d);
}

DEFUN ("tanh", Ftanh, Stanh, 1, 1, 0,
691 692
       doc: /* Return the hyperbolic tangent of ARG.  */)
     (arg)
693
     register Lisp_Object arg;
Mike Rowan's avatar
Mike Rowan committed
694
{
695 696
  double d = extract_float (arg);
  IN_FLOAT (d = tanh (d), "tanh", arg);
Mike Rowan's avatar
Mike Rowan committed
697 698
  return make_float (d);
}
699
#endif
Mike Rowan's avatar
Mike Rowan committed
700 701

DEFUN ("abs", Fabs, Sabs, 1, 1, 0,
702 703
       doc: /* Return the absolute value of ARG.  */)
     (arg)
704
     register Lisp_Object arg;
Mike Rowan's avatar
Mike Rowan committed
705
{
706
  CHECK_NUMBER_OR_FLOAT (arg);
Mike Rowan's avatar
Mike Rowan committed
707

708
  if (FLOATP (arg))
709
    IN_FLOAT (arg = make_float (fabs (XFLOAT_DATA (arg))), "abs", arg);
710
  else if (XINT (arg) < 0)
711
    XSETINT (arg, - XINT (arg));
Mike Rowan's avatar
Mike Rowan committed
712

713
  return arg;
Mike Rowan's avatar
Mike Rowan committed
714 715 716
}

DEFUN ("float", Ffloat, Sfloat, 1, 1, 0,
717 718
       doc: /* Return the floating point number equal to ARG.  */)
     (arg)
719
     register Lisp_Object arg;
Mike Rowan's avatar
Mike Rowan committed
720
{
721
  CHECK_NUMBER_OR_FLOAT (arg);
Mike Rowan's avatar
Mike Rowan committed
722

723
  if (INTEGERP (arg))
724
    return make_float ((double) XINT (arg));
Mike Rowan's avatar
Mike Rowan committed
725
  else				/* give 'em the same float back */
726
    return arg;
Mike Rowan's avatar
Mike Rowan committed
727 728 729
}

DEFUN ("logb", Flogb, Slogb, 1, 1, 0,
730 731
       doc: /* Returns largest integer <= the base 2 log of the magnitude of ARG.
This is the same as the exponent of a float.  */)
732 733
     (arg)
     Lisp_Object arg;
Mike Rowan's avatar
Mike Rowan committed
734
{
735
  Lisp_Object val;
736
  EMACS_INT value;
737
  double f = extract_float (arg);
738

739
  if (f == 0.0)
Stefan Monnier's avatar
Stefan Monnier committed
740
    value = MOST_NEGATIVE_FIXNUM;
741 742
  else
    {
743
#ifdef HAVE_LOGB
744
      IN_FLOAT (value = logb (f), "logb", arg);
745 746
#else
#ifdef HAVE_FREXP
747 748 749
      int ivalue;
      IN_FLOAT (frexp (f, &ivalue), "logb", arg);
      value = ivalue - 1;
750
#else
751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769
      int i;
      double d;
      if (f < 0.0)
	f = -f;
      value = -1;
      while (f < 0.5)
	{
	  for (i = 1, d = 0.5; d * d >= f; i += i)
	    d *= d;
	  f /= d;
	  value -= i;
	}
      while (f >= 1.0)
	{
	  for (i = 1, d = 2.0; d * d <= f; i += i)
	    d *= d;
	  f /= d;
	  value += i;
	}
770
#endif
771
#endif
772
    }
773
  XSETINT (val, value);
774
  return val;
Mike Rowan's avatar
Mike Rowan committed
775 776
}

777

778 779 780 781
/* the rounding functions  */

static Lisp_Object
rounding_driver (arg, divisor, double_round, int_round2, name)
782
     register Lisp_Object arg, divisor;
783 784 785
     double (*double_round) ();
     EMACS_INT (*int_round2) ();
     char *name;
Mike Rowan's avatar
Mike Rowan committed
786
{
787
  CHECK_NUMBER_OR_FLOAT (arg);
Mike Rowan's avatar
Mike Rowan committed
788

789 790
  if (! NILP (divisor))
    {
791
      EMACS_INT i1, i2;
792

793
      CHECK_NUMBER_OR_FLOAT (divisor);
794

795
      if (FLOATP (arg) || FLOATP (divisor))
796 797 798
	{
	  double f1, f2;

799 800
	  f1 = FLOATP (arg) ? XFLOAT_DATA (arg) : XINT (arg);
	  f2 = (FLOATP (divisor) ? XFLOAT_DATA (divisor) : XINT (divisor));
801
	  if (! IEEE_FLOATING_POINT && f2 == 0)
802
	    xsignal0 (Qarith_error);
803

804 805
	  IN_FLOAT2 (f1 = (*double_round) (f1 / f2), name, arg, divisor);
	  FLOAT_TO_INT2 (f1, arg, name, arg, divisor);
806 807 808 809 810 811 812
	  return arg;
	}

      i1 = XINT (arg);
      i2 = XINT (divisor);

      if (i2 == 0)
813
	xsignal0 (Qarith_error);
814

815
      XSETINT (arg, (*int_round2) (i1, i2));
816 817 818
      return arg;
    }

819
  if (FLOATP (arg))
820 821
    {
      double d;
822

823
      IN_FLOAT (d = (*double_round) (XFLOAT_DATA (arg)), name, arg);
824
      FLOAT_TO_INT (d, arg, name, arg);
825
    }
Mike Rowan's avatar
Mike Rowan committed
826

827
  return arg;
Mike Rowan's avatar
Mike Rowan committed
828 829
}

830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876
/* With C's /, the result is implementation-defined if either operand
   is negative, so take care with negative operands in the following
   integer functions.  */

static EMACS_INT
ceiling2 (i1, i2)
     EMACS_INT i1, i2;
{
  return (i2 < 0
	  ? (i1 < 0  ?  ((-1 - i1) / -i2) + 1  :  - (i1 / -i2))
	  : (i1 <= 0  ?  - (-i1 / i2)  :  ((i1 - 1) / i2) + 1));
}

static EMACS_INT
floor2 (i1, i2)
     EMACS_INT i1, i2;
{
  return (i2 < 0
	  ? (i1 <= 0  ?  -i1 / -i2  :  -1 - ((i1 - 1) / -i2))
	  : (i1 < 0  ?  -1 - ((-1 - i1) / i2)  :  i1 / i2));
}

static EMACS_INT
truncate2 (i1, i2)
     EMACS_INT i1, i2;
{
  return (i2 < 0
	  ? (i1 < 0  ?  -i1 / -i2  :  - (i1 / -i2))
	  : (i1 < 0  ?  - (-i1 / i2)  :  i1 / i2));
}

static EMACS_INT
round2 (i1, i2)
     EMACS_INT i1, i2;
{
  /* The C language's division operator gives us one remainder R, but
     we want the remainder R1 on the other side of 0 if R1 is closer
     to 0 than R is; because we want to round to even, we also want R1
     if R and R1 are the same distance from 0 and if C's quotient is
     odd.  */
  EMACS_INT q = i1 / i2;
  EMACS_INT r = i1 % i2;
  EMACS_INT abs_r = r < 0 ? -r : r;
  EMACS_INT abs_r1 = (i2 < 0 ? -i2 : i2) - abs_r;
  return q + (abs_r + (q & 1) <= abs_r1 ? 0 : (i2 ^ r) < 0 ? -1 : 1);
}

877 878 879 880 881
/* The code uses emacs_rint, so that it works to undefine HAVE_RINT
   if `rint' exists but does not work right.  */
#ifdef HAVE_RINT
#define emacs_rint rint
#else
882
static double
883
emacs_rint (d)
884 885
     double d;
{
Richard M. Stallman's avatar
Richard M. Stallman committed
886
  return floor (d + 0.5);
887 888 889
}
#endif

890 891 892 893 894 895 896 897
static double
double_identity (d)
     double d;
{
  return d;
}

DEFUN ("ceiling", Fceiling, Sceiling, 1, 2, 0,
898 899
       doc: /* Return the smallest integer no less than ARG.
This rounds the value towards +inf.
900 901
With optional DIVISOR, return the smallest integer no less than ARG/DIVISOR.  */)
     (arg, divisor)
902 903 904 905 906 907
     Lisp_Object arg, divisor;
{
  return rounding_driver (arg, divisor, ceil, ceiling2, "ceiling");
}

DEFUN ("floor", Ffloor, Sfloor, 1, 2, 0,
908
       doc: /* Return the largest integer no greater than ARG.
Lute Kamstra's avatar
Lute Kamstra committed
909
This rounds the value towards -inf.
910 911
With optional DIVISOR, return the largest integer no greater than ARG/DIVISOR.  */)
     (arg, divisor)
912 913 914 915 916 917
     Lisp_Object arg, divisor;
{
  return rounding_driver (arg, divisor, floor, floor2, "floor");
}

DEFUN ("round", Fround, Sround, 1, 2, 0,
918
       doc: /* Return the nearest integer to ARG.
Eli Zaretskii's avatar
Eli Zaretskii committed
919 920
With optional DIVISOR, return the nearest integer to ARG/DIVISOR.

Eli Zaretskii's avatar
Eli Zaretskii committed
921 922 923
Rounding a value equidistant between two integers may choose the
integer closer to zero, or it may prefer an even integer, depending on
your machine.  For example, \(round 2.5\) can return 3 on some
Eli Zaretskii's avatar
Eli Zaretskii committed
924
systems, but 2 on others.  */)
925
     (arg, divisor)
926 927
     Lisp_Object arg, divisor;
{
928
  return rounding_driver (arg, divisor, emacs_rint, round2, "round");
929 930 931
}

DEFUN ("truncate", Ftruncate, Struncate, 1, 2, 0,
932 933 934 935
       doc: /* Truncate a floating point number to an int.
Rounds ARG toward zero.
With optional DIVISOR, truncate ARG/DIVISOR.  */)
     (arg, divisor)
936 937 938 939 940 941
     Lisp_Object arg, divisor;
{
  return rounding_driver (arg, divisor, double_identity, truncate2,
			  "truncate");
}

942

943
Lisp_Object
944
fmod_float (Lisp_Object x, Lisp_Object y)
945 946 947
{
  double f1, f2;

948 949
  f1 = FLOATP (x) ? XFLOAT_DATA (x) : XINT (x);
  f2 = FLOATP (y) ? XFLOAT_DATA (y) : XINT (y);
950 951

  if (! IEEE_FLOATING_POINT && f2 == 0)
952
    xsignal0 (Qarith_error);
953 954 955 956 957 958 959

  /* If the "remainder" comes out with the wrong sign, fix it.  */
  IN_FLOAT2 ((f1 = fmod (f1, f2),
	      f1 = (f2 < 0 ? f1 > 0 : f1 < 0) ? f1 + f2 : f1),
	     "mod", x, y);
  return make_float (f1);
}
960 961 962 963

/* It's not clear these are worth adding.  */

DEFUN ("fceiling", Ffceiling, Sfceiling, 1, 1, 0,
964 965 966
       doc: /* Return the smallest integer no less than ARG, as a float.
\(Round toward +inf.\)  */)
     (arg)
967 968 969 970 971 972 973 974
     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,