floatfns.c 27.7 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 (double);
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

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

double
Andreas Schwab's avatar
Andreas Schwab committed
195
extract_float (Lisp_Object num)
Mike Rowan's avatar
Mike Rowan committed
196
{
197
  CHECK_NUMBER_OR_FLOAT (num);
Mike Rowan's avatar
Mike Rowan committed
198

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

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

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

232 233 234 235 236 237
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.  */)
Dan Nicolaescu's avatar
Dan Nicolaescu committed
238
  (register Lisp_Object y, Lisp_Object x)
Mike Rowan's avatar
Mike Rowan committed
239
{
240 241 242 243 244 245 246 247 248 249
  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
250 251 252
  return make_float (d);
}

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

262
DEFUN ("sin", Fsin, Ssin, 1, 1, 0,
263
       doc: /* Return the sine of ARG.  */)
Dan Nicolaescu's avatar
Dan Nicolaescu committed
264
  (register Lisp_Object arg)
Mike Rowan's avatar
Mike Rowan committed
265
{
266 267
  double d = extract_float (arg);
  IN_FLOAT (d = sin (d), "sin", arg);
Mike Rowan's avatar
Mike Rowan committed
268 269 270
  return make_float (d);
}

271
DEFUN ("tan", Ftan, Stan, 1, 1, 0,
272
       doc: /* Return the tangent of ARG.  */)
Dan Nicolaescu's avatar
Dan Nicolaescu committed
273
  (register Lisp_Object arg)
274 275 276 277 278 279 280 281
{
  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
282 283
  return make_float (d);
}
284 285 286 287

#if defined HAVE_ISNAN && defined HAVE_COPYSIGN
DEFUN ("isnan", Fisnan, Sisnan, 1, 1, 0,
       doc: /* Return non nil iff argument X is a NaN.  */)
Dan Nicolaescu's avatar
Dan Nicolaescu committed
288
  (Lisp_Object x)
289 290 291 292 293 294 295 296
{
  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.  */)
Dan Nicolaescu's avatar
Dan Nicolaescu committed
297
  (Lisp_Object x1, Lisp_Object x2)
298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319
{
  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.  */)
Dan Nicolaescu's avatar
Dan Nicolaescu committed
320
  (Lisp_Object x)
321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337
{
  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).   */)
Dan Nicolaescu's avatar
Dan Nicolaescu committed
338
  (Lisp_Object sgnfcand, Lisp_Object exp)
339 340 341 342 343
{
  CHECK_NUMBER (exp);
  return make_float (ldexp (XFLOATINT (sgnfcand), XINT (exp)));
}
#endif
Mike Rowan's avatar
Mike Rowan committed
344

345 346
#if 0 /* Leave these out unless we find there's a reason for them.  */

Mike Rowan's avatar
Mike Rowan committed
347
DEFUN ("bessel-j0", Fbessel_j0, Sbessel_j0, 1, 1, 0,
348
       doc: /* Return the bessel function j0 of ARG.  */)
Dan Nicolaescu's avatar
Dan Nicolaescu committed
349
  (register Lisp_Object arg)
Mike Rowan's avatar
Mike Rowan committed
350
{
351 352
  double d = extract_float (arg);
  IN_FLOAT (d = j0 (d), "bessel-j0", arg);
Mike Rowan's avatar
Mike Rowan committed
353 354 355 356
  return make_float (d);
}

DEFUN ("bessel-j1", Fbessel_j1, Sbessel_j1, 1, 1, 0,
357
       doc: /* Return the bessel function j1 of ARG.  */)
Dan Nicolaescu's avatar
Dan Nicolaescu committed
358
  (register Lisp_Object arg)
Mike Rowan's avatar
Mike Rowan committed
359
{
360 361
  double d = extract_float (arg);
  IN_FLOAT (d = j1 (d), "bessel-j1", arg);
Mike Rowan's avatar
Mike Rowan committed
362 363 364 365
  return make_float (d);
}

DEFUN ("bessel-jn", Fbessel_jn, Sbessel_jn, 2, 2, 0,
366 367
       doc: /* Return the order N bessel function output jn of ARG.
The first arg (the order) is truncated to an integer.  */)
Dan Nicolaescu's avatar
Dan Nicolaescu committed
368
  (register Lisp_Object n, Lisp_Object arg)
Mike Rowan's avatar
Mike Rowan committed
369
{
370 371
  int i1 = extract_float (n);
  double f2 = extract_float (arg);
Mike Rowan's avatar
Mike Rowan committed
372

373
  IN_FLOAT (f2 = jn (i1, f2), "bessel-jn", n);
Mike Rowan's avatar
Mike Rowan committed
374 375 376 377
  return make_float (f2);
}

DEFUN ("bessel-y0", Fbessel_y0, Sbessel_y0, 1, 1, 0,
378
       doc: /* Return the bessel function y0 of ARG.  */)
Dan Nicolaescu's avatar
Dan Nicolaescu committed
379
  (register Lisp_Object arg)
Mike Rowan's avatar
Mike Rowan committed
380
{
381 382
  double d = extract_float (arg);
  IN_FLOAT (d = y0 (d), "bessel-y0", arg);
Mike Rowan's avatar
Mike Rowan committed
383 384 385 386
  return make_float (d);
}

DEFUN ("bessel-y1", Fbessel_y1, Sbessel_y1, 1, 1, 0,
387
       doc: /* Return the bessel function y1 of ARG.  */)
Dan Nicolaescu's avatar
Dan Nicolaescu committed
388
  (register Lisp_Object arg)
Mike Rowan's avatar
Mike Rowan committed
389
{
390 391
  double d = extract_float (arg);
  IN_FLOAT (d = y1 (d), "bessel-y0", arg);
Mike Rowan's avatar
Mike Rowan committed
392 393 394 395
  return make_float (d);
}

DEFUN ("bessel-yn", Fbessel_yn, Sbessel_yn, 2, 2, 0,
396 397
       doc: /* Return the order N bessel function output yn of ARG.
The first arg (the order) is truncated to an integer.  */)
Dan Nicolaescu's avatar
Dan Nicolaescu committed
398
  (register Lisp_Object n, Lisp_Object arg)
Mike Rowan's avatar
Mike Rowan committed
399
{
400 401
  int i1 = extract_float (n);
  double f2 = extract_float (arg);
Mike Rowan's avatar
Mike Rowan committed
402

403
  IN_FLOAT (f2 = yn (i1, f2), "bessel-yn", n);
Mike Rowan's avatar
Mike Rowan committed
404 405 406
  return make_float (f2);
}

407 408 409
#endif

#if 0 /* Leave these out unless we see they are worth having.  */
Mike Rowan's avatar
Mike Rowan committed
410 411

DEFUN ("erf", Ferf, Serf, 1, 1, 0,
412
       doc: /* Return the mathematical error function of ARG.  */)
Dan Nicolaescu's avatar
Dan Nicolaescu committed
413
  (register Lisp_Object arg)
Mike Rowan's avatar
Mike Rowan committed
414
{
415 416
  double d = extract_float (arg);
  IN_FLOAT (d = erf (d), "erf", arg);
Mike Rowan's avatar
Mike Rowan committed
417 418 419 420
  return make_float (d);
}

DEFUN ("erfc", Ferfc, Serfc, 1, 1, 0,
421
       doc: /* Return the complementary error function of ARG.  */)
Dan Nicolaescu's avatar
Dan Nicolaescu committed
422
  (register Lisp_Object arg)
Mike Rowan's avatar
Mike Rowan committed
423
{
424 425
  double d = extract_float (arg);
  IN_FLOAT (d = erfc (d), "erfc", arg);
Mike Rowan's avatar
Mike Rowan committed
426 427 428 429
  return make_float (d);
}

DEFUN ("log-gamma", Flog_gamma, Slog_gamma, 1, 1, 0,
430
       doc: /* Return the log gamma of ARG.  */)
Dan Nicolaescu's avatar
Dan Nicolaescu committed
431
  (register Lisp_Object arg)
Mike Rowan's avatar
Mike Rowan committed
432
{
433 434
  double d = extract_float (arg);
  IN_FLOAT (d = lgamma (d), "log-gamma", arg);
Mike Rowan's avatar
Mike Rowan committed
435 436 437
  return make_float (d);
}

438
DEFUN ("cube-root", Fcube_root, Scube_root, 1, 1, 0,
439
       doc: /* Return the cube root of ARG.  */)
Dan Nicolaescu's avatar
Dan Nicolaescu committed
440
  (register Lisp_Object arg)
Mike Rowan's avatar
Mike Rowan committed
441
{
442 443 444 445 446 447 448 449 450
  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
451 452 453
  return make_float (d);
}

Richard M. Stallman's avatar
Richard M. Stallman committed
454 455
#endif

456
DEFUN ("exp", Fexp, Sexp, 1, 1, 0,
457
       doc: /* Return the exponential base e of ARG.  */)
Dan Nicolaescu's avatar
Dan Nicolaescu committed
458
  (register Lisp_Object arg)
459 460 461 462 463 464 465 466 467 468
{
  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
469 470 471 472
  return make_float (d);
}

DEFUN ("expt", Fexpt, Sexpt, 2, 2, 0,
473
       doc: /* Return the exponential ARG1 ** ARG2.  */)
Dan Nicolaescu's avatar
Dan Nicolaescu committed
474
  (register Lisp_Object arg1, Lisp_Object arg2)
Mike Rowan's avatar
Mike Rowan committed
475
{
476
  double f1, f2, f3;
Mike Rowan's avatar
Mike Rowan committed
477

478 479
  CHECK_NUMBER_OR_FLOAT (arg1);
  CHECK_NUMBER_OR_FLOAT (arg2);
480
  if (INTEGERP (arg1)     /* common lisp spec */
481 482
      && 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
483
    {				/* this can be improved by pre-calculating */
484
      EMACS_INT acc, x, y;	/* some binary powers of x then accumulating */
485 486
      Lisp_Object val;

487 488
      x = XINT (arg1);
      y = XINT (arg2);
Mike Rowan's avatar
Mike Rowan committed
489
      acc = 1;
490

Mike Rowan's avatar
Mike Rowan committed
491 492
      if (y < 0)
	{
493 494 495 496 497 498
	  if (x == 1)
	    acc = 1;
	  else if (x == -1)
	    acc = (y & 1) ? -1 : 1;
	  else
	    acc = 0;
Mike Rowan's avatar
Mike Rowan committed
499 500 501
	}
      else
	{
502 503 504 505 506 507 508
	  while (y > 0)
	    {
	      if (y & 1)
		acc *= x;
	      x *= x;
	      y = (unsigned)y >> 1;
	    }
Mike Rowan's avatar
Mike Rowan committed
509
	}
510
      XSETINT (val, acc);
511
      return val;
Mike Rowan's avatar
Mike Rowan committed
512
    }
513 514
  f1 = FLOATP (arg1) ? XFLOAT_DATA (arg1) : XINT (arg1);
  f2 = FLOATP (arg2) ? XFLOAT_DATA (arg2) : XINT (arg2);
515 516 517 518 519 520 521
  /* 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
522 523 524 525 526
  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
527
}
528

529
DEFUN ("log", Flog, Slog, 1, 2, 0,
530
       doc: /* Return the natural logarithm of ARG.
Richard M. Stallman's avatar
Richard M. Stallman committed
531
If the optional argument BASE is given, return log ARG using that base.  */)
Dan Nicolaescu's avatar
Dan Nicolaescu committed
532
  (register Lisp_Object arg, Lisp_Object base)
Mike Rowan's avatar
Mike Rowan committed
533
{
534
  double d = extract_float (arg);
535

536 537 538 539
#ifdef FLOAT_CHECK_DOMAIN
  if (d <= 0.0)
    domain_error2 ("log", arg, base);
#endif
540
  if (NILP (base))
541
    IN_FLOAT (d = log (d), "log", arg);
542 543 544 545
  else
    {
      double b = extract_float (base);

546 547 548 549 550 551 552
#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
553
	IN_FLOAT2 (d = log (d) / log (b), "log", arg, base);
554
    }
Mike Rowan's avatar
Mike Rowan committed
555 556 557
  return make_float (d);
}

558
DEFUN ("log10", Flog10, Slog10, 1, 1, 0,
559
       doc: /* Return the logarithm base 10 of ARG.  */)
Dan Nicolaescu's avatar
Dan Nicolaescu committed
560
  (register Lisp_Object arg)
Mike Rowan's avatar
Mike Rowan committed
561
{
562 563 564 565 566 567
  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);
568 569 570
  return make_float (d);
}

Mike Rowan's avatar
Mike Rowan committed
571
DEFUN ("sqrt", Fsqrt, Ssqrt, 1, 1, 0,
572
       doc: /* Return the square root of ARG.  */)
Dan Nicolaescu's avatar
Dan Nicolaescu committed
573
  (register Lisp_Object arg)
Mike Rowan's avatar
Mike Rowan committed
574
{
575 576 577 578 579 580
  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
581 582
  return make_float (d);
}
583

Richard M. Stallman's avatar
Richard M. Stallman committed
584
#if 0 /* Not clearly worth adding.  */
Mike Rowan's avatar
Mike Rowan committed
585

586
DEFUN ("acosh", Facosh, Sacosh, 1, 1, 0,
587
       doc: /* Return the inverse hyperbolic cosine of ARG.  */)
Dan Nicolaescu's avatar
Dan Nicolaescu committed
588
  (register Lisp_Object arg)
Mike Rowan's avatar
Mike Rowan committed
589
{
590 591 592 593 594 595 596 597 598 599
  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
600 601 602 603
  return make_float (d);
}

DEFUN ("asinh", Fasinh, Sasinh, 1, 1, 0,
604
       doc: /* Return the inverse hyperbolic sine of ARG.  */)
Dan Nicolaescu's avatar
Dan Nicolaescu committed
605
  (register Lisp_Object arg)
606
{
607 608 609 610 611 612
  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
613 614 615 616
  return make_float (d);
}

DEFUN ("atanh", Fatanh, Satanh, 1, 1, 0,
617
       doc: /* Return the inverse hyperbolic tangent of ARG.  */)
Dan Nicolaescu's avatar
Dan Nicolaescu committed
618
  (register Lisp_Object arg)
619
{
620 621 622 623 624 625 626 627 628 629
  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
630 631 632 633
  return make_float (d);
}

DEFUN ("cosh", Fcosh, Scosh, 1, 1, 0,
634
       doc: /* Return the hyperbolic cosine of ARG.  */)
Dan Nicolaescu's avatar
Dan Nicolaescu committed
635
  (register Lisp_Object arg)
636
{
637 638 639 640 641 642
  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);
643 644 645 646
  return make_float (d);
}

DEFUN ("sinh", Fsinh, Ssinh, 1, 1, 0,
647
       doc: /* Return the hyperbolic sine of ARG.  */)
Dan Nicolaescu's avatar
Dan Nicolaescu committed
648
  (register Lisp_Object arg)
649
{
650 651 652 653 654 655
  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
656 657 658 659
  return make_float (d);
}

DEFUN ("tanh", Ftanh, Stanh, 1, 1, 0,
660
       doc: /* Return the hyperbolic tangent of ARG.  */)
Dan Nicolaescu's avatar
Dan Nicolaescu committed
661
  (register Lisp_Object arg)
Mike Rowan's avatar
Mike Rowan committed
662
{
663 664
  double d = extract_float (arg);
  IN_FLOAT (d = tanh (d), "tanh", arg);
Mike Rowan's avatar
Mike Rowan committed
665 666
  return make_float (d);
}
667
#endif
Mike Rowan's avatar
Mike Rowan committed
668 669

DEFUN ("abs", Fabs, Sabs, 1, 1, 0,
670
       doc: /* Return the absolute value of ARG.  */)
Dan Nicolaescu's avatar
Dan Nicolaescu committed
671
  (register Lisp_Object arg)
Mike Rowan's avatar
Mike Rowan committed
672
{
673
  CHECK_NUMBER_OR_FLOAT (arg);
Mike Rowan's avatar
Mike Rowan committed
674

675
  if (FLOATP (arg))
676
    IN_FLOAT (arg = make_float (fabs (XFLOAT_DATA (arg))), "abs", arg);
677
  else if (XINT (arg) < 0)
678
    XSETINT (arg, - XINT (arg));
Mike Rowan's avatar
Mike Rowan committed
679

680
  return arg;
Mike Rowan's avatar
Mike Rowan committed
681 682 683
}

DEFUN ("float", Ffloat, Sfloat, 1, 1, 0,
684
       doc: /* Return the floating point number equal to ARG.  */)
Dan Nicolaescu's avatar
Dan Nicolaescu committed
685
  (register Lisp_Object arg)
Mike Rowan's avatar
Mike Rowan committed
686
{
687
  CHECK_NUMBER_OR_FLOAT (arg);
Mike Rowan's avatar
Mike Rowan committed
688

689
  if (INTEGERP (arg))
690
    return make_float ((double) XINT (arg));
Mike Rowan's avatar
Mike Rowan committed
691
  else				/* give 'em the same float back */
692
    return arg;
Mike Rowan's avatar
Mike Rowan committed
693 694 695
}

DEFUN ("logb", Flogb, Slogb, 1, 1, 0,
696 697
       doc: /* Returns largest integer <= the base 2 log of the magnitude of ARG.
This is the same as the exponent of a float.  */)
Dan Nicolaescu's avatar
Dan Nicolaescu committed
698
  (Lisp_Object arg)
Mike Rowan's avatar
Mike Rowan committed
699
{
700
  Lisp_Object val;
701
  EMACS_INT value;
702
  double f = extract_float (arg);
703

704
  if (f == 0.0)
Stefan Monnier's avatar
Stefan Monnier committed
705
    value = MOST_NEGATIVE_FIXNUM;
706 707
  else
    {
708
#ifdef HAVE_LOGB
709
      IN_FLOAT (value = logb (f), "logb", arg);
710 711
#else
#ifdef HAVE_FREXP
712 713 714
      int ivalue;
      IN_FLOAT (frexp (f, &ivalue), "logb", arg);
      value = ivalue - 1;
715
#else
716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734
      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;
	}
735
#endif
736
#endif
737
    }
738
  XSETINT (val, value);
739
  return val;
Mike Rowan's avatar
Mike Rowan committed
740 741
}

742

743 744 745
/* the rounding functions  */

static Lisp_Object
746 747 748 749
rounding_driver (Lisp_Object arg, Lisp_Object divisor,
		 double (*double_round) (double),
		 EMACS_INT (*int_round2) (EMACS_INT, EMACS_INT),
		 char *name)
Mike Rowan's avatar
Mike Rowan committed
750
{
751
  CHECK_NUMBER_OR_FLOAT (arg);
Mike Rowan's avatar
Mike Rowan committed
752

753 754
  if (! NILP (divisor))
    {
755
      EMACS_INT i1, i2;
756

757
      CHECK_NUMBER_OR_FLOAT (divisor);
758

759
      if (FLOATP (arg) || FLOATP (divisor))
760 761 762
	{
	  double f1, f2;

763 764
	  f1 = FLOATP (arg) ? XFLOAT_DATA (arg) : XINT (arg);
	  f2 = (FLOATP (divisor) ? XFLOAT_DATA (divisor) : XINT (divisor));
765
	  if (! IEEE_FLOATING_POINT && f2 == 0)
766
	    xsignal0 (Qarith_error);
767

768 769
	  IN_FLOAT2 (f1 = (*double_round) (f1 / f2), name, arg, divisor);
	  FLOAT_TO_INT2 (f1, arg, name, arg, divisor);
770 771 772 773 774 775 776
	  return arg;
	}

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

      if (i2 == 0)
777
	xsignal0 (Qarith_error);
778

779
      XSETINT (arg, (*int_round2) (i1, i2));
780 781 782
      return arg;
    }

783
  if (FLOATP (arg))
784 785
    {
      double d;
786

787
      IN_FLOAT (d = (*double_round) (XFLOAT_DATA (arg)), name, arg);
788
      FLOAT_TO_INT (d, arg, name, arg);
789
    }
Mike Rowan's avatar
Mike Rowan committed
790

791
  return arg;
Mike Rowan's avatar
Mike Rowan committed
792 793
}

794 795 796 797 798
/* 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
799
ceiling2 (EMACS_INT i1, EMACS_INT i2)
800 801 802 803 804 805 806
{
  return (i2 < 0
	  ? (i1 < 0  ?  ((-1 - i1) / -i2) + 1  :  - (i1 / -i2))
	  : (i1 <= 0  ?  - (-i1 / i2)  :  ((i1 - 1) / i2) + 1));
}

static EMACS_INT
807
floor2 (EMACS_INT i1, EMACS_INT i2)
808 809 810 811 812 813 814
{
  return (i2 < 0
	  ? (i1 <= 0  ?  -i1 / -i2  :  -1 - ((i1 - 1) / -i2))
	  : (i1 < 0  ?  -1 - ((-1 - i1) / i2)  :  i1 / i2));
}

static EMACS_INT
815
truncate2 (EMACS_INT i1, EMACS_INT i2)
816 817 818 819 820 821 822
{
  return (i2 < 0
	  ? (i1 < 0  ?  -i1 / -i2  :  - (i1 / -i2))
	  : (i1 < 0  ?  - (-i1 / i2)  :  i1 / i2));
}

static EMACS_INT
823
round2 (EMACS_INT i1, EMACS_INT i2)
824 825 826 827 828 829 830 831 832 833 834 835 836
{
  /* 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);
}

837 838 839 840 841
/* 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
842
static double
843
emacs_rint (double d)
844
{
Richard M. Stallman's avatar
Richard M. Stallman committed
845
  return floor (d + 0.5);
846 847 848
}
#endif

849
static double
850
double_identity (double d)
851 852 853 854 855
{
  return d;
}

DEFUN ("ceiling", Fceiling, Sceiling, 1, 2, 0,
856 857
       doc: /* Return the smallest integer no less than ARG.
This rounds the value towards +inf.
858
With optional DIVISOR, return the smallest integer no less than ARG/DIVISOR.  */)
Dan Nicolaescu's avatar
Dan Nicolaescu committed
859
  (Lisp_Object arg, Lisp_Object divisor)
860 861 862 863 864
{
  return rounding_driver (arg, divisor, ceil, ceiling2, "ceiling");
}

DEFUN ("floor", Ffloor, Sfloor, 1, 2, 0,
865
       doc: /* Return the largest integer no greater than ARG.
Lute Kamstra's avatar
Lute Kamstra committed
866
This rounds the value towards -inf.
867
With optional DIVISOR, return the largest integer no greater than ARG/DIVISOR.  */)
Dan Nicolaescu's avatar
Dan Nicolaescu committed
868
  (Lisp_Object arg, Lisp_Object divisor)
869 870 871 872 873
{
  return rounding_driver (arg, divisor, floor, floor2, "floor");
}

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

Eli Zaretskii's avatar
Eli Zaretskii committed
877 878 879
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
880
systems, but 2 on others.  */)
Dan Nicolaescu's avatar
Dan Nicolaescu committed
881
  (Lisp_Object arg, Lisp_Object divisor)
882
{
883
  return rounding_driver (arg, divisor, emacs_rint, round2, "round");
884 885 886
}

DEFUN ("truncate", Ftruncate, Struncate, 1, 2, 0,
887 888 889
       doc: /* Truncate a floating point number to an int.
Rounds ARG toward zero.
With optional DIVISOR, truncate ARG/DIVISOR.  */)
Dan Nicolaescu's avatar
Dan Nicolaescu committed
890
  (Lisp_Object arg, Lisp_Object divisor)
891 892 893 894 895
{
  return rounding_driver (arg, divisor, double_identity, truncate2,
			  "truncate");
}

896

897
Lisp_Object
898
fmod_float (Lisp_Object x, Lisp_Object y)
899 900 901
{
  double f1, f2;

902 903
  f1 = FLOATP (x) ? XFLOAT_DATA (x) : XINT (x);
  f2 = FLOATP (y) ? XFLOAT_DATA (y) : XINT (y);
904 905

  if (! IEEE_FLOATING_POINT && f2 == 0)
906
    xsignal0 (Qarith_error);