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