Ruby 1.9.3p327(2012-11-10revision37606)
|
00001 /* 00002 complex.c: Coded by Tadayoshi Funaba 2008-2011 00003 00004 This implementation is based on Keiju Ishitsuka's Complex library 00005 which is written in ruby. 00006 */ 00007 00008 #include "ruby.h" 00009 #include "internal.h" 00010 #include <math.h> 00011 00012 #define NDEBUG 00013 #include <assert.h> 00014 00015 #define ZERO INT2FIX(0) 00016 #define ONE INT2FIX(1) 00017 #define TWO INT2FIX(2) 00018 00019 VALUE rb_cComplex; 00020 00021 static ID id_abs, id_abs2, id_arg, id_cmp, id_conj, id_convert, 00022 id_denominator, id_divmod, id_eqeq_p, id_expt, id_fdiv, id_floor, 00023 id_idiv, id_imag, id_inspect, id_negate, id_numerator, id_quo, 00024 id_real, id_real_p, id_to_f, id_to_i, id_to_r, id_to_s; 00025 00026 #define f_boolcast(x) ((x) ? Qtrue : Qfalse) 00027 00028 #define binop(n,op) \ 00029 inline static VALUE \ 00030 f_##n(VALUE x, VALUE y)\ 00031 {\ 00032 return rb_funcall(x, (op), 1, y);\ 00033 } 00034 00035 #define fun1(n) \ 00036 inline static VALUE \ 00037 f_##n(VALUE x)\ 00038 {\ 00039 return rb_funcall(x, id_##n, 0);\ 00040 } 00041 00042 #define fun2(n) \ 00043 inline static VALUE \ 00044 f_##n(VALUE x, VALUE y)\ 00045 {\ 00046 return rb_funcall(x, id_##n, 1, y);\ 00047 } 00048 00049 #define math1(n) \ 00050 inline static VALUE \ 00051 m_##n(VALUE x)\ 00052 {\ 00053 return rb_funcall(rb_mMath, id_##n, 1, x);\ 00054 } 00055 00056 #define math2(n) \ 00057 inline static VALUE \ 00058 m_##n(VALUE x, VALUE y)\ 00059 {\ 00060 return rb_funcall(rb_mMath, id_##n, 2, x, y);\ 00061 } 00062 00063 #define PRESERVE_SIGNEDZERO 00064 00065 inline static VALUE 00066 f_add(VALUE x, VALUE y) 00067 { 00068 #ifndef PRESERVE_SIGNEDZERO 00069 if (FIXNUM_P(y) && FIX2LONG(y) == 0) 00070 return x; 00071 else if (FIXNUM_P(x) && FIX2LONG(x) == 0) 00072 return y; 00073 #endif 00074 return rb_funcall(x, '+', 1, y); 00075 } 00076 00077 inline static VALUE 00078 f_cmp(VALUE x, VALUE y) 00079 { 00080 if (FIXNUM_P(x) && FIXNUM_P(y)) { 00081 long c = FIX2LONG(x) - FIX2LONG(y); 00082 if (c > 0) 00083 c = 1; 00084 else if (c < 0) 00085 c = -1; 00086 return INT2FIX(c); 00087 } 00088 return rb_funcall(x, id_cmp, 1, y); 00089 } 00090 00091 inline static VALUE 00092 f_div(VALUE x, VALUE y) 00093 { 00094 if (FIXNUM_P(y) && FIX2LONG(y) == 1) 00095 return x; 00096 return rb_funcall(x, '/', 1, y); 00097 } 00098 00099 inline static VALUE 00100 f_gt_p(VALUE x, VALUE y) 00101 { 00102 if (FIXNUM_P(x) && FIXNUM_P(y)) 00103 return f_boolcast(FIX2LONG(x) > FIX2LONG(y)); 00104 return rb_funcall(x, '>', 1, y); 00105 } 00106 00107 inline static VALUE 00108 f_lt_p(VALUE x, VALUE y) 00109 { 00110 if (FIXNUM_P(x) && FIXNUM_P(y)) 00111 return f_boolcast(FIX2LONG(x) < FIX2LONG(y)); 00112 return rb_funcall(x, '<', 1, y); 00113 } 00114 00115 binop(mod, '%') 00116 00117 inline static VALUE 00118 f_mul(VALUE x, VALUE y) 00119 { 00120 #ifndef PRESERVE_SIGNEDZERO 00121 if (FIXNUM_P(y)) { 00122 long iy = FIX2LONG(y); 00123 if (iy == 0) { 00124 if (FIXNUM_P(x) || TYPE(x) == T_BIGNUM) 00125 return ZERO; 00126 } 00127 else if (iy == 1) 00128 return x; 00129 } 00130 else if (FIXNUM_P(x)) { 00131 long ix = FIX2LONG(x); 00132 if (ix == 0) { 00133 if (FIXNUM_P(y) || TYPE(y) == T_BIGNUM) 00134 return ZERO; 00135 } 00136 else if (ix == 1) 00137 return y; 00138 } 00139 #endif 00140 return rb_funcall(x, '*', 1, y); 00141 } 00142 00143 inline static VALUE 00144 f_sub(VALUE x, VALUE y) 00145 { 00146 #ifndef PRESERVE_SIGNEDZERO 00147 if (FIXNUM_P(y) && FIX2LONG(y) == 0) 00148 return x; 00149 #endif 00150 return rb_funcall(x, '-', 1, y); 00151 } 00152 00153 fun1(abs) 00154 fun1(abs2) 00155 fun1(arg) 00156 fun1(conj) 00157 fun1(denominator) 00158 fun1(floor) 00159 fun1(imag) 00160 fun1(inspect) 00161 fun1(negate) 00162 fun1(numerator) 00163 fun1(real) 00164 fun1(real_p) 00165 00166 inline static VALUE 00167 f_to_i(VALUE x) 00168 { 00169 if (TYPE(x) == T_STRING) 00170 return rb_str_to_inum(x, 10, 0); 00171 return rb_funcall(x, id_to_i, 0); 00172 } 00173 inline static VALUE 00174 f_to_f(VALUE x) 00175 { 00176 if (TYPE(x) == T_STRING) 00177 return DBL2NUM(rb_str_to_dbl(x, 0)); 00178 return rb_funcall(x, id_to_f, 0); 00179 } 00180 00181 fun1(to_r) 00182 fun1(to_s) 00183 00184 fun2(divmod) 00185 00186 inline static VALUE 00187 f_eqeq_p(VALUE x, VALUE y) 00188 { 00189 if (FIXNUM_P(x) && FIXNUM_P(y)) 00190 return f_boolcast(FIX2LONG(x) == FIX2LONG(y)); 00191 return rb_funcall(x, id_eqeq_p, 1, y); 00192 } 00193 00194 fun2(expt) 00195 fun2(fdiv) 00196 fun2(idiv) 00197 fun2(quo) 00198 00199 inline static VALUE 00200 f_negative_p(VALUE x) 00201 { 00202 if (FIXNUM_P(x)) 00203 return f_boolcast(FIX2LONG(x) < 0); 00204 return rb_funcall(x, '<', 1, ZERO); 00205 } 00206 00207 #define f_positive_p(x) (!f_negative_p(x)) 00208 00209 inline static VALUE 00210 f_zero_p(VALUE x) 00211 { 00212 switch (TYPE(x)) { 00213 case T_FIXNUM: 00214 return f_boolcast(FIX2LONG(x) == 0); 00215 case T_BIGNUM: 00216 return Qfalse; 00217 case T_RATIONAL: 00218 { 00219 VALUE num = RRATIONAL(x)->num; 00220 00221 return f_boolcast(FIXNUM_P(num) && FIX2LONG(num) == 0); 00222 } 00223 } 00224 return rb_funcall(x, id_eqeq_p, 1, ZERO); 00225 } 00226 00227 #define f_nonzero_p(x) (!f_zero_p(x)) 00228 00229 inline static VALUE 00230 f_one_p(VALUE x) 00231 { 00232 switch (TYPE(x)) { 00233 case T_FIXNUM: 00234 return f_boolcast(FIX2LONG(x) == 1); 00235 case T_BIGNUM: 00236 return Qfalse; 00237 case T_RATIONAL: 00238 { 00239 VALUE num = RRATIONAL(x)->num; 00240 VALUE den = RRATIONAL(x)->den; 00241 00242 return f_boolcast(FIXNUM_P(num) && FIX2LONG(num) == 1 && 00243 FIXNUM_P(den) && FIX2LONG(den) == 1); 00244 } 00245 } 00246 return rb_funcall(x, id_eqeq_p, 1, ONE); 00247 } 00248 00249 inline static VALUE 00250 f_kind_of_p(VALUE x, VALUE c) 00251 { 00252 return rb_obj_is_kind_of(x, c); 00253 } 00254 00255 inline static VALUE 00256 k_numeric_p(VALUE x) 00257 { 00258 return f_kind_of_p(x, rb_cNumeric); 00259 } 00260 00261 inline static VALUE 00262 k_integer_p(VALUE x) 00263 { 00264 return f_kind_of_p(x, rb_cInteger); 00265 } 00266 00267 inline static VALUE 00268 k_fixnum_p(VALUE x) 00269 { 00270 return f_kind_of_p(x, rb_cFixnum); 00271 } 00272 00273 inline static VALUE 00274 k_bignum_p(VALUE x) 00275 { 00276 return f_kind_of_p(x, rb_cBignum); 00277 } 00278 00279 inline static VALUE 00280 k_float_p(VALUE x) 00281 { 00282 return f_kind_of_p(x, rb_cFloat); 00283 } 00284 00285 inline static VALUE 00286 k_rational_p(VALUE x) 00287 { 00288 return f_kind_of_p(x, rb_cRational); 00289 } 00290 00291 inline static VALUE 00292 k_complex_p(VALUE x) 00293 { 00294 return f_kind_of_p(x, rb_cComplex); 00295 } 00296 00297 #define k_exact_p(x) (!k_float_p(x)) 00298 #define k_inexact_p(x) k_float_p(x) 00299 00300 #define k_exact_zero_p(x) (k_exact_p(x) && f_zero_p(x)) 00301 #define k_exact_one_p(x) (k_exact_p(x) && f_one_p(x)) 00302 00303 #define get_dat1(x) \ 00304 struct RComplex *dat;\ 00305 dat = ((struct RComplex *)(x)) 00306 00307 #define get_dat2(x,y) \ 00308 struct RComplex *adat, *bdat;\ 00309 adat = ((struct RComplex *)(x));\ 00310 bdat = ((struct RComplex *)(y)) 00311 00312 inline static VALUE 00313 nucomp_s_new_internal(VALUE klass, VALUE real, VALUE imag) 00314 { 00315 NEWOBJ(obj, struct RComplex); 00316 OBJSETUP(obj, klass, T_COMPLEX); 00317 00318 obj->real = real; 00319 obj->imag = imag; 00320 00321 return (VALUE)obj; 00322 } 00323 00324 static VALUE 00325 nucomp_s_alloc(VALUE klass) 00326 { 00327 return nucomp_s_new_internal(klass, ZERO, ZERO); 00328 } 00329 00330 #if 0 00331 static VALUE 00332 nucomp_s_new_bang(int argc, VALUE *argv, VALUE klass) 00333 { 00334 VALUE real, imag; 00335 00336 switch (rb_scan_args(argc, argv, "11", &real, &imag)) { 00337 case 1: 00338 if (!k_numeric_p(real)) 00339 real = f_to_i(real); 00340 imag = ZERO; 00341 break; 00342 default: 00343 if (!k_numeric_p(real)) 00344 real = f_to_i(real); 00345 if (!k_numeric_p(imag)) 00346 imag = f_to_i(imag); 00347 break; 00348 } 00349 00350 return nucomp_s_new_internal(klass, real, imag); 00351 } 00352 #endif 00353 00354 inline static VALUE 00355 f_complex_new_bang1(VALUE klass, VALUE x) 00356 { 00357 assert(!k_complex_p(x)); 00358 return nucomp_s_new_internal(klass, x, ZERO); 00359 } 00360 00361 inline static VALUE 00362 f_complex_new_bang2(VALUE klass, VALUE x, VALUE y) 00363 { 00364 assert(!k_complex_p(x)); 00365 assert(!k_complex_p(y)); 00366 return nucomp_s_new_internal(klass, x, y); 00367 } 00368 00369 #ifdef CANONICALIZATION_FOR_MATHN 00370 #define CANON 00371 #endif 00372 00373 #ifdef CANON 00374 static int canonicalization = 0; 00375 00376 RUBY_FUNC_EXPORTED void 00377 nucomp_canonicalization(int f) 00378 { 00379 canonicalization = f; 00380 } 00381 #endif 00382 00383 inline static void 00384 nucomp_real_check(VALUE num) 00385 { 00386 switch (TYPE(num)) { 00387 case T_FIXNUM: 00388 case T_BIGNUM: 00389 case T_FLOAT: 00390 case T_RATIONAL: 00391 break; 00392 default: 00393 if (!k_numeric_p(num) || !f_real_p(num)) 00394 rb_raise(rb_eTypeError, "not a real"); 00395 } 00396 } 00397 00398 inline static VALUE 00399 nucomp_s_canonicalize_internal(VALUE klass, VALUE real, VALUE imag) 00400 { 00401 #ifdef CANON 00402 #define CL_CANON 00403 #ifdef CL_CANON 00404 if (k_exact_zero_p(imag) && canonicalization) 00405 return real; 00406 #else 00407 if (f_zero_p(imag) && canonicalization) 00408 return real; 00409 #endif 00410 #endif 00411 if (f_real_p(real) && f_real_p(imag)) 00412 return nucomp_s_new_internal(klass, real, imag); 00413 else if (f_real_p(real)) { 00414 get_dat1(imag); 00415 00416 return nucomp_s_new_internal(klass, 00417 f_sub(real, dat->imag), 00418 f_add(ZERO, dat->real)); 00419 } 00420 else if (f_real_p(imag)) { 00421 get_dat1(real); 00422 00423 return nucomp_s_new_internal(klass, 00424 dat->real, 00425 f_add(dat->imag, imag)); 00426 } 00427 else { 00428 get_dat2(real, imag); 00429 00430 return nucomp_s_new_internal(klass, 00431 f_sub(adat->real, bdat->imag), 00432 f_add(adat->imag, bdat->real)); 00433 } 00434 } 00435 00436 /* 00437 * call-seq: 00438 * Complex.rect(real[, imag]) -> complex 00439 * Complex.rectangular(real[, imag]) -> complex 00440 * 00441 * Returns a complex object which denotes the given rectangular form. 00442 */ 00443 static VALUE 00444 nucomp_s_new(int argc, VALUE *argv, VALUE klass) 00445 { 00446 VALUE real, imag; 00447 00448 switch (rb_scan_args(argc, argv, "11", &real, &imag)) { 00449 case 1: 00450 nucomp_real_check(real); 00451 imag = ZERO; 00452 break; 00453 default: 00454 nucomp_real_check(real); 00455 nucomp_real_check(imag); 00456 break; 00457 } 00458 00459 return nucomp_s_canonicalize_internal(klass, real, imag); 00460 } 00461 00462 inline static VALUE 00463 f_complex_new1(VALUE klass, VALUE x) 00464 { 00465 assert(!k_complex_p(x)); 00466 return nucomp_s_canonicalize_internal(klass, x, ZERO); 00467 } 00468 00469 inline static VALUE 00470 f_complex_new2(VALUE klass, VALUE x, VALUE y) 00471 { 00472 assert(!k_complex_p(x)); 00473 return nucomp_s_canonicalize_internal(klass, x, y); 00474 } 00475 00476 /* 00477 * call-seq: 00478 * Complex(x[, y]) -> numeric 00479 * 00480 * Returns x+i*y; 00481 */ 00482 static VALUE 00483 nucomp_f_complex(int argc, VALUE *argv, VALUE klass) 00484 { 00485 return rb_funcall2(rb_cComplex, id_convert, argc, argv); 00486 } 00487 00488 #define imp1(n) \ 00489 inline static VALUE \ 00490 m_##n##_bang(VALUE x)\ 00491 {\ 00492 return rb_math_##n(x);\ 00493 } 00494 00495 #define imp2(n) \ 00496 inline static VALUE \ 00497 m_##n##_bang(VALUE x, VALUE y)\ 00498 {\ 00499 return rb_math_##n(x, y);\ 00500 } 00501 00502 imp2(atan2) 00503 imp1(cos) 00504 imp1(cosh) 00505 imp1(exp) 00506 imp2(hypot) 00507 00508 #define m_hypot(x,y) m_hypot_bang((x),(y)) 00509 00510 static VALUE 00511 m_log_bang(VALUE x) 00512 { 00513 return rb_math_log(1, &x); 00514 } 00515 00516 imp1(sin) 00517 imp1(sinh) 00518 imp1(sqrt) 00519 00520 static VALUE 00521 m_cos(VALUE x) 00522 { 00523 if (f_real_p(x)) 00524 return m_cos_bang(x); 00525 { 00526 get_dat1(x); 00527 return f_complex_new2(rb_cComplex, 00528 f_mul(m_cos_bang(dat->real), 00529 m_cosh_bang(dat->imag)), 00530 f_mul(f_negate(m_sin_bang(dat->real)), 00531 m_sinh_bang(dat->imag))); 00532 } 00533 } 00534 00535 static VALUE 00536 m_sin(VALUE x) 00537 { 00538 if (f_real_p(x)) 00539 return m_sin_bang(x); 00540 { 00541 get_dat1(x); 00542 return f_complex_new2(rb_cComplex, 00543 f_mul(m_sin_bang(dat->real), 00544 m_cosh_bang(dat->imag)), 00545 f_mul(m_cos_bang(dat->real), 00546 m_sinh_bang(dat->imag))); 00547 } 00548 } 00549 00550 #if 0 00551 static VALUE 00552 m_sqrt(VALUE x) 00553 { 00554 if (f_real_p(x)) { 00555 if (f_positive_p(x)) 00556 return m_sqrt_bang(x); 00557 return f_complex_new2(rb_cComplex, ZERO, m_sqrt_bang(f_negate(x))); 00558 } 00559 else { 00560 get_dat1(x); 00561 00562 if (f_negative_p(dat->imag)) 00563 return f_conj(m_sqrt(f_conj(x))); 00564 else { 00565 VALUE a = f_abs(x); 00566 return f_complex_new2(rb_cComplex, 00567 m_sqrt_bang(f_div(f_add(a, dat->real), TWO)), 00568 m_sqrt_bang(f_div(f_sub(a, dat->real), TWO))); 00569 } 00570 } 00571 } 00572 #endif 00573 00574 inline static VALUE 00575 f_complex_polar(VALUE klass, VALUE x, VALUE y) 00576 { 00577 assert(!k_complex_p(x)); 00578 assert(!k_complex_p(y)); 00579 return nucomp_s_canonicalize_internal(klass, 00580 f_mul(x, m_cos(y)), 00581 f_mul(x, m_sin(y))); 00582 } 00583 00584 /* 00585 * call-seq: 00586 * Complex.polar(abs[, arg]) -> complex 00587 * 00588 * Returns a complex object which denotes the given polar form. 00589 * 00590 * Complex.polar(3, 0) #=> (3.0+0.0i) 00591 * Complex.polar(3, Math::PI/2) #=> (1.836909530733566e-16+3.0i) 00592 * Complex.polar(3, Math::PI) #=> (-3.0+3.673819061467132e-16i) 00593 * Complex.polar(3, -Math::PI/2) #=> (1.836909530733566e-16-3.0i) 00594 */ 00595 static VALUE 00596 nucomp_s_polar(int argc, VALUE *argv, VALUE klass) 00597 { 00598 VALUE abs, arg; 00599 00600 switch (rb_scan_args(argc, argv, "11", &abs, &arg)) { 00601 case 1: 00602 nucomp_real_check(abs); 00603 arg = ZERO; 00604 break; 00605 default: 00606 nucomp_real_check(abs); 00607 nucomp_real_check(arg); 00608 break; 00609 } 00610 return f_complex_polar(klass, abs, arg); 00611 } 00612 00613 /* 00614 * call-seq: 00615 * cmp.real -> real 00616 * 00617 * Returns the real part. 00618 */ 00619 static VALUE 00620 nucomp_real(VALUE self) 00621 { 00622 get_dat1(self); 00623 return dat->real; 00624 } 00625 00626 /* 00627 * call-seq: 00628 * cmp.imag -> real 00629 * cmp.imaginary -> real 00630 * 00631 * Returns the imaginary part. 00632 */ 00633 static VALUE 00634 nucomp_imag(VALUE self) 00635 { 00636 get_dat1(self); 00637 return dat->imag; 00638 } 00639 00640 /* 00641 * call-seq: 00642 * -cmp -> complex 00643 * 00644 * Returns negation of the value. 00645 */ 00646 static VALUE 00647 nucomp_negate(VALUE self) 00648 { 00649 get_dat1(self); 00650 return f_complex_new2(CLASS_OF(self), 00651 f_negate(dat->real), f_negate(dat->imag)); 00652 } 00653 00654 inline static VALUE 00655 f_addsub(VALUE self, VALUE other, 00656 VALUE (*func)(VALUE, VALUE), ID id) 00657 { 00658 if (k_complex_p(other)) { 00659 VALUE real, imag; 00660 00661 get_dat2(self, other); 00662 00663 real = (*func)(adat->real, bdat->real); 00664 imag = (*func)(adat->imag, bdat->imag); 00665 00666 return f_complex_new2(CLASS_OF(self), real, imag); 00667 } 00668 if (k_numeric_p(other) && f_real_p(other)) { 00669 get_dat1(self); 00670 00671 return f_complex_new2(CLASS_OF(self), 00672 (*func)(dat->real, other), dat->imag); 00673 } 00674 return rb_num_coerce_bin(self, other, id); 00675 } 00676 00677 /* 00678 * call-seq: 00679 * cmp + numeric -> complex 00680 * 00681 * Performs addition. 00682 */ 00683 static VALUE 00684 nucomp_add(VALUE self, VALUE other) 00685 { 00686 return f_addsub(self, other, f_add, '+'); 00687 } 00688 00689 /* 00690 * call-seq: 00691 * cmp - numeric -> complex 00692 * 00693 * Performs subtraction. 00694 */ 00695 static VALUE 00696 nucomp_sub(VALUE self, VALUE other) 00697 { 00698 return f_addsub(self, other, f_sub, '-'); 00699 } 00700 00701 /* 00702 * call-seq: 00703 * cmp * numeric -> complex 00704 * 00705 * Performs multiplication. 00706 */ 00707 static VALUE 00708 nucomp_mul(VALUE self, VALUE other) 00709 { 00710 if (k_complex_p(other)) { 00711 VALUE real, imag; 00712 00713 get_dat2(self, other); 00714 00715 real = f_sub(f_mul(adat->real, bdat->real), 00716 f_mul(adat->imag, bdat->imag)); 00717 imag = f_add(f_mul(adat->real, bdat->imag), 00718 f_mul(adat->imag, bdat->real)); 00719 00720 return f_complex_new2(CLASS_OF(self), real, imag); 00721 } 00722 if (k_numeric_p(other) && f_real_p(other)) { 00723 get_dat1(self); 00724 00725 return f_complex_new2(CLASS_OF(self), 00726 f_mul(dat->real, other), 00727 f_mul(dat->imag, other)); 00728 } 00729 return rb_num_coerce_bin(self, other, '*'); 00730 } 00731 00732 inline static VALUE 00733 f_divide(VALUE self, VALUE other, 00734 VALUE (*func)(VALUE, VALUE), ID id) 00735 { 00736 if (k_complex_p(other)) { 00737 int flo; 00738 get_dat2(self, other); 00739 00740 flo = (k_float_p(adat->real) || k_float_p(adat->imag) || 00741 k_float_p(bdat->real) || k_float_p(bdat->imag)); 00742 00743 if (f_gt_p(f_abs(bdat->real), f_abs(bdat->imag))) { 00744 VALUE r, n; 00745 00746 r = (*func)(bdat->imag, bdat->real); 00747 n = f_mul(bdat->real, f_add(ONE, f_mul(r, r))); 00748 if (flo) 00749 return f_complex_new2(CLASS_OF(self), 00750 (*func)(self, n), 00751 (*func)(f_negate(f_mul(self, r)), n)); 00752 return f_complex_new2(CLASS_OF(self), 00753 (*func)(f_add(adat->real, 00754 f_mul(adat->imag, r)), n), 00755 (*func)(f_sub(adat->imag, 00756 f_mul(adat->real, r)), n)); 00757 } 00758 else { 00759 VALUE r, n; 00760 00761 r = (*func)(bdat->real, bdat->imag); 00762 n = f_mul(bdat->imag, f_add(ONE, f_mul(r, r))); 00763 if (flo) 00764 return f_complex_new2(CLASS_OF(self), 00765 (*func)(f_mul(self, r), n), 00766 (*func)(f_negate(self), n)); 00767 return f_complex_new2(CLASS_OF(self), 00768 (*func)(f_add(f_mul(adat->real, r), 00769 adat->imag), n), 00770 (*func)(f_sub(f_mul(adat->imag, r), 00771 adat->real), n)); 00772 } 00773 } 00774 if (k_numeric_p(other) && f_real_p(other)) { 00775 get_dat1(self); 00776 00777 return f_complex_new2(CLASS_OF(self), 00778 (*func)(dat->real, other), 00779 (*func)(dat->imag, other)); 00780 } 00781 return rb_num_coerce_bin(self, other, id); 00782 } 00783 00784 #define rb_raise_zerodiv() rb_raise(rb_eZeroDivError, "divided by 0") 00785 00786 /* 00787 * call-seq: 00788 * cmp / numeric -> complex 00789 * cmp.quo(numeric) -> complex 00790 * 00791 * Performs division. 00792 * 00793 * For example: 00794 * 00795 * Complex(10.0) / 3 #=> (3.3333333333333335+(0/1)*i) 00796 * Complex(10) / 3 #=> ((10/3)+(0/1)*i) # not (3+0i) 00797 */ 00798 static VALUE 00799 nucomp_div(VALUE self, VALUE other) 00800 { 00801 return f_divide(self, other, f_quo, id_quo); 00802 } 00803 00804 #define nucomp_quo nucomp_div 00805 00806 /* 00807 * call-seq: 00808 * cmp.fdiv(numeric) -> complex 00809 * 00810 * Performs division as each part is a float, never returns a float. 00811 * 00812 * For example: 00813 * 00814 * Complex(11,22).fdiv(3) #=> (3.6666666666666665+7.333333333333333i) 00815 */ 00816 static VALUE 00817 nucomp_fdiv(VALUE self, VALUE other) 00818 { 00819 return f_divide(self, other, f_fdiv, id_fdiv); 00820 } 00821 00822 inline static VALUE 00823 f_reciprocal(VALUE x) 00824 { 00825 return f_quo(ONE, x); 00826 } 00827 00828 /* 00829 * call-seq: 00830 * cmp ** numeric -> complex 00831 * 00832 * Performs exponentiation. 00833 * 00834 * For example: 00835 * 00836 * Complex('i') ** 2 #=> (-1+0i) 00837 * Complex(-8) ** Rational(1,3) #=> (1.0000000000000002+1.7320508075688772i) 00838 */ 00839 static VALUE 00840 nucomp_expt(VALUE self, VALUE other) 00841 { 00842 if (k_numeric_p(other) && k_exact_zero_p(other)) 00843 return f_complex_new_bang1(CLASS_OF(self), ONE); 00844 00845 if (k_rational_p(other) && f_one_p(f_denominator(other))) 00846 other = f_numerator(other); /* c14n */ 00847 00848 if (k_complex_p(other)) { 00849 get_dat1(other); 00850 00851 if (k_exact_zero_p(dat->imag)) 00852 other = dat->real; /* c14n */ 00853 } 00854 00855 if (k_complex_p(other)) { 00856 VALUE r, theta, nr, ntheta; 00857 00858 get_dat1(other); 00859 00860 r = f_abs(self); 00861 theta = f_arg(self); 00862 00863 nr = m_exp_bang(f_sub(f_mul(dat->real, m_log_bang(r)), 00864 f_mul(dat->imag, theta))); 00865 ntheta = f_add(f_mul(theta, dat->real), 00866 f_mul(dat->imag, m_log_bang(r))); 00867 return f_complex_polar(CLASS_OF(self), nr, ntheta); 00868 } 00869 if (k_fixnum_p(other)) { 00870 if (f_gt_p(other, ZERO)) { 00871 VALUE x, z; 00872 long n; 00873 00874 x = self; 00875 z = x; 00876 n = FIX2LONG(other) - 1; 00877 00878 while (n) { 00879 long q, r; 00880 00881 while (1) { 00882 get_dat1(x); 00883 00884 q = n / 2; 00885 r = n % 2; 00886 00887 if (r) 00888 break; 00889 00890 x = f_complex_new2(CLASS_OF(self), 00891 f_sub(f_mul(dat->real, dat->real), 00892 f_mul(dat->imag, dat->imag)), 00893 f_mul(f_mul(TWO, dat->real), dat->imag)); 00894 n = q; 00895 } 00896 z = f_mul(z, x); 00897 n--; 00898 } 00899 return z; 00900 } 00901 return f_expt(f_reciprocal(self), f_negate(other)); 00902 } 00903 if (k_numeric_p(other) && f_real_p(other)) { 00904 VALUE r, theta; 00905 00906 if (k_bignum_p(other)) 00907 rb_warn("in a**b, b may be too big"); 00908 00909 r = f_abs(self); 00910 theta = f_arg(self); 00911 00912 return f_complex_polar(CLASS_OF(self), f_expt(r, other), 00913 f_mul(theta, other)); 00914 } 00915 return rb_num_coerce_bin(self, other, id_expt); 00916 } 00917 00918 /* 00919 * call-seq: 00920 * cmp == object -> true or false 00921 * 00922 * Returns true if cmp equals object numerically. 00923 */ 00924 static VALUE 00925 nucomp_eqeq_p(VALUE self, VALUE other) 00926 { 00927 if (k_complex_p(other)) { 00928 get_dat2(self, other); 00929 00930 return f_boolcast(f_eqeq_p(adat->real, bdat->real) && 00931 f_eqeq_p(adat->imag, bdat->imag)); 00932 } 00933 if (k_numeric_p(other) && f_real_p(other)) { 00934 get_dat1(self); 00935 00936 return f_boolcast(f_eqeq_p(dat->real, other) && f_zero_p(dat->imag)); 00937 } 00938 return f_eqeq_p(other, self); 00939 } 00940 00941 /* :nodoc: */ 00942 static VALUE 00943 nucomp_coerce(VALUE self, VALUE other) 00944 { 00945 if (k_numeric_p(other) && f_real_p(other)) 00946 return rb_assoc_new(f_complex_new_bang1(CLASS_OF(self), other), self); 00947 if (TYPE(other) == T_COMPLEX) 00948 return rb_assoc_new(other, self); 00949 00950 rb_raise(rb_eTypeError, "%s can't be coerced into %s", 00951 rb_obj_classname(other), rb_obj_classname(self)); 00952 return Qnil; 00953 } 00954 00955 /* 00956 * call-seq: 00957 * cmp.abs -> real 00958 * cmp.magnitude -> real 00959 * 00960 * Returns the absolute part of its polar form. 00961 */ 00962 static VALUE 00963 nucomp_abs(VALUE self) 00964 { 00965 get_dat1(self); 00966 00967 if (f_zero_p(dat->real)) { 00968 VALUE a = f_abs(dat->imag); 00969 if (k_float_p(dat->real) && !k_float_p(dat->imag)) 00970 a = f_to_f(a); 00971 return a; 00972 } 00973 if (f_zero_p(dat->imag)) { 00974 VALUE a = f_abs(dat->real); 00975 if (!k_float_p(dat->real) && k_float_p(dat->imag)) 00976 a = f_to_f(a); 00977 return a; 00978 } 00979 return m_hypot(dat->real, dat->imag); 00980 } 00981 00982 /* 00983 * call-seq: 00984 * cmp.abs2 -> real 00985 * 00986 * Returns square of the absolute value. 00987 */ 00988 static VALUE 00989 nucomp_abs2(VALUE self) 00990 { 00991 get_dat1(self); 00992 return f_add(f_mul(dat->real, dat->real), 00993 f_mul(dat->imag, dat->imag)); 00994 } 00995 00996 /* 00997 * call-seq: 00998 * cmp.arg -> float 00999 * cmp.angle -> float 01000 * cmp.phase -> float 01001 * 01002 * Returns the angle part of its polar form. 01003 * 01004 * Complex.polar(3, Math::PI/2).arg #=> 1.5707963267948966 01005 * 01006 */ 01007 static VALUE 01008 nucomp_arg(VALUE self) 01009 { 01010 get_dat1(self); 01011 return m_atan2_bang(dat->imag, dat->real); 01012 } 01013 01014 /* 01015 * call-seq: 01016 * cmp.rect -> array 01017 * cmp.rectangular -> array 01018 * 01019 * Returns an array; [cmp.real, cmp.imag]. 01020 */ 01021 static VALUE 01022 nucomp_rect(VALUE self) 01023 { 01024 get_dat1(self); 01025 return rb_assoc_new(dat->real, dat->imag); 01026 } 01027 01028 /* 01029 * call-seq: 01030 * cmp.polar -> array 01031 * 01032 * Returns an array; [cmp.abs, cmp.arg]. 01033 */ 01034 static VALUE 01035 nucomp_polar(VALUE self) 01036 { 01037 return rb_assoc_new(f_abs(self), f_arg(self)); 01038 } 01039 01040 /* 01041 * call-seq: 01042 * cmp.conj -> complex 01043 * cmp.conjugate -> complex 01044 * 01045 * Returns the complex conjugate. 01046 */ 01047 static VALUE 01048 nucomp_conj(VALUE self) 01049 { 01050 get_dat1(self); 01051 return f_complex_new2(CLASS_OF(self), dat->real, f_negate(dat->imag)); 01052 } 01053 01054 #if 0 01055 /* :nodoc: */ 01056 static VALUE 01057 nucomp_true(VALUE self) 01058 { 01059 return Qtrue; 01060 } 01061 #endif 01062 01063 /* 01064 * call-seq: 01065 * cmp.real? -> false 01066 * 01067 * Returns false. 01068 */ 01069 static VALUE 01070 nucomp_false(VALUE self) 01071 { 01072 return Qfalse; 01073 } 01074 01075 #if 0 01076 /* :nodoc: */ 01077 static VALUE 01078 nucomp_exact_p(VALUE self) 01079 { 01080 get_dat1(self); 01081 return f_boolcast(k_exact_p(dat->real) && k_exact_p(dat->imag)); 01082 } 01083 01084 /* :nodoc: */ 01085 static VALUE 01086 nucomp_inexact_p(VALUE self) 01087 { 01088 return f_boolcast(!nucomp_exact_p(self)); 01089 } 01090 #endif 01091 01092 /* 01093 * call-seq: 01094 * cmp.denominator -> integer 01095 * 01096 * Returns the denominator (lcm of both denominator - real and imag). 01097 * 01098 * See numerator. 01099 */ 01100 static VALUE 01101 nucomp_denominator(VALUE self) 01102 { 01103 get_dat1(self); 01104 return rb_lcm(f_denominator(dat->real), f_denominator(dat->imag)); 01105 } 01106 01107 /* 01108 * call-seq: 01109 * cmp.numerator -> numeric 01110 * 01111 * Returns the numerator. 01112 * 01113 * For example: 01114 * 01115 * 1 2 3+4i <- numerator 01116 * - + -i -> ---- 01117 * 2 3 6 <- denominator 01118 * 01119 * c = Complex('1/2+2/3i') #=> ((1/2)+(2/3)*i) 01120 * n = c.numerator #=> (3+4i) 01121 * d = c.denominator #=> 6 01122 * n / d #=> ((1/2)+(2/3)*i) 01123 * Complex(Rational(n.real, d), Rational(n.imag, d)) 01124 * #=> ((1/2)+(2/3)*i) 01125 * See denominator. 01126 */ 01127 static VALUE 01128 nucomp_numerator(VALUE self) 01129 { 01130 VALUE cd; 01131 01132 get_dat1(self); 01133 01134 cd = f_denominator(self); 01135 return f_complex_new2(CLASS_OF(self), 01136 f_mul(f_numerator(dat->real), 01137 f_div(cd, f_denominator(dat->real))), 01138 f_mul(f_numerator(dat->imag), 01139 f_div(cd, f_denominator(dat->imag)))); 01140 } 01141 01142 /* :nodoc: */ 01143 static VALUE 01144 nucomp_hash(VALUE self) 01145 { 01146 st_index_t v, h[2]; 01147 VALUE n; 01148 01149 get_dat1(self); 01150 n = rb_hash(dat->real); 01151 h[0] = NUM2LONG(n); 01152 n = rb_hash(dat->imag); 01153 h[1] = NUM2LONG(n); 01154 v = rb_memhash(h, sizeof(h)); 01155 return LONG2FIX(v); 01156 } 01157 01158 /* :nodoc: */ 01159 static VALUE 01160 nucomp_eql_p(VALUE self, VALUE other) 01161 { 01162 if (k_complex_p(other)) { 01163 get_dat2(self, other); 01164 01165 return f_boolcast((CLASS_OF(adat->real) == CLASS_OF(bdat->real)) && 01166 (CLASS_OF(adat->imag) == CLASS_OF(bdat->imag)) && 01167 f_eqeq_p(self, other)); 01168 01169 } 01170 return Qfalse; 01171 } 01172 01173 inline static VALUE 01174 f_signbit(VALUE x) 01175 { 01176 #if defined(HAVE_SIGNBIT) && defined(__GNUC__) && defined(__sun__) && \ 01177 !defined(signbit) 01178 extern int signbit(double); 01179 #endif 01180 switch (TYPE(x)) { 01181 case T_FLOAT: { 01182 double f = RFLOAT_VALUE(x); 01183 return f_boolcast(!isnan(f) && signbit(f)); 01184 } 01185 } 01186 return f_negative_p(x); 01187 } 01188 01189 inline static VALUE 01190 f_tpositive_p(VALUE x) 01191 { 01192 return f_boolcast(!f_signbit(x)); 01193 } 01194 01195 static VALUE 01196 f_format(VALUE self, VALUE (*func)(VALUE)) 01197 { 01198 VALUE s, impos; 01199 01200 get_dat1(self); 01201 01202 impos = f_tpositive_p(dat->imag); 01203 01204 s = (*func)(dat->real); 01205 rb_str_cat2(s, !impos ? "-" : "+"); 01206 01207 rb_str_concat(s, (*func)(f_abs(dat->imag))); 01208 if (!rb_isdigit(RSTRING_PTR(s)[RSTRING_LEN(s) - 1])) 01209 rb_str_cat2(s, "*"); 01210 rb_str_cat2(s, "i"); 01211 01212 return s; 01213 } 01214 01215 /* 01216 * call-seq: 01217 * cmp.to_s -> string 01218 * 01219 * Returns the value as a string. 01220 */ 01221 static VALUE 01222 nucomp_to_s(VALUE self) 01223 { 01224 return f_format(self, f_to_s); 01225 } 01226 01227 /* 01228 * call-seq: 01229 * cmp.inspect -> string 01230 * 01231 * Returns the value as a string for inspection. 01232 */ 01233 static VALUE 01234 nucomp_inspect(VALUE self) 01235 { 01236 VALUE s; 01237 01238 s = rb_usascii_str_new2("("); 01239 rb_str_concat(s, f_format(self, f_inspect)); 01240 rb_str_cat2(s, ")"); 01241 01242 return s; 01243 } 01244 01245 /* :nodoc: */ 01246 static VALUE 01247 nucomp_marshal_dump(VALUE self) 01248 { 01249 VALUE a; 01250 get_dat1(self); 01251 01252 a = rb_assoc_new(dat->real, dat->imag); 01253 rb_copy_generic_ivar(a, self); 01254 return a; 01255 } 01256 01257 /* :nodoc: */ 01258 static VALUE 01259 nucomp_marshal_load(VALUE self, VALUE a) 01260 { 01261 get_dat1(self); 01262 Check_Type(a, T_ARRAY); 01263 if (RARRAY_LEN(a) != 2) 01264 rb_raise(rb_eArgError, "marshaled complex must have an array whose length is 2 but %ld", RARRAY_LEN(a)); 01265 dat->real = RARRAY_PTR(a)[0]; 01266 dat->imag = RARRAY_PTR(a)[1]; 01267 rb_copy_generic_ivar(self, a); 01268 return self; 01269 } 01270 01271 /* --- */ 01272 01273 VALUE 01274 rb_complex_raw(VALUE x, VALUE y) 01275 { 01276 return nucomp_s_new_internal(rb_cComplex, x, y); 01277 } 01278 01279 VALUE 01280 rb_complex_new(VALUE x, VALUE y) 01281 { 01282 return nucomp_s_canonicalize_internal(rb_cComplex, x, y); 01283 } 01284 01285 VALUE 01286 rb_complex_polar(VALUE x, VALUE y) 01287 { 01288 return f_complex_polar(rb_cComplex, x, y); 01289 } 01290 01291 static VALUE nucomp_s_convert(int argc, VALUE *argv, VALUE klass); 01292 01293 VALUE 01294 rb_Complex(VALUE x, VALUE y) 01295 { 01296 VALUE a[2]; 01297 a[0] = x; 01298 a[1] = y; 01299 return nucomp_s_convert(2, a, rb_cComplex); 01300 } 01301 01302 /* 01303 * call-seq: 01304 * cmp.to_i -> integer 01305 * 01306 * Returns the value as an integer if possible. 01307 */ 01308 static VALUE 01309 nucomp_to_i(VALUE self) 01310 { 01311 get_dat1(self); 01312 01313 if (k_inexact_p(dat->imag) || f_nonzero_p(dat->imag)) { 01314 VALUE s = f_to_s(self); 01315 rb_raise(rb_eRangeError, "can't convert %s into Integer", 01316 StringValuePtr(s)); 01317 } 01318 return f_to_i(dat->real); 01319 } 01320 01321 /* 01322 * call-seq: 01323 * cmp.to_f -> float 01324 * 01325 * Returns the value as a float if possible. 01326 */ 01327 static VALUE 01328 nucomp_to_f(VALUE self) 01329 { 01330 get_dat1(self); 01331 01332 if (k_inexact_p(dat->imag) || f_nonzero_p(dat->imag)) { 01333 VALUE s = f_to_s(self); 01334 rb_raise(rb_eRangeError, "can't convert %s into Float", 01335 StringValuePtr(s)); 01336 } 01337 return f_to_f(dat->real); 01338 } 01339 01340 /* 01341 * call-seq: 01342 * cmp.to_r -> rational 01343 * 01344 * If the imaginary part is exactly 0, returns the real part as a Rational, 01345 * otherwise a RangeError is raised. 01346 */ 01347 static VALUE 01348 nucomp_to_r(VALUE self) 01349 { 01350 get_dat1(self); 01351 01352 if (k_inexact_p(dat->imag) || f_nonzero_p(dat->imag)) { 01353 VALUE s = f_to_s(self); 01354 rb_raise(rb_eRangeError, "can't convert %s into Rational", 01355 StringValuePtr(s)); 01356 } 01357 return f_to_r(dat->real); 01358 } 01359 01360 /* 01361 * call-seq: 01362 * cmp.rationalize([eps]) -> rational 01363 * 01364 * If the imaginary part is exactly 0, returns the real part as a Rational, 01365 * otherwise a RangeError is raised. 01366 */ 01367 static VALUE 01368 nucomp_rationalize(int argc, VALUE *argv, VALUE self) 01369 { 01370 get_dat1(self); 01371 01372 rb_scan_args(argc, argv, "01", NULL); 01373 01374 if (k_inexact_p(dat->imag) || f_nonzero_p(dat->imag)) { 01375 VALUE s = f_to_s(self); 01376 rb_raise(rb_eRangeError, "can't convert %s into Rational", 01377 StringValuePtr(s)); 01378 } 01379 return rb_funcall2(dat->real, rb_intern("rationalize"), argc, argv); 01380 } 01381 01382 /* 01383 * call-seq: 01384 * nil.to_c -> (0+0i) 01385 * 01386 * Returns zero as a complex. 01387 */ 01388 static VALUE 01389 nilclass_to_c(VALUE self) 01390 { 01391 return rb_complex_new1(INT2FIX(0)); 01392 } 01393 01394 /* 01395 * call-seq: 01396 * num.to_c -> complex 01397 * 01398 * Returns the value as a complex. 01399 */ 01400 static VALUE 01401 numeric_to_c(VALUE self) 01402 { 01403 return rb_complex_new1(self); 01404 } 01405 01406 static VALUE comp_pat0, comp_pat1, comp_pat2, a_slash, a_dot_and_an_e, 01407 null_string, underscores_pat, an_underscore; 01408 01409 #define WS "\\s*" 01410 #define DIGITS "(?:[0-9](?:_[0-9]|[0-9])*)" 01411 #define NUMERATOR "(?:" DIGITS "?\\.)?" DIGITS "(?:[eE][-+]?" DIGITS ")?" 01412 #define DENOMINATOR DIGITS 01413 #define NUMBER "[-+]?" NUMERATOR "(?:\\/" DENOMINATOR ")?" 01414 #define NUMBERNOS NUMERATOR "(?:\\/" DENOMINATOR ")?" 01415 #define PATTERN0 "\\A" WS "(" NUMBER ")@(" NUMBER ")" WS 01416 #define PATTERN1 "\\A" WS "([-+])?(" NUMBER ")?[iIjJ]" WS 01417 #define PATTERN2 "\\A" WS "(" NUMBER ")(([-+])(" NUMBERNOS ")?[iIjJ])?" WS 01418 01419 static void 01420 make_patterns(void) 01421 { 01422 static const char comp_pat0_source[] = PATTERN0; 01423 static const char comp_pat1_source[] = PATTERN1; 01424 static const char comp_pat2_source[] = PATTERN2; 01425 static const char underscores_pat_source[] = "_+"; 01426 01427 if (comp_pat0) return; 01428 01429 comp_pat0 = rb_reg_new(comp_pat0_source, sizeof comp_pat0_source - 1, 0); 01430 rb_gc_register_mark_object(comp_pat0); 01431 01432 comp_pat1 = rb_reg_new(comp_pat1_source, sizeof comp_pat1_source - 1, 0); 01433 rb_gc_register_mark_object(comp_pat1); 01434 01435 comp_pat2 = rb_reg_new(comp_pat2_source, sizeof comp_pat2_source - 1, 0); 01436 rb_gc_register_mark_object(comp_pat2); 01437 01438 a_slash = rb_usascii_str_new2("/"); 01439 rb_gc_register_mark_object(a_slash); 01440 01441 a_dot_and_an_e = rb_usascii_str_new2(".eE"); 01442 rb_gc_register_mark_object(a_dot_and_an_e); 01443 01444 null_string = rb_usascii_str_new2(""); 01445 rb_gc_register_mark_object(null_string); 01446 01447 underscores_pat = rb_reg_new(underscores_pat_source, 01448 sizeof underscores_pat_source - 1, 0); 01449 rb_gc_register_mark_object(underscores_pat); 01450 01451 an_underscore = rb_usascii_str_new2("_"); 01452 rb_gc_register_mark_object(an_underscore); 01453 } 01454 01455 #define id_match rb_intern("match") 01456 #define f_match(x,y) rb_funcall((x), id_match, 1, (y)) 01457 01458 #define id_gsub_bang rb_intern("gsub!") 01459 #define f_gsub_bang(x,y,z) rb_funcall((x), id_gsub_bang, 2, (y), (z)) 01460 01461 static VALUE 01462 string_to_c_internal(VALUE self) 01463 { 01464 VALUE s; 01465 01466 s = self; 01467 01468 if (RSTRING_LEN(s) == 0) 01469 return rb_assoc_new(Qnil, self); 01470 01471 { 01472 VALUE m, sr, si, re, r, i; 01473 int po; 01474 01475 m = f_match(comp_pat0, s); 01476 if (!NIL_P(m)) { 01477 sr = rb_reg_nth_match(1, m); 01478 si = rb_reg_nth_match(2, m); 01479 re = rb_reg_match_post(m); 01480 po = 1; 01481 } 01482 if (NIL_P(m)) { 01483 m = f_match(comp_pat1, s); 01484 if (!NIL_P(m)) { 01485 sr = Qnil; 01486 si = rb_reg_nth_match(1, m); 01487 if (NIL_P(si)) 01488 si = rb_usascii_str_new2(""); 01489 { 01490 VALUE t; 01491 01492 t = rb_reg_nth_match(2, m); 01493 if (NIL_P(t)) 01494 t = rb_usascii_str_new2("1"); 01495 rb_str_concat(si, t); 01496 } 01497 re = rb_reg_match_post(m); 01498 po = 0; 01499 } 01500 } 01501 if (NIL_P(m)) { 01502 m = f_match(comp_pat2, s); 01503 if (NIL_P(m)) 01504 return rb_assoc_new(Qnil, self); 01505 sr = rb_reg_nth_match(1, m); 01506 if (NIL_P(rb_reg_nth_match(2, m))) 01507 si = Qnil; 01508 else { 01509 VALUE t; 01510 01511 si = rb_reg_nth_match(3, m); 01512 t = rb_reg_nth_match(4, m); 01513 if (NIL_P(t)) 01514 t = rb_usascii_str_new2("1"); 01515 rb_str_concat(si, t); 01516 } 01517 re = rb_reg_match_post(m); 01518 po = 0; 01519 } 01520 r = INT2FIX(0); 01521 i = INT2FIX(0); 01522 if (!NIL_P(sr)) { 01523 if (strchr(RSTRING_PTR(sr), '/')) 01524 r = f_to_r(sr); 01525 else if (strpbrk(RSTRING_PTR(sr), ".eE")) 01526 r = f_to_f(sr); 01527 else 01528 r = f_to_i(sr); 01529 } 01530 if (!NIL_P(si)) { 01531 if (strchr(RSTRING_PTR(si), '/')) 01532 i = f_to_r(si); 01533 else if (strpbrk(RSTRING_PTR(si), ".eE")) 01534 i = f_to_f(si); 01535 else 01536 i = f_to_i(si); 01537 } 01538 if (po) 01539 return rb_assoc_new(rb_complex_polar(r, i), re); 01540 else 01541 return rb_assoc_new(rb_complex_new2(r, i), re); 01542 } 01543 } 01544 01545 static VALUE 01546 string_to_c_strict(VALUE self) 01547 { 01548 VALUE a = string_to_c_internal(self); 01549 if (NIL_P(RARRAY_PTR(a)[0]) || RSTRING_LEN(RARRAY_PTR(a)[1]) > 0) { 01550 VALUE s = f_inspect(self); 01551 rb_raise(rb_eArgError, "invalid value for convert(): %s", 01552 StringValuePtr(s)); 01553 } 01554 return RARRAY_PTR(a)[0]; 01555 } 01556 01557 #define id_gsub rb_intern("gsub") 01558 #define f_gsub(x,y,z) rb_funcall((x), id_gsub, 2, (y), (z)) 01559 01560 /* 01561 * call-seq: 01562 * str.to_c -> complex 01563 * 01564 * Returns a complex which denotes the string form. The parser 01565 * ignores leading whitespaces and trailing garbage. Any digit 01566 * sequences can be separated by an underscore. Returns zero for null 01567 * or garbage string. 01568 * 01569 * For example: 01570 * 01571 * '9'.to_c #=> (9+0i) 01572 * '2.5'.to_c #=> (2.5+0i) 01573 * '2.5/1'.to_c #=> ((5/2)+0i) 01574 * '-3/2'.to_c #=> ((-3/2)+0i) 01575 * '-i'.to_c #=> (0-1i) 01576 * '45i'.to_c #=> (0+45i) 01577 * '3-4i'.to_c #=> (3-4i) 01578 * '-4e2-4e-2i'.to_c #=> (-400.0-0.04i) 01579 * '-0.0-0.0i'.to_c #=> (-0.0-0.0i) 01580 * '1/2+3/4i'.to_c #=> ((1/2)+(3/4)*i) 01581 * 'ruby'.to_c #=> (0+0i) 01582 */ 01583 static VALUE 01584 string_to_c(VALUE self) 01585 { 01586 VALUE s, a, backref; 01587 01588 backref = rb_backref_get(); 01589 rb_match_busy(backref); 01590 01591 s = f_gsub(self, underscores_pat, an_underscore); 01592 a = string_to_c_internal(s); 01593 01594 rb_backref_set(backref); 01595 01596 if (!NIL_P(RARRAY_PTR(a)[0])) 01597 return RARRAY_PTR(a)[0]; 01598 return rb_complex_new1(INT2FIX(0)); 01599 } 01600 01601 static VALUE 01602 nucomp_s_convert(int argc, VALUE *argv, VALUE klass) 01603 { 01604 VALUE a1, a2, backref; 01605 01606 rb_scan_args(argc, argv, "11", &a1, &a2); 01607 01608 if (NIL_P(a1) || (argc == 2 && NIL_P(a2))) 01609 rb_raise(rb_eTypeError, "can't convert nil into Complex"); 01610 01611 backref = rb_backref_get(); 01612 rb_match_busy(backref); 01613 01614 switch (TYPE(a1)) { 01615 case T_FIXNUM: 01616 case T_BIGNUM: 01617 case T_FLOAT: 01618 break; 01619 case T_STRING: 01620 a1 = string_to_c_strict(a1); 01621 break; 01622 } 01623 01624 switch (TYPE(a2)) { 01625 case T_FIXNUM: 01626 case T_BIGNUM: 01627 case T_FLOAT: 01628 break; 01629 case T_STRING: 01630 a2 = string_to_c_strict(a2); 01631 break; 01632 } 01633 01634 rb_backref_set(backref); 01635 01636 switch (TYPE(a1)) { 01637 case T_COMPLEX: 01638 { 01639 get_dat1(a1); 01640 01641 if (k_exact_zero_p(dat->imag)) 01642 a1 = dat->real; 01643 } 01644 } 01645 01646 switch (TYPE(a2)) { 01647 case T_COMPLEX: 01648 { 01649 get_dat1(a2); 01650 01651 if (k_exact_zero_p(dat->imag)) 01652 a2 = dat->real; 01653 } 01654 } 01655 01656 switch (TYPE(a1)) { 01657 case T_COMPLEX: 01658 if (argc == 1 || (k_exact_zero_p(a2))) 01659 return a1; 01660 } 01661 01662 if (argc == 1) { 01663 if (k_numeric_p(a1) && !f_real_p(a1)) 01664 return a1; 01665 /* should raise exception for consistency */ 01666 if (!k_numeric_p(a1)) 01667 return rb_convert_type(a1, T_COMPLEX, "Complex", "to_c"); 01668 } 01669 else { 01670 if ((k_numeric_p(a1) && k_numeric_p(a2)) && 01671 (!f_real_p(a1) || !f_real_p(a2))) 01672 return f_add(a1, 01673 f_mul(a2, 01674 f_complex_new_bang2(rb_cComplex, ZERO, ONE))); 01675 } 01676 01677 { 01678 VALUE argv2[2]; 01679 argv2[0] = a1; 01680 argv2[1] = a2; 01681 return nucomp_s_new(argc, argv2, klass); 01682 } 01683 } 01684 01685 /* --- */ 01686 01687 /* 01688 * call-seq: 01689 * num.real -> self 01690 * 01691 * Returns self. 01692 */ 01693 static VALUE 01694 numeric_real(VALUE self) 01695 { 01696 return self; 01697 } 01698 01699 /* 01700 * call-seq: 01701 * num.imag -> 0 01702 * num.imaginary -> 0 01703 * 01704 * Returns zero. 01705 */ 01706 static VALUE 01707 numeric_imag(VALUE self) 01708 { 01709 return INT2FIX(0); 01710 } 01711 01712 /* 01713 * call-seq: 01714 * num.abs2 -> real 01715 * 01716 * Returns square of self. 01717 */ 01718 static VALUE 01719 numeric_abs2(VALUE self) 01720 { 01721 return f_mul(self, self); 01722 } 01723 01724 #define id_PI rb_intern("PI") 01725 01726 /* 01727 * call-seq: 01728 * num.arg -> 0 or float 01729 * num.angle -> 0 or float 01730 * num.phase -> 0 or float 01731 * 01732 * Returns 0 if the value is positive, pi otherwise. 01733 */ 01734 static VALUE 01735 numeric_arg(VALUE self) 01736 { 01737 if (f_positive_p(self)) 01738 return INT2FIX(0); 01739 return rb_const_get(rb_mMath, id_PI); 01740 } 01741 01742 /* 01743 * call-seq: 01744 * num.rect -> array 01745 * 01746 * Returns an array; [num, 0]. 01747 */ 01748 static VALUE 01749 numeric_rect(VALUE self) 01750 { 01751 return rb_assoc_new(self, INT2FIX(0)); 01752 } 01753 01754 /* 01755 * call-seq: 01756 * num.polar -> array 01757 * 01758 * Returns an array; [num.abs, num.arg]. 01759 */ 01760 static VALUE 01761 numeric_polar(VALUE self) 01762 { 01763 return rb_assoc_new(f_abs(self), f_arg(self)); 01764 } 01765 01766 /* 01767 * call-seq: 01768 * num.conj -> self 01769 * num.conjugate -> self 01770 * 01771 * Returns self. 01772 */ 01773 static VALUE 01774 numeric_conj(VALUE self) 01775 { 01776 return self; 01777 } 01778 01779 /* 01780 * call-seq: 01781 * flo.arg -> 0 or float 01782 * flo.angle -> 0 or float 01783 * flo.phase -> 0 or float 01784 * 01785 * Returns 0 if the value is positive, pi otherwise. 01786 */ 01787 static VALUE 01788 float_arg(VALUE self) 01789 { 01790 if (isnan(RFLOAT_VALUE(self))) 01791 return self; 01792 if (f_tpositive_p(self)) 01793 return INT2FIX(0); 01794 return rb_const_get(rb_mMath, id_PI); 01795 } 01796 01797 /* 01798 * A complex number can be represented as a paired real number with 01799 * imaginary unit; a+bi. Where a is real part, b is imaginary part 01800 * and i is imaginary unit. Real a equals complex a+0i 01801 * mathematically. 01802 * 01803 * In ruby, you can create complex object with Complex, Complex::rect, 01804 * Complex::polar or to_c method. 01805 * 01806 * Complex(1) #=> (1+0i) 01807 * Complex(2, 3) #=> (2+3i) 01808 * Complex.polar(2, 3) #=> (-1.9799849932008908+0.2822400161197344i) 01809 * 3.to_c #=> (3+0i) 01810 * 01811 * You can also create complex object from floating-point numbers or 01812 * strings. 01813 * 01814 * Complex(0.3) #=> (0.3+0i) 01815 * Complex('0.3-0.5i') #=> (0.3-0.5i) 01816 * Complex('2/3+3/4i') #=> ((2/3)+(3/4)*i) 01817 * Complex('1@2') #=> (-0.4161468365471424+0.9092974268256817i) 01818 * 01819 * 0.3.to_c #=> (0.3+0i) 01820 * '0.3-0.5i'.to_c #=> (0.3-0.5i) 01821 * '2/3+3/4i'.to_c #=> ((2/3)+(3/4)*i) 01822 * '1@2'.to_c #=> (-0.4161468365471424+0.9092974268256817i) 01823 * 01824 * A complex object is either an exact or an inexact number. 01825 * 01826 * Complex(1, 1) / 2 #=> ((1/2)+(1/2)*i) 01827 * Complex(1, 1) / 2.0 #=> (0.5+0.5i) 01828 */ 01829 void 01830 Init_Complex(void) 01831 { 01832 #undef rb_intern 01833 #define rb_intern(str) rb_intern_const(str) 01834 01835 assert(fprintf(stderr, "assert() is now active\n")); 01836 01837 id_abs = rb_intern("abs"); 01838 id_abs2 = rb_intern("abs2"); 01839 id_arg = rb_intern("arg"); 01840 id_cmp = rb_intern("<=>"); 01841 id_conj = rb_intern("conj"); 01842 id_convert = rb_intern("convert"); 01843 id_denominator = rb_intern("denominator"); 01844 id_divmod = rb_intern("divmod"); 01845 id_eqeq_p = rb_intern("=="); 01846 id_expt = rb_intern("**"); 01847 id_fdiv = rb_intern("fdiv"); 01848 id_floor = rb_intern("floor"); 01849 id_idiv = rb_intern("div"); 01850 id_imag = rb_intern("imag"); 01851 id_inspect = rb_intern("inspect"); 01852 id_negate = rb_intern("-@"); 01853 id_numerator = rb_intern("numerator"); 01854 id_quo = rb_intern("quo"); 01855 id_real = rb_intern("real"); 01856 id_real_p = rb_intern("real?"); 01857 id_to_f = rb_intern("to_f"); 01858 id_to_i = rb_intern("to_i"); 01859 id_to_r = rb_intern("to_r"); 01860 id_to_s = rb_intern("to_s"); 01861 01862 rb_cComplex = rb_define_class("Complex", rb_cNumeric); 01863 01864 rb_define_alloc_func(rb_cComplex, nucomp_s_alloc); 01865 rb_undef_method(CLASS_OF(rb_cComplex), "allocate"); 01866 01867 #if 0 01868 rb_define_private_method(CLASS_OF(rb_cComplex), "new!", nucomp_s_new_bang, -1); 01869 rb_define_private_method(CLASS_OF(rb_cComplex), "new", nucomp_s_new, -1); 01870 #else 01871 rb_undef_method(CLASS_OF(rb_cComplex), "new"); 01872 #endif 01873 01874 rb_define_singleton_method(rb_cComplex, "rectangular", nucomp_s_new, -1); 01875 rb_define_singleton_method(rb_cComplex, "rect", nucomp_s_new, -1); 01876 rb_define_singleton_method(rb_cComplex, "polar", nucomp_s_polar, -1); 01877 01878 rb_define_global_function("Complex", nucomp_f_complex, -1); 01879 01880 rb_undef_method(rb_cComplex, "%"); 01881 rb_undef_method(rb_cComplex, "<"); 01882 rb_undef_method(rb_cComplex, "<="); 01883 rb_undef_method(rb_cComplex, "<=>"); 01884 rb_undef_method(rb_cComplex, ">"); 01885 rb_undef_method(rb_cComplex, ">="); 01886 rb_undef_method(rb_cComplex, "between?"); 01887 rb_undef_method(rb_cComplex, "div"); 01888 rb_undef_method(rb_cComplex, "divmod"); 01889 rb_undef_method(rb_cComplex, "floor"); 01890 rb_undef_method(rb_cComplex, "ceil"); 01891 rb_undef_method(rb_cComplex, "modulo"); 01892 rb_undef_method(rb_cComplex, "remainder"); 01893 rb_undef_method(rb_cComplex, "round"); 01894 rb_undef_method(rb_cComplex, "step"); 01895 rb_undef_method(rb_cComplex, "truncate"); 01896 rb_undef_method(rb_cComplex, "i"); 01897 01898 #if 0 /* NUBY */ 01899 rb_undef_method(rb_cComplex, "//"); 01900 #endif 01901 01902 rb_define_method(rb_cComplex, "real", nucomp_real, 0); 01903 rb_define_method(rb_cComplex, "imaginary", nucomp_imag, 0); 01904 rb_define_method(rb_cComplex, "imag", nucomp_imag, 0); 01905 01906 rb_define_method(rb_cComplex, "-@", nucomp_negate, 0); 01907 rb_define_method(rb_cComplex, "+", nucomp_add, 1); 01908 rb_define_method(rb_cComplex, "-", nucomp_sub, 1); 01909 rb_define_method(rb_cComplex, "*", nucomp_mul, 1); 01910 rb_define_method(rb_cComplex, "/", nucomp_div, 1); 01911 rb_define_method(rb_cComplex, "quo", nucomp_quo, 1); 01912 rb_define_method(rb_cComplex, "fdiv", nucomp_fdiv, 1); 01913 rb_define_method(rb_cComplex, "**", nucomp_expt, 1); 01914 01915 rb_define_method(rb_cComplex, "==", nucomp_eqeq_p, 1); 01916 rb_define_method(rb_cComplex, "coerce", nucomp_coerce, 1); 01917 01918 rb_define_method(rb_cComplex, "abs", nucomp_abs, 0); 01919 rb_define_method(rb_cComplex, "magnitude", nucomp_abs, 0); 01920 rb_define_method(rb_cComplex, "abs2", nucomp_abs2, 0); 01921 rb_define_method(rb_cComplex, "arg", nucomp_arg, 0); 01922 rb_define_method(rb_cComplex, "angle", nucomp_arg, 0); 01923 rb_define_method(rb_cComplex, "phase", nucomp_arg, 0); 01924 rb_define_method(rb_cComplex, "rectangular", nucomp_rect, 0); 01925 rb_define_method(rb_cComplex, "rect", nucomp_rect, 0); 01926 rb_define_method(rb_cComplex, "polar", nucomp_polar, 0); 01927 rb_define_method(rb_cComplex, "conjugate", nucomp_conj, 0); 01928 rb_define_method(rb_cComplex, "conj", nucomp_conj, 0); 01929 #if 0 01930 rb_define_method(rb_cComplex, "~", nucomp_conj, 0); /* gcc */ 01931 #endif 01932 01933 rb_define_method(rb_cComplex, "real?", nucomp_false, 0); 01934 #if 0 01935 rb_define_method(rb_cComplex, "complex?", nucomp_true, 0); 01936 rb_define_method(rb_cComplex, "exact?", nucomp_exact_p, 0); 01937 rb_define_method(rb_cComplex, "inexact?", nucomp_inexact_p, 0); 01938 #endif 01939 01940 rb_define_method(rb_cComplex, "numerator", nucomp_numerator, 0); 01941 rb_define_method(rb_cComplex, "denominator", nucomp_denominator, 0); 01942 01943 rb_define_method(rb_cComplex, "hash", nucomp_hash, 0); 01944 rb_define_method(rb_cComplex, "eql?", nucomp_eql_p, 1); 01945 01946 rb_define_method(rb_cComplex, "to_s", nucomp_to_s, 0); 01947 rb_define_method(rb_cComplex, "inspect", nucomp_inspect, 0); 01948 01949 rb_define_method(rb_cComplex, "marshal_dump", nucomp_marshal_dump, 0); 01950 rb_define_method(rb_cComplex, "marshal_load", nucomp_marshal_load, 1); 01951 01952 /* --- */ 01953 01954 rb_define_method(rb_cComplex, "to_i", nucomp_to_i, 0); 01955 rb_define_method(rb_cComplex, "to_f", nucomp_to_f, 0); 01956 rb_define_method(rb_cComplex, "to_r", nucomp_to_r, 0); 01957 rb_define_method(rb_cComplex, "rationalize", nucomp_rationalize, -1); 01958 rb_define_method(rb_cNilClass, "to_c", nilclass_to_c, 0); 01959 rb_define_method(rb_cNumeric, "to_c", numeric_to_c, 0); 01960 01961 make_patterns(); 01962 01963 rb_define_method(rb_cString, "to_c", string_to_c, 0); 01964 01965 rb_define_private_method(CLASS_OF(rb_cComplex), "convert", nucomp_s_convert, -1); 01966 01967 /* --- */ 01968 01969 rb_define_method(rb_cNumeric, "real", numeric_real, 0); 01970 rb_define_method(rb_cNumeric, "imaginary", numeric_imag, 0); 01971 rb_define_method(rb_cNumeric, "imag", numeric_imag, 0); 01972 rb_define_method(rb_cNumeric, "abs2", numeric_abs2, 0); 01973 rb_define_method(rb_cNumeric, "arg", numeric_arg, 0); 01974 rb_define_method(rb_cNumeric, "angle", numeric_arg, 0); 01975 rb_define_method(rb_cNumeric, "phase", numeric_arg, 0); 01976 rb_define_method(rb_cNumeric, "rectangular", numeric_rect, 0); 01977 rb_define_method(rb_cNumeric, "rect", numeric_rect, 0); 01978 rb_define_method(rb_cNumeric, "polar", numeric_polar, 0); 01979 rb_define_method(rb_cNumeric, "conjugate", numeric_conj, 0); 01980 rb_define_method(rb_cNumeric, "conj", numeric_conj, 0); 01981 01982 rb_define_method(rb_cFloat, "arg", float_arg, 0); 01983 rb_define_method(rb_cFloat, "angle", float_arg, 0); 01984 rb_define_method(rb_cFloat, "phase", float_arg, 0); 01985 01986 rb_define_const(rb_cComplex, "I", 01987 f_complex_new_bang2(rb_cComplex, ZERO, ONE)); 01988 } 01989 01990 /* 01991 Local variables: 01992 c-file-style: "ruby" 01993 End: 01994 */ 01995