00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #include "ruby/ruby.h"
00013 #include <math.h>
00014 #include <errno.h>
00015
00016 #define numberof(array) (int)(sizeof(array) / sizeof((array)[0]))
00017
00018 VALUE rb_mMath;
00019 VALUE rb_eMathDomainError;
00020
00021 extern VALUE rb_to_float(VALUE val);
00022 #define Need_Float(x) do {if (TYPE(x) != T_FLOAT) {(x) = rb_to_float(x);}} while(0)
00023 #define Need_Float2(x,y) do {\
00024 Need_Float(x);\
00025 Need_Float(y);\
00026 } while (0)
00027
00028 #define domain_error(msg) \
00029 rb_raise(rb_eMathDomainError, "Numerical argument is out of domain - " #msg);
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051 static VALUE
00052 math_atan2(VALUE obj, VALUE y, VALUE x)
00053 {
00054 double dx, dy;
00055 Need_Float2(y, x);
00056 dx = RFLOAT_VALUE(x);
00057 dy = RFLOAT_VALUE(y);
00058 if (dx == 0.0 && dy == 0.0) domain_error("atan2");
00059 if (isinf(dx) && isinf(dy)) domain_error("atan2");
00060 return DBL2NUM(atan2(dy, dx));
00061 }
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072 static VALUE
00073 math_cos(VALUE obj, VALUE x)
00074 {
00075 Need_Float(x);
00076 return DBL2NUM(cos(RFLOAT_VALUE(x)));
00077 }
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087 static VALUE
00088 math_sin(VALUE obj, VALUE x)
00089 {
00090 Need_Float(x);
00091
00092 return DBL2NUM(sin(RFLOAT_VALUE(x)));
00093 }
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103 static VALUE
00104 math_tan(VALUE obj, VALUE x)
00105 {
00106 Need_Float(x);
00107
00108 return DBL2NUM(tan(RFLOAT_VALUE(x)));
00109 }
00110
00111
00112
00113
00114
00115
00116
00117
00118 static VALUE
00119 math_acos(VALUE obj, VALUE x)
00120 {
00121 double d0, d;
00122
00123 Need_Float(x);
00124 d0 = RFLOAT_VALUE(x);
00125
00126 if (d0 < -1.0 || 1.0 < d0) domain_error("acos");
00127 d = acos(d0);
00128 return DBL2NUM(d);
00129 }
00130
00131
00132
00133
00134
00135
00136
00137
00138 static VALUE
00139 math_asin(VALUE obj, VALUE x)
00140 {
00141 double d0, d;
00142
00143 Need_Float(x);
00144 d0 = RFLOAT_VALUE(x);
00145
00146 if (d0 < -1.0 || 1.0 < d0) domain_error("asin");
00147 d = asin(d0);
00148 return DBL2NUM(d);
00149 }
00150
00151
00152
00153
00154
00155
00156
00157
00158 static VALUE
00159 math_atan(VALUE obj, VALUE x)
00160 {
00161 Need_Float(x);
00162 return DBL2NUM(atan(RFLOAT_VALUE(x)));
00163 }
00164
00165 #ifndef HAVE_COSH
00166 double
00167 cosh(double x)
00168 {
00169 return (exp(x) + exp(-x)) / 2;
00170 }
00171 #endif
00172
00173
00174
00175
00176
00177
00178
00179
00180 static VALUE
00181 math_cosh(VALUE obj, VALUE x)
00182 {
00183 Need_Float(x);
00184
00185 return DBL2NUM(cosh(RFLOAT_VALUE(x)));
00186 }
00187
00188 #ifndef HAVE_SINH
00189 double
00190 sinh(double x)
00191 {
00192 return (exp(x) - exp(-x)) / 2;
00193 }
00194 #endif
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204 static VALUE
00205 math_sinh(VALUE obj, VALUE x)
00206 {
00207 Need_Float(x);
00208 return DBL2NUM(sinh(RFLOAT_VALUE(x)));
00209 }
00210
00211 #ifndef HAVE_TANH
00212 double
00213 tanh(double x)
00214 {
00215 return sinh(x) / cosh(x);
00216 }
00217 #endif
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227 static VALUE
00228 math_tanh(VALUE obj, VALUE x)
00229 {
00230 Need_Float(x);
00231 return DBL2NUM(tanh(RFLOAT_VALUE(x)));
00232 }
00233
00234
00235
00236
00237
00238
00239
00240
00241 static VALUE
00242 math_acosh(VALUE obj, VALUE x)
00243 {
00244 double d0, d;
00245
00246 Need_Float(x);
00247 d0 = RFLOAT_VALUE(x);
00248
00249 if (d0 < 1.0) domain_error("acosh");
00250 d = acosh(d0);
00251 return DBL2NUM(d);
00252 }
00253
00254
00255
00256
00257
00258
00259
00260
00261 static VALUE
00262 math_asinh(VALUE obj, VALUE x)
00263 {
00264 Need_Float(x);
00265 return DBL2NUM(asinh(RFLOAT_VALUE(x)));
00266 }
00267
00268
00269
00270
00271
00272
00273
00274
00275 static VALUE
00276 math_atanh(VALUE obj, VALUE x)
00277 {
00278 double d0, d;
00279
00280 Need_Float(x);
00281 d0 = RFLOAT_VALUE(x);
00282
00283 if (d0 < -1.0 || +1.0 < d0) domain_error("atanh");
00284
00285 if (d0 == -1.0) return DBL2NUM(-INFINITY);
00286 if (d0 == +1.0) return DBL2NUM(+INFINITY);
00287 d = atanh(d0);
00288 return DBL2NUM(d);
00289 }
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303 static VALUE
00304 math_exp(VALUE obj, VALUE x)
00305 {
00306 Need_Float(x);
00307 return DBL2NUM(exp(RFLOAT_VALUE(x)));
00308 }
00309
00310 #if defined __CYGWIN__
00311 # include <cygwin/version.h>
00312 # if CYGWIN_VERSION_DLL_MAJOR < 1005
00313 # define nan(x) nan()
00314 # endif
00315 # define log(x) ((x) < 0.0 ? nan("") : log(x))
00316 # define log10(x) ((x) < 0.0 ? nan("") : log10(x))
00317 #endif
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335 static VALUE
00336 math_log(int argc, VALUE *argv)
00337 {
00338 VALUE x, base;
00339 double d0, d;
00340
00341 rb_scan_args(argc, argv, "11", &x, &base);
00342 Need_Float(x);
00343 d0 = RFLOAT_VALUE(x);
00344
00345 if (d0 < 0.0) domain_error("log");
00346
00347 if (d0 == 0.0) return DBL2NUM(-INFINITY);
00348 d = log(d0);
00349 if (argc == 2) {
00350 Need_Float(base);
00351 d /= log(RFLOAT_VALUE(base));
00352 }
00353 return DBL2NUM(d);
00354 }
00355
00356 #ifndef log2
00357 #ifndef HAVE_LOG2
00358 double
00359 log2(double x)
00360 {
00361 return log10(x)/log10(2.0);
00362 }
00363 #else
00364 extern double log2(double);
00365 #endif
00366 #endif
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381 static VALUE
00382 math_log2(VALUE obj, VALUE x)
00383 {
00384 double d0, d;
00385
00386 Need_Float(x);
00387 d0 = RFLOAT_VALUE(x);
00388
00389 if (d0 < 0.0) domain_error("log2");
00390
00391 if (d0 == 0.0) return DBL2NUM(-INFINITY);
00392 d = log2(d0);
00393 return DBL2NUM(d);
00394 }
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408 static VALUE
00409 math_log10(VALUE obj, VALUE x)
00410 {
00411 double d0, d;
00412
00413 Need_Float(x);
00414 d0 = RFLOAT_VALUE(x);
00415
00416 if (d0 < 0.0) domain_error("log10");
00417
00418 if (d0 == 0.0) return DBL2NUM(-INFINITY);
00419 d = log10(d0);
00420 return DBL2NUM(d);
00421 }
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444
00445
00446
00447 static VALUE
00448 math_sqrt(VALUE obj, VALUE x)
00449 {
00450 double d0, d;
00451
00452 Need_Float(x);
00453 d0 = RFLOAT_VALUE(x);
00454
00455 if (d0 < 0.0) domain_error("sqrt");
00456 if (d0 == 0.0) return DBL2NUM(0.0);
00457 d = sqrt(d0);
00458 return DBL2NUM(d);
00459 }
00460
00461
00462
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493 static VALUE
00494 math_cbrt(VALUE obj, VALUE x)
00495 {
00496 Need_Float(x);
00497 return DBL2NUM(cbrt(RFLOAT_VALUE(x)));
00498 }
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511
00512 static VALUE
00513 math_frexp(VALUE obj, VALUE x)
00514 {
00515 double d;
00516 int exp;
00517
00518 Need_Float(x);
00519
00520 d = frexp(RFLOAT_VALUE(x), &exp);
00521 return rb_assoc_new(DBL2NUM(d), INT2NUM(exp));
00522 }
00523
00524
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534 static VALUE
00535 math_ldexp(VALUE obj, VALUE x, VALUE n)
00536 {
00537 Need_Float(x);
00538 return DBL2NUM(ldexp(RFLOAT_VALUE(x), NUM2INT(n)));
00539 }
00540
00541
00542
00543
00544
00545
00546
00547
00548
00549
00550
00551 static VALUE
00552 math_hypot(VALUE obj, VALUE x, VALUE y)
00553 {
00554 Need_Float2(x, y);
00555 return DBL2NUM(hypot(RFLOAT_VALUE(x), RFLOAT_VALUE(y)));
00556 }
00557
00558
00559
00560
00561
00562
00563
00564
00565 static VALUE
00566 math_erf(VALUE obj, VALUE x)
00567 {
00568 Need_Float(x);
00569 return DBL2NUM(erf(RFLOAT_VALUE(x)));
00570 }
00571
00572
00573
00574
00575
00576
00577
00578
00579 static VALUE
00580 math_erfc(VALUE obj, VALUE x)
00581 {
00582 Need_Float(x);
00583 return DBL2NUM(erfc(RFLOAT_VALUE(x)));
00584 }
00585
00586
00587
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616
00617
00618
00619
00620
00621
00622
00623
00624
00625
00626 static VALUE
00627 math_gamma(VALUE obj, VALUE x)
00628 {
00629 static const double fact_table[] = {
00630 1.0,
00631 1.0,
00632 2.0,
00633 6.0,
00634 24.0,
00635 120.0,
00636 720.0,
00637 5040.0,
00638 40320.0,
00639 362880.0,
00640 3628800.0,
00641 39916800.0,
00642 479001600.0,
00643 6227020800.0,
00644 87178291200.0,
00645 1307674368000.0,
00646 20922789888000.0,
00647 355687428096000.0,
00648 6402373705728000.0,
00649 121645100408832000.0,
00650 2432902008176640000.0,
00651 51090942171709440000.0,
00652 1124000727777607680000.0,
00653
00654
00655
00656 };
00657 double d0, d;
00658 double intpart, fracpart;
00659 Need_Float(x);
00660 d0 = RFLOAT_VALUE(x);
00661
00662 if (isinf(d0) && signbit(d0)) domain_error("gamma");
00663 fracpart = modf(d0, &intpart);
00664 if (fracpart == 0.0) {
00665 if (intpart < 0) domain_error("gamma");
00666 if (0 < intpart &&
00667 intpart - 1 < (double)numberof(fact_table)) {
00668 return DBL2NUM(fact_table[(int)intpart - 1]);
00669 }
00670 }
00671 d = tgamma(d0);
00672 return DBL2NUM(d);
00673 }
00674
00675
00676
00677
00678
00679
00680
00681
00682
00683
00684
00685
00686
00687 static VALUE
00688 math_lgamma(VALUE obj, VALUE x)
00689 {
00690 double d0, d;
00691 int sign=1;
00692 VALUE v;
00693 Need_Float(x);
00694 d0 = RFLOAT_VALUE(x);
00695
00696 if (isinf(d0)) {
00697 if (signbit(d0)) domain_error("lgamma");
00698 return rb_assoc_new(DBL2NUM(INFINITY), INT2FIX(1));
00699 }
00700 d = lgamma_r(d0, &sign);
00701 v = DBL2NUM(d);
00702 return rb_assoc_new(v, INT2FIX(sign));
00703 }
00704
00705
00706 #define exp1(n) \
00707 VALUE \
00708 rb_math_##n(VALUE x)\
00709 {\
00710 return math_##n(rb_mMath, x);\
00711 }
00712
00713 #define exp2(n) \
00714 VALUE \
00715 rb_math_##n(VALUE x, VALUE y)\
00716 {\
00717 return math_##n(rb_mMath, x, y);\
00718 }
00719
00720 exp2(atan2)
00721 exp1(cos)
00722 exp1(cosh)
00723 exp1(exp)
00724 exp2(hypot)
00725
00726 VALUE
00727 rb_math_log(int argc, VALUE *argv)
00728 {
00729 return math_log(argc, argv);
00730 }
00731
00732 exp1(sin)
00733 exp1(sinh)
00734 exp1(sqrt)
00735
00736
00737
00738
00739
00740
00741
00742
00743
00744
00745
00746
00747
00748
00749
00750
00751
00752
00753
00754
00755
00756
00757
00758
00759
00760
00761 void
00762 Init_Math(void)
00763 {
00764 rb_mMath = rb_define_module("Math");
00765 rb_eMathDomainError = rb_define_class_under(rb_mMath, "DomainError", rb_eStandardError);
00766
00767 #ifdef M_PI
00768 rb_define_const(rb_mMath, "PI", DBL2NUM(M_PI));
00769 #else
00770 rb_define_const(rb_mMath, "PI", DBL2NUM(atan(1.0)*4.0));
00771 #endif
00772
00773 #ifdef M_E
00774 rb_define_const(rb_mMath, "E", DBL2NUM(M_E));
00775 #else
00776 rb_define_const(rb_mMath, "E", DBL2NUM(exp(1.0)));
00777 #endif
00778
00779 rb_define_module_function(rb_mMath, "atan2", math_atan2, 2);
00780 rb_define_module_function(rb_mMath, "cos", math_cos, 1);
00781 rb_define_module_function(rb_mMath, "sin", math_sin, 1);
00782 rb_define_module_function(rb_mMath, "tan", math_tan, 1);
00783
00784 rb_define_module_function(rb_mMath, "acos", math_acos, 1);
00785 rb_define_module_function(rb_mMath, "asin", math_asin, 1);
00786 rb_define_module_function(rb_mMath, "atan", math_atan, 1);
00787
00788 rb_define_module_function(rb_mMath, "cosh", math_cosh, 1);
00789 rb_define_module_function(rb_mMath, "sinh", math_sinh, 1);
00790 rb_define_module_function(rb_mMath, "tanh", math_tanh, 1);
00791
00792 rb_define_module_function(rb_mMath, "acosh", math_acosh, 1);
00793 rb_define_module_function(rb_mMath, "asinh", math_asinh, 1);
00794 rb_define_module_function(rb_mMath, "atanh", math_atanh, 1);
00795
00796 rb_define_module_function(rb_mMath, "exp", math_exp, 1);
00797 rb_define_module_function(rb_mMath, "log", math_log, -1);
00798 rb_define_module_function(rb_mMath, "log2", math_log2, 1);
00799 rb_define_module_function(rb_mMath, "log10", math_log10, 1);
00800 rb_define_module_function(rb_mMath, "sqrt", math_sqrt, 1);
00801 rb_define_module_function(rb_mMath, "cbrt", math_cbrt, 1);
00802
00803 rb_define_module_function(rb_mMath, "frexp", math_frexp, 1);
00804 rb_define_module_function(rb_mMath, "ldexp", math_ldexp, 2);
00805
00806 rb_define_module_function(rb_mMath, "hypot", math_hypot, 2);
00807
00808 rb_define_module_function(rb_mMath, "erf", math_erf, 1);
00809 rb_define_module_function(rb_mMath, "erfc", math_erfc, 1);
00810
00811 rb_define_module_function(rb_mMath, "gamma", math_gamma, 1);
00812 rb_define_module_function(rb_mMath, "lgamma", math_lgamma, 1);
00813 }
00814