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