Ruby 1.9.2p136(2010-12-25revision30365)

- Main Page
- Modules
- Data Structures
- Files

- File List
- Globals

# complex.c

Go to the documentation of this file.

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