math.c

   1 /**********************************************************************
   2
   3   math.c -
   4
   5   $Author$
   6   created at: Tue Jan 25 14:12:56 JST 1994
   7
   8   Copyright (C) 1993-2007 Yukihiro Matsumoto
   9
  10 **********************************************************************/
  11
  12 #include "ruby/internal/config.h"
  13
  14 #ifdef _MSC_VER
  15 # define _USE_MATH_DEFINES 1
  16 #endif
  17
  18 #include <errno.h>
  19 #include <float.h>
  20 #include <math.h>
  21
  22 #include "internal.h"
  23 #include "internal/bignum.h"
  24 #include "internal/complex.h"
  25 #include "internal/math.h"
  26 #include "internal/object.h"
  27 #include "internal/vm.h"
  28
  29 VALUE rb_mMath;
  30 VALUE rb_eMathDomainError;
  31
  32 #define Get_Double(x) rb_num_to_dbl(x)
  33
  34 #define domain_error(msg) \
  35     rb_raise(rb_eMathDomainError, "Numerical argument is out of domain - " msg)
  36 #define domain_check_min(val, min, msg) \
  37     ((val) < (min) ? domain_error(msg) : (void)0)
  38 #define domain_check_range(val, min, max, msg) \
  39     ((val) < (min) || (max) < (val) ? domain_error(msg) : (void)0)
  40
  41 /*
  42  *  call-seq:
  43  *     Math.atan2(y, x) -> float
  44  *
  45  *  Returns the {arc tangent}[https://en.wikipedia.org/wiki/Atan2] of +y+ and +x+
  46  *  in {radians}[https://en.wikipedia.org/wiki/Trigonometric_functions#Radians_versus_degrees].
  47  *
  48  *  - Domain of +y+: <tt>[-INFINITY, INFINITY]</tt>.
  49  *  - Domain of +x+: <tt>[-INFINITY, INFINITY]</tt>.
  50  *  - Range: <tt>[-PI, PI]</tt>.
  51  *
  52  *  Examples:
  53  *
  54  *    atan2(-1.0, -1.0) # => -2.356194490192345  # -3*PI/4
  55  *    atan2(-1.0, 0.0)  # => -1.5707963267948966 # -PI/2
  56  *    atan2(-1.0, 1.0)  # => -0.7853981633974483 # -PI/4
  57  *    atan2(0.0, -1.0)  # => 3.141592653589793   # PI
  58  *
  59  */
  60
  61 static VALUE
  62 math_atan2(VALUE unused_obj, VALUE y, VALUE x)
  63 {
  64     double dx, dy;
  65     dx = Get_Double(x);
  66     dy = Get_Double(y);
  67     if (dx == 0.0 && dy == 0.0) {
  68         if (!signbit(dx))
  69             return DBL2NUM(dy);
  70         if (!signbit(dy))
  71             return DBL2NUM(M_PI);
  72         return DBL2NUM(-M_PI);
  73     }
  74 #ifndef ATAN2_INF_C99
  75     if (isinf(dx) && isinf(dy)) {
  76         /* optimization for FLONUM */
  77         if (dx < 0.0) {
  78             const double dz = (3.0 * M_PI / 4.0);
  79             return (dy < 0.0) ? DBL2NUM(-dz) : DBL2NUM(dz);
  80         }
  81         else {
  82             const double dz = (M_PI / 4.0);
  83             return (dy < 0.0) ? DBL2NUM(-dz) : DBL2NUM(dz);
  84         }
  85     }
  86 #endif
  87     return DBL2NUM(atan2(dy, dx));
  88 }
  89
  90
  91 /*
  92  *  call-seq:
  93  *    Math.cos(x) -> float
  94  *
  95  *  Returns the
  96  *  {cosine}[https://en.wikipedia.org/wiki/Sine_and_cosine] of +x+
  97  *  in {radians}[https://en.wikipedia.org/wiki/Trigonometric_functions#Radians_versus_degrees].
  98  *
  99  *  - Domain: <tt>(-INFINITY, INFINITY)</tt>.
 100  *  - Range: <tt>[-1.0, 1.0]</tt>.
 101  *
 102  *  Examples:
 103  *
 104  *    cos(-PI)   # => -1.0
 105  *    cos(-PI/2) # => 6.123031769111886e-17 # 0.0000000000000001
 106  *    cos(0.0)   # => 1.0
 107  *    cos(PI/2)  # => 6.123031769111886e-17 # 0.0000000000000001
 108  *    cos(PI)    # => -1.0
 109  *
 110  */
 111
 112 static VALUE
 113 math_cos(VALUE unused_obj, VALUE x)
 114 {
 115     return DBL2NUM(cos(Get_Double(x)));
 116 }
 117
 118 /*
 119  *  call-seq:
 120  *    Math.sin(x) -> float
 121  *
 122  *  Returns the
 123  *  {sine}[https://en.wikipedia.org/wiki/Sine_and_cosine] of +x+
 124  *  in {radians}[https://en.wikipedia.org/wiki/Trigonometric_functions#Radians_versus_degrees].
 125  *
 126  *  - Domain: <tt>(-INFINITY, INFINITY)</tt>.
 127  *  - Range: <tt>[-1.0, 1.0]</tt>.
 128  *
 129  *  Examples:
 130  *
 131  *    sin(-PI)   # => -1.2246063538223773e-16 # -0.0000000000000001
 132  *    sin(-PI/2) # => -1.0
 133  *    sin(0.0)   # => 0.0
 134  *    sin(PI/2)  # => 1.0
 135  *    sin(PI)    # => 1.2246063538223773e-16  # 0.0000000000000001
 136  *
 137  */
 138
 139 static VALUE
 140 math_sin(VALUE unused_obj, VALUE x)
 141 {
 142     return DBL2NUM(sin(Get_Double(x)));
 143 }
 144
 145
 146 /*
 147  *  call-seq:
 148  *    Math.tan(x) -> float
 149  *
 150  *  Returns the
 151  *  {tangent}[https://en.wikipedia.org/wiki/Trigonometric_functions] of +x+
 152  *  in {radians}[https://en.wikipedia.org/wiki/Trigonometric_functions#Radians_versus_degrees].
 153  *
 154  *  - Domain: <tt>(-INFINITY, INFINITY)</tt>.
 155  *  - Range: <tt>(-INFINITY, INFINITY)</tt>.
 156  *
 157  *  Examples:
 158  *
 159  *    tan(-PI)   # => 1.2246467991473532e-16  # -0.0000000000000001
 160  *    tan(-PI/2) # => -1.633123935319537e+16  # -16331239353195370.0
 161  *    tan(0.0)   # => 0.0
 162  *    tan(PI/2)  # => 1.633123935319537e+16   # 16331239353195370.0
 163  *    tan(PI)    # => -1.2246467991473532e-16 # -0.0000000000000001
 164  *
 165  */
 166
 167 static VALUE
 168 math_tan(VALUE unused_obj, VALUE x)
 169 {
 170     return DBL2NUM(tan(Get_Double(x)));
 171 }
 172
 173 #define math_arc(num, func) \
 174     double d; \
 175     d = Get_Double((num)); \
 176     domain_check_range(d, -1.0, 1.0, #func); \
 177     return DBL2NUM(func(d));
 178
 179 /*
 180  *  call-seq:
 181  *     Math.acos(x) -> float
 182  *
 183  *  Returns the {arc cosine}[https://en.wikipedia.org/wiki/Inverse_trigonometric_functions] of +x+.
 184  *
 185  *  - Domain: <tt>[-1, 1]</tt>.
 186  *  - Range: <tt>[0, PI]</tt>.
 187  *
 188  *  Examples:
 189  *
 190  *    acos(-1.0) # => 3.141592653589793  # PI
 191  *    acos(0.0)  # => 1.5707963267948966 # PI/2
 192  *    acos(1.0)  # => 0.0
 193  *
 194  */
 195
 196 static VALUE
 197 math_acos(VALUE unused_obj, VALUE x)
 198 {
 199     math_arc(x, acos)
 200 }
 201
 202 /*
 203  *  call-seq:
 204  *     Math.asin(x) -> float
 205  *
 206  *  Returns the {arc sine}[https://en.wikipedia.org/wiki/Inverse_trigonometric_functions] of +x+.
 207  *
 208  *  - Domain: <tt>[-1, -1]</tt>.
 209  *  - Range: <tt>[-PI/2, PI/2]</tt>.
 210  *
 211  *  Examples:
 212  *
 213  *    asin(-1.0) # => -1.5707963267948966 # -PI/2
 214  *    asin(0.0)  # => 0.0
 215  *    asin(1.0)  # => 1.5707963267948966  # PI/2
 216  *
 217  */
 218
 219 static VALUE
 220 math_asin(VALUE unused_obj, VALUE x)
 221 {
 222     math_arc(x, asin)
 223 }
 224
 225 /*
 226  *  call-seq:
 227  *     Math.atan(x)    -> Float
 228  *
 229  *  Returns the {arc tangent}[https://en.wikipedia.org/wiki/Inverse_trigonometric_functions] of +x+.
 230  *
 231  *  - Domain: <tt>[-INFINITY, INFINITY]</tt>.
 232  *  - Range: <tt>[-PI/2, PI/2]  </tt>.
 233  *
 234  *  Examples:
 235  *
 236  *    atan(-INFINITY) # => -1.5707963267948966 # -PI2
 237  *    atan(-PI)       # => -1.2626272556789115
 238  *    atan(-PI/2)     # => -1.0038848218538872
 239  *    atan(0.0)       # => 0.0
 240  *    atan(PI/2)      # => 1.0038848218538872
 241  *    atan(PI)        # => 1.2626272556789115
 242  *    atan(INFINITY)  # => 1.5707963267948966  # PI/2
 243  *
 244  */
 245
 246 static VALUE
 247 math_atan(VALUE unused_obj, VALUE x)
 248 {
 249     return DBL2NUM(atan(Get_Double(x)));
 250 }
 251
 252 #ifndef HAVE_COSH
 253 double
 254 cosh(double x)
 255 {
 256     return (exp(x) + exp(-x)) / 2;
 257 }
 258 #endif
 259
 260 /*
 261  *  call-seq:
 262  *    Math.cosh(x) -> float
 263  *
 264  *  Returns the {hyperbolic cosine}[https://en.wikipedia.org/wiki/Hyperbolic_functions] of +x+
 265  *  in {radians}[https://en.wikipedia.org/wiki/Trigonometric_functions#Radians_versus_degrees].
 266  *
 267  *  - Domain: <tt>[-INFINITY, INFINITY]</tt>.
 268  *  - Range: <tt>[1, INFINITY]</tt>.
 269  *
 270  *  Examples:
 271  *
 272  *    cosh(-INFINITY) # => Infinity
 273  *    cosh(0.0)       # => 1.0
 274  *    cosh(INFINITY)  # => Infinity
 275  *
 276  */
 277
 278 static VALUE
 279 math_cosh(VALUE unused_obj, VALUE x)
 280 {
 281     return DBL2NUM(cosh(Get_Double(x)));
 282 }
 283
 284 #ifndef HAVE_SINH
 285 double
 286 sinh(double x)
 287 {
 288     return (exp(x) - exp(-x)) / 2;
 289 }
 290 #endif
 291
 292 /*
 293  *  call-seq:
 294  *    Math.sinh(x) -> float
 295  *
 296  *  Returns the {hyperbolic sine}[https://en.wikipedia.org/wiki/Hyperbolic_functions] of +x+
 297  *  in {radians}[https://en.wikipedia.org/wiki/Trigonometric_functions#Radians_versus_degrees].
 298  *
 299  *  - Domain: <tt>[-INFINITY, INFINITY]</tt>.
 300  *  - Range: <tt>[-INFINITY, INFINITY]</tt>.
 301  *
 302  *  Examples:
 303  *
 304  *    sinh(-INFINITY) # => -Infinity
 305  *    sinh(0.0)       # => 0.0
 306  *    sinh(INFINITY)  # => Infinity
 307  *
 308  */
 309
 310 static VALUE
 311 math_sinh(VALUE unused_obj, VALUE x)
 312 {
 313     return DBL2NUM(sinh(Get_Double(x)));
 314 }
 315
 316 #ifndef HAVE_TANH
 317 double
 318 tanh(double x)
 319 {
 320 # if defined(HAVE_SINH) && defined(HAVE_COSH)
 321     const double c = cosh(x);
 322     if (!isinf(c)) return sinh(x) / c;
 323 # else
 324     const double e = exp(x+x);
 325     if (!isinf(e)) return (e - 1) / (e + 1);
 326 # endif
 327     return x > 0 ? 1.0 : -1.0;
 328 }
 329 #endif
 330
 331 /*
 332  *  call-seq:
 333  *    Math.tanh(x) -> float
 334  *
 335  *  Returns the {hyperbolic tangent}[https://en.wikipedia.org/wiki/Hyperbolic_functions] of +x+
 336  *  in {radians}[https://en.wikipedia.org/wiki/Trigonometric_functions#Radians_versus_degrees].
 337  *
 338  *  - Domain: <tt>[-INFINITY, INFINITY]</tt>.
 339  *  - Range: <tt>[-1, 1]</tt>.
 340  *
 341  *  Examples:
 342  *
 343  *    tanh(-INFINITY) # => -1.0
 344  *    tanh(0.0)       # => 0.0
 345  *    tanh(INFINITY)  # => 1.0
 346  *
 347  */
 348
 349 static VALUE
 350 math_tanh(VALUE unused_obj, VALUE x)
 351 {
 352     return DBL2NUM(tanh(Get_Double(x)));
 353 }
 354
 355 /*
 356  *  call-seq:
 357  *    Math.acosh(x) -> float
 358  *
 359  *  Returns the {inverse hyperbolic cosine}[https://en.wikipedia.org/wiki/Inverse_hyperbolic_functions] of +x+.
 360  *
 361  *  - Domain: <tt>[1, INFINITY]</tt>.
 362  *  - Range: <tt>[0, INFINITY]</tt>.
 363  *
 364  *  Examples:
 365  *
 366  *    acosh(1.0)      # => 0.0
 367  *    acosh(INFINITY) # => Infinity
 368  *
 369  */
 370
 371 static VALUE
 372 math_acosh(VALUE unused_obj, VALUE x)
 373 {
 374     double d;
 375
 376     d = Get_Double(x);
 377     domain_check_min(d, 1.0, "acosh");
 378     return DBL2NUM(acosh(d));
 379 }
 380
 381 /*
 382  *  call-seq:
 383  *    Math.asinh(x) -> float
 384  *
 385  *  Returns the {inverse hyperbolic sine}[https://en.wikipedia.org/wiki/Inverse_hyperbolic_functions] of +x+.
 386  *
 387  *  - Domain: <tt>[-INFINITY, INFINITY]</tt>.
 388  *  - Range: <tt>[-INFINITY, INFINITY]</tt>.
 389  *
 390  *  Examples:
 391  *
 392  *    asinh(-INFINITY) # => -Infinity
 393  *    asinh(0.0)       # => 0.0
 394  *    asinh(INFINITY)  # => Infinity
 395  *
 396  */
 397
 398 static VALUE
 399 math_asinh(VALUE unused_obj, VALUE x)
 400 {
 401     return DBL2NUM(asinh(Get_Double(x)));
 402 }
 403
 404 /*
 405  *  call-seq:
 406  *    Math.atanh(x) -> float
 407  *
 408  *  Returns the {inverse hyperbolic tangent}[https://en.wikipedia.org/wiki/Inverse_hyperbolic_functions] of +x+.
 409  *
 410  *  - Domain: <tt>[-1, 1]</tt>.
 411  *  - Range: <tt>[-INFINITY, INFINITY]</tt>.
 412  *
 413  *  Examples:
 414  *
 415  *    atanh(-1.0) # => -Infinity
 416  *    atanh(0.0)  # => 0.0
 417  *    atanh(1.0)  # => Infinity
 418  *
 419  */
 420
 421 static VALUE
 422 math_atanh(VALUE unused_obj, VALUE x)
 423 {
 424     double d;
 425
 426     d = Get_Double(x);
 427     domain_check_range(d, -1.0, +1.0, "atanh");
 428     /* check for pole error */
 429     if (d == -1.0) return DBL2NUM(-HUGE_VAL);
 430     if (d == +1.0) return DBL2NUM(+HUGE_VAL);
 431     return DBL2NUM(atanh(d));
 432 }
 433
 434 /*
 435  *  call-seq:
 436  *    Math.exp(x) -> float
 437  *
 438  *  Returns +e+ raised to the +x+ power.
 439  *
 440  *  - Domain: <tt>[-INFINITY, INFINITY]</tt>.
 441  *  - Range: <tt>[0, INFINITY]</tt>.
 442  *
 443  *  Examples:
 444  *
 445  *    exp(-INFINITY) # => 0.0
 446  *    exp(-1.0)      # => 0.36787944117144233 # 1.0/E
 447  *    exp(0.0)       # => 1.0
 448  *    exp(0.5)       # => 1.6487212707001282  # sqrt(E)
 449  *    exp(1.0)       # => 2.718281828459045   # E
 450  *    exp(2.0)       # => 7.38905609893065    # E**2
 451  *    exp(INFINITY)  # => Infinity
 452  *
 453  */
 454
 455 static VALUE
 456 math_exp(VALUE unused_obj, VALUE x)
 457 {
 458     return DBL2NUM(exp(Get_Double(x)));
 459 }
 460
 461 #if defined __CYGWIN__
 462 # include <cygwin/version.h>
 463 # if CYGWIN_VERSION_DLL_MAJOR < 1005
 464 #  define nan(x) nan()
 465 # endif
 466 # define log(x) ((x) < 0.0 ? nan("") : log(x))
 467 # define log10(x) ((x) < 0.0 ? nan("") : log10(x))
 468 #endif
 469
 470 #ifndef M_LN2
 471 # define M_LN2 0.693147180559945309417232121458176568
 472 #endif
 473 #ifndef M_LN10
 474 # define M_LN10 2.30258509299404568401799145468436421
 475 #endif
 476
 477 FUNC_MINIMIZED(static VALUE math_log(int, const VALUE *, VALUE));
 478
 479 /*
 480  *  call-seq:
 481  *    Math.log(x, base = Math::E) -> Float
 482  *
 483  *  Returns the base +base+ {logarithm}[https://en.wikipedia.org/wiki/Logarithm] of +x+.
 484  *
 485  *  - Domain: <tt>[0, INFINITY]</tt>.
 486  *  - Range: <tt>[-INFINITY, INFINITY)]</tt>.
 487  *
 488  *  Examples:
 489  *
 490  *    log(0.0)        # => -Infinity
 491  *    log(1.0)        # => 0.0
 492  *    log(E)          # => 1.0
 493  *    log(INFINITY)   # => Infinity
 494  *
 495  *    log(0.0, 2.0)   # => -Infinity
 496  *    log(1.0, 2.0)   # => 0.0
 497  *    log(2.0, 2.0)   # => 1.0
 498  *
 499  *    log(0.0, 10.0)  # => -Infinity
 500  *    log(1.0, 10.0)  # => 0.0
 501  *    log(10.0, 10.0) # => 1.0
 502  *
 503  */
 504
 505 static VALUE
 506 math_log(int argc, const VALUE *argv, VALUE unused_obj)
 507 {
 508     return rb_math_log(argc, argv);
 509 }
 510
 511 static double
 512 get_double_rshift(VALUE x, size_t *pnumbits)
 513 {
 514     size_t numbits;
 515
 516     if (RB_BIGNUM_TYPE_P(x) && BIGNUM_POSITIVE_P(x) &&
 517             DBL_MAX_EXP <= (numbits = rb_absint_numwords(x, 1, NULL))) {
 518         numbits -= DBL_MANT_DIG;
 519         x = rb_big_rshift(x, SIZET2NUM(numbits));
 520     }
 521     else {
 522         numbits = 0;
 523     }
 524     *pnumbits = numbits;
 525     return Get_Double(x);
 526 }
 527
 528 static double
 529 math_log_split(VALUE x, size_t *numbits)
 530 {
 531     double d = get_double_rshift(x, numbits);
 532
 533     domain_check_min(d, 0.0, "log");
 534     return d;
 535 }
 536
 537 #if defined(log2) || defined(HAVE_LOG2)
 538 # define log_intermediate log2
 539 #else
 540 # define log_intermediate log10
 541 double log2(double x);
 542 #endif
 543
 544 VALUE
 545 rb_math_log(int argc, const VALUE *argv)
 546 {
 547     VALUE x, base;
 548     double d;
 549     size_t numbits;
 550
 551     argc = rb_scan_args(argc, argv, "11", &x, &base);
 552     d = math_log_split(x, &numbits);
 553     if (argc == 2) {
 554         size_t numbits_2;
 555         double b = math_log_split(base, &numbits_2);
 556         /* check for pole error */
 557         if (d == 0.0) {
 558             // Already DomainError if b < 0.0
 559             return b ? DBL2NUM(-HUGE_VAL) : DBL2NUM(NAN);
 560         }
 561         else if (b == 0.0) {
 562             return DBL2NUM(-0.0);
 563         }
 564         d = log_intermediate(d) / log_intermediate(b);
 565         d += (numbits - numbits_2) / log2(b);
 566     }
 567     else {
 568         /* check for pole error */
 569         if (d == 0.0) return DBL2NUM(-HUGE_VAL);
 570         d = log(d);
 571         d += numbits * M_LN2;
 572     }
 573     return DBL2NUM(d);
 574 }
 575
 576 #ifndef log2
 577 #ifndef HAVE_LOG2
 578 double
 579 log2(double x)
 580 {
 581     return log10(x)/log10(2.0);
 582 }
 583 #else
 584 extern double log2(double);
 585 #endif
 586 #endif
 587
 588 /*
 589  *  call-seq:
 590  *    Math.log2(x) -> float
 591  *
 592  *  Returns the base 2 {logarithm}[https://en.wikipedia.org/wiki/Logarithm] of +x+.
 593  *
 594  *  - Domain: <tt>[0, INFINITY]</tt>.
 595  *  - Range: <tt>[-INFINITY, INFINITY]</tt>.
 596  *
 597  *  Examples:
 598  *
 599  *    log2(0.0)      # => -Infinity
 600  *    log2(1.0)      # => 0.0
 601  *    log2(2.0)      # => 1.0
 602  *    log2(INFINITY) # => Infinity
 603  *
 604  */
 605
 606 static VALUE
 607 math_log2(VALUE unused_obj, VALUE x)
 608 {
 609     size_t numbits;
 610     double d = get_double_rshift(x, &numbits);
 611
 612     domain_check_min(d, 0.0, "log2");
 613     /* check for pole error */
 614     if (d == 0.0) return DBL2NUM(-HUGE_VAL);
 615
 616     return DBL2NUM(log2(d) + numbits); /* log2(d * 2 ** numbits) */
 617 }
 618
 619 /*
 620  *  call-seq:
 621  *    Math.log10(x) -> float
 622  *
 623  *  Returns the base 10 {logarithm}[https://en.wikipedia.org/wiki/Logarithm] of +x+.
 624  *
 625  *  - Domain: <tt>[0, INFINITY]</tt>.
 626  *  - Range: <tt>[-INFINITY, INFINITY]</tt>.
 627  *
 628  *  Examples:
 629  *
 630  *    log10(0.0)      # => -Infinity
 631  *    log10(1.0)      # => 0.0
 632  *    log10(10.0)     # => 1.0
 633  *    log10(INFINITY) # => Infinity
 634  *
 635  */
 636
 637 static VALUE
 638 math_log10(VALUE unused_obj, VALUE x)
 639 {
 640     size_t numbits;
 641     double d = get_double_rshift(x, &numbits);
 642
 643     domain_check_min(d, 0.0, "log10");
 644     /* check for pole error */
 645     if (d == 0.0) return DBL2NUM(-HUGE_VAL);
 646
 647     return DBL2NUM(log10(d) + numbits * log10(2)); /* log10(d * 2 ** numbits) */
 648 }
 649
 650 static VALUE rb_math_sqrt(VALUE x);
 651
 652 /*
 653  *  call-seq:
 654  *    Math.sqrt(x) -> float
 655  *
 656  *  Returns the principal (non-negative) {square root}[https://en.wikipedia.org/wiki/Square_root] of +x+.
 657  *
 658  *  - Domain: <tt>[0, INFINITY]</tt>.
 659  *  - Range: <tt>[0, INFINITY]</tt>.
 660  *
 661  *  Examples:
 662  *
 663  *    sqrt(0.0)      # => 0.0
 664  *    sqrt(0.5)      # => 0.7071067811865476
 665  *    sqrt(1.0)      # => 1.0
 666  *    sqrt(2.0)      # => 1.4142135623730951
 667  *    sqrt(4.0)      # => 2.0
 668  *    sqrt(9.0)      # => 3.0
 669  *    sqrt(INFINITY) # => Infinity
 670  *
 671  */
 672
 673 static VALUE
 674 math_sqrt(VALUE unused_obj, VALUE x)
 675 {
 676     return rb_math_sqrt(x);
 677 }
 678
 679 inline static VALUE
 680 f_negative_p(VALUE x)
 681 {
 682     if (FIXNUM_P(x))
 683         return RBOOL(FIX2LONG(x) < 0);
 684     return rb_funcall(x, '<', 1, INT2FIX(0));
 685 }
 686 inline static VALUE
 687 f_signbit(VALUE x)
 688 {
 689     if (RB_FLOAT_TYPE_P(x)) {
 690         double f = RFLOAT_VALUE(x);
 691         return RBOOL(!isnan(f) && signbit(f));
 692     }
 693     return f_negative_p(x);
 694 }
 695
 696 static VALUE
 697 rb_math_sqrt(VALUE x)
 698 {
 699     double d;
 700
 701     if (RB_TYPE_P(x, T_COMPLEX)) {
 702         VALUE neg = f_signbit(RCOMPLEX(x)->imag);
 703         double re = Get_Double(RCOMPLEX(x)->real), im;
 704         d = Get_Double(rb_complex_abs(x));
 705         im = sqrt((d - re) / 2.0);
 706         re = sqrt((d + re) / 2.0);
 707         if (neg) im = -im;
 708         return rb_complex_new(DBL2NUM(re), DBL2NUM(im));
 709     }
 710     d = Get_Double(x);
 711     domain_check_min(d, 0.0, "sqrt");
 712     if (d == 0.0) return DBL2NUM(0.0);
 713     return DBL2NUM(sqrt(d));
 714 }
 715
 716 /*
 717  *  call-seq:
 718  *    Math.cbrt(x) -> float
 719  *
 720  *  Returns the {cube root}[https://en.wikipedia.org/wiki/Cube_root] of +x+.
 721  *
 722  *  - Domain: <tt>[-INFINITY, INFINITY]</tt>.
 723  *  - Range: <tt>[-INFINITY, INFINITY]</tt>.
 724  *
 725  *  Examples:
 726  *
 727  *    cbrt(-INFINITY) # => -Infinity
 728  *    cbrt(-27.0)     # => -3.0
 729  *    cbrt(-8.0)      # => -2.0
 730  *    cbrt(-2.0)      # => -1.2599210498948732
 731  *    cbrt(1.0)       # => 1.0
 732  *    cbrt(0.0)       # => 0.0
 733  *    cbrt(1.0)       # => 1.0
 734  *    cbrt(2.0)       # => 1.2599210498948732
 735  *    cbrt(8.0)       # => 2.0
 736  *    cbrt(27.0)      # => 3.0
 737  *    cbrt(INFINITY)  # => Infinity
 738  *
 739  */
 740
 741 static VALUE
 742 math_cbrt(VALUE unused_obj, VALUE x)
 743 {
 744     double f = Get_Double(x);
 745     double r = cbrt(f);
 746 #if defined __GLIBC__
 747     if (isfinite(r) && !(f == 0.0 && r == 0.0)) {
 748         r = (2.0 * r + (f / r / r)) / 3.0;
 749     }
 750 #endif
 751     return DBL2NUM(r);
 752 }
 753
 754 /*
 755  *  call-seq:
 756  *    Math.frexp(x) -> [fraction, exponent]
 757  *
 758  *  Returns a 2-element array containing the normalized signed float +fraction+
 759  *  and integer +exponent+ of +x+ such that:
 760  *
 761  *    x = fraction * 2**exponent
 762  *
 763  *  See {IEEE 754 double-precision binary floating-point format: binary64}[https://en.wikipedia.org/wiki/Double-precision_floating-point_format#IEEE_754_double-precision_binary_floating-point_format:_binary64].
 764  *
 765  *  - Domain: <tt>[-INFINITY, INFINITY]</tt>.
 766  *  - Range <tt>[-INFINITY, INFINITY]</tt>.
 767  *
 768  *  Examples:
 769  *
 770  *    frexp(-INFINITY) # => [-Infinity, -1]
 771  *    frexp(-2.0)      # => [-0.5, 2]
 772  *    frexp(-1.0)      # => [-0.5, 1]
 773  *    frexp(0.0)       # => [0.0, 0]
 774  *    frexp(1.0)       # => [0.5, 1]
 775  *    frexp(2.0)       # => [0.5, 2]
 776  *    frexp(INFINITY)  # => [Infinity, -1]
 777  *
 778  *  Related: Math.ldexp (inverse of Math.frexp).
 779  *
 780  */
 781
 782 static VALUE
 783 math_frexp(VALUE unused_obj, VALUE x)
 784 {
 785     double d;
 786     int exp;
 787
 788     d = frexp(Get_Double(x), &exp);
 789     return rb_assoc_new(DBL2NUM(d), INT2NUM(exp));
 790 }
 791
 792 /*
 793  *  call-seq:
 794  *    Math.ldexp(fraction, exponent) -> float
 795  *
 796  *  Returns the value of <tt>fraction * 2**exponent</tt>.
 797  *
 798  *  - Domain of +fraction+: <tt>[0.0, 1.0)</tt>.
 799  *  - Domain of +exponent+: <tt>[0, 1024]</tt>
 800  *    (larger values are equivalent to 1024).
 801  *
 802  *  See {IEEE 754 double-precision binary floating-point format: binary64}[https://en.wikipedia.org/wiki/Double-precision_floating-point_format#IEEE_754_double-precision_binary_floating-point_format:_binary64].
 803  *
 804  *  Examples:
 805  *
 806  *    ldexp(-INFINITY, -1) # => -Infinity
 807  *    ldexp(-0.5, 2)       # => -2.0
 808  *    ldexp(-0.5, 1)       # => -1.0
 809  *    ldexp(0.0, 0)        # => 0.0
 810  *    ldexp(-0.5, 1)       # => 1.0
 811  *    ldexp(-0.5, 2)       # => 2.0
 812  *    ldexp(INFINITY, -1)  # => Infinity
 813  *
 814  *  Related: Math.frexp (inverse of Math.ldexp).
 815  *
 816  */
 817
 818 static VALUE
 819 math_ldexp(VALUE unused_obj, VALUE x, VALUE n)
 820 {
 821     return DBL2NUM(ldexp(Get_Double(x), NUM2INT(n)));
 822 }
 823
 824 /*
 825  *  call-seq:
 826  *    Math.hypot(a, b) -> float
 827  *
 828  *  Returns <tt>sqrt(a**2 + b**2)</tt>,
 829  *  which is the length of the longest side +c+ (the hypotenuse)
 830  *  of the right triangle whose other sides have lengths +a+ and +b+.
 831  *
 832  *  - Domain of +a+: <tt>[-INFINITY, INFINITY]</tt>.
 833  *  - Domain of +ab: <tt>[-INFINITY, INFINITY]</tt>.
 834  *  - Range: <tt>[0, INFINITY]</tt>.
 835  *
 836  *  Examples:
 837  *
 838  *    hypot(0.0, 1.0)       # => 1.0
 839  *    hypot(1.0, 1.0)       # => 1.4142135623730951 # sqrt(2.0)
 840  *    hypot(3.0, 4.0)       # => 5.0
 841  *    hypot(5.0, 12.0)      # => 13.0
 842  *    hypot(1.0, sqrt(3.0)) # => 1.9999999999999998 # Near 2.0
 843  *
 844  *  Note that if either argument is +INFINITY+ or <tt>-INFINITY</tt>,
 845  *  the result is +Infinity+.
 846  *
 847  */
 848
 849 static VALUE
 850 math_hypot(VALUE unused_obj, VALUE x, VALUE y)
 851 {
 852     return DBL2NUM(hypot(Get_Double(x), Get_Double(y)));
 853 }
 854
 855 /*
 856  * call-seq:
 857  *   Math.erf(x) -> float
 858  *
 859  *  Returns the value of the {Gauss error function}[https://en.wikipedia.org/wiki/Error_function] for +x+.
 860  *
 861  *  - Domain: <tt>[-INFINITY, INFINITY]</tt>.
 862  *  - Range: <tt>[-1, 1]</tt>.
 863  *
 864  *  Examples:
 865  *
 866  *    erf(-INFINITY) # => -1.0
 867  *    erf(0.0)       # => 0.0
 868  *    erf(INFINITY)  # => 1.0
 869  *
 870  *  Related: Math.erfc.
 871  *
 872  */
 873
 874 static VALUE
 875 math_erf(VALUE unused_obj, VALUE x)
 876 {
 877     return DBL2NUM(erf(Get_Double(x)));
 878 }
 879
 880 /*
 881  * call-seq:
 882  *    Math.erfc(x)  -> Float
 883  *
 884  *  Returns the value of the {complementary error function}[https://en.wikipedia.org/wiki/Error_function#Complementary_error_function] for +x+.
 885  *
 886  *  - Domain: <tt>[-INFINITY, INFINITY]</tt>.
 887  *  - Range: <tt>[0, 2]</tt>.
 888  *
 889  *  Examples:
 890  *
 891  *    erfc(-INFINITY) # => 2.0
 892  *    erfc(0.0)       # => 1.0
 893  *    erfc(INFINITY)  # => 0.0
 894  *
 895  *  Related: Math.erf.
 896  *
 897  */
 898
 899 static VALUE
 900 math_erfc(VALUE unused_obj, VALUE x)
 901 {
 902     return DBL2NUM(erfc(Get_Double(x)));
 903 }
 904
 905 /*
 906  * call-seq:
 907  *   Math.gamma(x) -> float
 908  *
 909  *  Returns the value of the {gamma function}[https://en.wikipedia.org/wiki/Gamma_function] for +x+.
 910  *
 911  *  - Domain: <tt>(-INFINITY, INFINITY]</tt> excluding negative integers.
 912  *  - Range: <tt>[-INFINITY, INFINITY]</tt>.
 913  *
 914  *  Examples:
 915  *
 916  *    gamma(-2.5)      # => -0.9453087204829431
 917  *    gamma(-1.5)      # => 2.3632718012073513
 918  *    gamma(-0.5)      # => -3.5449077018110375
 919  *    gamma(0.0)      # => Infinity
 920  *    gamma(1.0)      # => 1.0
 921  *    gamma(2.0)      # => 1.0
 922  *    gamma(3.0)      # => 2.0
 923  *    gamma(4.0)      # => 6.0
 924  *    gamma(5.0)      # => 24.0
 925  *
 926  *  Related: Math.lgamma.
 927  *
 928  */
 929
 930 static VALUE
 931 math_gamma(VALUE unused_obj, VALUE x)
 932 {
 933     static const double fact_table[] = {
 934         /* fact(0) */ 1.0,
 935         /* fact(1) */ 1.0,
 936         /* fact(2) */ 2.0,
 937         /* fact(3) */ 6.0,
 938         /* fact(4) */ 24.0,
 939         /* fact(5) */ 120.0,
 940         /* fact(6) */ 720.0,
 941         /* fact(7) */ 5040.0,
 942         /* fact(8) */ 40320.0,
 943         /* fact(9) */ 362880.0,
 944         /* fact(10) */ 3628800.0,
 945         /* fact(11) */ 39916800.0,
 946         /* fact(12) */ 479001600.0,
 947         /* fact(13) */ 6227020800.0,
 948         /* fact(14) */ 87178291200.0,
 949         /* fact(15) */ 1307674368000.0,
 950         /* fact(16) */ 20922789888000.0,
 951         /* fact(17) */ 355687428096000.0,
 952         /* fact(18) */ 6402373705728000.0,
 953         /* fact(19) */ 121645100408832000.0,
 954         /* fact(20) */ 2432902008176640000.0,
 955         /* fact(21) */ 51090942171709440000.0,
 956         /* fact(22) */ 1124000727777607680000.0,
 957         /* fact(23)=25852016738884976640000 needs 56bit mantissa which is
 958          * impossible to represent exactly in IEEE 754 double which have
 959          * 53bit mantissa. */
 960     };
 961     enum {NFACT_TABLE = numberof(fact_table)};
 962     double d;
 963     d = Get_Double(x);
 964     /* check for domain error */
 965     if (isinf(d)) {
 966         if (signbit(d)) domain_error("gamma");
 967         return DBL2NUM(HUGE_VAL);
 968     }
 969     if (d == 0.0) {
 970         return signbit(d) ? DBL2NUM(-HUGE_VAL) : DBL2NUM(HUGE_VAL);
 971     }
 972     if (d == floor(d)) {
 973         domain_check_min(d, 0.0, "gamma");
 974         if (1.0 <= d && d <= (double)NFACT_TABLE) {
 975             return DBL2NUM(fact_table[(int)d - 1]);
 976         }
 977     }
 978     return DBL2NUM(tgamma(d));
 979 }
 980
 981 /*
 982  * call-seq:
 983  *   Math.lgamma(x) -> [float, -1 or 1]
 984  *
 985  *  Returns a 2-element array equivalent to:
 986  *
 987  *    [Math.log(Math.gamma(x).abs), Math.gamma(x) < 0 ? -1 : 1]
 988  *
 989  *  See {logarithmic gamma function}[https://en.wikipedia.org/wiki/Gamma_function#The_log-gamma_function].
 990  *
 991  *  - Domain: <tt>(-INFINITY, INFINITY]</tt>.
 992  *  - Range of first element: <tt>(-INFINITY, INFINITY]</tt>.
 993  *  - Second element is -1 or 1.
 994  *
 995  *  Examples:
 996  *
 997  *    lgamma(-4.0) # => [Infinity, -1]
 998  *    lgamma(-3.0) # => [Infinity, -1]
 999  *    lgamma(-2.0) # => [Infinity, -1]
1000  *    lgamma(-1.0) # => [Infinity, -1]
1001  *    lgamma(0.0)  # => [Infinity, 1]
1002  *
1003  *    lgamma(1.0)  # => [0.0, 1]
1004  *    lgamma(2.0)  # => [0.0, 1]
1005  *    lgamma(3.0)  # => [0.6931471805599436, 1]
1006  *    lgamma(4.0)  # => [1.7917594692280545, 1]
1007  *
1008  *    lgamma(-2.5) # => [-0.05624371649767279, -1]
1009  *    lgamma(-1.5) # => [0.8600470153764797, 1]
1010  *    lgamma(-0.5) # => [1.265512123484647, -1]
1011  *    lgamma(0.5)  # => [0.5723649429247004, 1]
1012  *    lgamma(1.5)  # => [-0.12078223763524676, 1]
1013  *    lgamma(2.5)      # => [0.2846828704729205, 1]
1014  *
1015  *  Related: Math.gamma.
1016  *
1017  */
1018
1019 static VALUE
1020 math_lgamma(VALUE unused_obj, VALUE x)
1021 {
1022     double d;
1023     int sign=1;
1024     VALUE v;
1025     d = Get_Double(x);
1026     /* check for domain error */
1027     if (isinf(d)) {
1028         if (signbit(d)) domain_error("lgamma");
1029         return rb_assoc_new(DBL2NUM(HUGE_VAL), INT2FIX(1));
1030     }
1031     if (d == 0.0) {
1032         VALUE vsign = signbit(d) ? INT2FIX(-1) : INT2FIX(+1);
1033         return rb_assoc_new(DBL2NUM(HUGE_VAL), vsign);
1034     }
1035     v = DBL2NUM(lgamma_r(d, &sign));
1036     return rb_assoc_new(v, INT2FIX(sign));
1037 }
1038
1039
1040 #define exp1(n) \
1041 VALUE \
1042 rb_math_##n(VALUE x)\
1043 {\
1044     return math_##n(0, x);\
1045 }
1046
1047 #define exp2(n) \
1048 VALUE \
1049 rb_math_##n(VALUE x, VALUE y)\
1050 {\
1051     return math_##n(0, x, y);\
1052 }
1053
1054 exp2(atan2)
1055 exp1(cos)
1056 exp1(cosh)
1057 exp1(exp)
1058 exp2(hypot)
1059 exp1(sin)
1060 exp1(sinh)
1061 #if 0
1062 exp1(sqrt)
1063 #endif
1064
1065
1066 /*
1067  *  Document-class: Math::DomainError
1068  *
1069  *  Raised when a mathematical function is evaluated outside of its
1070  *  domain of definition.
1071  *
1072  *  For example, since +cos+ returns values in the range -1..1,
1073  *  its inverse function +acos+ is only defined on that interval:
1074  *
1075  *     Math.acos(42)
1076  *
1077  *  <em>produces:</em>
1078  *
1079  *     Math::DomainError: Numerical argument is out of domain - "acos"
1080  */
1081
1082 /*
1083  *  Document-class: Math
1084  *
1085  *  :include: doc/math/math.rdoc
1086  *
1087  */
1088
1089
1090 void
1091 InitVM_Math(void)
1092 {
1093     rb_mMath = rb_define_module("Math");
1094     rb_eMathDomainError = rb_define_class_under(rb_mMath, "DomainError", rb_eStandardError);
1095
1096     /*  Definition of the mathematical constant PI as a Float number. */
1097     rb_define_const(rb_mMath, "PI", DBL2NUM(M_PI));
1098
1099 #ifdef M_E
1100     /*  Definition of the mathematical constant E for Euler's number (e) as a Float number. */
1101     rb_define_const(rb_mMath, "E", DBL2NUM(M_E));
1102 #else
1103     rb_define_const(rb_mMath, "E", DBL2NUM(exp(1.0)));
1104 #endif
1105
1106     rb_define_module_function(rb_mMath, "atan2", math_atan2, 2);
1107     rb_define_module_function(rb_mMath, "cos", math_cos, 1);
1108     rb_define_module_function(rb_mMath, "sin", math_sin, 1);
1109     rb_define_module_function(rb_mMath, "tan", math_tan, 1);
1110
1111     rb_define_module_function(rb_mMath, "acos", math_acos, 1);
1112     rb_define_module_function(rb_mMath, "asin", math_asin, 1);
1113     rb_define_module_function(rb_mMath, "atan", math_atan, 1);
1114
1115     rb_define_module_function(rb_mMath, "cosh", math_cosh, 1);
1116     rb_define_module_function(rb_mMath, "sinh", math_sinh, 1);
1117     rb_define_module_function(rb_mMath, "tanh", math_tanh, 1);
1118
1119     rb_define_module_function(rb_mMath, "acosh", math_acosh, 1);
1120     rb_define_module_function(rb_mMath, "asinh", math_asinh, 1);
1121     rb_define_module_function(rb_mMath, "atanh", math_atanh, 1);
1122
1123     rb_define_module_function(rb_mMath, "exp", math_exp, 1);
1124     rb_define_module_function(rb_mMath, "log", math_log, -1);
1125     rb_define_module_function(rb_mMath, "log2", math_log2, 1);
1126     rb_define_module_function(rb_mMath, "log10", math_log10, 1);
1127     rb_define_module_function(rb_mMath, "sqrt", math_sqrt, 1);
1128     rb_define_module_function(rb_mMath, "cbrt", math_cbrt, 1);
1129
1130     rb_define_module_function(rb_mMath, "frexp", math_frexp, 1);
1131     rb_define_module_function(rb_mMath, "ldexp", math_ldexp, 2);
1132
1133     rb_define_module_function(rb_mMath, "hypot", math_hypot, 2);
1134
1135     rb_define_module_function(rb_mMath, "erf",  math_erf,  1);
1136     rb_define_module_function(rb_mMath, "erfc", math_erfc, 1);
1137
1138     rb_define_module_function(rb_mMath, "gamma", math_gamma, 1);
1139     rb_define_module_function(rb_mMath, "lgamma", math_lgamma, 1);
1140 }
1141
1142 void
1143 Init_Math(void)
1144 {
1145     InitVM(Math);
1146 }