floatfns.c 26.4 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  Free Software Foundation, Inc.
Mike Rowan's avatar
Mike Rowan committed
4 5 6

This file is part of GNU Emacs.

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

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
18
along with GNU Emacs.  If not, see <http://www.gnu.org/licenses/>.  */
Mike Rowan's avatar
Mike Rowan committed
19 20


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

46
#include <config.h>
47
#include <signal.h>
48 49 50
#include "lisp.h"
#include "syssignal.h"

51 52 53 54
#if STDC_HEADERS
#include <float.h>
#endif

55 56 57 58 59 60 61 62 63 64
/* 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
65
#include <math.h>
66

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

72 73 74 75 76 77 78
#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
79
#ifdef NO_MATHERR
80 81 82
#undef HAVE_MATHERR
#endif

83 84 85 86 87 88 89 90 91 92 93 94 95 96 97
#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>
Jim Blandy's avatar
Jim Blandy committed
98

99
#ifndef USE_CRT_DLL
Jim Blandy's avatar
Jim Blandy committed
100
extern int errno;
101
#endif
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 291
  return make_float (d);
}

292 293
#if 0 /* Leave these out unless we find there's a reason for them.  */

Mike Rowan's avatar
Mike Rowan committed
294
DEFUN ("bessel-j0", Fbessel_j0, Sbessel_j0, 1, 1, 0,
295 296
       doc: /* Return the bessel function j0 of ARG.  */)
     (arg)
297
     register Lisp_Object arg;
Mike Rowan's avatar
Mike Rowan committed
298
{
299 300
  double d = extract_float (arg);
  IN_FLOAT (d = j0 (d), "bessel-j0", arg);
Mike Rowan's avatar
Mike Rowan committed
301 302 303 304
  return make_float (d);
}

DEFUN ("bessel-j1", Fbessel_j1, Sbessel_j1, 1, 1, 0,
305 306
       doc: /* Return the bessel function j1 of ARG.  */)
     (arg)
307
     register Lisp_Object arg;
Mike Rowan's avatar
Mike Rowan committed
308
{
309 310
  double d = extract_float (arg);
  IN_FLOAT (d = j1 (d), "bessel-j1", arg);
Mike Rowan's avatar
Mike Rowan committed
311 312 313 314
  return make_float (d);
}

DEFUN ("bessel-jn", Fbessel_jn, Sbessel_jn, 2, 2, 0,
315 316 317
       doc: /* Return the order N bessel function output jn of ARG.
The first arg (the order) is truncated to an integer.  */)
     (n, arg)
318
     register Lisp_Object n, arg;
Mike Rowan's avatar
Mike Rowan committed
319
{
320 321
  int i1 = extract_float (n);
  double f2 = extract_float (arg);
Mike Rowan's avatar
Mike Rowan committed
322

323
  IN_FLOAT (f2 = jn (i1, f2), "bessel-jn", n);
Mike Rowan's avatar
Mike Rowan committed
324 325 326 327
  return make_float (f2);
}

DEFUN ("bessel-y0", Fbessel_y0, Sbessel_y0, 1, 1, 0,
328 329
       doc: /* Return the bessel function y0 of ARG.  */)
     (arg)
330
     register Lisp_Object arg;
Mike Rowan's avatar
Mike Rowan committed
331
{
332 333
  double d = extract_float (arg);
  IN_FLOAT (d = y0 (d), "bessel-y0", arg);
Mike Rowan's avatar
Mike Rowan committed
334 335 336 337
  return make_float (d);
}

DEFUN ("bessel-y1", Fbessel_y1, Sbessel_y1, 1, 1, 0,
338 339
       doc: /* Return the bessel function y1 of ARG.  */)
     (arg)
340
     register Lisp_Object arg;
Mike Rowan's avatar
Mike Rowan committed
341
{
342 343
  double d = extract_float (arg);
  IN_FLOAT (d = y1 (d), "bessel-y0", arg);
Mike Rowan's avatar
Mike Rowan committed
344 345 346 347
  return make_float (d);
}

DEFUN ("bessel-yn", Fbessel_yn, Sbessel_yn, 2, 2, 0,
348 349 350
       doc: /* Return the order N bessel function output yn of ARG.
The first arg (the order) is truncated to an integer.  */)
     (n, arg)
351
     register Lisp_Object n, arg;
Mike Rowan's avatar
Mike Rowan committed
352
{
353 354
  int i1 = extract_float (n);
  double f2 = extract_float (arg);
Mike Rowan's avatar
Mike Rowan committed
355

356
  IN_FLOAT (f2 = yn (i1, f2), "bessel-yn", n);
Mike Rowan's avatar
Mike Rowan committed
357 358 359
  return make_float (f2);
}

360 361 362
#endif

#if 0 /* Leave these out unless we see they are worth having.  */
Mike Rowan's avatar
Mike Rowan committed
363 364

DEFUN ("erf", Ferf, Serf, 1, 1, 0,
365 366
       doc: /* Return the mathematical error function of ARG.  */)
     (arg)
367
     register Lisp_Object arg;
Mike Rowan's avatar
Mike Rowan committed
368
{
369 370
  double d = extract_float (arg);
  IN_FLOAT (d = erf (d), "erf", arg);
Mike Rowan's avatar
Mike Rowan committed
371 372 373 374
  return make_float (d);
}

DEFUN ("erfc", Ferfc, Serfc, 1, 1, 0,
375 376
       doc: /* Return the complementary error function of ARG.  */)
     (arg)
377
     register Lisp_Object arg;
Mike Rowan's avatar
Mike Rowan committed
378
{
379 380
  double d = extract_float (arg);
  IN_FLOAT (d = erfc (d), "erfc", arg);
Mike Rowan's avatar
Mike Rowan committed
381 382 383 384
  return make_float (d);
}

DEFUN ("log-gamma", Flog_gamma, Slog_gamma, 1, 1, 0,
385 386
       doc: /* Return the log gamma of ARG.  */)
     (arg)
387
     register Lisp_Object arg;
Mike Rowan's avatar
Mike Rowan committed
388
{
389 390
  double d = extract_float (arg);
  IN_FLOAT (d = lgamma (d), "log-gamma", arg);
Mike Rowan's avatar
Mike Rowan committed
391 392 393
  return make_float (d);
}

394
DEFUN ("cube-root", Fcube_root, Scube_root, 1, 1, 0,
395 396
       doc: /* Return the cube root of ARG.  */)
     (arg)
397
     register Lisp_Object arg;
Mike Rowan's avatar
Mike Rowan committed
398
{
399 400 401 402 403 404 405 406 407
  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
408 409 410
  return make_float (d);
}

Richard M. Stallman's avatar
Richard M. Stallman committed
411 412
#endif

413
DEFUN ("exp", Fexp, Sexp, 1, 1, 0,
414 415
       doc: /* Return the exponential base e of ARG.  */)
     (arg)
416 417 418 419 420 421 422 423 424 425 426
     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
427 428 429 430
  return make_float (d);
}

DEFUN ("expt", Fexpt, Sexpt, 2, 2, 0,
431 432
       doc: /* Return the exponential ARG1 ** ARG2.  */)
     (arg1, arg2)
433
     register Lisp_Object arg1, arg2;
Mike Rowan's avatar
Mike Rowan committed
434
{
435
  double f1, f2, f3;
Mike Rowan's avatar
Mike Rowan committed
436

437 438
  CHECK_NUMBER_OR_FLOAT (arg1);
  CHECK_NUMBER_OR_FLOAT (arg2);
439
  if (INTEGERP (arg1)     /* common lisp spec */
440 441
      && 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
442
    {				/* this can be improved by pre-calculating */
443
      EMACS_INT acc, x, y;	/* some binary powers of x then accumulating */
444 445
      Lisp_Object val;

446 447
      x = XINT (arg1);
      y = XINT (arg2);
Mike Rowan's avatar
Mike Rowan committed
448
      acc = 1;
449

Mike Rowan's avatar
Mike Rowan committed
450 451
      if (y < 0)
	{
452 453 454 455 456 457
	  if (x == 1)
	    acc = 1;
	  else if (x == -1)
	    acc = (y & 1) ? -1 : 1;
	  else
	    acc = 0;
Mike Rowan's avatar
Mike Rowan committed
458 459 460
	}
      else
	{
461 462 463 464 465 466 467
	  while (y > 0)
	    {
	      if (y & 1)
		acc *= x;
	      x *= x;
	      y = (unsigned)y >> 1;
	    }
Mike Rowan's avatar
Mike Rowan committed
468
	}
469
      XSETINT (val, acc);
470
      return val;
Mike Rowan's avatar
Mike Rowan committed
471
    }
472 473
  f1 = FLOATP (arg1) ? XFLOAT_DATA (arg1) : XINT (arg1);
  f2 = FLOATP (arg2) ? XFLOAT_DATA (arg2) : XINT (arg2);
474 475 476 477 478 479 480
  /* 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
481 482 483 484 485
  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
486
}
487

488
DEFUN ("log", Flog, Slog, 1, 2, 0,
489
       doc: /* Return the natural logarithm of ARG.
Richard M. Stallman's avatar
Richard M. Stallman committed
490
If the optional argument BASE is given, return log ARG using that base.  */)
491
     (arg, base)
492
     register Lisp_Object arg, base;
Mike Rowan's avatar
Mike Rowan committed
493
{
494
  double d = extract_float (arg);
495

496 497 498 499
#ifdef FLOAT_CHECK_DOMAIN
  if (d <= 0.0)
    domain_error2 ("log", arg, base);
#endif
500
  if (NILP (base))
501
    IN_FLOAT (d = log (d), "log", arg);
502 503 504 505
  else
    {
      double b = extract_float (base);

506 507 508 509 510 511 512
#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
513
	IN_FLOAT2 (d = log (d) / log (b), "log", arg, base);
514
    }
Mike Rowan's avatar
Mike Rowan committed
515 516 517
  return make_float (d);
}

518
DEFUN ("log10", Flog10, Slog10, 1, 1, 0,
519 520
       doc: /* Return the logarithm base 10 of ARG.  */)
     (arg)
521
     register Lisp_Object arg;
Mike Rowan's avatar
Mike Rowan committed
522
{
523 524 525 526 527 528
  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);
529 530 531
  return make_float (d);
}

Mike Rowan's avatar
Mike Rowan committed
532
DEFUN ("sqrt", Fsqrt, Ssqrt, 1, 1, 0,
533 534
       doc: /* Return the square root of ARG.  */)
     (arg)
535
     register Lisp_Object arg;
Mike Rowan's avatar
Mike Rowan committed
536
{
537 538 539 540 541 542
  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
543 544
  return make_float (d);
}
545

Richard M. Stallman's avatar
Richard M. Stallman committed
546
#if 0 /* Not clearly worth adding.  */
Mike Rowan's avatar
Mike Rowan committed
547

548
DEFUN ("acosh", Facosh, Sacosh, 1, 1, 0,
549 550
       doc: /* Return the inverse hyperbolic cosine of ARG.  */)
     (arg)
551
     register Lisp_Object arg;
Mike Rowan's avatar
Mike Rowan committed
552
{
553 554 555 556 557 558 559 560 561 562
  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
563 564 565 566
  return make_float (d);
}

DEFUN ("asinh", Fasinh, Sasinh, 1, 1, 0,
567 568
       doc: /* Return the inverse hyperbolic sine of ARG.  */)
     (arg)
569
     register Lisp_Object arg;
570
{
571 572 573 574 575 576
  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
577 578 579 580
  return make_float (d);
}

DEFUN ("atanh", Fatanh, Satanh, 1, 1, 0,
581 582
       doc: /* Return the inverse hyperbolic tangent of ARG.  */)
     (arg)
583
     register Lisp_Object arg;
584
{
585 586 587 588 589 590 591 592 593 594
  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
595 596 597 598
  return make_float (d);
}

DEFUN ("cosh", Fcosh, Scosh, 1, 1, 0,
599 600
       doc: /* Return the hyperbolic cosine of ARG.  */)
     (arg)
601
     register Lisp_Object arg;
602
{
603 604 605 606 607 608
  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);
609 610 611 612
  return make_float (d);
}

DEFUN ("sinh", Fsinh, Ssinh, 1, 1, 0,
613 614
       doc: /* Return the hyperbolic sine of ARG.  */)
     (arg)
615
     register Lisp_Object arg;
616
{
617 618 619 620 621 622
  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
623 624 625 626
  return make_float (d);
}

DEFUN ("tanh", Ftanh, Stanh, 1, 1, 0,
627 628
       doc: /* Return the hyperbolic tangent of ARG.  */)
     (arg)
629
     register Lisp_Object arg;
Mike Rowan's avatar
Mike Rowan committed
630
{
631 632
  double d = extract_float (arg);
  IN_FLOAT (d = tanh (d), "tanh", arg);
Mike Rowan's avatar
Mike Rowan committed
633 634
  return make_float (d);
}
635
#endif
Mike Rowan's avatar
Mike Rowan committed
636 637

DEFUN ("abs", Fabs, Sabs, 1, 1, 0,
638 639
       doc: /* Return the absolute value of ARG.  */)
     (arg)
640
     register Lisp_Object arg;
Mike Rowan's avatar
Mike Rowan committed
641
{
642
  CHECK_NUMBER_OR_FLOAT (arg);
Mike Rowan's avatar
Mike Rowan committed
643

644
  if (FLOATP (arg))
645
    IN_FLOAT (arg = make_float (fabs (XFLOAT_DATA (arg))), "abs", arg);
646
  else if (XINT (arg) < 0)
647
    XSETINT (arg, - XINT (arg));
Mike Rowan's avatar
Mike Rowan committed
648

649
  return arg;
Mike Rowan's avatar
Mike Rowan committed
650 651 652
}

DEFUN ("float", Ffloat, Sfloat, 1, 1, 0,
653 654
       doc: /* Return the floating point number equal to ARG.  */)
     (arg)
655
     register Lisp_Object arg;
Mike Rowan's avatar
Mike Rowan committed
656
{
657
  CHECK_NUMBER_OR_FLOAT (arg);
Mike Rowan's avatar
Mike Rowan committed
658

659
  if (INTEGERP (arg))
660
    return make_float ((double) XINT (arg));
Mike Rowan's avatar
Mike Rowan committed
661
  else				/* give 'em the same float back */
662
    return arg;
Mike Rowan's avatar
Mike Rowan committed
663 664 665
}

DEFUN ("logb", Flogb, Slogb, 1, 1, 0,
666 667
       doc: /* Returns largest integer <= the base 2 log of the magnitude of ARG.
This is the same as the exponent of a float.  */)
668 669
     (arg)
     Lisp_Object arg;
Mike Rowan's avatar
Mike Rowan committed
670
{
671
  Lisp_Object val;
672
  EMACS_INT value;
673
  double f = extract_float (arg);
674

675
  if (f == 0.0)
Stefan Monnier's avatar
Stefan Monnier committed
676
    value = MOST_NEGATIVE_FIXNUM;
677 678
  else
    {
679
#ifdef HAVE_LOGB
680
      IN_FLOAT (value = logb (f), "logb", arg);
681 682
#else
#ifdef HAVE_FREXP
683 684 685
      int ivalue;
      IN_FLOAT (frexp (f, &ivalue), "logb", arg);
      value = ivalue - 1;
686
#else
687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705
      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;
	}
706
#endif
707
#endif
708
    }
709
  XSETINT (val, value);
710
  return val;
Mike Rowan's avatar
Mike Rowan committed
711 712
}

713

714 715 716 717
/* the rounding functions  */

static Lisp_Object
rounding_driver (arg, divisor, double_round, int_round2, name)
718
     register Lisp_Object arg, divisor;
719 720 721
     double (*double_round) ();
     EMACS_INT (*int_round2) ();
     char *name;
Mike Rowan's avatar
Mike Rowan committed
722
{
723
  CHECK_NUMBER_OR_FLOAT (arg);
Mike Rowan's avatar
Mike Rowan committed
724

725 726
  if (! NILP (divisor))
    {
727
      EMACS_INT i1, i2;
728

729
      CHECK_NUMBER_OR_FLOAT (divisor);
730

731
      if (FLOATP (arg) || FLOATP (divisor))
732 733 734
	{
	  double f1, f2;

735 736
	  f1 = FLOATP (arg) ? XFLOAT_DATA (arg) : XINT (arg);
	  f2 = (FLOATP (divisor) ? XFLOAT_DATA (divisor) : XINT (divisor));
737
	  if (! IEEE_FLOATING_POINT && f2 == 0)
738
	    xsignal0 (Qarith_error);
739

740 741
	  IN_FLOAT2 (f1 = (*double_round) (f1 / f2), name, arg, divisor);
	  FLOAT_TO_INT2 (f1, arg, name, arg, divisor);
742 743 744 745 746 747 748
	  return arg;
	}

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

      if (i2 == 0)
749
	xsignal0 (Qarith_error);
750

751
      XSETINT (arg, (*int_round2) (i1, i2));
752 753 754
      return arg;
    }

755
  if (FLOATP (arg))
756 757
    {
      double d;
758

759
      IN_FLOAT (d = (*double_round) (XFLOAT_DATA (arg)), name, arg);
760
      FLOAT_TO_INT (d, arg, name, arg);
761
    }
Mike Rowan's avatar
Mike Rowan committed
762

763
  return arg;
Mike Rowan's avatar
Mike Rowan committed
764 765
}

766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812
/* 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);
}

813 814 815 816 817
/* 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
818
static double
819
emacs_rint (d)
820 821
     double d;
{
Richard M. Stallman's avatar
Richard M. Stallman committed
822
  return floor (d + 0.5);
823 824 825
}
#endif

826 827 828 829 830 831 832 833
static double
double_identity (d)
     double d;
{
  return d;
}

DEFUN ("ceiling", Fceiling, Sceiling, 1, 2, 0,
834 835
       doc: /* Return the smallest integer no less than ARG.
This rounds the value towards +inf.
836 837
With optional DIVISOR, return the smallest integer no less than ARG/DIVISOR.  */)
     (arg, divisor)
838 839 840 841 842 843
     Lisp_Object arg, divisor;
{
  return rounding_driver (arg, divisor, ceil, ceiling2, "ceiling");
}

DEFUN ("floor", Ffloor, Sfloor, 1, 2, 0,
844
       doc: /* Return the largest integer no greater than ARG.
Lute Kamstra's avatar
Lute Kamstra committed
845
This rounds the value towards -inf.
846 847
With optional DIVISOR, return the largest integer no greater than ARG/DIVISOR.  */)
     (arg, divisor)
848 849 850 851 852 853
     Lisp_Object arg, divisor;
{
  return rounding_driver (arg, divisor, floor, floor2, "floor");
}

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

Eli Zaretskii's avatar
Eli Zaretskii committed
857 858 859
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
860
systems, but 2 on others.  */)
861
     (arg, divisor)
862 863
     Lisp_Object arg, divisor;
{
864
  return rounding_driver (arg, divisor, emacs_rint, round2, "round");
865 866 867
}

DEFUN ("truncate", Ftruncate, Struncate, 1, 2, 0,
868 869 870 871
       doc: /* Truncate a floating point number to an int.
Rounds ARG toward zero.
With optional DIVISOR, truncate ARG/DIVISOR.  */)
     (arg, divisor)
872 873 874 875 876 877
     Lisp_Object arg, divisor;
{
  return rounding_driver (arg, divisor, double_identity, truncate2,
			  "truncate");
}

878

879 880 881 882 883 884
Lisp_Object
fmod_float (x, y)
     register Lisp_Object x, y;
{
  double f1, f2;

885 886
  f1 = FLOATP (x) ? XFLOAT_DATA (x) : XINT (x);
  f2 = FLOATP (y) ? XFLOAT_DATA (y) : XINT (y);
887 888

  if (! IEEE_FLOATING_POINT && f2 == 0)
889
    xsignal0 (Qarith_error);
890 891 892 893 894 895 896

  /* 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);
}
897 898 899 900

/* It's not clear these are worth adding.  */

DEFUN ("fceiling", Ffceiling, Sfceiling, 1, 1, 0,
901 902 903
       doc: /* Return the smallest integer no less than ARG, as a float.
\(Round toward +inf.\)  */)
     (arg)
904 905 906 907 908 909 910 911
     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,
912 913 914
       doc: /* Return the largest integer no greater than ARG, as a float.
\(Round towards -inf.\)  */)
     (arg)
915 916 917 918 919 920
     register Lisp_Object arg;
{
  double d = extract_float (arg);
  IN_FLOAT (d = floor (d), "ffloor", arg);
  return make_float (d);
}
Mike Rowan's avatar
Mike Rowan committed
921

922
DEFUN ("fround", Ffround, Sfround, 1, 1, 0,
923 924
       doc: /* Return the nearest integer to ARG, as a float.  */)
     (arg)
925 926 927
     register Lisp_Object arg;
{
  double d = extract_float (arg);
928
  IN_FLOAT (d = emacs_rint (d), "fround", arg);
929 930 931 932
  return make_float (d);
}

DEFUN ("ftruncate", Fftruncate, Sftruncate, 1, 1, 0,
933 934 935
       doc: /* Truncate a floating point number to an integral float value.
Rounds the value toward zero.  */)
     (arg)
936 937 938 939 940 941
     register Lisp_Object arg;
{
  double d = extract_float (arg);
  if (d >= 0.0)
    IN_FLOAT (d = floor (d), "ftruncate", arg);
  else
942
    IN_FLOAT (d = ceil (d), "ftruncate", arg);
943
  return m