Ruby 1.9.3p327(2012-11-10revision37606)
complex.c
Go to the documentation of this file.
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