Ruby  3.4.0dev (2024-11-22 revision 37a72b0150ec36b4ea27175039afc28c62207b0c)
math.c (37a72b0150ec36b4ea27175039afc28c62207b0c)
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 
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");
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 }
VALUE rb_define_class_under(VALUE outer, const char *name, VALUE super)
Defines a class under the namespace of outer.
Definition: class.c:1012
VALUE rb_define_module(const char *name)
Defines a top-level module.
Definition: class.c:1095
void rb_define_module_function(VALUE module, const char *name, VALUE(*func)(ANYARGS), int argc)
Defines a module function for a module.
Definition: class.c:2329
int rb_scan_args(int argc, const VALUE *argv, const char *fmt,...)
Retrieves argument from argc and argv to given VALUE references according to the format string.
Definition: class.c:2635
#define T_COMPLEX
Old name of RUBY_T_COMPLEX.
Definition: value_type.h:59
#define RFLOAT_VALUE
Old name of rb_float_value.
Definition: double.h:28
#define INT2FIX
Old name of RB_INT2FIX.
Definition: long.h:48
#define SIZET2NUM
Old name of RB_SIZE2NUM.
Definition: size_t.h:62
#define NUM2INT
Old name of RB_NUM2INT.
Definition: int.h:44
#define INT2NUM
Old name of RB_INT2NUM.
Definition: int.h:43
#define FIX2LONG
Old name of RB_FIX2LONG.
Definition: long.h:46
#define DBL2NUM
Old name of rb_float_new.
Definition: double.h:29
#define FIXNUM_P
Old name of RB_FIXNUM_P.
VALUE rb_eStandardError
StandardError exception.
Definition: error.c:1405
VALUE rb_eMathDomainError
Math::DomainError exception.
Definition: math.c:30
VALUE rb_mMath
Math module.
Definition: math.c:29
VALUE rb_funcall(VALUE recv, ID mid, int n,...)
Calls a method.
Definition: vm_eval.c:1099
VALUE rb_assoc_new(VALUE car, VALUE cdr)
Identical to rb_ary_new_from_values(), except it expects exactly two parameters.
Definition: array.c:1001
size_t rb_absint_numwords(VALUE val, size_t word_numbits, size_t *nlz_bits_ret)
Calculates the number of words needed represent the absolute value of the passed integer.
Definition: bignum.c:3424
VALUE rb_big_rshift(VALUE x, VALUE y)
Performs shift right.
Definition: bignum.c:6686
VALUE rb_complex_new(VALUE real, VALUE imag)
Constructs a Complex, by first multiplying the imaginary part with 1i then adds it to the real part.
Definition: complex.c:1755
VALUE rb_complex_abs(VALUE z)
Queries the absolute (or the magnitude) of the passed object.
Definition: complex.c:1331
void rb_define_const(VALUE klass, const char *name, VALUE val)
Defines a Ruby level constant under a namespace.
Definition: variable.c:3714
#define InitVM(ext)
This macro is for internal use.
Definition: ruby.h:231
uintptr_t VALUE
Type that represents a Ruby object.
Definition: value.h:40
static bool RB_FLOAT_TYPE_P(VALUE obj)
Queries if the object is an instance of rb_cFloat.
Definition: value_type.h:264
static bool RB_TYPE_P(VALUE obj, enum ruby_value_type t)
Queries if the given object is of given type.
Definition: value_type.h:376