Ruby: math.c Source File

00001 /**********************************************************************
00002 
00003   math.c -
00004 
00005   $Author: nobu $
00006   created at: Tue Jan 25 14:12:56 JST 1994
00007 
00008   Copyright (C) 1993-2007 Yukihiro Matsumoto
00009 
00010 **********************************************************************/
00011 
00012 #include "ruby/ruby.h"
00013 #include "internal.h"
00014 #include <float.h>
00015 #include <math.h>
00016 #include <errno.h>
00017 
00018 #if defined(HAVE_SIGNBIT) && defined(__GNUC__) && defined(__sun) && \
00019     !defined(signbit)
00020     extern int signbit(double);
00021 #endif
00022 
00023 #define RB_BIGNUM_TYPE_P(x) RB_TYPE_P((x), T_BIGNUM)
00024 
00025 VALUE rb_mMath;
00026 VALUE rb_eMathDomainError;
00027 
00028 #define Need_Float(x) do {if (!RB_TYPE_P(x, T_FLOAT)) {(x) = rb_to_float(x);}} while(0)
00029 #define Need_Float2(x,y) do {\
00030     Need_Float(x);\
00031     Need_Float(y);\
00032 } while (0)
00033 
00034 #define domain_error(msg) \
00035     rb_raise(rb_eMathDomainError, "Numerical argument is out of domain - " #msg)
00036 
00037 /*
00038  *  call-seq:
00039  *     Math.atan2(y, x)  -> Float
00040  *
00041  *  Computes the arc tangent given +y+ and +x+.
00042  *  Returns a Float in the range -PI..PI.
00043  *
00044  *  Domain: (-INFINITY, INFINITY)
00045  *
00046  *  Codomain: [-PI, PI]
00047  *
00048  *    Math.atan2(-0.0, -1.0) #=> -3.141592653589793
00049  *    Math.atan2(-1.0, -1.0) #=> -2.356194490192345
00050  *    Math.atan2(-1.0, 0.0)  #=> -1.5707963267948966
00051  *    Math.atan2(-1.0, 1.0)  #=> -0.7853981633974483
00052  *    Math.atan2(-0.0, 1.0)  #=> -0.0
00053  *    Math.atan2(0.0, 1.0)   #=> 0.0
00054  *    Math.atan2(1.0, 1.0)   #=> 0.7853981633974483
00055  *    Math.atan2(1.0, 0.0)   #=> 1.5707963267948966
00056  *    Math.atan2(1.0, -1.0)  #=> 2.356194490192345
00057  *    Math.atan2(0.0, -1.0)  #=> 3.141592653589793
00058  *
00059  */
00060 
00061 static VALUE
00062 math_atan2(VALUE obj, VALUE y, VALUE x)
00063 {
00064 #ifndef M_PI
00065 # define M_PI 3.14159265358979323846
00066 #endif
00067     double dx, dy;
00068     Need_Float2(y, x);
00069     dx = RFLOAT_VALUE(x);
00070     dy = RFLOAT_VALUE(y);
00071     if (dx == 0.0 && dy == 0.0) {
00072         if (!signbit(dx))
00073             return DBL2NUM(dy);
00074         if (!signbit(dy))
00075             return DBL2NUM(M_PI);
00076         return DBL2NUM(-M_PI);
00077     }
00078     if (isinf(dx) && isinf(dy)) domain_error("atan2");
00079     return DBL2NUM(atan2(dy, dx));
00080 }
00081 
00082 
00083 /*
00084  *  call-seq:
00085  *     Math.cos(x)    -> Float
00086  *
00087  *  Computes the cosine of +x+ (expressed in radians).
00088  *  Returns a Float in the range -1.0..1.0.
00089  *
00090  *  Domain: (-INFINITY, INFINITY)
00091  *
00092  *  Codomain: [-1, 1]
00093  *
00094  *    Math.cos(Math::PI) #=> -1.0
00095  *
00096  */
00097 
00098 static VALUE
00099 math_cos(VALUE obj, VALUE x)
00100 {
00101     Need_Float(x);
00102     return DBL2NUM(cos(RFLOAT_VALUE(x)));
00103 }
00104 
00105 /*
00106  *  call-seq:
00107  *     Math.sin(x)    -> Float
00108  *
00109  *  Computes the sine of +x+ (expressed in radians).
00110  *  Returns a Float in the range -1.0..1.0.
00111  *
00112  *  Domain: (-INFINITY, INFINITY)
00113  *
00114  *  Codomain: [-1, 1]
00115  *
00116  *    Math.sin(Math::PI/2) #=> 1.0
00117  *
00118  */
00119 
00120 static VALUE
00121 math_sin(VALUE obj, VALUE x)
00122 {
00123     Need_Float(x);
00124     return DBL2NUM(sin(RFLOAT_VALUE(x)));
00125 }
00126 
00127 
00128 /*
00129  *  call-seq:
00130  *     Math.tan(x)    -> Float
00131  *
00132  *  Computes the tangent of +x+ (expressed in radians).
00133  *
00134  *  Domain: (-INFINITY, INFINITY)
00135  *
00136  *  Codomain: (-INFINITY, INFINITY)
00137  *
00138  *    Math.tan(0) #=> 0.0
00139  *
00140  */
00141 
00142 static VALUE
00143 math_tan(VALUE obj, VALUE x)
00144 {
00145     Need_Float(x);
00146     return DBL2NUM(tan(RFLOAT_VALUE(x)));
00147 }
00148 
00149 /*
00150  *  call-seq:
00151  *     Math.acos(x)    -> Float
00152  *
00153  *  Computes the arc cosine of +x+. Returns 0..PI.
00154  *
00155  *  Domain: [-1, 1]
00156  *
00157  *  Codomain: [0, PI]
00158  *
00159  *    Math.acos(0) == Math::PI/2  #=> true
00160  *
00161  */
00162 
00163 static VALUE
00164 math_acos(VALUE obj, VALUE x)
00165 {
00166     double d0, d;
00167 
00168     Need_Float(x);
00169     d0 = RFLOAT_VALUE(x);
00170     /* check for domain error */
00171     if (d0 < -1.0 || 1.0 < d0) domain_error("acos");
00172     d = acos(d0);
00173     return DBL2NUM(d);
00174 }
00175 
00176 /*
00177  *  call-seq:
00178  *     Math.asin(x)    -> Float
00179  *
00180  *  Computes the arc sine of +x+. Returns -PI/2..PI/2.
00181  *
00182  *  Domain: [-1, -1]
00183  *
00184  *  Codomain: [-PI/2, PI/2]
00185  *
00186  *    Math.asin(1) == Math::PI/2  #=> true
00187  */
00188 
00189 static VALUE
00190 math_asin(VALUE obj, VALUE x)
00191 {
00192     double d0, d;
00193 
00194     Need_Float(x);
00195     d0 = RFLOAT_VALUE(x);
00196     /* check for domain error */
00197     if (d0 < -1.0 || 1.0 < d0) domain_error("asin");
00198     d = asin(d0);
00199     return DBL2NUM(d);
00200 }
00201 
00202 /*
00203  *  call-seq:
00204  *     Math.atan(x)    -> Float
00205  *
00206  *  Computes the arc tangent of +x+. Returns -PI/2..PI/2.
00207  *
00208  *  Domain: (-INFINITY, INFINITY)
00209  *
00210  *  Codomain: (-PI/2, PI/2)
00211  *
00212  *    Math.atan(0) #=> 0.0
00213  */
00214 
00215 static VALUE
00216 math_atan(VALUE obj, VALUE x)
00217 {
00218     Need_Float(x);
00219     return DBL2NUM(atan(RFLOAT_VALUE(x)));
00220 }
00221 
00222 #ifndef HAVE_COSH
00223 double
00224 cosh(double x)
00225 {
00226     return (exp(x) + exp(-x)) / 2;
00227 }
00228 #endif
00229 
00230 /*
00231  *  call-seq:
00232  *     Math.cosh(x)    -> Float
00233  *
00234  *  Computes the hyperbolic cosine of +x+ (expressed in radians).
00235  *
00236  *  Domain: (-INFINITY, INFINITY)
00237  *
00238  *  Codomain: [1, INFINITY)
00239  *
00240  *    Math.cosh(0) #=> 1.0
00241  *
00242  */
00243 
00244 static VALUE
00245 math_cosh(VALUE obj, VALUE x)
00246 {
00247     Need_Float(x);
00248     return DBL2NUM(cosh(RFLOAT_VALUE(x)));
00249 }
00250 
00251 #ifndef HAVE_SINH
00252 double
00253 sinh(double x)
00254 {
00255     return (exp(x) - exp(-x)) / 2;
00256 }
00257 #endif
00258 
00259 /*
00260  *  call-seq:
00261  *     Math.sinh(x)    -> Float
00262  *
00263  *  Computes the hyperbolic sine of +x+ (expressed in radians).
00264  *
00265  *  Domain: (-INFINITY, INFINITY)
00266  *
00267  *  Codomain: (-INFINITY, INFINITY)
00268  *
00269  *    Math.sinh(0) #=> 0.0
00270  *
00271  */
00272 
00273 static VALUE
00274 math_sinh(VALUE obj, VALUE x)
00275 {
00276     Need_Float(x);
00277     return DBL2NUM(sinh(RFLOAT_VALUE(x)));
00278 }
00279 
00280 #ifndef HAVE_TANH
00281 double
00282 tanh(double x)
00283 {
00284     return sinh(x) / cosh(x);
00285 }
00286 #endif
00287 
00288 /*
00289  *  call-seq:
00290  *     Math.tanh(x)    -> Float
00291  *
00292  *  Computes the hyperbolic tangent of +x+ (expressed in radians).
00293  *
00294  *  Domain: (-INFINITY, INFINITY)
00295  *
00296  *  Codomain: (-1, 1)
00297  *
00298  *    Math.tanh(0) #=> 0.0
00299  *
00300  */
00301 
00302 static VALUE
00303 math_tanh(VALUE obj, VALUE x)
00304 {
00305     Need_Float(x);
00306     return DBL2NUM(tanh(RFLOAT_VALUE(x)));
00307 }
00308 
00309 /*
00310  *  call-seq:
00311  *     Math.acosh(x)    -> Float
00312  *
00313  *  Computes the inverse hyperbolic cosine of +x+.
00314  *
00315  *  Domain: [1, INFINITY)
00316  *
00317  *  Codomain: [0, INFINITY)
00318  *
00319  *    Math.acosh(1) #=> 0.0
00320  *
00321  */
00322 
00323 static VALUE
00324 math_acosh(VALUE obj, VALUE x)
00325 {
00326     double d0, d;
00327 
00328     Need_Float(x);
00329     d0 = RFLOAT_VALUE(x);
00330     /* check for domain error */
00331     if (d0 < 1.0) domain_error("acosh");
00332     d = acosh(d0);
00333     return DBL2NUM(d);
00334 }
00335 
00336 /*
00337  *  call-seq:
00338  *     Math.asinh(x)    -> Float
00339  *
00340  *  Computes the inverse hyperbolic sine of +x+.
00341  *
00342  *  Domain: (-INFINITY, INFINITY)
00343  *
00344  *  Codomain: (-INFINITY, INFINITY)
00345  *
00346  *    Math.asinh(1) #=> 0.881373587019543
00347  *
00348  */
00349 
00350 static VALUE
00351 math_asinh(VALUE obj, VALUE x)
00352 {
00353     Need_Float(x);
00354     return DBL2NUM(asinh(RFLOAT_VALUE(x)));
00355 }
00356 
00357 /*
00358  *  call-seq:
00359  *     Math.atanh(x)    -> Float
00360  *
00361  *  Computes the inverse hyperbolic tangent of +x+.
00362  *
00363  *  Domain: (-1, 1)
00364  *
00365  *  Codomain: (-INFINITY, INFINITY)
00366  *
00367  *   Math.atanh(1) #=> Infinity
00368  *
00369  */
00370 
00371 static VALUE
00372 math_atanh(VALUE obj, VALUE x)
00373 {
00374     double d0, d;
00375 
00376     Need_Float(x);
00377     d0 = RFLOAT_VALUE(x);
00378     /* check for domain error */
00379     if (d0 <  -1.0 || +1.0 <  d0) domain_error("atanh");
00380     /* check for pole error */
00381     if (d0 == -1.0) return DBL2NUM(-INFINITY);
00382     if (d0 == +1.0) return DBL2NUM(+INFINITY);
00383     d = atanh(d0);
00384     return DBL2NUM(d);
00385 }
00386 
00387 /*
00388  *  call-seq:
00389  *     Math.exp(x)    -> Float
00390  *
00391  *  Returns e**x.
00392  *
00393  *  Domain: (-INFINITY, INFINITY)
00394  *
00395  *  Codomain: (0, INFINITY)
00396  *
00397  *    Math.exp(0)       #=> 1.0
00398  *    Math.exp(1)       #=> 2.718281828459045
00399  *    Math.exp(1.5)     #=> 4.4816890703380645
00400  *
00401  */
00402 
00403 static VALUE
00404 math_exp(VALUE obj, VALUE x)
00405 {
00406     Need_Float(x);
00407     return DBL2NUM(exp(RFLOAT_VALUE(x)));
00408 }
00409 
00410 #if defined __CYGWIN__
00411 # include <cygwin/version.h>
00412 # if CYGWIN_VERSION_DLL_MAJOR < 1005
00413 #  define nan(x) nan()
00414 # endif
00415 # define log(x) ((x) < 0.0 ? nan("") : log(x))
00416 # define log10(x) ((x) < 0.0 ? nan("") : log10(x))
00417 #endif
00418 
00419 /*
00420  *  call-seq:
00421  *     Math.log(x)          -> Float
00422  *     Math.log(x, base)    -> Float
00423  *
00424  *  Returns the logarithm of +x+.
00425  *  If additional second argument is given, it will be the base
00426  *  of logarithm. Otherwise it is +e+ (for the natural logarithm).
00427  *
00428  *  Domain: (0, INFINITY)
00429  *
00430  *  Codomain: (-INFINITY, INFINITY)
00431  *
00432  *    Math.log(0)          #=> -Infinity
00433  *    Math.log(1)          #=> 0.0
00434  *    Math.log(Math::E)    #=> 1.0
00435  *    Math.log(Math::E**3) #=> 3.0
00436  *    Math.log(12, 3)      #=> 2.2618595071429146
00437  *
00438  */
00439 
00440 static VALUE
00441 math_log(int argc, VALUE *argv)
00442 {
00443     VALUE x, base;
00444     double d0, d;
00445     size_t numbits;
00446 
00447     rb_scan_args(argc, argv, "11", &x, &base);
00448 
00449     if (RB_BIGNUM_TYPE_P(x) && RBIGNUM_POSITIVE_P(x) &&
00450             DBL_MAX_EXP <= (numbits = rb_absint_numwords(x, 1, NULL))) {
00451         numbits -= DBL_MANT_DIG;
00452         x = rb_big_rshift(x, SIZET2NUM(numbits));
00453     }
00454     else {
00455         numbits = 0;
00456     }
00457 
00458     Need_Float(x);
00459     d0 = RFLOAT_VALUE(x);
00460     /* check for domain error */
00461     if (d0 < 0.0) domain_error("log");
00462     /* check for pole error */
00463     if (d0 == 0.0) return DBL2NUM(-INFINITY);
00464     d = log(d0);
00465     if (numbits)
00466         d += numbits * log(2); /* log(2**numbits) */
00467     if (argc == 2) {
00468         Need_Float(base);
00469         d /= log(RFLOAT_VALUE(base));
00470     }
00471     return DBL2NUM(d);
00472 }
00473 
00474 #ifndef log2
00475 #ifndef HAVE_LOG2
00476 double
00477 log2(double x)
00478 {
00479     return log10(x)/log10(2.0);
00480 }
00481 #else
00482 extern double log2(double);
00483 #endif
00484 #endif
00485 
00486 /*
00487  *  call-seq:
00488  *     Math.log2(x)    -> Float
00489  *
00490  *  Returns the base 2 logarithm of +x+.
00491  *
00492  *  Domain: (0, INFINITY)
00493  *
00494  *  Codomain: (-INFINITY, INFINITY)
00495  *
00496  *    Math.log2(1)      #=> 0.0
00497  *    Math.log2(2)      #=> 1.0
00498  *    Math.log2(32768)  #=> 15.0
00499  *    Math.log2(65536)  #=> 16.0
00500  *
00501  */
00502 
00503 static VALUE
00504 math_log2(VALUE obj, VALUE x)
00505 {
00506     double d0, d;
00507     size_t numbits;
00508 
00509     if (RB_BIGNUM_TYPE_P(x) && RBIGNUM_POSITIVE_P(x) &&
00510             DBL_MAX_EXP <= (numbits = rb_absint_numwords(x, 1, NULL))) {
00511         numbits -= DBL_MANT_DIG;
00512         x = rb_big_rshift(x, SIZET2NUM(numbits));
00513     }
00514     else {
00515         numbits = 0;
00516     }
00517 
00518     Need_Float(x);
00519     d0 = RFLOAT_VALUE(x);
00520     /* check for domain error */
00521     if (d0 < 0.0) domain_error("log2");
00522     /* check for pole error */
00523     if (d0 == 0.0) return DBL2NUM(-INFINITY);
00524     d = log2(d0);
00525     d += numbits;
00526     return DBL2NUM(d);
00527 }
00528 
00529 /*
00530  *  call-seq:
00531  *     Math.log10(x)    -> Float
00532  *
00533  *  Returns the base 10 logarithm of +x+.
00534  *
00535  *  Domain: (0, INFINITY)
00536  *
00537  *  Codomain: (-INFINITY, INFINITY)
00538  *
00539  *    Math.log10(1)       #=> 0.0
00540  *    Math.log10(10)      #=> 1.0
00541  *    Math.log10(10**100) #=> 100.0
00542  *
00543  */
00544 
00545 static VALUE
00546 math_log10(VALUE obj, VALUE x)
00547 {
00548     double d0, d;
00549     size_t numbits;
00550 
00551     if (RB_BIGNUM_TYPE_P(x) && RBIGNUM_POSITIVE_P(x) &&
00552             DBL_MAX_EXP <= (numbits = rb_absint_numwords(x, 1, NULL))) {
00553         numbits -= DBL_MANT_DIG;
00554         x = rb_big_rshift(x, SIZET2NUM(numbits));
00555     }
00556     else {
00557         numbits = 0;
00558     }
00559 
00560     Need_Float(x);
00561     d0 = RFLOAT_VALUE(x);
00562     /* check for domain error */
00563     if (d0 < 0.0) domain_error("log10");
00564     /* check for pole error */
00565     if (d0 == 0.0) return DBL2NUM(-INFINITY);
00566     d = log10(d0);
00567     if (numbits)
00568         d += numbits * log10(2); /* log10(2**numbits) */
00569     return DBL2NUM(d);
00570 }
00571 
00572 /*
00573  *  call-seq:
00574  *     Math.sqrt(x)    -> Float
00575  *
00576  *  Returns the non-negative square root of +x+.
00577  *
00578  *  Domain: [0, INFINITY)
00579  *
00580  *  Codomain:[0, INFINITY)
00581  *
00582  *    0.upto(10) {|x|
00583  *      p [x, Math.sqrt(x), Math.sqrt(x)**2]
00584  *    }
00585  *    #=> [0, 0.0, 0.0]
00586  *    #   [1, 1.0, 1.0]
00587  *    #   [2, 1.4142135623731, 2.0]
00588  *    #   [3, 1.73205080756888, 3.0]
00589  *    #   [4, 2.0, 4.0]
00590  *    #   [5, 2.23606797749979, 5.0]
00591  *    #   [6, 2.44948974278318, 6.0]
00592  *    #   [7, 2.64575131106459, 7.0]
00593  *    #   [8, 2.82842712474619, 8.0]
00594  *    #   [9, 3.0, 9.0]
00595  *    #   [10, 3.16227766016838, 10.0]
00596  */
00597 
00598 static VALUE
00599 math_sqrt(VALUE obj, VALUE x)
00600 {
00601     double d0, d;
00602 
00603     Need_Float(x);
00604     d0 = RFLOAT_VALUE(x);
00605     /* check for domain error */
00606     if (d0 < 0.0) domain_error("sqrt");
00607     if (d0 == 0.0) return DBL2NUM(0.0);
00608     d = sqrt(d0);
00609     return DBL2NUM(d);
00610 }
00611 
00612 /*
00613  *  call-seq:
00614  *     Math.cbrt(x)    -> Float
00615  *
00616  *  Returns the cube root of +x+.
00617  *
00618  *  Domain: [0, INFINITY)
00619  *
00620  *  Codomain:[0, INFINITY)
00621  *
00622  *    -9.upto(9) {|x|
00623  *      p [x, Math.cbrt(x), Math.cbrt(x)**3]
00624  *    }
00625  *    #=> [-9, -2.0800838230519, -9.0]
00626  *    #   [-8, -2.0, -8.0]
00627  *    #   [-7, -1.91293118277239, -7.0]
00628  *    #   [-6, -1.81712059283214, -6.0]
00629  *    #   [-5, -1.7099759466767, -5.0]
00630  *    #   [-4, -1.5874010519682, -4.0]
00631  *    #   [-3, -1.44224957030741, -3.0]
00632  *    #   [-2, -1.25992104989487, -2.0]
00633  *    #   [-1, -1.0, -1.0]
00634  *    #   [0, 0.0, 0.0]
00635  *    #   [1, 1.0, 1.0]
00636  *    #   [2, 1.25992104989487, 2.0]
00637  *    #   [3, 1.44224957030741, 3.0]
00638  *    #   [4, 1.5874010519682, 4.0]
00639  *    #   [5, 1.7099759466767, 5.0]
00640  *    #   [6, 1.81712059283214, 6.0]
00641  *    #   [7, 1.91293118277239, 7.0]
00642  *    #   [8, 2.0, 8.0]
00643  *    #   [9, 2.0800838230519, 9.0]
00644  *
00645  */
00646 
00647 static VALUE
00648 math_cbrt(VALUE obj, VALUE x)
00649 {
00650     Need_Float(x);
00651     return DBL2NUM(cbrt(RFLOAT_VALUE(x)));
00652 }
00653 
00654 /*
00655  *  call-seq:
00656  *     Math.frexp(x)    -> [fraction, exponent]
00657  *
00658  *  Returns a two-element array containing the normalized fraction (a Float)
00659  *  and exponent (a Fixnum) of +x+.
00660  *
00661  *     fraction, exponent = Math.frexp(1234)   #=> [0.6025390625, 11]
00662  *     fraction * 2**exponent                  #=> 1234.0
00663  */
00664 
00665 static VALUE
00666 math_frexp(VALUE obj, VALUE x)
00667 {
00668     double d;
00669     int exp;
00670 
00671     Need_Float(x);
00672 
00673     d = frexp(RFLOAT_VALUE(x), &exp);
00674     return rb_assoc_new(DBL2NUM(d), INT2NUM(exp));
00675 }
00676 
00677 /*
00678  *  call-seq:
00679  *     Math.ldexp(fraction, exponent) -> float
00680  *
00681  *  Returns the value of +fraction+*(2**+exponent+).
00682  *
00683  *     fraction, exponent = Math.frexp(1234)
00684  *     Math.ldexp(fraction, exponent)   #=> 1234.0
00685  */
00686 
00687 static VALUE
00688 math_ldexp(VALUE obj, VALUE x, VALUE n)
00689 {
00690     Need_Float(x);
00691     return DBL2NUM(ldexp(RFLOAT_VALUE(x), NUM2INT(n)));
00692 }
00693 
00694 /*
00695  *  call-seq:
00696  *     Math.hypot(x, y)    -> Float
00697  *
00698  *  Returns sqrt(x**2 + y**2), the hypotenuse of a right-angled triangle with
00699  *  sides +x+ and +y+.
00700  *
00701  *     Math.hypot(3, 4)   #=> 5.0
00702  */
00703 
00704 static VALUE
00705 math_hypot(VALUE obj, VALUE x, VALUE y)
00706 {
00707     Need_Float2(x, y);
00708     return DBL2NUM(hypot(RFLOAT_VALUE(x), RFLOAT_VALUE(y)));
00709 }
00710 
00711 /*
00712  * call-seq:
00713  *    Math.erf(x)  -> Float
00714  *
00715  *  Calculates the error function of +x+.
00716  *
00717  *  Domain: (-INFINITY, INFINITY)
00718  *
00719  *  Codomain: (-1, 1)
00720  *
00721  *    Math.erf(0) #=> 0.0
00722  *
00723  */
00724 
00725 static VALUE
00726 math_erf(VALUE obj, VALUE x)
00727 {
00728     Need_Float(x);
00729     return DBL2NUM(erf(RFLOAT_VALUE(x)));
00730 }
00731 
00732 /*
00733  * call-seq:
00734  *    Math.erfc(x)  -> Float
00735  *
00736  *  Calculates the complementary error function of x.
00737  *
00738  *  Domain: (-INFINITY, INFINITY)
00739  *
00740  *  Codomain: (0, 2)
00741  *
00742  *    Math.erfc(0) #=> 1.0
00743  *
00744  */
00745 
00746 static VALUE
00747 math_erfc(VALUE obj, VALUE x)
00748 {
00749     Need_Float(x);
00750     return DBL2NUM(erfc(RFLOAT_VALUE(x)));
00751 }
00752 
00753 /*
00754  * call-seq:
00755  *    Math.gamma(x)  -> Float
00756  *
00757  *  Calculates the gamma function of x.
00758  *
00759  *  Note that gamma(n) is same as fact(n-1) for integer n > 0.
00760  *  However gamma(n) returns float and can be an approximation.
00761  *
00762  *   def fact(n) (1..n).inject(1) {|r,i| r*i } end
00763  *   1.upto(26) {|i| p [i, Math.gamma(i), fact(i-1)] }
00764  *   #=> [1, 1.0, 1]
00765  *   #   [2, 1.0, 1]
00766  *   #   [3, 2.0, 2]
00767  *   #   [4, 6.0, 6]
00768  *   #   [5, 24.0, 24]
00769  *   #   [6, 120.0, 120]
00770  *   #   [7, 720.0, 720]
00771  *   #   [8, 5040.0, 5040]
00772  *   #   [9, 40320.0, 40320]
00773  *   #   [10, 362880.0, 362880]
00774  *   #   [11, 3628800.0, 3628800]
00775  *   #   [12, 39916800.0, 39916800]
00776  *   #   [13, 479001600.0, 479001600]
00777  *   #   [14, 6227020800.0, 6227020800]
00778  *   #   [15, 87178291200.0, 87178291200]
00779  *   #   [16, 1307674368000.0, 1307674368000]
00780  *   #   [17, 20922789888000.0, 20922789888000]
00781  *   #   [18, 355687428096000.0, 355687428096000]
00782  *   #   [19, 6.402373705728e+15, 6402373705728000]
00783  *   #   [20, 1.21645100408832e+17, 121645100408832000]
00784  *   #   [21, 2.43290200817664e+18, 2432902008176640000]
00785  *   #   [22, 5.109094217170944e+19, 51090942171709440000]
00786  *   #   [23, 1.1240007277776077e+21, 1124000727777607680000]
00787  *   #   [24, 2.5852016738885062e+22, 25852016738884976640000]
00788  *   #   [25, 6.204484017332391e+23, 620448401733239439360000]
00789  *   #   [26, 1.5511210043330954e+25, 15511210043330985984000000]
00790  *
00791  */
00792 
00793 static VALUE
00794 math_gamma(VALUE obj, VALUE x)
00795 {
00796     static const double fact_table[] = {
00797         /* fact(0) */ 1.0,
00798         /* fact(1) */ 1.0,
00799         /* fact(2) */ 2.0,
00800         /* fact(3) */ 6.0,
00801         /* fact(4) */ 24.0,
00802         /* fact(5) */ 120.0,
00803         /* fact(6) */ 720.0,
00804         /* fact(7) */ 5040.0,
00805         /* fact(8) */ 40320.0,
00806         /* fact(9) */ 362880.0,
00807         /* fact(10) */ 3628800.0,
00808         /* fact(11) */ 39916800.0,
00809         /* fact(12) */ 479001600.0,
00810         /* fact(13) */ 6227020800.0,
00811         /* fact(14) */ 87178291200.0,
00812         /* fact(15) */ 1307674368000.0,
00813         /* fact(16) */ 20922789888000.0,
00814         /* fact(17) */ 355687428096000.0,
00815         /* fact(18) */ 6402373705728000.0,
00816         /* fact(19) */ 121645100408832000.0,
00817         /* fact(20) */ 2432902008176640000.0,
00818         /* fact(21) */ 51090942171709440000.0,
00819         /* fact(22) */ 1124000727777607680000.0,
00820         /* fact(23)=25852016738884976640000 needs 56bit mantissa which is
00821          * impossible to represent exactly in IEEE 754 double which have
00822          * 53bit mantissa. */
00823     };
00824     double d0, d;
00825     double intpart, fracpart;
00826     Need_Float(x);
00827     d0 = RFLOAT_VALUE(x);
00828     /* check for domain error */
00829     if (isinf(d0) && signbit(d0)) domain_error("gamma");
00830     fracpart = modf(d0, &intpart);
00831     if (fracpart == 0.0) {
00832         if (intpart < 0) domain_error("gamma");
00833         if (0 < intpart &&
00834             intpart - 1 < (double)numberof(fact_table)) {
00835             return DBL2NUM(fact_table[(int)intpart - 1]);
00836         }
00837     }
00838     d = tgamma(d0);
00839     return DBL2NUM(d);
00840 }
00841 
00842 /*
00843  * call-seq:
00844  *    Math.lgamma(x)  -> [float, -1 or 1]
00845  *
00846  *  Calculates the logarithmic gamma of +x+ and the sign of gamma of +x+.
00847  *
00848  *  Math.lgamma(x) is same as
00849  *   [Math.log(Math.gamma(x).abs), Math.gamma(x) < 0 ? -1 : 1]
00850  *  but avoid overflow by Math.gamma(x) for large x.
00851  *
00852  *    Math.lgamma(0) #=> [Infinity, 1]
00853  *
00854  */
00855 
00856 static VALUE
00857 math_lgamma(VALUE obj, VALUE x)
00858 {
00859     double d0, d;
00860     int sign=1;
00861     VALUE v;
00862     Need_Float(x);
00863     d0 = RFLOAT_VALUE(x);
00864     /* check for domain error */
00865     if (isinf(d0)) {
00866         if (signbit(d0)) domain_error("lgamma");
00867         return rb_assoc_new(DBL2NUM(INFINITY), INT2FIX(1));
00868     }
00869     d = lgamma_r(d0, &sign);
00870     v = DBL2NUM(d);
00871     return rb_assoc_new(v, INT2FIX(sign));
00872 }
00873 
00874 
00875 #define exp1(n) \
00876 VALUE \
00877 rb_math_##n(VALUE x)\
00878 {\
00879     return math_##n(rb_mMath, x);\
00880 }
00881 
00882 #define exp2(n) \
00883 VALUE \
00884 rb_math_##n(VALUE x, VALUE y)\
00885 {\
00886     return math_##n(rb_mMath, x, y);\
00887 }
00888 
00889 exp2(atan2)
00890 exp1(cos)
00891 exp1(cosh)
00892 exp1(exp)
00893 exp2(hypot)
00894 
00895 VALUE
00896 rb_math_log(int argc, VALUE *argv)
00897 {
00898     return math_log(argc, argv);
00899 }
00900 
00901 exp1(sin)
00902 exp1(sinh)
00903 exp1(sqrt)
00904 
00905 
00906 /*
00907  *  Document-class: Math::DomainError
00908  *
00909  *  Raised when a mathematical function is evaluated outside of its
00910  *  domain of definition.
00911  *
00912  *  For example, since +cos+ returns values in the range -1..1,
00913  *  its inverse function +acos+ is only defined on that interval:
00914  *
00915  *     Math.acos(42)
00916  *
00917  *  <em>produces:</em>
00918  *
00919  *     Math::DomainError: Numerical argument is out of domain - "acos"
00920  */
00921 
00922 /*
00923  *  Document-class: Math
00924  *
00925  *  The Math module contains module functions for basic
00926  *  trigonometric and transcendental functions. See class
00927  *  Float for a list of constants that
00928  *  define Ruby's floating point accuracy.
00929  *
00930  *  Domains and codomains are given only for real (not complex) numbers.
00931  */
00932 
00933 
00934 void
00935 Init_Math(void)
00936 {
00937     rb_mMath = rb_define_module("Math");
00938     rb_eMathDomainError = rb_define_class_under(rb_mMath, "DomainError", rb_eStandardError);
00939 
00940 #ifdef M_PI
00941     /*  Definition of the mathematical constant PI as a Float number. */
00942     rb_define_const(rb_mMath, "PI", DBL2NUM(M_PI));
00943 #else
00944     rb_define_const(rb_mMath, "PI", DBL2NUM(atan(1.0)*4.0));
00945 #endif
00946 
00947 #ifdef M_E
00948     /*  Definition of the mathematical constant E (e) as a Float number. */
00949     rb_define_const(rb_mMath, "E", DBL2NUM(M_E));
00950 #else
00951     rb_define_const(rb_mMath, "E", DBL2NUM(exp(1.0)));
00952 #endif
00953 
00954     rb_define_module_function(rb_mMath, "atan2", math_atan2, 2);
00955     rb_define_module_function(rb_mMath, "cos", math_cos, 1);
00956     rb_define_module_function(rb_mMath, "sin", math_sin, 1);
00957     rb_define_module_function(rb_mMath, "tan", math_tan, 1);
00958 
00959     rb_define_module_function(rb_mMath, "acos", math_acos, 1);
00960     rb_define_module_function(rb_mMath, "asin", math_asin, 1);
00961     rb_define_module_function(rb_mMath, "atan", math_atan, 1);
00962 
00963     rb_define_module_function(rb_mMath, "cosh", math_cosh, 1);
00964     rb_define_module_function(rb_mMath, "sinh", math_sinh, 1);
00965     rb_define_module_function(rb_mMath, "tanh", math_tanh, 1);
00966 
00967     rb_define_module_function(rb_mMath, "acosh", math_acosh, 1);
00968     rb_define_module_function(rb_mMath, "asinh", math_asinh, 1);
00969     rb_define_module_function(rb_mMath, "atanh", math_atanh, 1);
00970 
00971     rb_define_module_function(rb_mMath, "exp", math_exp, 1);
00972     rb_define_module_function(rb_mMath, "log", math_log, -1);
00973     rb_define_module_function(rb_mMath, "log2", math_log2, 1);
00974     rb_define_module_function(rb_mMath, "log10", math_log10, 1);
00975     rb_define_module_function(rb_mMath, "sqrt", math_sqrt, 1);
00976     rb_define_module_function(rb_mMath, "cbrt", math_cbrt, 1);
00977 
00978     rb_define_module_function(rb_mMath, "frexp", math_frexp, 1);
00979     rb_define_module_function(rb_mMath, "ldexp", math_ldexp, 2);
00980 
00981     rb_define_module_function(rb_mMath, "hypot", math_hypot, 2);
00982 
00983     rb_define_module_function(rb_mMath, "erf",  math_erf,  1);
00984     rb_define_module_function(rb_mMath, "erfc", math_erfc, 1);
00985 
00986     rb_define_module_function(rb_mMath, "gamma", math_gamma, 1);
00987     rb_define_module_function(rb_mMath, "lgamma", math_lgamma, 1);
00988 }
00989