complex.c

   1 /*
   2   complex.c: Coded by Tadayoshi Funaba 2008-2012
   3
   4   This implementation is based on Keiju Ishitsuka's Complex library
   5   which is written in ruby.
   6 */
   7
   8 #include "ruby/internal/config.h"
   9
  10 #if defined _MSC_VER
  11 /* Microsoft Visual C does not define M_PI and others by default */
  12 # define _USE_MATH_DEFINES 1
  13 #endif
  14
  15 #include <ctype.h>
  16 #include <math.h>
  17
  18 #include "id.h"
  19 #include "internal.h"
  20 #include "internal/array.h"
  21 #include "internal/class.h"
  22 #include "internal/complex.h"
  23 #include "internal/math.h"
  24 #include "internal/numeric.h"
  25 #include "internal/object.h"
  26 #include "internal/rational.h"
  27 #include "ruby_assert.h"
  28
  29 #define ZERO INT2FIX(0)
  30 #define ONE INT2FIX(1)
  31 #define TWO INT2FIX(2)
  32 #if USE_FLONUM
  33 #define RFLOAT_0 DBL2NUM(0)
  34 #else
  35 static VALUE RFLOAT_0;
  36 #endif
  37
  38 VALUE rb_cComplex;
  39
  40 static ID id_abs, id_arg,
  41     id_denominator, id_numerator,
  42     id_real_p, id_i_real, id_i_imag,
  43     id_finite_p, id_infinite_p, id_rationalize,
  44     id_PI;
  45 #define id_to_i idTo_i
  46 #define id_to_r idTo_r
  47 #define id_negate idUMinus
  48 #define id_expt idPow
  49 #define id_to_f idTo_f
  50 #define id_quo idQuo
  51 #define id_fdiv idFdiv
  52
  53 #define fun1(n) \
  54 inline static VALUE \
  55 f_##n(VALUE x)\
  56 {\
  57     return rb_funcall(x, id_##n, 0);\
  58 }
  59
  60 #define fun2(n) \
  61 inline static VALUE \
  62 f_##n(VALUE x, VALUE y)\
  63 {\
  64     return rb_funcall(x, id_##n, 1, y);\
  65 }
  66
  67 #define PRESERVE_SIGNEDZERO
  68
  69 inline static VALUE
  70 f_add(VALUE x, VALUE y)
  71 {
  72     if (RB_INTEGER_TYPE_P(x) &&
  73         LIKELY(rb_method_basic_definition_p(rb_cInteger, idPLUS))) {
  74         if (FIXNUM_ZERO_P(x))
  75             return y;
  76         if (FIXNUM_ZERO_P(y))
  77             return x;
  78         return rb_int_plus(x, y);
  79     }
  80     else if (RB_FLOAT_TYPE_P(x) &&
  81              LIKELY(rb_method_basic_definition_p(rb_cFloat, idPLUS))) {
  82         if (FIXNUM_ZERO_P(y))
  83             return x;
  84         return rb_float_plus(x, y);
  85     }
  86     else if (RB_TYPE_P(x, T_RATIONAL) &&
  87              LIKELY(rb_method_basic_definition_p(rb_cRational, idPLUS))) {
  88         if (FIXNUM_ZERO_P(y))
  89             return x;
  90         return rb_rational_plus(x, y);
  91     }
  92
  93     return rb_funcall(x, '+', 1, y);
  94 }
  95
  96 inline static VALUE
  97 f_div(VALUE x, VALUE y)
  98 {
  99     if (FIXNUM_P(y) && FIX2LONG(y) == 1)
 100         return x;
 101     return rb_funcall(x, '/', 1, y);
 102 }
 103
 104 inline static int
 105 f_gt_p(VALUE x, VALUE y)
 106 {
 107     if (RB_INTEGER_TYPE_P(x)) {
 108         if (FIXNUM_P(x) && FIXNUM_P(y))
 109             return (SIGNED_VALUE)x > (SIGNED_VALUE)y;
 110         return RTEST(rb_int_gt(x, y));
 111     }
 112     else if (RB_FLOAT_TYPE_P(x))
 113         return RTEST(rb_float_gt(x, y));
 114     else if (RB_TYPE_P(x, T_RATIONAL)) {
 115         int const cmp = rb_cmpint(rb_rational_cmp(x, y), x, y);
 116         return cmp > 0;
 117     }
 118     return RTEST(rb_funcall(x, '>', 1, y));
 119 }
 120
 121 inline static VALUE
 122 f_mul(VALUE x, VALUE y)
 123 {
 124     if (RB_INTEGER_TYPE_P(x) &&
 125         LIKELY(rb_method_basic_definition_p(rb_cInteger, idMULT))) {
 126         if (FIXNUM_ZERO_P(y))
 127             return ZERO;
 128         if (FIXNUM_ZERO_P(x) && RB_INTEGER_TYPE_P(y))
 129             return ZERO;
 130         if (x == ONE) return y;
 131         if (y == ONE) return x;
 132         return rb_int_mul(x, y);
 133     }
 134     else if (RB_FLOAT_TYPE_P(x) &&
 135              LIKELY(rb_method_basic_definition_p(rb_cFloat, idMULT))) {
 136         if (y == ONE) return x;
 137         return rb_float_mul(x, y);
 138     }
 139     else if (RB_TYPE_P(x, T_RATIONAL) &&
 140              LIKELY(rb_method_basic_definition_p(rb_cRational, idMULT))) {
 141         if (y == ONE) return x;
 142         return rb_rational_mul(x, y);
 143     }
 144     else if (LIKELY(rb_method_basic_definition_p(CLASS_OF(x), idMULT))) {
 145         if (y == ONE) return x;
 146     }
 147     return rb_funcall(x, '*', 1, y);
 148 }
 149
 150 inline static VALUE
 151 f_sub(VALUE x, VALUE y)
 152 {
 153     if (FIXNUM_ZERO_P(y) &&
 154         LIKELY(rb_method_basic_definition_p(CLASS_OF(x), idMINUS))) {
 155         return x;
 156     }
 157     return rb_funcall(x, '-', 1, y);
 158 }
 159
 160 inline static VALUE
 161 f_abs(VALUE x)
 162 {
 163     if (RB_INTEGER_TYPE_P(x)) {
 164         return rb_int_abs(x);
 165     }
 166     else if (RB_FLOAT_TYPE_P(x)) {
 167         return rb_float_abs(x);
 168     }
 169     else if (RB_TYPE_P(x, T_RATIONAL)) {
 170         return rb_rational_abs(x);
 171     }
 172     else if (RB_TYPE_P(x, T_COMPLEX)) {
 173         return rb_complex_abs(x);
 174     }
 175     return rb_funcall(x, id_abs, 0);
 176 }
 177
 178 static VALUE numeric_arg(VALUE self);
 179 static VALUE float_arg(VALUE self);
 180
 181 inline static VALUE
 182 f_arg(VALUE x)
 183 {
 184     if (RB_INTEGER_TYPE_P(x)) {
 185         return numeric_arg(x);
 186     }
 187     else if (RB_FLOAT_TYPE_P(x)) {
 188         return float_arg(x);
 189     }
 190     else if (RB_TYPE_P(x, T_RATIONAL)) {
 191         return numeric_arg(x);
 192     }
 193     else if (RB_TYPE_P(x, T_COMPLEX)) {
 194         return rb_complex_arg(x);
 195     }
 196     return rb_funcall(x, id_arg, 0);
 197 }
 198
 199 inline static VALUE
 200 f_numerator(VALUE x)
 201 {
 202     if (RB_TYPE_P(x, T_RATIONAL)) {
 203         return RRATIONAL(x)->num;
 204     }
 205     if (RB_FLOAT_TYPE_P(x)) {
 206         return rb_float_numerator(x);
 207     }
 208     return x;
 209 }
 210
 211 inline static VALUE
 212 f_denominator(VALUE x)
 213 {
 214     if (RB_TYPE_P(x, T_RATIONAL)) {
 215         return RRATIONAL(x)->den;
 216     }
 217     if (RB_FLOAT_TYPE_P(x)) {
 218         return rb_float_denominator(x);
 219     }
 220     return INT2FIX(1);
 221 }
 222
 223 inline static VALUE
 224 f_negate(VALUE x)
 225 {
 226     if (RB_INTEGER_TYPE_P(x)) {
 227         return rb_int_uminus(x);
 228     }
 229     else if (RB_FLOAT_TYPE_P(x)) {
 230         return rb_float_uminus(x);
 231     }
 232     else if (RB_TYPE_P(x, T_RATIONAL)) {
 233         return rb_rational_uminus(x);
 234     }
 235     else if (RB_TYPE_P(x, T_COMPLEX)) {
 236         return rb_complex_uminus(x);
 237     }
 238     return rb_funcall(x, id_negate, 0);
 239 }
 240
 241 static bool nucomp_real_p(VALUE self);
 242
 243 static inline bool
 244 f_real_p(VALUE x)
 245 {
 246     if (RB_INTEGER_TYPE_P(x)) {
 247         return true;
 248     }
 249     else if (RB_FLOAT_TYPE_P(x)) {
 250         return true;
 251     }
 252     else if (RB_TYPE_P(x, T_RATIONAL)) {
 253         return true;
 254     }
 255     else if (RB_TYPE_P(x, T_COMPLEX)) {
 256         return nucomp_real_p(x);
 257     }
 258     return rb_funcall(x, id_real_p, 0);
 259 }
 260
 261 inline static VALUE
 262 f_to_i(VALUE x)
 263 {
 264     if (RB_TYPE_P(x, T_STRING))
 265         return rb_str_to_inum(x, 10, 0);
 266     return rb_funcall(x, id_to_i, 0);
 267 }
 268
 269 inline static VALUE
 270 f_to_f(VALUE x)
 271 {
 272     if (RB_TYPE_P(x, T_STRING))
 273         return DBL2NUM(rb_str_to_dbl(x, 0));
 274     return rb_funcall(x, id_to_f, 0);
 275 }
 276
 277 fun1(to_r)
 278
 279 inline static int
 280 f_eqeq_p(VALUE x, VALUE y)
 281 {
 282     if (FIXNUM_P(x) && FIXNUM_P(y))
 283         return x == y;
 284     else if (RB_FLOAT_TYPE_P(x) || RB_FLOAT_TYPE_P(y))
 285         return NUM2DBL(x) == NUM2DBL(y);
 286     return (int)rb_equal(x, y);
 287 }
 288
 289 fun2(expt)
 290 fun2(fdiv)
 291
 292 static VALUE
 293 f_quo(VALUE x, VALUE y)
 294 {
 295     if (RB_INTEGER_TYPE_P(x))
 296         return rb_numeric_quo(x, y);
 297     if (RB_FLOAT_TYPE_P(x))
 298         return rb_float_div(x, y);
 299     if (RB_TYPE_P(x, T_RATIONAL))
 300         return rb_numeric_quo(x, y);
 301
 302     return rb_funcallv(x, id_quo, 1, &y);
 303 }
 304
 305 inline static int
 306 f_negative_p(VALUE x)
 307 {
 308     if (RB_INTEGER_TYPE_P(x))
 309         return INT_NEGATIVE_P(x);
 310     else if (RB_FLOAT_TYPE_P(x))
 311         return RFLOAT_VALUE(x) < 0.0;
 312     else if (RB_TYPE_P(x, T_RATIONAL))
 313         return INT_NEGATIVE_P(RRATIONAL(x)->num);
 314     return rb_num_negative_p(x);
 315 }
 316
 317 #define f_positive_p(x) (!f_negative_p(x))
 318
 319 inline static int
 320 f_zero_p(VALUE x)
 321 {
 322     if (RB_FLOAT_TYPE_P(x)) {
 323         return FLOAT_ZERO_P(x);
 324     }
 325     else if (RB_INTEGER_TYPE_P(x)) {
 326         return FIXNUM_ZERO_P(x);
 327     }
 328     else if (RB_TYPE_P(x, T_RATIONAL)) {
 329         const VALUE num = RRATIONAL(x)->num;
 330         return FIXNUM_ZERO_P(num);
 331     }
 332     return (int)rb_equal(x, ZERO);
 333 }
 334
 335 #define f_nonzero_p(x) (!f_zero_p(x))
 336
 337 static inline bool
 338 always_finite_type_p(VALUE x)
 339 {
 340     if (FIXNUM_P(x)) return true;
 341     if (FLONUM_P(x)) return true; /* Infinity can't be a flonum */
 342     return (RB_INTEGER_TYPE_P(x) || RB_TYPE_P(x, T_RATIONAL));
 343 }
 344
 345 inline static int
 346 f_finite_p(VALUE x)
 347 {
 348     if (always_finite_type_p(x)) {
 349         return TRUE;
 350     }
 351     else if (RB_FLOAT_TYPE_P(x)) {
 352         return isfinite(RFLOAT_VALUE(x));
 353     }
 354     return RTEST(rb_funcallv(x, id_finite_p, 0, 0));
 355 }
 356
 357 inline static int
 358 f_infinite_p(VALUE x)
 359 {
 360     if (always_finite_type_p(x)) {
 361         return FALSE;
 362     }
 363     else if (RB_FLOAT_TYPE_P(x)) {
 364         return isinf(RFLOAT_VALUE(x));
 365     }
 366     return RTEST(rb_funcallv(x, id_infinite_p, 0, 0));
 367 }
 368
 369 inline static int
 370 f_kind_of_p(VALUE x, VALUE c)
 371 {
 372     return (int)rb_obj_is_kind_of(x, c);
 373 }
 374
 375 inline static int
 376 k_numeric_p(VALUE x)
 377 {
 378     return f_kind_of_p(x, rb_cNumeric);
 379 }
 380
 381 #define k_exact_p(x) (!RB_FLOAT_TYPE_P(x))
 382
 383 #define k_exact_zero_p(x) (k_exact_p(x) && f_zero_p(x))
 384
 385 #define get_dat1(x) \
 386     struct RComplex *dat = RCOMPLEX(x)
 387
 388 #define get_dat2(x,y) \
 389     struct RComplex *adat = RCOMPLEX(x), *bdat = RCOMPLEX(y)
 390
 391 inline static VALUE
 392 nucomp_s_new_internal(VALUE klass, VALUE real, VALUE imag)
 393 {
 394     NEWOBJ_OF(obj, struct RComplex, klass, T_COMPLEX | (RGENGC_WB_PROTECTED_COMPLEX ? FL_WB_PROTECTED : 0));
 395
 396     RCOMPLEX_SET_REAL(obj, real);
 397     RCOMPLEX_SET_IMAG(obj, imag);
 398     OBJ_FREEZE_RAW((VALUE)obj);
 399
 400     return (VALUE)obj;
 401 }
 402
 403 static VALUE
 404 nucomp_s_alloc(VALUE klass)
 405 {
 406     return nucomp_s_new_internal(klass, ZERO, ZERO);
 407 }
 408
 409 inline static VALUE
 410 f_complex_new_bang1(VALUE klass, VALUE x)
 411 {
 412     assert(!RB_TYPE_P(x, T_COMPLEX));
 413     return nucomp_s_new_internal(klass, x, ZERO);
 414 }
 415
 416 inline static VALUE
 417 f_complex_new_bang2(VALUE klass, VALUE x, VALUE y)
 418 {
 419     assert(!RB_TYPE_P(x, T_COMPLEX));
 420     assert(!RB_TYPE_P(y, T_COMPLEX));
 421     return nucomp_s_new_internal(klass, x, y);
 422 }
 423
 424 inline static void
 425 nucomp_real_check(VALUE num)
 426 {
 427     if (!RB_INTEGER_TYPE_P(num) &&
 428         !RB_FLOAT_TYPE_P(num) &&
 429         !RB_TYPE_P(num, T_RATIONAL)) {
 430         if (!k_numeric_p(num) || !f_real_p(num))
 431             rb_raise(rb_eTypeError, "not a real");
 432     }
 433 }
 434
 435 inline static VALUE
 436 nucomp_s_canonicalize_internal(VALUE klass, VALUE real, VALUE imag)
 437 {
 438     int complex_r, complex_i;
 439     complex_r = RB_TYPE_P(real, T_COMPLEX);
 440     complex_i = RB_TYPE_P(imag, T_COMPLEX);
 441     if (!complex_r && !complex_i) {
 442         return nucomp_s_new_internal(klass, real, imag);
 443     }
 444     else if (!complex_r) {
 445         get_dat1(imag);
 446
 447         return nucomp_s_new_internal(klass,
 448                                      f_sub(real, dat->imag),
 449                                      f_add(ZERO, dat->real));
 450     }
 451     else if (!complex_i) {
 452         get_dat1(real);
 453
 454         return nucomp_s_new_internal(klass,
 455                                      dat->real,
 456                                      f_add(dat->imag, imag));
 457     }
 458     else {
 459         get_dat2(real, imag);
 460
 461         return nucomp_s_new_internal(klass,
 462                                      f_sub(adat->real, bdat->imag),
 463                                      f_add(adat->imag, bdat->real));
 464     }
 465 }
 466
 467 /*
 468  * call-seq:
 469  *    Complex.rect(real[, imag])         ->  complex
 470  *    Complex.rectangular(real[, imag])  ->  complex
 471  *
 472  * Returns a complex object which denotes the given rectangular form.
 473  *
 474  *    Complex.rectangular(1, 2)  #=> (1+2i)
 475  */
 476 static VALUE
 477 nucomp_s_new(int argc, VALUE *argv, VALUE klass)
 478 {
 479     VALUE real, imag;
 480
 481     switch (rb_scan_args(argc, argv, "11", &real, &imag)) {
 482       case 1:
 483         nucomp_real_check(real);
 484         imag = ZERO;
 485         break;
 486       default:
 487         nucomp_real_check(real);
 488         nucomp_real_check(imag);
 489         break;
 490     }
 491
 492     return nucomp_s_canonicalize_internal(klass, real, imag);
 493 }
 494
 495 inline static VALUE
 496 f_complex_new2(VALUE klass, VALUE x, VALUE y)
 497 {
 498     assert(!RB_TYPE_P(x, T_COMPLEX));
 499     return nucomp_s_canonicalize_internal(klass, x, y);
 500 }
 501
 502 static VALUE nucomp_convert(VALUE klass, VALUE a1, VALUE a2, int raise);
 503 static VALUE nucomp_s_convert(int argc, VALUE *argv, VALUE klass);
 504
 505 /*
 506  * call-seq:
 507  *    Complex(x[, y], exception: true)  ->  numeric or nil
 508  *
 509  * Returns x+i*y;
 510  *
 511  *    Complex(1, 2)    #=> (1+2i)
 512  *    Complex('1+2i')  #=> (1+2i)
 513  *    Complex(nil)     #=> TypeError
 514  *    Complex(1, nil)  #=> TypeError
 515  *
 516  *    Complex(1, nil, exception: false)  #=> nil
 517  *    Complex('1+2', exception: false)   #=> nil
 518  *
 519  * Syntax of string form:
 520  *
 521  *   string form = extra spaces , complex , extra spaces ;
 522  *   complex = real part | [ sign ] , imaginary part
 523  *           | real part , sign , imaginary part
 524  *           | rational , "@" , rational ;
 525  *   real part = rational ;
 526  *   imaginary part = imaginary unit | unsigned rational , imaginary unit ;
 527  *   rational = [ sign ] , unsigned rational ;
 528  *   unsigned rational = numerator | numerator , "/" , denominator ;
 529  *   numerator = integer part | fractional part | integer part , fractional part ;
 530  *   denominator = digits ;
 531  *   integer part = digits ;
 532  *   fractional part = "." , digits , [ ( "e" | "E" ) , [ sign ] , digits ] ;
 533  *   imaginary unit = "i" | "I" | "j" | "J" ;
 534  *   sign = "-" | "+" ;
 535  *   digits = digit , { digit | "_" , digit };
 536  *   digit = "0" | "1" | "2" | "3" | "4" | "5" | "6" | "7" | "8" | "9" ;
 537  *   extra spaces = ? \s* ? ;
 538  *
 539  * See String#to_c.
 540  */
 541 static VALUE
 542 nucomp_f_complex(int argc, VALUE *argv, VALUE klass)
 543 {
 544     VALUE a1, a2, opts = Qnil;
 545     int raise = TRUE;
 546
 547     if (rb_scan_args(argc, argv, "11:", &a1, &a2, &opts) == 1) {
 548         a2 = Qundef;
 549     }
 550     if (!NIL_P(opts)) {
 551         raise = rb_opts_exception_p(opts, raise);
 552     }
 553     if (argc > 0 && CLASS_OF(a1) == rb_cComplex && a2 == Qundef) {
 554         return a1;
 555     }
 556     return nucomp_convert(rb_cComplex, a1, a2, raise);
 557 }
 558
 559 #define imp1(n) \
 560 inline static VALUE \
 561 m_##n##_bang(VALUE x)\
 562 {\
 563     return rb_math_##n(x);\
 564 }
 565
 566 imp1(cos)
 567 imp1(cosh)
 568 imp1(exp)
 569
 570 static VALUE
 571 m_log_bang(VALUE x)
 572 {
 573     return rb_math_log(1, &x);
 574 }
 575
 576 imp1(sin)
 577 imp1(sinh)
 578
 579 static VALUE
 580 m_cos(VALUE x)
 581 {
 582     if (!RB_TYPE_P(x, T_COMPLEX))
 583         return m_cos_bang(x);
 584     {
 585         get_dat1(x);
 586         return f_complex_new2(rb_cComplex,
 587                               f_mul(m_cos_bang(dat->real),
 588                                     m_cosh_bang(dat->imag)),
 589                               f_mul(f_negate(m_sin_bang(dat->real)),
 590                                     m_sinh_bang(dat->imag)));
 591     }
 592 }
 593
 594 static VALUE
 595 m_sin(VALUE x)
 596 {
 597     if (!RB_TYPE_P(x, T_COMPLEX))
 598         return m_sin_bang(x);
 599     {
 600         get_dat1(x);
 601         return f_complex_new2(rb_cComplex,
 602                               f_mul(m_sin_bang(dat->real),
 603                                     m_cosh_bang(dat->imag)),
 604                               f_mul(m_cos_bang(dat->real),
 605                                     m_sinh_bang(dat->imag)));
 606     }
 607 }
 608
 609 static VALUE
 610 f_complex_polar(VALUE klass, VALUE x, VALUE y)
 611 {
 612     assert(!RB_TYPE_P(x, T_COMPLEX));
 613     assert(!RB_TYPE_P(y, T_COMPLEX));
 614     if (f_zero_p(x) || f_zero_p(y)) {
 615         return nucomp_s_new_internal(klass, x, RFLOAT_0);
 616     }
 617     if (RB_FLOAT_TYPE_P(y)) {
 618         const double arg = RFLOAT_VALUE(y);
 619         if (arg == M_PI) {
 620             x = f_negate(x);
 621             y = RFLOAT_0;
 622         }
 623         else if (arg == M_PI_2) {
 624             y = x;
 625             x = RFLOAT_0;
 626         }
 627         else if (arg == M_PI_2+M_PI) {
 628             y = f_negate(x);
 629             x = RFLOAT_0;
 630         }
 631         else if (RB_FLOAT_TYPE_P(x)) {
 632             const double abs = RFLOAT_VALUE(x);
 633             const double real = abs * cos(arg), imag = abs * sin(arg);
 634             x = DBL2NUM(real);
 635             y = DBL2NUM(imag);
 636         }
 637         else {
 638             const double ax = sin(arg), ay = cos(arg);
 639             y = f_mul(x, DBL2NUM(ax));
 640             x = f_mul(x, DBL2NUM(ay));
 641         }
 642         return nucomp_s_new_internal(klass, x, y);
 643     }
 644     return nucomp_s_canonicalize_internal(klass,
 645                                           f_mul(x, m_cos(y)),
 646                                           f_mul(x, m_sin(y)));
 647 }
 648
 649 #ifdef HAVE___COSPI
 650 # define cospi(x) __cospi(x)
 651 #else
 652 # define cospi(x) cos((x) * M_PI)
 653 #endif
 654 #ifdef HAVE___SINPI
 655 # define sinpi(x) __sinpi(x)
 656 #else
 657 # define sinpi(x) sin((x) * M_PI)
 658 #endif
 659 /* returns a Complex or Float of ang*PI-rotated abs */
 660 VALUE
 661 rb_dbl_complex_new_polar_pi(double abs, double ang)
 662 {
 663     double fi;
 664     const double fr = modf(ang, &fi);
 665     int pos = fr == +0.5;
 666
 667     if (pos || fr == -0.5) {
 668         if ((modf(fi / 2.0, &fi) != fr) ^ pos) abs = -abs;
 669         return rb_complex_new(RFLOAT_0, DBL2NUM(abs));
 670     }
 671     else if (fr == 0.0) {
 672         if (modf(fi / 2.0, &fi) != 0.0) abs = -abs;
 673         return DBL2NUM(abs);
 674     }
 675     else {
 676         const double real = abs * cospi(ang), imag = abs * sinpi(ang);
 677         return rb_complex_new(DBL2NUM(real), DBL2NUM(imag));
 678     }
 679 }
 680
 681 /*
 682  * call-seq:
 683  *    Complex.polar(abs[, arg])  ->  complex
 684  *
 685  * Returns a complex object which denotes the given polar form.
 686  *
 687  *    Complex.polar(3, 0)            #=> (3.0+0.0i)
 688  *    Complex.polar(3, Math::PI/2)   #=> (1.836909530733566e-16+3.0i)
 689  *    Complex.polar(3, Math::PI)     #=> (-3.0+3.673819061467132e-16i)
 690  *    Complex.polar(3, -Math::PI/2)  #=> (1.836909530733566e-16-3.0i)
 691  */
 692 static VALUE
 693 nucomp_s_polar(int argc, VALUE *argv, VALUE klass)
 694 {
 695     VALUE abs, arg;
 696
 697     switch (rb_scan_args(argc, argv, "11", &abs, &arg)) {
 698       case 1:
 699         nucomp_real_check(abs);
 700         return nucomp_s_new_internal(klass, abs, ZERO);
 701       default:
 702         nucomp_real_check(abs);
 703         nucomp_real_check(arg);
 704         break;
 705     }
 706     if (RB_TYPE_P(abs, T_COMPLEX)) {
 707         get_dat1(abs);
 708         abs = dat->real;
 709     }
 710     if (RB_TYPE_P(arg, T_COMPLEX)) {
 711         get_dat1(arg);
 712         arg = dat->real;
 713     }
 714     return f_complex_polar(klass, abs, arg);
 715 }
 716
 717 /*
 718  * call-seq:
 719  *    cmp.real  ->  real
 720  *
 721  * Returns the real part.
 722  *
 723  *    Complex(7).real      #=> 7
 724  *    Complex(9, -4).real  #=> 9
 725  */
 726 VALUE
 727 rb_complex_real(VALUE self)
 728 {
 729     get_dat1(self);
 730     return dat->real;
 731 }
 732
 733 /*
 734  * call-seq:
 735  *    cmp.imag       ->  real
 736  *    cmp.imaginary  ->  real
 737  *
 738  * Returns the imaginary part.
 739  *
 740  *    Complex(7).imaginary      #=> 0
 741  *    Complex(9, -4).imaginary  #=> -4
 742  */
 743 VALUE
 744 rb_complex_imag(VALUE self)
 745 {
 746     get_dat1(self);
 747     return dat->imag;
 748 }
 749
 750 /*
 751  * call-seq:
 752  *    -cmp  ->  complex
 753  *
 754  * Returns negation of the value.
 755  *
 756  *    -Complex(1, 2)  #=> (-1-2i)
 757  */
 758 VALUE
 759 rb_complex_uminus(VALUE self)
 760 {
 761     get_dat1(self);
 762     return f_complex_new2(CLASS_OF(self),
 763                           f_negate(dat->real), f_negate(dat->imag));
 764 }
 765
 766 /*
 767  * call-seq:
 768  *    cmp + numeric  ->  complex
 769  *
 770  * Performs addition.
 771  *
 772  *    Complex(2, 3)  + Complex(2, 3)   #=> (4+6i)
 773  *    Complex(900)   + Complex(1)      #=> (901+0i)
 774  *    Complex(-2, 9) + Complex(-9, 2)  #=> (-11+11i)
 775  *    Complex(9, 8)  + 4               #=> (13+8i)
 776  *    Complex(20, 9) + 9.8             #=> (29.8+9i)
 777  */
 778 VALUE
 779 rb_complex_plus(VALUE self, VALUE other)
 780 {
 781     if (RB_TYPE_P(other, T_COMPLEX)) {
 782         VALUE real, imag;
 783
 784         get_dat2(self, other);
 785
 786         real = f_add(adat->real, bdat->real);
 787         imag = f_add(adat->imag, bdat->imag);
 788
 789         return f_complex_new2(CLASS_OF(self), real, imag);
 790     }
 791     if (k_numeric_p(other) && f_real_p(other)) {
 792         get_dat1(self);
 793
 794         return f_complex_new2(CLASS_OF(self),
 795                               f_add(dat->real, other), dat->imag);
 796     }
 797     return rb_num_coerce_bin(self, other, '+');
 798 }
 799
 800 /*
 801  * call-seq:
 802  *    cmp - numeric  ->  complex
 803  *
 804  * Performs subtraction.
 805  *
 806  *    Complex(2, 3)  - Complex(2, 3)   #=> (0+0i)
 807  *    Complex(900)   - Complex(1)      #=> (899+0i)
 808  *    Complex(-2, 9) - Complex(-9, 2)  #=> (7+7i)
 809  *    Complex(9, 8)  - 4               #=> (5+8i)
 810  *    Complex(20, 9) - 9.8             #=> (10.2+9i)
 811  */
 812 VALUE
 813 rb_complex_minus(VALUE self, VALUE other)
 814 {
 815     if (RB_TYPE_P(other, T_COMPLEX)) {
 816         VALUE real, imag;
 817
 818         get_dat2(self, other);
 819
 820         real = f_sub(adat->real, bdat->real);
 821         imag = f_sub(adat->imag, bdat->imag);
 822
 823         return f_complex_new2(CLASS_OF(self), real, imag);
 824     }
 825     if (k_numeric_p(other) && f_real_p(other)) {
 826         get_dat1(self);
 827
 828         return f_complex_new2(CLASS_OF(self),
 829                               f_sub(dat->real, other), dat->imag);
 830     }
 831     return rb_num_coerce_bin(self, other, '-');
 832 }
 833
 834 static VALUE
 835 safe_mul(VALUE a, VALUE b, int az, int bz)
 836 {
 837     double v;
 838     if (!az && bz && RB_FLOAT_TYPE_P(a) && (v = RFLOAT_VALUE(a), !isnan(v))) {
 839         a = signbit(v) ? DBL2NUM(-1.0) : DBL2NUM(1.0);
 840     }
 841     if (!bz && az && RB_FLOAT_TYPE_P(b) && (v = RFLOAT_VALUE(b), !isnan(v))) {
 842         b = signbit(v) ? DBL2NUM(-1.0) : DBL2NUM(1.0);
 843     }
 844     return f_mul(a, b);
 845 }
 846
 847 static void
 848 comp_mul(VALUE areal, VALUE aimag, VALUE breal, VALUE bimag, VALUE *real, VALUE *imag)
 849 {
 850     int arzero = f_zero_p(areal);
 851     int aizero = f_zero_p(aimag);
 852     int brzero = f_zero_p(breal);
 853     int bizero = f_zero_p(bimag);
 854     *real = f_sub(safe_mul(areal, breal, arzero, brzero),
 855                   safe_mul(aimag, bimag, aizero, bizero));
 856     *imag = f_add(safe_mul(areal, bimag, arzero, bizero),
 857                   safe_mul(aimag, breal, aizero, brzero));
 858 }
 859
 860 /*
 861  * call-seq:
 862  *    cmp * numeric  ->  complex
 863  *
 864  * Performs multiplication.
 865  *
 866  *    Complex(2, 3)  * Complex(2, 3)   #=> (-5+12i)
 867  *    Complex(900)   * Complex(1)      #=> (900+0i)
 868  *    Complex(-2, 9) * Complex(-9, 2)  #=> (0-85i)
 869  *    Complex(9, 8)  * 4               #=> (36+32i)
 870  *    Complex(20, 9) * 9.8             #=> (196.0+88.2i)
 871  */
 872 VALUE
 873 rb_complex_mul(VALUE self, VALUE other)
 874 {
 875     if (RB_TYPE_P(other, T_COMPLEX)) {
 876         VALUE real, imag;
 877         get_dat2(self, other);
 878
 879         comp_mul(adat->real, adat->imag, bdat->real, bdat->imag, &real, &imag);
 880
 881         return f_complex_new2(CLASS_OF(self), real, imag);
 882     }
 883     if (k_numeric_p(other) && f_real_p(other)) {
 884         get_dat1(self);
 885
 886         return f_complex_new2(CLASS_OF(self),
 887                               f_mul(dat->real, other),
 888                               f_mul(dat->imag, other));
 889     }
 890     return rb_num_coerce_bin(self, other, '*');
 891 }
 892
 893 inline static VALUE
 894 f_divide(VALUE self, VALUE other,
 895          VALUE (*func)(VALUE, VALUE), ID id)
 896 {
 897     if (RB_TYPE_P(other, T_COMPLEX)) {
 898         VALUE r, n, x, y;
 899         int flo;
 900         get_dat2(self, other);
 901
 902         flo = (RB_FLOAT_TYPE_P(adat->real) || RB_FLOAT_TYPE_P(adat->imag) ||
 903                RB_FLOAT_TYPE_P(bdat->real) || RB_FLOAT_TYPE_P(bdat->imag));
 904
 905         if (f_gt_p(f_abs(bdat->real), f_abs(bdat->imag))) {
 906             r = (*func)(bdat->imag, bdat->real);
 907             n = f_mul(bdat->real, f_add(ONE, f_mul(r, r)));
 908             x = (*func)(f_add(adat->real, f_mul(adat->imag, r)), n);
 909             y = (*func)(f_sub(adat->imag, f_mul(adat->real, r)), n);
 910         }
 911         else {
 912             r = (*func)(bdat->real, bdat->imag);
 913             n = f_mul(bdat->imag, f_add(ONE, f_mul(r, r)));
 914             x = (*func)(f_add(f_mul(adat->real, r), adat->imag), n);
 915             y = (*func)(f_sub(f_mul(adat->imag, r), adat->real), n);
 916         }
 917         if (!flo) {
 918             x = rb_rational_canonicalize(x);
 919             y = rb_rational_canonicalize(y);
 920         }
 921         return f_complex_new2(CLASS_OF(self), x, y);
 922     }
 923     if (k_numeric_p(other) && f_real_p(other)) {
 924         VALUE x, y;
 925         get_dat1(self);
 926         x = rb_rational_canonicalize((*func)(dat->real, other));
 927         y = rb_rational_canonicalize((*func)(dat->imag, other));
 928         return f_complex_new2(CLASS_OF(self), x, y);
 929     }
 930     return rb_num_coerce_bin(self, other, id);
 931 }
 932
 933 #define rb_raise_zerodiv() rb_raise(rb_eZeroDivError, "divided by 0")
 934
 935 /*
 936  * call-seq:
 937  *    cmp / numeric     ->  complex
 938  *    cmp.quo(numeric)  ->  complex
 939  *
 940  * Performs division.
 941  *
 942  *    Complex(2, 3)  / Complex(2, 3)   #=> ((1/1)+(0/1)*i)
 943  *    Complex(900)   / Complex(1)      #=> ((900/1)+(0/1)*i)
 944  *    Complex(-2, 9) / Complex(-9, 2)  #=> ((36/85)-(77/85)*i)
 945  *    Complex(9, 8)  / 4               #=> ((9/4)+(2/1)*i)
 946  *    Complex(20, 9) / 9.8             #=> (2.0408163265306123+0.9183673469387754i)
 947  */
 948 VALUE
 949 rb_complex_div(VALUE self, VALUE other)
 950 {
 951     return f_divide(self, other, f_quo, id_quo);
 952 }
 953
 954 #define nucomp_quo rb_complex_div
 955
 956 /*
 957  * call-seq:
 958  *    cmp.fdiv(numeric)  ->  complex
 959  *
 960  * Performs division as each part is a float, never returns a float.
 961  *
 962  *    Complex(11, 22).fdiv(3)  #=> (3.6666666666666665+7.333333333333333i)
 963  */
 964 static VALUE
 965 nucomp_fdiv(VALUE self, VALUE other)
 966 {
 967     return f_divide(self, other, f_fdiv, id_fdiv);
 968 }
 969
 970 inline static VALUE
 971 f_reciprocal(VALUE x)
 972 {
 973     return f_quo(ONE, x);
 974 }
 975
 976 /*
 977  * call-seq:
 978  *    cmp ** numeric  ->  complex
 979  *
 980  * Performs exponentiation.
 981  *
 982  *    Complex('i') ** 2              #=> (-1+0i)
 983  *    Complex(-8) ** Rational(1, 3)  #=> (1.0000000000000002+1.7320508075688772i)
 984  */
 985 VALUE
 986 rb_complex_pow(VALUE self, VALUE other)
 987 {
 988     if (k_numeric_p(other) && k_exact_zero_p(other))
 989         return f_complex_new_bang1(CLASS_OF(self), ONE);
 990
 991     if (RB_TYPE_P(other, T_RATIONAL) && RRATIONAL(other)->den == LONG2FIX(1))
 992         other = RRATIONAL(other)->num; /* c14n */
 993
 994     if (RB_TYPE_P(other, T_COMPLEX)) {
 995         get_dat1(other);
 996
 997         if (k_exact_zero_p(dat->imag))
 998             other = dat->real; /* c14n */
 999     }
1000
1001     if (RB_TYPE_P(other, T_COMPLEX)) {
1002         VALUE r, theta, nr, ntheta;
1003
1004         get_dat1(other);
1005
1006         r = f_abs(self);
1007         theta = f_arg(self);
1008
1009         nr = m_exp_bang(f_sub(f_mul(dat->real, m_log_bang(r)),
1010                               f_mul(dat->imag, theta)));
1011         ntheta = f_add(f_mul(theta, dat->real),
1012                        f_mul(dat->imag, m_log_bang(r)));
1013         return f_complex_polar(CLASS_OF(self), nr, ntheta);
1014     }
1015     if (FIXNUM_P(other)) {
1016         long n = FIX2LONG(other);
1017         if (n == 0) {
1018             return nucomp_s_new_internal(CLASS_OF(self), ONE, ZERO);
1019         }
1020         if (n < 0) {
1021             self = f_reciprocal(self);
1022             other = rb_int_uminus(other);
1023             n = -n;
1024         }
1025         {
1026             get_dat1(self);
1027             VALUE xr = dat->real, xi = dat->imag, zr = xr, zi = xi;
1028
1029             if (f_zero_p(xi)) {
1030                 zr = rb_num_pow(zr, other);
1031             }
1032             else if (f_zero_p(xr)) {
1033                 zi = rb_num_pow(zi, other);
1034                 if (n & 2) zi = f_negate(zi);
1035                 if (!(n & 1)) {
1036                     VALUE tmp = zr;
1037                     zr = zi;
1038                     zi = tmp;
1039                 }
1040             }
1041             else {
1042                 while (--n) {
1043                     long q, r;
1044
1045                     for (; q = n / 2, r = n % 2, r == 0; n = q) {
1046                         VALUE tmp = f_sub(f_mul(xr, xr), f_mul(xi, xi));
1047                         xi = f_mul(f_mul(TWO, xr), xi);
1048                         xr = tmp;
1049                     }
1050                     comp_mul(zr, zi, xr, xi, &zr, &zi);
1051                 }
1052             }
1053             return nucomp_s_new_internal(CLASS_OF(self), zr, zi);
1054         }
1055     }
1056     if (k_numeric_p(other) && f_real_p(other)) {
1057         VALUE r, theta;
1058
1059         if (RB_BIGNUM_TYPE_P(other))
1060             rb_warn("in a**b, b may be too big");
1061
1062         r = f_abs(self);
1063         theta = f_arg(self);
1064
1065         return f_complex_polar(CLASS_OF(self), f_expt(r, other),
1066                                f_mul(theta, other));
1067     }
1068     return rb_num_coerce_bin(self, other, id_expt);
1069 }
1070
1071 /*
1072  * call-seq:
1073  *    cmp == object  ->  true or false
1074  *
1075  * Returns true if cmp equals object numerically.
1076  *
1077  *    Complex(2, 3)  == Complex(2, 3)   #=> true
1078  *    Complex(5)     == 5               #=> true
1079  *    Complex(0)     == 0.0             #=> true
1080  *    Complex('1/3') == 0.33            #=> false
1081  *    Complex('1/2') == '1/2'           #=> false
1082  */
1083 static VALUE
1084 nucomp_eqeq_p(VALUE self, VALUE other)
1085 {
1086     if (RB_TYPE_P(other, T_COMPLEX)) {
1087         get_dat2(self, other);
1088
1089         return RBOOL(f_eqeq_p(adat->real, bdat->real) &&
1090                           f_eqeq_p(adat->imag, bdat->imag));
1091     }
1092     if (k_numeric_p(other) && f_real_p(other)) {
1093         get_dat1(self);
1094
1095         return RBOOL(f_eqeq_p(dat->real, other) && f_zero_p(dat->imag));
1096     }
1097     return RBOOL(f_eqeq_p(other, self));
1098 }
1099
1100 static bool
1101 nucomp_real_p(VALUE self)
1102 {
1103     get_dat1(self);
1104     return(f_zero_p(dat->imag) ? true : false);
1105 }
1106
1107 /*
1108  * call-seq:
1109  *    cmp <=> object  ->  0, 1, -1, or nil
1110  *
1111  * If +cmp+'s imaginary part is zero, and +object+ is also a
1112  * real number (or a Complex number where the imaginary part is zero),
1113  * compare the real part of +cmp+ to object.  Otherwise, return nil.
1114  *
1115  *    Complex(2, 3)  <=> Complex(2, 3)   #=> nil
1116  *    Complex(2, 3)  <=> 1               #=> nil
1117  *    Complex(2)     <=> 1               #=> 1
1118  *    Complex(2)     <=> 2               #=> 0
1119  *    Complex(2)     <=> 3               #=> -1
1120  */
1121 static VALUE
1122 nucomp_cmp(VALUE self, VALUE other)
1123 {
1124     if (nucomp_real_p(self) && k_numeric_p(other)) {
1125         if (RB_TYPE_P(other, T_COMPLEX) && nucomp_real_p(other)) {
1126             get_dat2(self, other);
1127             return rb_funcall(adat->real, idCmp, 1, bdat->real);
1128         }
1129         else if (f_real_p(other)) {
1130             get_dat1(self);
1131             return rb_funcall(dat->real, idCmp, 1, other);
1132         }
1133     }
1134     return Qnil;
1135 }
1136
1137 /* :nodoc: */
1138 static VALUE
1139 nucomp_coerce(VALUE self, VALUE other)
1140 {
1141     if (RB_TYPE_P(other, T_COMPLEX))
1142         return rb_assoc_new(other, self);
1143     if (k_numeric_p(other) && f_real_p(other))
1144         return rb_assoc_new(f_complex_new_bang1(CLASS_OF(self), other), self);
1145
1146     rb_raise(rb_eTypeError, "%"PRIsVALUE" can't be coerced into %"PRIsVALUE,
1147              rb_obj_class(other), rb_obj_class(self));
1148     return Qnil;
1149 }
1150
1151 /*
1152  * call-seq:
1153  *    cmp.abs        ->  real
1154  *    cmp.magnitude  ->  real
1155  *
1156  * Returns the absolute part of its polar form.
1157  *
1158  *    Complex(-1).abs         #=> 1
1159  *    Complex(3.0, -4.0).abs  #=> 5.0
1160  */
1161 VALUE
1162 rb_complex_abs(VALUE self)
1163 {
1164     get_dat1(self);
1165
1166     if (f_zero_p(dat->real)) {
1167         VALUE a = f_abs(dat->imag);
1168         if (RB_FLOAT_TYPE_P(dat->real) && !RB_FLOAT_TYPE_P(dat->imag))
1169             a = f_to_f(a);
1170         return a;
1171     }
1172     if (f_zero_p(dat->imag)) {
1173         VALUE a = f_abs(dat->real);
1174         if (!RB_FLOAT_TYPE_P(dat->real) && RB_FLOAT_TYPE_P(dat->imag))
1175             a = f_to_f(a);
1176         return a;
1177     }
1178     return rb_math_hypot(dat->real, dat->imag);
1179 }
1180
1181 /*
1182  * call-seq:
1183  *    cmp.abs2  ->  real
1184  *
1185  * Returns square of the absolute value.
1186  *
1187  *    Complex(-1).abs2         #=> 1
1188  *    Complex(3.0, -4.0).abs2  #=> 25.0
1189  */
1190 static VALUE
1191 nucomp_abs2(VALUE self)
1192 {
1193     get_dat1(self);
1194     return f_add(f_mul(dat->real, dat->real),
1195                  f_mul(dat->imag, dat->imag));
1196 }
1197
1198 /*
1199  * call-seq:
1200  *    cmp.arg    ->  float
1201  *    cmp.angle  ->  float
1202  *    cmp.phase  ->  float
1203  *
1204  * Returns the angle part of its polar form.
1205  *
1206  *    Complex.polar(3, Math::PI/2).arg  #=> 1.5707963267948966
1207  */
1208 VALUE
1209 rb_complex_arg(VALUE self)
1210 {
1211     get_dat1(self);
1212     return rb_math_atan2(dat->imag, dat->real);
1213 }
1214
1215 /*
1216  * call-seq:
1217  *    cmp.rect         ->  array
1218  *    cmp.rectangular  ->  array
1219  *
1220  * Returns an array; [cmp.real, cmp.imag].
1221  *
1222  *    Complex(1, 2).rectangular  #=> [1, 2]
1223  */
1224 static VALUE
1225 nucomp_rect(VALUE self)
1226 {
1227     get_dat1(self);
1228     return rb_assoc_new(dat->real, dat->imag);
1229 }
1230
1231 /*
1232  * call-seq:
1233  *    cmp.polar  ->  array
1234  *
1235  * Returns an array; [cmp.abs, cmp.arg].
1236  *
1237  *    Complex(1, 2).polar  #=> [2.23606797749979, 1.1071487177940904]
1238  */
1239 static VALUE
1240 nucomp_polar(VALUE self)
1241 {
1242     return rb_assoc_new(f_abs(self), f_arg(self));
1243 }
1244
1245 /*
1246  * call-seq:
1247  *    cmp.conj       ->  complex
1248  *    cmp.conjugate  ->  complex
1249  *
1250  * Returns the complex conjugate.
1251  *
1252  *    Complex(1, 2).conjugate  #=> (1-2i)
1253  */
1254 VALUE
1255 rb_complex_conjugate(VALUE self)
1256 {
1257     get_dat1(self);
1258     return f_complex_new2(CLASS_OF(self), dat->real, f_negate(dat->imag));
1259 }
1260
1261 /*
1262  * call-seq:
1263  *    Complex(1).real?     ->  false
1264  *    Complex(1, 2).real?  ->  false
1265  *
1266  * Returns false, even if the complex number has no imaginary part.
1267  */
1268 static VALUE
1269 nucomp_real_p_m(VALUE self)
1270 {
1271     return Qfalse;
1272 }
1273
1274 /*
1275  * call-seq:
1276  *    cmp.denominator  ->  integer
1277  *
1278  * Returns the denominator (lcm of both denominator - real and imag).
1279  *
1280  * See numerator.
1281  */
1282 static VALUE
1283 nucomp_denominator(VALUE self)
1284 {
1285     get_dat1(self);
1286     return rb_lcm(f_denominator(dat->real), f_denominator(dat->imag));
1287 }
1288
1289 /*
1290  * call-seq:
1291  *    cmp.numerator  ->  numeric
1292  *
1293  * Returns the numerator.
1294  *
1295  *        1   2       3+4i  <-  numerator
1296  *        - + -i  ->  ----
1297  *        2   3        6    <-  denominator
1298  *
1299  *    c = Complex('1/2+2/3i')  #=> ((1/2)+(2/3)*i)
1300  *    n = c.numerator          #=> (3+4i)
1301  *    d = c.denominator        #=> 6
1302  *    n / d                    #=> ((1/2)+(2/3)*i)
1303  *    Complex(Rational(n.real, d), Rational(n.imag, d))
1304  *                             #=> ((1/2)+(2/3)*i)
1305  * See denominator.
1306  */
1307 static VALUE
1308 nucomp_numerator(VALUE self)
1309 {
1310     VALUE cd;
1311
1312     get_dat1(self);
1313
1314     cd = nucomp_denominator(self);
1315     return f_complex_new2(CLASS_OF(self),
1316                           f_mul(f_numerator(dat->real),
1317                                 f_div(cd, f_denominator(dat->real))),
1318                           f_mul(f_numerator(dat->imag),
1319                                 f_div(cd, f_denominator(dat->imag))));
1320 }
1321
1322 /* :nodoc: */
1323 st_index_t
1324 rb_complex_hash(VALUE self)
1325 {
1326     st_index_t v, h[2];
1327     VALUE n;
1328
1329     get_dat1(self);
1330     n = rb_hash(dat->real);
1331     h[0] = NUM2LONG(n);
1332     n = rb_hash(dat->imag);
1333     h[1] = NUM2LONG(n);
1334     v = rb_memhash(h, sizeof(h));
1335     return v;
1336 }
1337
1338 static VALUE
1339 nucomp_hash(VALUE self)
1340 {
1341     return ST2FIX(rb_complex_hash(self));
1342 }
1343
1344 /* :nodoc: */
1345 static VALUE
1346 nucomp_eql_p(VALUE self, VALUE other)
1347 {
1348     if (RB_TYPE_P(other, T_COMPLEX)) {
1349         get_dat2(self, other);
1350
1351         return RBOOL((CLASS_OF(adat->real) == CLASS_OF(bdat->real)) &&
1352                           (CLASS_OF(adat->imag) == CLASS_OF(bdat->imag)) &&
1353                           f_eqeq_p(self, other));
1354
1355     }
1356     return Qfalse;
1357 }
1358
1359 inline static int
1360 f_signbit(VALUE x)
1361 {
1362     if (RB_FLOAT_TYPE_P(x)) {
1363         double f = RFLOAT_VALUE(x);
1364         return !isnan(f) && signbit(f);
1365     }
1366     return f_negative_p(x);
1367 }
1368
1369 inline static int
1370 f_tpositive_p(VALUE x)
1371 {
1372     return !f_signbit(x);
1373 }
1374
1375 static VALUE
1376 f_format(VALUE self, VALUE (*func)(VALUE))
1377 {
1378     VALUE s;
1379     int impos;
1380
1381     get_dat1(self);
1382
1383     impos = f_tpositive_p(dat->imag);
1384
1385     s = (*func)(dat->real);
1386     rb_str_cat2(s, !impos ? "-" : "+");
1387
1388     rb_str_concat(s, (*func)(f_abs(dat->imag)));
1389     if (!rb_isdigit(RSTRING_PTR(s)[RSTRING_LEN(s) - 1]))
1390         rb_str_cat2(s, "*");
1391     rb_str_cat2(s, "i");
1392
1393     return s;
1394 }
1395
1396 /*
1397  * call-seq:
1398  *    cmp.to_s  ->  string
1399  *
1400  * Returns the value as a string.
1401  *
1402  *    Complex(2).to_s                       #=> "2+0i"
1403  *    Complex('-8/6').to_s                  #=> "-4/3+0i"
1404  *    Complex('1/2i').to_s                  #=> "0+1/2i"
1405  *    Complex(0, Float::INFINITY).to_s      #=> "0+Infinity*i"
1406  *    Complex(Float::NAN, Float::NAN).to_s  #=> "NaN+NaN*i"
1407  */
1408 static VALUE
1409 nucomp_to_s(VALUE self)
1410 {
1411     return f_format(self, rb_String);
1412 }
1413
1414 /*
1415  * call-seq:
1416  *    cmp.inspect  ->  string
1417  *
1418  * Returns the value as a string for inspection.
1419  *
1420  *    Complex(2).inspect                       #=> "(2+0i)"
1421  *    Complex('-8/6').inspect                  #=> "((-4/3)+0i)"
1422  *    Complex('1/2i').inspect                  #=> "(0+(1/2)*i)"
1423  *    Complex(0, Float::INFINITY).inspect      #=> "(0+Infinity*i)"
1424  *    Complex(Float::NAN, Float::NAN).inspect  #=> "(NaN+NaN*i)"
1425  */
1426 static VALUE
1427 nucomp_inspect(VALUE self)
1428 {
1429     VALUE s;
1430
1431     s = rb_usascii_str_new2("(");
1432     rb_str_concat(s, f_format(self, rb_inspect));
1433     rb_str_cat2(s, ")");
1434
1435     return s;
1436 }
1437
1438 #define FINITE_TYPE_P(v) (RB_INTEGER_TYPE_P(v) || RB_TYPE_P(v, T_RATIONAL))
1439
1440 /*
1441  * call-seq:
1442  *    cmp.finite?  ->  true or false
1443  *
1444  * Returns +true+ if +cmp+'s real and imaginary parts are both finite numbers,
1445  * otherwise returns +false+.
1446  */
1447 static VALUE
1448 rb_complex_finite_p(VALUE self)
1449 {
1450     get_dat1(self);
1451
1452     return RBOOL(f_finite_p(dat->real) && f_finite_p(dat->imag));
1453 }
1454
1455 /*
1456  * call-seq:
1457  *    cmp.infinite?  ->  nil or 1
1458  *
1459  * Returns +1+ if +cmp+'s real or imaginary part is an infinite number,
1460  * otherwise returns +nil+.
1461  *
1462  *  For example:
1463  *
1464  *     (1+1i).infinite?                   #=> nil
1465  *     (Float::INFINITY + 1i).infinite?   #=> 1
1466  */
1467 static VALUE
1468 rb_complex_infinite_p(VALUE self)
1469 {
1470     get_dat1(self);
1471
1472     if (!f_infinite_p(dat->real) && !f_infinite_p(dat->imag)) {
1473         return Qnil;
1474     }
1475     return ONE;
1476 }
1477
1478 /* :nodoc: */
1479 static VALUE
1480 nucomp_dumper(VALUE self)
1481 {
1482     return self;
1483 }
1484
1485 /* :nodoc: */
1486 static VALUE
1487 nucomp_loader(VALUE self, VALUE a)
1488 {
1489     get_dat1(self);
1490
1491     RCOMPLEX_SET_REAL(dat, rb_ivar_get(a, id_i_real));
1492     RCOMPLEX_SET_IMAG(dat, rb_ivar_get(a, id_i_imag));
1493     OBJ_FREEZE_RAW(self);
1494
1495     return self;
1496 }
1497
1498 /* :nodoc: */
1499 static VALUE
1500 nucomp_marshal_dump(VALUE self)
1501 {
1502     VALUE a;
1503     get_dat1(self);
1504
1505     a = rb_assoc_new(dat->real, dat->imag);
1506     rb_copy_generic_ivar(a, self);
1507     return a;
1508 }
1509
1510 /* :nodoc: */
1511 static VALUE
1512 nucomp_marshal_load(VALUE self, VALUE a)
1513 {
1514     Check_Type(a, T_ARRAY);
1515     if (RARRAY_LEN(a) != 2)
1516         rb_raise(rb_eArgError, "marshaled complex must have an array whose length is 2 but %ld", RARRAY_LEN(a));
1517     rb_ivar_set(self, id_i_real, RARRAY_AREF(a, 0));
1518     rb_ivar_set(self, id_i_imag, RARRAY_AREF(a, 1));
1519     return self;
1520 }
1521
1522 VALUE
1523 rb_complex_raw(VALUE x, VALUE y)
1524 {
1525     return nucomp_s_new_internal(rb_cComplex, x, y);
1526 }
1527
1528 VALUE
1529 rb_complex_new(VALUE x, VALUE y)
1530 {
1531     return nucomp_s_canonicalize_internal(rb_cComplex, x, y);
1532 }
1533
1534 VALUE
1535 rb_complex_new_polar(VALUE x, VALUE y)
1536 {
1537     return f_complex_polar(rb_cComplex, x, y);
1538 }
1539
1540 VALUE
1541 rb_complex_polar(VALUE x, VALUE y)
1542 {
1543     return rb_complex_new_polar(x, y);
1544 }
1545
1546 VALUE
1547 rb_Complex(VALUE x, VALUE y)
1548 {
1549     VALUE a[2];
1550     a[0] = x;
1551     a[1] = y;
1552     return nucomp_s_convert(2, a, rb_cComplex);
1553 }
1554
1555 VALUE
1556 rb_dbl_complex_new(double real, double imag)
1557 {
1558     return rb_complex_raw(DBL2NUM(real), DBL2NUM(imag));
1559 }
1560
1561 /*
1562  * call-seq:
1563  *    cmp.to_i  ->  integer
1564  *
1565  * Returns the value as an integer if possible (the imaginary part
1566  * should be exactly zero).
1567  *
1568  *    Complex(1, 0).to_i    #=> 1
1569  *    Complex(1, 0.0).to_i  # RangeError
1570  *    Complex(1, 2).to_i    # RangeError
1571  */
1572 static VALUE
1573 nucomp_to_i(VALUE self)
1574 {
1575     get_dat1(self);
1576
1577     if (!k_exact_zero_p(dat->imag)) {
1578         rb_raise(rb_eRangeError, "can't convert %"PRIsVALUE" into Integer",
1579                  self);
1580     }
1581     return f_to_i(dat->real);
1582 }
1583
1584 /*
1585  * call-seq:
1586  *    cmp.to_f  ->  float
1587  *
1588  * Returns the value as a float if possible (the imaginary part should
1589  * be exactly zero).
1590  *
1591  *    Complex(1, 0).to_f    #=> 1.0
1592  *    Complex(1, 0.0).to_f  # RangeError
1593  *    Complex(1, 2).to_f    # RangeError
1594  */
1595 static VALUE
1596 nucomp_to_f(VALUE self)
1597 {
1598     get_dat1(self);
1599
1600     if (!k_exact_zero_p(dat->imag)) {
1601         rb_raise(rb_eRangeError, "can't convert %"PRIsVALUE" into Float",
1602                  self);
1603     }
1604     return f_to_f(dat->real);
1605 }
1606
1607 /*
1608  * call-seq:
1609  *    cmp.to_r  ->  rational
1610  *
1611  * Returns the value as a rational if possible (the imaginary part
1612  * should be exactly zero).
1613  *
1614  *    Complex(1, 0).to_r    #=> (1/1)
1615  *    Complex(1, 0.0).to_r  # RangeError
1616  *    Complex(1, 2).to_r    # RangeError
1617  *
1618  * See rationalize.
1619  */
1620 static VALUE
1621 nucomp_to_r(VALUE self)
1622 {
1623     get_dat1(self);
1624
1625     if (!k_exact_zero_p(dat->imag)) {
1626         rb_raise(rb_eRangeError, "can't convert %"PRIsVALUE" into Rational",
1627                  self);
1628     }
1629     return f_to_r(dat->real);
1630 }
1631
1632 /*
1633  * call-seq:
1634  *    cmp.rationalize([eps])  ->  rational
1635  *
1636  * Returns the value as a rational if possible (the imaginary part
1637  * should be exactly zero).
1638  *
1639  *    Complex(1.0/3, 0).rationalize  #=> (1/3)
1640  *    Complex(1, 0.0).rationalize    # RangeError
1641  *    Complex(1, 2).rationalize      # RangeError
1642  *
1643  * See to_r.
1644  */
1645 static VALUE
1646 nucomp_rationalize(int argc, VALUE *argv, VALUE self)
1647 {
1648     get_dat1(self);
1649
1650     rb_check_arity(argc, 0, 1);
1651
1652     if (!k_exact_zero_p(dat->imag)) {
1653        rb_raise(rb_eRangeError, "can't convert %"PRIsVALUE" into Rational",
1654                 self);
1655     }
1656     return rb_funcallv(dat->real, id_rationalize, argc, argv);
1657 }
1658
1659 /*
1660  * call-seq:
1661  *    complex.to_c  ->  self
1662  *
1663  * Returns self.
1664  *
1665  *    Complex(2).to_c      #=> (2+0i)
1666  *    Complex(-8, 6).to_c  #=> (-8+6i)
1667  */
1668 static VALUE
1669 nucomp_to_c(VALUE self)
1670 {
1671     return self;
1672 }
1673
1674 /*
1675  * call-seq:
1676  *    nil.to_c  ->  (0+0i)
1677  *
1678  * Returns zero as a complex.
1679  */
1680 static VALUE
1681 nilclass_to_c(VALUE self)
1682 {
1683     return rb_complex_new1(INT2FIX(0));
1684 }
1685
1686 /*
1687  * call-seq:
1688  *    num.to_c  ->  complex
1689  *
1690  * Returns the value as a complex.
1691  */
1692 static VALUE
1693 numeric_to_c(VALUE self)
1694 {
1695     return rb_complex_new1(self);
1696 }
1697
1698 inline static int
1699 issign(int c)
1700 {
1701     return (c == '-' || c == '+');
1702 }
1703
1704 static int
1705 read_sign(const char **s,
1706           char **b)
1707 {
1708     int sign = '?';
1709
1710     if (issign(**s)) {
1711         sign = **b = **s;
1712         (*s)++;
1713         (*b)++;
1714     }
1715     return sign;
1716 }
1717
1718 inline static int
1719 isdecimal(int c)
1720 {
1721     return isdigit((unsigned char)c);
1722 }
1723
1724 static int
1725 read_digits(const char **s, int strict,
1726             char **b)
1727 {
1728     int us = 1;
1729
1730     if (!isdecimal(**s))
1731         return 0;
1732
1733     while (isdecimal(**s) || **s == '_') {
1734         if (**s == '_') {
1735             if (strict) {
1736                 if (us)
1737                     return 0;
1738             }
1739             us = 1;
1740         }
1741         else {
1742             **b = **s;
1743             (*b)++;
1744             us = 0;
1745         }
1746         (*s)++;
1747     }
1748     if (us)
1749         do {
1750             (*s)--;
1751         } while (**s == '_');
1752     return 1;
1753 }
1754
1755 inline static int
1756 islettere(int c)
1757 {
1758     return (c == 'e' || c == 'E');
1759 }
1760
1761 static int
1762 read_num(const char **s, int strict,
1763          char **b)
1764 {
1765     if (**s != '.') {
1766         if (!read_digits(s, strict, b))
1767             return 0;
1768     }
1769
1770     if (**s == '.') {
1771         **b = **s;
1772         (*s)++;
1773         (*b)++;
1774         if (!read_digits(s, strict, b)) {
1775             (*b)--;
1776             return 0;
1777         }
1778     }
1779
1780     if (islettere(**s)) {
1781         **b = **s;
1782         (*s)++;
1783         (*b)++;
1784         read_sign(s, b);
1785         if (!read_digits(s, strict, b)) {
1786             (*b)--;
1787             return 0;
1788         }
1789     }
1790     return 1;
1791 }
1792
1793 inline static int
1794 read_den(const char **s, int strict,
1795          char **b)
1796 {
1797     if (!read_digits(s, strict, b))
1798         return 0;
1799     return 1;
1800 }
1801
1802 static int
1803 read_rat_nos(const char **s, int strict,
1804              char **b)
1805 {
1806     if (!read_num(s, strict, b))
1807         return 0;
1808     if (**s == '/') {
1809         **b = **s;
1810         (*s)++;
1811         (*b)++;
1812         if (!read_den(s, strict, b)) {
1813             (*b)--;
1814             return 0;
1815         }
1816     }
1817     return 1;
1818 }
1819
1820 static int
1821 read_rat(const char **s, int strict,
1822          char **b)
1823 {
1824     read_sign(s, b);
1825     if (!read_rat_nos(s, strict, b))
1826         return 0;
1827     return 1;
1828 }
1829
1830 inline static int
1831 isimagunit(int c)
1832 {
1833     return (c == 'i' || c == 'I' ||
1834             c == 'j' || c == 'J');
1835 }
1836
1837 static VALUE
1838 str2num(char *s)
1839 {
1840     if (strchr(s, '/'))
1841         return rb_cstr_to_rat(s, 0);
1842     if (strpbrk(s, ".eE"))
1843         return DBL2NUM(rb_cstr_to_dbl(s, 0));
1844     return rb_cstr_to_inum(s, 10, 0);
1845 }
1846
1847 static int
1848 read_comp(const char **s, int strict,
1849           VALUE *ret, char **b)
1850 {
1851     char *bb;
1852     int sign;
1853     VALUE num, num2;
1854
1855     bb = *b;
1856
1857     sign = read_sign(s, b);
1858
1859     if (isimagunit(**s)) {
1860         (*s)++;
1861         num = INT2FIX((sign == '-') ? -1 : + 1);
1862         *ret = rb_complex_new2(ZERO, num);
1863         return 1; /* e.g. "i" */
1864     }
1865
1866     if (!read_rat_nos(s, strict, b)) {
1867         **b = '\0';
1868         num = str2num(bb);
1869         *ret = rb_complex_new2(num, ZERO);
1870         return 0; /* e.g. "-" */
1871     }
1872     **b = '\0';
1873     num = str2num(bb);
1874
1875     if (isimagunit(**s)) {
1876         (*s)++;
1877         *ret = rb_complex_new2(ZERO, num);
1878         return 1; /* e.g. "3i" */
1879     }
1880
1881     if (**s == '@') {
1882         int st;
1883
1884         (*s)++;
1885         bb = *b;
1886         st = read_rat(s, strict, b);
1887         **b = '\0';
1888         if (strlen(bb) < 1 ||
1889             !isdecimal(*(bb + strlen(bb) - 1))) {
1890             *ret = rb_complex_new2(num, ZERO);
1891             return 0; /* e.g. "1@-" */
1892         }
1893         num2 = str2num(bb);
1894         *ret = rb_complex_new_polar(num, num2);
1895         if (!st)
1896             return 0; /* e.g. "1@2." */
1897         else
1898             return 1; /* e.g. "1@2" */
1899     }
1900
1901     if (issign(**s)) {
1902         bb = *b;
1903         sign = read_sign(s, b);
1904         if (isimagunit(**s))
1905             num2 = INT2FIX((sign == '-') ? -1 : + 1);
1906         else {
1907             if (!read_rat_nos(s, strict, b)) {
1908                 *ret = rb_complex_new2(num, ZERO);
1909                 return 0; /* e.g. "1+xi" */
1910             }
1911             **b = '\0';
1912             num2 = str2num(bb);
1913         }
1914         if (!isimagunit(**s)) {
1915             *ret = rb_complex_new2(num, ZERO);
1916             return 0; /* e.g. "1+3x" */
1917         }
1918         (*s)++;
1919         *ret = rb_complex_new2(num, num2);
1920         return 1; /* e.g. "1+2i" */
1921     }
1922     /* !(@, - or +) */
1923     {
1924         *ret = rb_complex_new2(num, ZERO);
1925         return 1; /* e.g. "3" */
1926     }
1927 }
1928
1929 inline static void
1930 skip_ws(const char **s)
1931 {
1932     while (isspace((unsigned char)**s))
1933         (*s)++;
1934 }
1935
1936 static int
1937 parse_comp(const char *s, int strict, VALUE *num)
1938 {
1939     char *buf, *b;
1940     VALUE tmp;
1941     int ret = 1;
1942
1943     buf = ALLOCV_N(char, tmp, strlen(s) + 1);
1944     b = buf;
1945
1946     skip_ws(&s);
1947     if (!read_comp(&s, strict, num, &b)) {
1948         ret = 0;
1949     }
1950     else {
1951         skip_ws(&s);
1952
1953         if (strict)
1954             if (*s != '\0')
1955                 ret = 0;
1956     }
1957     ALLOCV_END(tmp);
1958
1959     return ret;
1960 }
1961
1962 static VALUE
1963 string_to_c_strict(VALUE self, int raise)
1964 {
1965     char *s;
1966     VALUE num;
1967
1968     rb_must_asciicompat(self);
1969
1970     s = RSTRING_PTR(self);
1971
1972     if (!s || memchr(s, '\0', RSTRING_LEN(self))) {
1973         if (!raise) return Qnil;
1974         rb_raise(rb_eArgError, "string contains null byte");
1975     }
1976
1977     if (s && s[RSTRING_LEN(self)]) {
1978         rb_str_modify(self);
1979         s = RSTRING_PTR(self);
1980         s[RSTRING_LEN(self)] = '\0';
1981     }
1982
1983     if (!s)
1984         s = (char *)"";
1985
1986     if (!parse_comp(s, 1, &num)) {
1987         if (!raise) return Qnil;
1988         rb_raise(rb_eArgError, "invalid value for convert(): %+"PRIsVALUE,
1989                  self);
1990     }
1991
1992     return num;
1993 }
1994
1995 /*
1996  * call-seq:
1997  *    str.to_c  ->  complex
1998  *
1999  * Returns a complex which denotes the string form.  The parser
2000  * ignores leading whitespaces and trailing garbage.  Any digit
2001  * sequences can be separated by an underscore.  Returns zero for null
2002  * or garbage string.
2003  *
2004  *    '9'.to_c           #=> (9+0i)
2005  *    '2.5'.to_c         #=> (2.5+0i)
2006  *    '2.5/1'.to_c       #=> ((5/2)+0i)
2007  *    '-3/2'.to_c        #=> ((-3/2)+0i)
2008  *    '-i'.to_c          #=> (0-1i)
2009  *    '45i'.to_c         #=> (0+45i)
2010  *    '3-4i'.to_c        #=> (3-4i)
2011  *    '-4e2-4e-2i'.to_c  #=> (-400.0-0.04i)
2012  *    '-0.0-0.0i'.to_c   #=> (-0.0-0.0i)
2013  *    '1/2+3/4i'.to_c    #=> ((1/2)+(3/4)*i)
2014  *    'ruby'.to_c        #=> (0+0i)
2015  *
2016  * See Kernel.Complex.
2017  */
2018 static VALUE
2019 string_to_c(VALUE self)
2020 {
2021     char *s;
2022     VALUE num;
2023
2024     rb_must_asciicompat(self);
2025
2026     s = RSTRING_PTR(self);
2027
2028     if (s && s[RSTRING_LEN(self)]) {
2029         rb_str_modify(self);
2030         s = RSTRING_PTR(self);
2031         s[RSTRING_LEN(self)] = '\0';
2032     }
2033
2034     if (!s)
2035         s = (char *)"";
2036
2037     (void)parse_comp(s, 0, &num);
2038
2039     return num;
2040 }
2041
2042 static VALUE
2043 to_complex(VALUE val)
2044 {
2045     return rb_convert_type(val, T_COMPLEX, "Complex", "to_c");
2046 }
2047
2048 static VALUE
2049 nucomp_convert(VALUE klass, VALUE a1, VALUE a2, int raise)
2050 {
2051     if (NIL_P(a1) || NIL_P(a2)) {
2052         if (!raise) return Qnil;
2053         rb_raise(rb_eTypeError, "can't convert nil into Complex");
2054     }
2055
2056     if (RB_TYPE_P(a1, T_STRING)) {
2057         a1 = string_to_c_strict(a1, raise);
2058         if (NIL_P(a1)) return Qnil;
2059     }
2060
2061     if (RB_TYPE_P(a2, T_STRING)) {
2062         a2 = string_to_c_strict(a2, raise);
2063         if (NIL_P(a2)) return Qnil;
2064     }
2065
2066     if (RB_TYPE_P(a1, T_COMPLEX)) {
2067         {
2068             get_dat1(a1);
2069
2070             if (k_exact_zero_p(dat->imag))
2071                 a1 = dat->real;
2072         }
2073     }
2074
2075     if (RB_TYPE_P(a2, T_COMPLEX)) {
2076         {
2077             get_dat1(a2);
2078
2079             if (k_exact_zero_p(dat->imag))
2080                 a2 = dat->real;
2081         }
2082     }
2083
2084     if (RB_TYPE_P(a1, T_COMPLEX)) {
2085         if (a2 == Qundef || (k_exact_zero_p(a2)))
2086             return a1;
2087     }
2088
2089     if (a2 == Qundef) {
2090         if (k_numeric_p(a1) && !f_real_p(a1))
2091             return a1;
2092         /* should raise exception for consistency */
2093         if (!k_numeric_p(a1)) {
2094             if (!raise)
2095                 return rb_protect(to_complex, a1, NULL);
2096             return to_complex(a1);
2097         }
2098     }
2099     else {
2100         if ((k_numeric_p(a1) && k_numeric_p(a2)) &&
2101             (!f_real_p(a1) || !f_real_p(a2)))
2102             return f_add(a1,
2103                          f_mul(a2,
2104                                f_complex_new_bang2(rb_cComplex, ZERO, ONE)));
2105     }
2106
2107     {
2108         int argc;
2109         VALUE argv2[2];
2110         argv2[0] = a1;
2111         if (a2 == Qundef) {
2112             argv2[1] = Qnil;
2113             argc = 1;
2114         }
2115         else {
2116             if (!raise && !RB_INTEGER_TYPE_P(a2) && !RB_FLOAT_TYPE_P(a2) && !RB_TYPE_P(a2, T_RATIONAL))
2117                 return Qnil;
2118             argv2[1] = a2;
2119             argc = 2;
2120         }
2121         return nucomp_s_new(argc, argv2, klass);
2122     }
2123 }
2124
2125 static VALUE
2126 nucomp_s_convert(int argc, VALUE *argv, VALUE klass)
2127 {
2128     VALUE a1, a2;
2129
2130     if (rb_scan_args(argc, argv, "11", &a1, &a2) == 1) {
2131         a2 = Qundef;
2132     }
2133
2134     return nucomp_convert(klass, a1, a2, TRUE);
2135 }
2136
2137 /*
2138  * call-seq:
2139  *    num.real  ->  self
2140  *
2141  * Returns self.
2142  */
2143 static VALUE
2144 numeric_real(VALUE self)
2145 {
2146     return self;
2147 }
2148
2149 /*
2150  * call-seq:
2151  *    num.imag       ->  0
2152  *    num.imaginary  ->  0
2153  *
2154  * Returns zero.
2155  */
2156 static VALUE
2157 numeric_imag(VALUE self)
2158 {
2159     return INT2FIX(0);
2160 }
2161
2162 /*
2163  * call-seq:
2164  *    num.abs2  ->  real
2165  *
2166  * Returns square of self.
2167  */
2168 static VALUE
2169 numeric_abs2(VALUE self)
2170 {
2171     return f_mul(self, self);
2172 }
2173
2174 /*
2175  * call-seq:
2176  *    num.arg    ->  0 or float
2177  *    num.angle  ->  0 or float
2178  *    num.phase  ->  0 or float
2179  *
2180  * Returns 0 if the value is positive, pi otherwise.
2181  */
2182 static VALUE
2183 numeric_arg(VALUE self)
2184 {
2185     if (f_positive_p(self))
2186         return INT2FIX(0);
2187     return DBL2NUM(M_PI);
2188 }
2189
2190 /*
2191  * call-seq:
2192  *    num.rect  ->  array
2193  *    num.rectangular  ->  array
2194  *
2195  * Returns an array; [num, 0].
2196  */
2197 static VALUE
2198 numeric_rect(VALUE self)
2199 {
2200     return rb_assoc_new(self, INT2FIX(0));
2201 }
2202
2203 /*
2204  * call-seq:
2205  *    num.polar  ->  array
2206  *
2207  * Returns an array; [num.abs, num.arg].
2208  */
2209 static VALUE
2210 numeric_polar(VALUE self)
2211 {
2212     VALUE abs, arg;
2213
2214     if (RB_INTEGER_TYPE_P(self)) {
2215         abs = rb_int_abs(self);
2216         arg = numeric_arg(self);
2217     }
2218     else if (RB_FLOAT_TYPE_P(self)) {
2219         abs = rb_float_abs(self);
2220         arg = float_arg(self);
2221     }
2222     else if (RB_TYPE_P(self, T_RATIONAL)) {
2223         abs = rb_rational_abs(self);
2224         arg = numeric_arg(self);
2225     }
2226     else {
2227         abs = f_abs(self);
2228         arg = f_arg(self);
2229     }
2230     return rb_assoc_new(abs, arg);
2231 }
2232
2233 /*
2234  * call-seq:
2235  *    num.conj       ->  self
2236  *    num.conjugate  ->  self
2237  *
2238  * Returns self.
2239  */
2240 static VALUE
2241 numeric_conj(VALUE self)
2242 {
2243     return self;
2244 }
2245
2246 /*
2247  * call-seq:
2248  *    flo.arg    ->  0 or float
2249  *    flo.angle  ->  0 or float
2250  *    flo.phase  ->  0 or float
2251  *
2252  * Returns 0 if the value is positive, pi otherwise.
2253  */
2254 static VALUE
2255 float_arg(VALUE self)
2256 {
2257     if (isnan(RFLOAT_VALUE(self)))
2258         return self;
2259     if (f_tpositive_p(self))
2260         return INT2FIX(0);
2261     return rb_const_get(rb_mMath, id_PI);
2262 }
2263
2264 /*
2265  * A complex number can be represented as a paired real number with
2266  * imaginary unit; a+bi.  Where a is real part, b is imaginary part
2267  * and i is imaginary unit.  Real a equals complex a+0i
2268  * mathematically.
2269  *
2270  * You can create a \Complex object explicitly with:
2271  *
2272  * - A {complex literal}[doc/syntax/literals_rdoc.html#label-Complex+Literals].
2273  *
2274  * You can convert certain objects to \Complex objects with:
2275  *
2276  * - \Method {Complex}[Kernel.html#method-i-Complex].
2277  *
2278  * Complex object can be created as literal, and also by using
2279  * Kernel#Complex, Complex::rect, Complex::polar or to_c method.
2280  *
2281  *    2+1i                 #=> (2+1i)
2282  *    Complex(1)           #=> (1+0i)
2283  *    Complex(2, 3)        #=> (2+3i)
2284  *    Complex.polar(2, 3)  #=> (-1.9799849932008908+0.2822400161197344i)
2285  *    3.to_c               #=> (3+0i)
2286  *
2287  * You can also create complex object from floating-point numbers or
2288  * strings.
2289  *
2290  *    Complex(0.3)         #=> (0.3+0i)
2291  *    Complex('0.3-0.5i')  #=> (0.3-0.5i)
2292  *    Complex('2/3+3/4i')  #=> ((2/3)+(3/4)*i)
2293  *    Complex('1@2')       #=> (-0.4161468365471424+0.9092974268256817i)
2294  *
2295  *    0.3.to_c             #=> (0.3+0i)
2296  *    '0.3-0.5i'.to_c      #=> (0.3-0.5i)
2297  *    '2/3+3/4i'.to_c      #=> ((2/3)+(3/4)*i)
2298  *    '1@2'.to_c           #=> (-0.4161468365471424+0.9092974268256817i)
2299  *
2300  * A complex object is either an exact or an inexact number.
2301  *
2302  *    Complex(1, 1) / 2    #=> ((1/2)+(1/2)*i)
2303  *    Complex(1, 1) / 2.0  #=> (0.5+0.5i)
2304  */
2305 void
2306 Init_Complex(void)
2307 {
2308     VALUE compat;
2309     id_abs = rb_intern_const("abs");
2310     id_arg = rb_intern_const("arg");
2311     id_denominator = rb_intern_const("denominator");
2312     id_numerator = rb_intern_const("numerator");
2313     id_real_p = rb_intern_const("real?");
2314     id_i_real = rb_intern_const("@real");
2315     id_i_imag = rb_intern_const("@image"); /* @image, not @imag */
2316     id_finite_p = rb_intern_const("finite?");
2317     id_infinite_p = rb_intern_const("infinite?");
2318     id_rationalize = rb_intern_const("rationalize");
2319     id_PI = rb_intern_const("PI");
2320
2321     rb_cComplex = rb_define_class("Complex", rb_cNumeric);
2322
2323     rb_define_alloc_func(rb_cComplex, nucomp_s_alloc);
2324     rb_undef_method(CLASS_OF(rb_cComplex), "allocate");
2325
2326     rb_undef_method(CLASS_OF(rb_cComplex), "new");
2327
2328     rb_define_singleton_method(rb_cComplex, "rectangular", nucomp_s_new, -1);
2329     rb_define_singleton_method(rb_cComplex, "rect", nucomp_s_new, -1);
2330     rb_define_singleton_method(rb_cComplex, "polar", nucomp_s_polar, -1);
2331
2332     rb_define_global_function("Complex", nucomp_f_complex, -1);
2333
2334     rb_undef_methods_from(rb_cComplex, RCLASS_ORIGIN(rb_mComparable));
2335     rb_undef_method(rb_cComplex, "%");
2336     rb_undef_method(rb_cComplex, "div");
2337     rb_undef_method(rb_cComplex, "divmod");
2338     rb_undef_method(rb_cComplex, "floor");
2339     rb_undef_method(rb_cComplex, "ceil");
2340     rb_undef_method(rb_cComplex, "modulo");
2341     rb_undef_method(rb_cComplex, "remainder");
2342     rb_undef_method(rb_cComplex, "round");
2343     rb_undef_method(rb_cComplex, "step");
2344     rb_undef_method(rb_cComplex, "truncate");
2345     rb_undef_method(rb_cComplex, "i");
2346
2347     rb_define_method(rb_cComplex, "real", rb_complex_real, 0);
2348     rb_define_method(rb_cComplex, "imaginary", rb_complex_imag, 0);
2349     rb_define_method(rb_cComplex, "imag", rb_complex_imag, 0);
2350
2351     rb_define_method(rb_cComplex, "-@", rb_complex_uminus, 0);
2352     rb_define_method(rb_cComplex, "+", rb_complex_plus, 1);
2353     rb_define_method(rb_cComplex, "-", rb_complex_minus, 1);
2354     rb_define_method(rb_cComplex, "*", rb_complex_mul, 1);
2355     rb_define_method(rb_cComplex, "/", rb_complex_div, 1);
2356     rb_define_method(rb_cComplex, "quo", nucomp_quo, 1);
2357     rb_define_method(rb_cComplex, "fdiv", nucomp_fdiv, 1);
2358     rb_define_method(rb_cComplex, "**", rb_complex_pow, 1);
2359
2360     rb_define_method(rb_cComplex, "==", nucomp_eqeq_p, 1);
2361     rb_define_method(rb_cComplex, "<=>", nucomp_cmp, 1);
2362     rb_define_method(rb_cComplex, "coerce", nucomp_coerce, 1);
2363
2364     rb_define_method(rb_cComplex, "abs", rb_complex_abs, 0);
2365     rb_define_method(rb_cComplex, "magnitude", rb_complex_abs, 0);
2366     rb_define_method(rb_cComplex, "abs2", nucomp_abs2, 0);
2367     rb_define_method(rb_cComplex, "arg", rb_complex_arg, 0);
2368     rb_define_method(rb_cComplex, "angle", rb_complex_arg, 0);
2369     rb_define_method(rb_cComplex, "phase", rb_complex_arg, 0);
2370     rb_define_method(rb_cComplex, "rectangular", nucomp_rect, 0);
2371     rb_define_method(rb_cComplex, "rect", nucomp_rect, 0);
2372     rb_define_method(rb_cComplex, "polar", nucomp_polar, 0);
2373     rb_define_method(rb_cComplex, "conjugate", rb_complex_conjugate, 0);
2374     rb_define_method(rb_cComplex, "conj", rb_complex_conjugate, 0);
2375
2376     rb_define_method(rb_cComplex, "real?", nucomp_real_p_m, 0);
2377
2378     rb_define_method(rb_cComplex, "numerator", nucomp_numerator, 0);
2379     rb_define_method(rb_cComplex, "denominator", nucomp_denominator, 0);
2380
2381     rb_define_method(rb_cComplex, "hash", nucomp_hash, 0);
2382     rb_define_method(rb_cComplex, "eql?", nucomp_eql_p, 1);
2383
2384     rb_define_method(rb_cComplex, "to_s", nucomp_to_s, 0);
2385     rb_define_method(rb_cComplex, "inspect", nucomp_inspect, 0);
2386
2387     rb_undef_method(rb_cComplex, "positive?");
2388     rb_undef_method(rb_cComplex, "negative?");
2389
2390     rb_define_method(rb_cComplex, "finite?", rb_complex_finite_p, 0);
2391     rb_define_method(rb_cComplex, "infinite?", rb_complex_infinite_p, 0);
2392
2393     rb_define_private_method(rb_cComplex, "marshal_dump", nucomp_marshal_dump, 0);
2394     /* :nodoc: */
2395     compat = rb_define_class_under(rb_cComplex, "compatible", rb_cObject);
2396     rb_define_private_method(compat, "marshal_load", nucomp_marshal_load, 1);
2397     rb_marshal_define_compat(rb_cComplex, compat, nucomp_dumper, nucomp_loader);
2398
2399     rb_define_method(rb_cComplex, "to_i", nucomp_to_i, 0);
2400     rb_define_method(rb_cComplex, "to_f", nucomp_to_f, 0);
2401     rb_define_method(rb_cComplex, "to_r", nucomp_to_r, 0);
2402     rb_define_method(rb_cComplex, "rationalize", nucomp_rationalize, -1);
2403     rb_define_method(rb_cComplex, "to_c", nucomp_to_c, 0);
2404     rb_define_method(rb_cNilClass, "to_c", nilclass_to_c, 0);
2405     rb_define_method(rb_cNumeric, "to_c", numeric_to_c, 0);
2406
2407     rb_define_method(rb_cString, "to_c", string_to_c, 0);
2408
2409     rb_define_private_method(CLASS_OF(rb_cComplex), "convert", nucomp_s_convert, -1);
2410
2411     rb_define_method(rb_cNumeric, "real", numeric_real, 0);
2412     rb_define_method(rb_cNumeric, "imaginary", numeric_imag, 0);
2413     rb_define_method(rb_cNumeric, "imag", numeric_imag, 0);
2414     rb_define_method(rb_cNumeric, "abs2", numeric_abs2, 0);
2415     rb_define_method(rb_cNumeric, "arg", numeric_arg, 0);
2416     rb_define_method(rb_cNumeric, "angle", numeric_arg, 0);
2417     rb_define_method(rb_cNumeric, "phase", numeric_arg, 0);
2418     rb_define_method(rb_cNumeric, "rectangular", numeric_rect, 0);
2419     rb_define_method(rb_cNumeric, "rect", numeric_rect, 0);
2420     rb_define_method(rb_cNumeric, "polar", numeric_polar, 0);
2421     rb_define_method(rb_cNumeric, "conjugate", numeric_conj, 0);
2422     rb_define_method(rb_cNumeric, "conj", numeric_conj, 0);
2423
2424     rb_define_method(rb_cFloat, "arg", float_arg, 0);
2425     rb_define_method(rb_cFloat, "angle", float_arg, 0);
2426     rb_define_method(rb_cFloat, "phase", float_arg, 0);
2427
2428     /*
2429      * The imaginary unit.
2430      */
2431     rb_define_const(rb_cComplex, "I",
2432                     f_complex_new_bang2(rb_cComplex, ZERO, ONE));
2433
2434 #if !USE_FLONUM
2435     rb_gc_register_mark_object(RFLOAT_0 = DBL2NUM(0.0));
2436 #endif
2437
2438     rb_provide("complex.so");   /* for backward compatibility */
2439 }