Ruby  2.1.10p492(2016-04-01revision54464)
math.c
Go to the documentation of this file.
1 /**********************************************************************
2 
3  math.c -
4 
5  $Author: nobu $
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/ruby.h"
13 #include "internal.h"
14 #include <float.h>
15 #include <math.h>
16 #include <errno.h>
17 
18 #if defined(HAVE_SIGNBIT) && defined(__GNUC__) && defined(__sun) && \
19  !defined(signbit)
20  extern int signbit(double);
21 #endif
22 
23 #define RB_BIGNUM_TYPE_P(x) RB_TYPE_P((x), T_BIGNUM)
24 
27 
28 #define Need_Float(x) do {if (!RB_TYPE_P(x, T_FLOAT)) {(x) = rb_to_float(x);}} while(0)
29 #define Need_Float2(x,y) do {\
30  Need_Float(x);\
31  Need_Float(y);\
32 } while (0)
33 
34 #define domain_error(msg) \
35  rb_raise(rb_eMathDomainError, "Numerical argument is out of domain - " #msg)
36 
37 /*
38  * call-seq:
39  * Math.atan2(y, x) -> Float
40  *
41  * Computes the arc tangent given +y+ and +x+.
42  * Returns a Float in the range -PI..PI.
43  *
44  * Domain: (-INFINITY, INFINITY)
45  *
46  * Codomain: [-PI, PI]
47  *
48  * Math.atan2(-0.0, -1.0) #=> -3.141592653589793
49  * Math.atan2(-1.0, -1.0) #=> -2.356194490192345
50  * Math.atan2(-1.0, 0.0) #=> -1.5707963267948966
51  * Math.atan2(-1.0, 1.0) #=> -0.7853981633974483
52  * Math.atan2(-0.0, 1.0) #=> -0.0
53  * Math.atan2(0.0, 1.0) #=> 0.0
54  * Math.atan2(1.0, 1.0) #=> 0.7853981633974483
55  * Math.atan2(1.0, 0.0) #=> 1.5707963267948966
56  * Math.atan2(1.0, -1.0) #=> 2.356194490192345
57  * Math.atan2(0.0, -1.0) #=> 3.141592653589793
58  *
59  */
60 
61 static VALUE
63 {
64 #ifndef M_PI
65 # define M_PI 3.14159265358979323846
66 #endif
67  double dx, dy;
68  Need_Float2(y, x);
69  dx = RFLOAT_VALUE(x);
70  dy = RFLOAT_VALUE(y);
71  if (dx == 0.0 && dy == 0.0) {
72  if (!signbit(dx))
73  return DBL2NUM(dy);
74  if (!signbit(dy))
75  return DBL2NUM(M_PI);
76  return DBL2NUM(-M_PI);
77  }
78  if (isinf(dx) && isinf(dy)) domain_error("atan2");
79  return DBL2NUM(atan2(dy, dx));
80 }
81 
82 
83 /*
84  * call-seq:
85  * Math.cos(x) -> Float
86  *
87  * Computes the cosine of +x+ (expressed in radians).
88  * Returns a Float in the range -1.0..1.0.
89  *
90  * Domain: (-INFINITY, INFINITY)
91  *
92  * Codomain: [-1, 1]
93  *
94  * Math.cos(Math::PI) #=> -1.0
95  *
96  */
97 
98 static VALUE
100 {
101  Need_Float(x);
102  return DBL2NUM(cos(RFLOAT_VALUE(x)));
103 }
104 
105 /*
106  * call-seq:
107  * Math.sin(x) -> Float
108  *
109  * Computes the sine of +x+ (expressed in radians).
110  * Returns a Float in the range -1.0..1.0.
111  *
112  * Domain: (-INFINITY, INFINITY)
113  *
114  * Codomain: [-1, 1]
115  *
116  * Math.sin(Math::PI/2) #=> 1.0
117  *
118  */
119 
120 static VALUE
122 {
123  Need_Float(x);
124  return DBL2NUM(sin(RFLOAT_VALUE(x)));
125 }
126 
127 
128 /*
129  * call-seq:
130  * Math.tan(x) -> Float
131  *
132  * Computes the tangent of +x+ (expressed in radians).
133  *
134  * Domain: (-INFINITY, INFINITY)
135  *
136  * Codomain: (-INFINITY, INFINITY)
137  *
138  * Math.tan(0) #=> 0.0
139  *
140  */
141 
142 static VALUE
144 {
145  Need_Float(x);
146  return DBL2NUM(tan(RFLOAT_VALUE(x)));
147 }
148 
149 /*
150  * call-seq:
151  * Math.acos(x) -> Float
152  *
153  * Computes the arc cosine of +x+. Returns 0..PI.
154  *
155  * Domain: [-1, 1]
156  *
157  * Codomain: [0, PI]
158  *
159  * Math.acos(0) == Math::PI/2 #=> true
160  *
161  */
162 
163 static VALUE
165 {
166  double d0, d;
167 
168  Need_Float(x);
169  d0 = RFLOAT_VALUE(x);
170  /* check for domain error */
171  if (d0 < -1.0 || 1.0 < d0) domain_error("acos");
172  d = acos(d0);
173  return DBL2NUM(d);
174 }
175 
176 /*
177  * call-seq:
178  * Math.asin(x) -> Float
179  *
180  * Computes the arc sine of +x+. Returns -PI/2..PI/2.
181  *
182  * Domain: [-1, -1]
183  *
184  * Codomain: [-PI/2, PI/2]
185  *
186  * Math.asin(1) == Math::PI/2 #=> true
187  */
188 
189 static VALUE
191 {
192  double d0, d;
193 
194  Need_Float(x);
195  d0 = RFLOAT_VALUE(x);
196  /* check for domain error */
197  if (d0 < -1.0 || 1.0 < d0) domain_error("asin");
198  d = asin(d0);
199  return DBL2NUM(d);
200 }
201 
202 /*
203  * call-seq:
204  * Math.atan(x) -> Float
205  *
206  * Computes the arc tangent of +x+. Returns -PI/2..PI/2.
207  *
208  * Domain: (-INFINITY, INFINITY)
209  *
210  * Codomain: (-PI/2, PI/2)
211  *
212  * Math.atan(0) #=> 0.0
213  */
214 
215 static VALUE
217 {
218  Need_Float(x);
219  return DBL2NUM(atan(RFLOAT_VALUE(x)));
220 }
221 
222 #ifndef HAVE_COSH
223 double
224 cosh(double x)
225 {
226  return (exp(x) + exp(-x)) / 2;
227 }
228 #endif
229 
230 /*
231  * call-seq:
232  * Math.cosh(x) -> Float
233  *
234  * Computes the hyperbolic cosine of +x+ (expressed in radians).
235  *
236  * Domain: (-INFINITY, INFINITY)
237  *
238  * Codomain: [1, INFINITY)
239  *
240  * Math.cosh(0) #=> 1.0
241  *
242  */
243 
244 static VALUE
246 {
247  Need_Float(x);
248  return DBL2NUM(cosh(RFLOAT_VALUE(x)));
249 }
250 
251 #ifndef HAVE_SINH
252 double
253 sinh(double x)
254 {
255  return (exp(x) - exp(-x)) / 2;
256 }
257 #endif
258 
259 /*
260  * call-seq:
261  * Math.sinh(x) -> Float
262  *
263  * Computes the hyperbolic sine of +x+ (expressed in radians).
264  *
265  * Domain: (-INFINITY, INFINITY)
266  *
267  * Codomain: (-INFINITY, INFINITY)
268  *
269  * Math.sinh(0) #=> 0.0
270  *
271  */
272 
273 static VALUE
275 {
276  Need_Float(x);
277  return DBL2NUM(sinh(RFLOAT_VALUE(x)));
278 }
279 
280 #ifndef HAVE_TANH
281 double
282 tanh(double x)
283 {
284  return sinh(x) / cosh(x);
285 }
286 #endif
287 
288 /*
289  * call-seq:
290  * Math.tanh(x) -> Float
291  *
292  * Computes the hyperbolic tangent of +x+ (expressed in radians).
293  *
294  * Domain: (-INFINITY, INFINITY)
295  *
296  * Codomain: (-1, 1)
297  *
298  * Math.tanh(0) #=> 0.0
299  *
300  */
301 
302 static VALUE
304 {
305  Need_Float(x);
306  return DBL2NUM(tanh(RFLOAT_VALUE(x)));
307 }
308 
309 /*
310  * call-seq:
311  * Math.acosh(x) -> Float
312  *
313  * Computes the inverse hyperbolic cosine of +x+.
314  *
315  * Domain: [1, INFINITY)
316  *
317  * Codomain: [0, INFINITY)
318  *
319  * Math.acosh(1) #=> 0.0
320  *
321  */
322 
323 static VALUE
325 {
326  double d0, d;
327 
328  Need_Float(x);
329  d0 = RFLOAT_VALUE(x);
330  /* check for domain error */
331  if (d0 < 1.0) domain_error("acosh");
332  d = acosh(d0);
333  return DBL2NUM(d);
334 }
335 
336 /*
337  * call-seq:
338  * Math.asinh(x) -> Float
339  *
340  * Computes the inverse hyperbolic sine of +x+.
341  *
342  * Domain: (-INFINITY, INFINITY)
343  *
344  * Codomain: (-INFINITY, INFINITY)
345  *
346  * Math.asinh(1) #=> 0.881373587019543
347  *
348  */
349 
350 static VALUE
352 {
353  Need_Float(x);
354  return DBL2NUM(asinh(RFLOAT_VALUE(x)));
355 }
356 
357 /*
358  * call-seq:
359  * Math.atanh(x) -> Float
360  *
361  * Computes the inverse hyperbolic tangent of +x+.
362  *
363  * Domain: (-1, 1)
364  *
365  * Codomain: (-INFINITY, INFINITY)
366  *
367  * Math.atanh(1) #=> Infinity
368  *
369  */
370 
371 static VALUE
373 {
374  double d0, d;
375 
376  Need_Float(x);
377  d0 = RFLOAT_VALUE(x);
378  /* check for domain error */
379  if (d0 < -1.0 || +1.0 < d0) domain_error("atanh");
380  /* check for pole error */
381  if (d0 == -1.0) return DBL2NUM(-INFINITY);
382  if (d0 == +1.0) return DBL2NUM(+INFINITY);
383  d = atanh(d0);
384  return DBL2NUM(d);
385 }
386 
387 /*
388  * call-seq:
389  * Math.exp(x) -> Float
390  *
391  * Returns e**x.
392  *
393  * Domain: (-INFINITY, INFINITY)
394  *
395  * Codomain: (0, INFINITY)
396  *
397  * Math.exp(0) #=> 1.0
398  * Math.exp(1) #=> 2.718281828459045
399  * Math.exp(1.5) #=> 4.4816890703380645
400  *
401  */
402 
403 static VALUE
405 {
406  Need_Float(x);
407  return DBL2NUM(exp(RFLOAT_VALUE(x)));
408 }
409 
410 #if defined __CYGWIN__
411 # include <cygwin/version.h>
412 # if CYGWIN_VERSION_DLL_MAJOR < 1005
413 # define nan(x) nan()
414 # endif
415 # define log(x) ((x) < 0.0 ? nan("") : log(x))
416 # define log10(x) ((x) < 0.0 ? nan("") : log10(x))
417 #endif
418 
419 /*
420  * call-seq:
421  * Math.log(x) -> Float
422  * Math.log(x, base) -> Float
423  *
424  * Returns the logarithm of +x+.
425  * If additional second argument is given, it will be the base
426  * of logarithm. Otherwise it is +e+ (for the natural logarithm).
427  *
428  * Domain: (0, INFINITY)
429  *
430  * Codomain: (-INFINITY, INFINITY)
431  *
432  * Math.log(0) #=> -Infinity
433  * Math.log(1) #=> 0.0
434  * Math.log(Math::E) #=> 1.0
435  * Math.log(Math::E**3) #=> 3.0
436  * Math.log(12, 3) #=> 2.2618595071429146
437  *
438  */
439 
440 static VALUE
442 {
443  VALUE x, base;
444  double d0, d;
445  size_t numbits;
446 
447  rb_scan_args(argc, argv, "11", &x, &base);
448 
449  if (RB_BIGNUM_TYPE_P(x) && RBIGNUM_POSITIVE_P(x) &&
450  DBL_MAX_EXP <= (numbits = rb_absint_numwords(x, 1, NULL))) {
451  numbits -= DBL_MANT_DIG;
452  x = rb_big_rshift(x, SIZET2NUM(numbits));
453  }
454  else {
455  numbits = 0;
456  }
457 
458  Need_Float(x);
459  d0 = RFLOAT_VALUE(x);
460  /* check for domain error */
461  if (d0 < 0.0) domain_error("log");
462  /* check for pole error */
463  if (d0 == 0.0) return DBL2NUM(-INFINITY);
464  d = log(d0);
465  if (numbits)
466  d += numbits * log(2); /* log(2**numbits) */
467  if (argc == 2) {
468  Need_Float(base);
469  d /= log(RFLOAT_VALUE(base));
470  }
471  return DBL2NUM(d);
472 }
473 
474 #ifndef log2
475 #ifndef HAVE_LOG2
476 double
477 log2(double x)
478 {
479  return log10(x)/log10(2.0);
480 }
481 #else
482 extern double log2(double);
483 #endif
484 #endif
485 
486 /*
487  * call-seq:
488  * Math.log2(x) -> Float
489  *
490  * Returns the base 2 logarithm of +x+.
491  *
492  * Domain: (0, INFINITY)
493  *
494  * Codomain: (-INFINITY, INFINITY)
495  *
496  * Math.log2(1) #=> 0.0
497  * Math.log2(2) #=> 1.0
498  * Math.log2(32768) #=> 15.0
499  * Math.log2(65536) #=> 16.0
500  *
501  */
502 
503 static VALUE
505 {
506  double d0, d;
507  size_t numbits;
508 
509  if (RB_BIGNUM_TYPE_P(x) && RBIGNUM_POSITIVE_P(x) &&
510  DBL_MAX_EXP <= (numbits = rb_absint_numwords(x, 1, NULL))) {
511  numbits -= DBL_MANT_DIG;
512  x = rb_big_rshift(x, SIZET2NUM(numbits));
513  }
514  else {
515  numbits = 0;
516  }
517 
518  Need_Float(x);
519  d0 = RFLOAT_VALUE(x);
520  /* check for domain error */
521  if (d0 < 0.0) domain_error("log2");
522  /* check for pole error */
523  if (d0 == 0.0) return DBL2NUM(-INFINITY);
524  d = log2(d0);
525  d += numbits;
526  return DBL2NUM(d);
527 }
528 
529 /*
530  * call-seq:
531  * Math.log10(x) -> Float
532  *
533  * Returns the base 10 logarithm of +x+.
534  *
535  * Domain: (0, INFINITY)
536  *
537  * Codomain: (-INFINITY, INFINITY)
538  *
539  * Math.log10(1) #=> 0.0
540  * Math.log10(10) #=> 1.0
541  * Math.log10(10**100) #=> 100.0
542  *
543  */
544 
545 static VALUE
547 {
548  double d0, d;
549  size_t numbits;
550 
551  if (RB_BIGNUM_TYPE_P(x) && RBIGNUM_POSITIVE_P(x) &&
552  DBL_MAX_EXP <= (numbits = rb_absint_numwords(x, 1, NULL))) {
553  numbits -= DBL_MANT_DIG;
554  x = rb_big_rshift(x, SIZET2NUM(numbits));
555  }
556  else {
557  numbits = 0;
558  }
559 
560  Need_Float(x);
561  d0 = RFLOAT_VALUE(x);
562  /* check for domain error */
563  if (d0 < 0.0) domain_error("log10");
564  /* check for pole error */
565  if (d0 == 0.0) return DBL2NUM(-INFINITY);
566  d = log10(d0);
567  if (numbits)
568  d += numbits * log10(2); /* log10(2**numbits) */
569  return DBL2NUM(d);
570 }
571 
572 /*
573  * call-seq:
574  * Math.sqrt(x) -> Float
575  *
576  * Returns the non-negative square root of +x+.
577  *
578  * Domain: [0, INFINITY)
579  *
580  * Codomain:[0, INFINITY)
581  *
582  * 0.upto(10) {|x|
583  * p [x, Math.sqrt(x), Math.sqrt(x)**2]
584  * }
585  * #=> [0, 0.0, 0.0]
586  * # [1, 1.0, 1.0]
587  * # [2, 1.4142135623731, 2.0]
588  * # [3, 1.73205080756888, 3.0]
589  * # [4, 2.0, 4.0]
590  * # [5, 2.23606797749979, 5.0]
591  * # [6, 2.44948974278318, 6.0]
592  * # [7, 2.64575131106459, 7.0]
593  * # [8, 2.82842712474619, 8.0]
594  * # [9, 3.0, 9.0]
595  * # [10, 3.16227766016838, 10.0]
596  */
597 
598 static VALUE
600 {
601  double d0, d;
602 
603  Need_Float(x);
604  d0 = RFLOAT_VALUE(x);
605  /* check for domain error */
606  if (d0 < 0.0) domain_error("sqrt");
607  if (d0 == 0.0) return DBL2NUM(0.0);
608  d = sqrt(d0);
609  return DBL2NUM(d);
610 }
611 
612 /*
613  * call-seq:
614  * Math.cbrt(x) -> Float
615  *
616  * Returns the cube root of +x+.
617  *
618  * Domain: [0, INFINITY)
619  *
620  * Codomain:[0, INFINITY)
621  *
622  * -9.upto(9) {|x|
623  * p [x, Math.cbrt(x), Math.cbrt(x)**3]
624  * }
625  * #=> [-9, -2.0800838230519, -9.0]
626  * # [-8, -2.0, -8.0]
627  * # [-7, -1.91293118277239, -7.0]
628  * # [-6, -1.81712059283214, -6.0]
629  * # [-5, -1.7099759466767, -5.0]
630  * # [-4, -1.5874010519682, -4.0]
631  * # [-3, -1.44224957030741, -3.0]
632  * # [-2, -1.25992104989487, -2.0]
633  * # [-1, -1.0, -1.0]
634  * # [0, 0.0, 0.0]
635  * # [1, 1.0, 1.0]
636  * # [2, 1.25992104989487, 2.0]
637  * # [3, 1.44224957030741, 3.0]
638  * # [4, 1.5874010519682, 4.0]
639  * # [5, 1.7099759466767, 5.0]
640  * # [6, 1.81712059283214, 6.0]
641  * # [7, 1.91293118277239, 7.0]
642  * # [8, 2.0, 8.0]
643  * # [9, 2.0800838230519, 9.0]
644  *
645  */
646 
647 static VALUE
649 {
650  Need_Float(x);
651  return DBL2NUM(cbrt(RFLOAT_VALUE(x)));
652 }
653 
654 /*
655  * call-seq:
656  * Math.frexp(x) -> [fraction, exponent]
657  *
658  * Returns a two-element array containing the normalized fraction (a Float)
659  * and exponent (a Fixnum) of +x+.
660  *
661  * fraction, exponent = Math.frexp(1234) #=> [0.6025390625, 11]
662  * fraction * 2**exponent #=> 1234.0
663  */
664 
665 static VALUE
667 {
668  double d;
669  int exp;
670 
671  Need_Float(x);
672 
673  d = frexp(RFLOAT_VALUE(x), &exp);
674  return rb_assoc_new(DBL2NUM(d), INT2NUM(exp));
675 }
676 
677 /*
678  * call-seq:
679  * Math.ldexp(fraction, exponent) -> float
680  *
681  * Returns the value of +fraction+*(2**+exponent+).
682  *
683  * fraction, exponent = Math.frexp(1234)
684  * Math.ldexp(fraction, exponent) #=> 1234.0
685  */
686 
687 static VALUE
689 {
690  Need_Float(x);
691  return DBL2NUM(ldexp(RFLOAT_VALUE(x), NUM2INT(n)));
692 }
693 
694 /*
695  * call-seq:
696  * Math.hypot(x, y) -> Float
697  *
698  * Returns sqrt(x**2 + y**2), the hypotenuse of a right-angled triangle with
699  * sides +x+ and +y+.
700  *
701  * Math.hypot(3, 4) #=> 5.0
702  */
703 
704 static VALUE
706 {
707  Need_Float2(x, y);
708  return DBL2NUM(hypot(RFLOAT_VALUE(x), RFLOAT_VALUE(y)));
709 }
710 
711 /*
712  * call-seq:
713  * Math.erf(x) -> Float
714  *
715  * Calculates the error function of +x+.
716  *
717  * Domain: (-INFINITY, INFINITY)
718  *
719  * Codomain: (-1, 1)
720  *
721  * Math.erf(0) #=> 0.0
722  *
723  */
724 
725 static VALUE
727 {
728  Need_Float(x);
729  return DBL2NUM(erf(RFLOAT_VALUE(x)));
730 }
731 
732 /*
733  * call-seq:
734  * Math.erfc(x) -> Float
735  *
736  * Calculates the complementary error function of x.
737  *
738  * Domain: (-INFINITY, INFINITY)
739  *
740  * Codomain: (0, 2)
741  *
742  * Math.erfc(0) #=> 1.0
743  *
744  */
745 
746 static VALUE
748 {
749  Need_Float(x);
750  return DBL2NUM(erfc(RFLOAT_VALUE(x)));
751 }
752 
753 /*
754  * call-seq:
755  * Math.gamma(x) -> Float
756  *
757  * Calculates the gamma function of x.
758  *
759  * Note that gamma(n) is same as fact(n-1) for integer n > 0.
760  * However gamma(n) returns float and can be an approximation.
761  *
762  * def fact(n) (1..n).inject(1) {|r,i| r*i } end
763  * 1.upto(26) {|i| p [i, Math.gamma(i), fact(i-1)] }
764  * #=> [1, 1.0, 1]
765  * # [2, 1.0, 1]
766  * # [3, 2.0, 2]
767  * # [4, 6.0, 6]
768  * # [5, 24.0, 24]
769  * # [6, 120.0, 120]
770  * # [7, 720.0, 720]
771  * # [8, 5040.0, 5040]
772  * # [9, 40320.0, 40320]
773  * # [10, 362880.0, 362880]
774  * # [11, 3628800.0, 3628800]
775  * # [12, 39916800.0, 39916800]
776  * # [13, 479001600.0, 479001600]
777  * # [14, 6227020800.0, 6227020800]
778  * # [15, 87178291200.0, 87178291200]
779  * # [16, 1307674368000.0, 1307674368000]
780  * # [17, 20922789888000.0, 20922789888000]
781  * # [18, 355687428096000.0, 355687428096000]
782  * # [19, 6.402373705728e+15, 6402373705728000]
783  * # [20, 1.21645100408832e+17, 121645100408832000]
784  * # [21, 2.43290200817664e+18, 2432902008176640000]
785  * # [22, 5.109094217170944e+19, 51090942171709440000]
786  * # [23, 1.1240007277776077e+21, 1124000727777607680000]
787  * # [24, 2.5852016738885062e+22, 25852016738884976640000]
788  * # [25, 6.204484017332391e+23, 620448401733239439360000]
789  * # [26, 1.5511210043330954e+25, 15511210043330985984000000]
790  *
791  */
792 
793 static VALUE
795 {
796  static const double fact_table[] = {
797  /* fact(0) */ 1.0,
798  /* fact(1) */ 1.0,
799  /* fact(2) */ 2.0,
800  /* fact(3) */ 6.0,
801  /* fact(4) */ 24.0,
802  /* fact(5) */ 120.0,
803  /* fact(6) */ 720.0,
804  /* fact(7) */ 5040.0,
805  /* fact(8) */ 40320.0,
806  /* fact(9) */ 362880.0,
807  /* fact(10) */ 3628800.0,
808  /* fact(11) */ 39916800.0,
809  /* fact(12) */ 479001600.0,
810  /* fact(13) */ 6227020800.0,
811  /* fact(14) */ 87178291200.0,
812  /* fact(15) */ 1307674368000.0,
813  /* fact(16) */ 20922789888000.0,
814  /* fact(17) */ 355687428096000.0,
815  /* fact(18) */ 6402373705728000.0,
816  /* fact(19) */ 121645100408832000.0,
817  /* fact(20) */ 2432902008176640000.0,
818  /* fact(21) */ 51090942171709440000.0,
819  /* fact(22) */ 1124000727777607680000.0,
820  /* fact(23)=25852016738884976640000 needs 56bit mantissa which is
821  * impossible to represent exactly in IEEE 754 double which have
822  * 53bit mantissa. */
823  };
824  double d0, d;
825  double intpart, fracpart;
826  Need_Float(x);
827  d0 = RFLOAT_VALUE(x);
828  /* check for domain error */
829  if (isinf(d0) && signbit(d0)) domain_error("gamma");
830  fracpart = modf(d0, &intpart);
831  if (fracpart == 0.0) {
832  if (intpart < 0) domain_error("gamma");
833  if (0 < intpart &&
834  intpart - 1 < (double)numberof(fact_table)) {
835  return DBL2NUM(fact_table[(int)intpart - 1]);
836  }
837  }
838  d = tgamma(d0);
839  return DBL2NUM(d);
840 }
841 
842 /*
843  * call-seq:
844  * Math.lgamma(x) -> [float, -1 or 1]
845  *
846  * Calculates the logarithmic gamma of +x+ and the sign of gamma of +x+.
847  *
848  * Math.lgamma(x) is same as
849  * [Math.log(Math.gamma(x).abs), Math.gamma(x) < 0 ? -1 : 1]
850  * but avoid overflow by Math.gamma(x) for large x.
851  *
852  * Math.lgamma(0) #=> [Infinity, 1]
853  *
854  */
855 
856 static VALUE
858 {
859  double d0, d;
860  int sign=1;
861  VALUE v;
862  Need_Float(x);
863  d0 = RFLOAT_VALUE(x);
864  /* check for domain error */
865  if (isinf(d0)) {
866  if (signbit(d0)) domain_error("lgamma");
867  return rb_assoc_new(DBL2NUM(INFINITY), INT2FIX(1));
868  }
869  d = lgamma_r(d0, &sign);
870  v = DBL2NUM(d);
871  return rb_assoc_new(v, INT2FIX(sign));
872 }
873 
874 
875 #define exp1(n) \
876 VALUE \
877 rb_math_##n(VALUE x)\
878 {\
879  return math_##n(rb_mMath, x);\
880 }
881 
882 #define exp2(n) \
883 VALUE \
884 rb_math_##n(VALUE x, VALUE y)\
885 {\
886  return math_##n(rb_mMath, x, y);\
887 }
888 
889 exp2(atan2)
890 exp1(cos)
891 exp1(cosh)
892 exp1(exp)
893 exp2(hypot)
894 
895 VALUE
896 rb_math_log(int argc, VALUE *argv)
897 {
898  return math_log(argc, argv);
899 }
900 
901 exp1(sin)
902 exp1(sinh)
903 exp1(sqrt)
904 
905 
906 /*
907  * Document-class: Math::DomainError
908  *
909  * Raised when a mathematical function is evaluated outside of its
910  * domain of definition.
911  *
912  * For example, since +cos+ returns values in the range -1..1,
913  * its inverse function +acos+ is only defined on that interval:
914  *
915  * Math.acos(42)
916  *
917  * <em>produces:</em>
918  *
919  * Math::DomainError: Numerical argument is out of domain - "acos"
920  */
921 
922 /*
923  * Document-class: Math
924  *
925  * The Math module contains module functions for basic
926  * trigonometric and transcendental functions. See class
927  * Float for a list of constants that
928  * define Ruby's floating point accuracy.
929  *
930  * Domains and codomains are given only for real (not complex) numbers.
931  */
932 
933 
934 void
935 Init_Math(void)
936 {
937  rb_mMath = rb_define_module("Math");
939 
940 #ifdef M_PI
941  /* Definition of the mathematical constant PI as a Float number. */
943 #else
944  rb_define_const(rb_mMath, "PI", DBL2NUM(atan(1.0)*4.0));
945 #endif
946 
947 #ifdef M_E
948  /* Definition of the mathematical constant E (e) as a Float number. */
949  rb_define_const(rb_mMath, "E", DBL2NUM(M_E));
950 #else
951  rb_define_const(rb_mMath, "E", DBL2NUM(exp(1.0)));
952 #endif
953 
958 
962 
966 
970 
977 
980 
982 
985 
988 }
#define d0
static VALUE math_asinh(VALUE obj, VALUE x)
Definition: math.c:351
VALUE rb_eStandardError
Definition: error.c:546
RUBY_EXTERN double cbrt(double)
Definition: cbrt.c:4
#define exp1(n)
Definition: math.c:875
double sinh(double x)
Definition: math.c:253
#define INT2NUM(x)
Definition: ruby.h:1296
#define NUM2INT(x)
Definition: ruby.h:630
RUBY_EXTERN int signbit(double x)
Definition: signbit.c:5
double log2(double x)
Definition: math.c:477
static VALUE math_atan2(VALUE obj, VALUE y, VALUE x)
Definition: math.c:62
static VALUE math_atanh(VALUE obj, VALUE x)
Definition: math.c:372
VALUE rb_eMathDomainError
Definition: math.c:26
VALUE rb_mMath
Definition: math.c:25
VALUE rb_define_class_under(VALUE outer, const char *name, VALUE super)
Defines a class under the namespace of outer.
Definition: class.c:657
double cosh(double x)
Definition: math.c:224
RUBY_EXTERN double tgamma(double)
Definition: tgamma.c:72
#define RB_BIGNUM_TYPE_P(x)
Definition: math.c:23
static VALUE math_ldexp(VALUE obj, VALUE x, VALUE n)
Definition: math.c:688
RUBY_EXTERN double lgamma_r(double, int *)
Definition: lgamma_r.c:63
static VALUE math_exp(VALUE obj, VALUE x)
Definition: math.c:404
VALUE rb_math_log(int argc, VALUE *argv)
#define RBIGNUM_POSITIVE_P(b)
Definition: ruby.h:1097
static VALUE math_log10(VALUE obj, VALUE x)
Definition: math.c:546
double tanh(double x)
Definition: math.c:282
void rb_define_const(VALUE, const char *, VALUE)
Definition: variable.c:2228
RUBY_SYMBOL_EXPORT_BEGIN RUBY_EXTERN double acosh(double)
Definition: acosh.c:36
int argc
Definition: ruby.c:131
RUBY_EXTERN double erfc(double)
Definition: erf.c:81
static VALUE math_tanh(VALUE obj, VALUE x)
Definition: math.c:303
RUBY_EXTERN int isinf(double)
Definition: isinf.c:56
#define numberof(array)
Definition: etc.c:602
void rb_define_module_function(VALUE module, const char *name, VALUE(*func)(ANYARGS), int argc)
Defines a module function for module.
Definition: class.c:1661
static VALUE math_sqrt(VALUE obj, VALUE x)
Definition: math.c:599
RUBY_EXTERN double atanh(double)
Definition: acosh.c:75
#define Need_Float2(x, y)
Definition: math.c:29
RUBY_EXTERN double hypot(double, double)
Definition: hypot.c:6
int rb_scan_args(int argc, const VALUE *argv, const char *fmt,...)
Definition: class.c:1719
static VALUE math_erf(VALUE obj, VALUE x)
Definition: math.c:726
VALUE rb_assoc_new(VALUE car, VALUE cdr)
Definition: array.c:620
static VALUE math_log(int argc, VALUE *argv)
Definition: math.c:441
#define domain_error(msg)
Definition: math.c:34
unsigned long VALUE
Definition: ruby.h:88
RUBY_EXTERN double asinh(double)
Definition: acosh.c:52
#define INFINITY
Definition: missing.h:141
#define DBL_MANT_DIG
Definition: acosh.c:19
#define DBL_MAX_EXP
Definition: numeric.c:58
#define RFLOAT_VALUE(v)
Definition: ruby.h:814
#define M_PI
#define INT2FIX(i)
Definition: ruby.h:231
static VALUE math_cos(VALUE obj, VALUE x)
Definition: math.c:99
static VALUE math_acos(VALUE obj, VALUE x)
Definition: math.c:164
static VALUE math_atan(VALUE obj, VALUE x)
Definition: math.c:216
static VALUE math_erfc(VALUE obj, VALUE x)
Definition: math.c:747
static VALUE math_acosh(VALUE obj, VALUE x)
Definition: math.c:324
static VALUE math_cbrt(VALUE obj, VALUE x)
Definition: math.c:648
static VALUE math_log2(VALUE obj, VALUE x)
Definition: math.c:504
static VALUE math_frexp(VALUE obj, VALUE x)
Definition: math.c:666
size_t rb_absint_numwords(VALUE val, size_t word_numbits, size_t *nlz_bits_ret)
Definition: bignum.c:3366
static VALUE math_sin(VALUE obj, VALUE x)
Definition: math.c:121
static VALUE math_asin(VALUE obj, VALUE x)
Definition: math.c:190
VALUE rb_big_rshift(VALUE x, VALUE y)
Definition: bignum.c:6798
static VALUE math_tan(VALUE obj, VALUE x)
Definition: math.c:143
VALUE rb_define_module(const char *name)
Definition: class.c:727
#define NULL
Definition: _sdbm.c:102
static VALUE math_gamma(VALUE obj, VALUE x)
Definition: math.c:794
#define exp2(n)
Definition: math.c:882
static VALUE math_hypot(VALUE obj, VALUE x, VALUE y)
Definition: math.c:705
RUBY_EXTERN double erf(double)
Definition: erf.c:71
#define SIZET2NUM(v)
Definition: ruby.h:262
static VALUE math_cosh(VALUE obj, VALUE x)
Definition: math.c:245
char ** argv
Definition: ruby.c:132
#define DBL2NUM(dbl)
Definition: ruby.h:815
static VALUE math_lgamma(VALUE obj, VALUE x)
Definition: math.c:857
static VALUE math_sinh(VALUE obj, VALUE x)
Definition: math.c:274
#define Need_Float(x)
Definition: math.c:28