src/pmc/complex.pmc

   1 /*
   2 Copyright (C) 2004-2009, Parrot Foundation.
   3 $Id$
   4
   5 =head1 NAME
   6
   7 src/pmc/complex.pmc - Complex Numbers PMC Class
   8
   9 =head1 DESCRIPTION
  10
  11 C<Complex> provides a representation of complex numbers. It handles
  12 string parsing/generating and basic mathematical operations.
  13
  14 =head2 Functions
  15
  16 =over 4
  17
  18 =cut
  19
  20 Equations used are sometimes listed.  At times, multiple equations are given,
  21 but those starting with => are the ones used
  22
  23 */
  24
  25 /*
  26
  27 =item C<static void
  28 complex_parse_string(PARROT_INTERP, FLOATVAL *re, FLOATVAL *im, STRING *value)>
  29
  30 Parses the string in C<value> to produce a complex number, represented
  31 by the real (C<*re>) and imaginary (C<*im>) parts. Raises an exception
  32 if it cannot understand the string.  The string should be of the form
  33 C<a+bi> with optional spaces around C<+> and before C<i>. You can also
  34 use C<j> instead of C<i>.
  35
  36 =cut
  37
  38 */
  39
  40 static void
  41 complex_parse_string(PARROT_INTERP, FLOATVAL *re, FLOATVAL *im, STRING *value) {
  42     char   * const str        = Parrot_str_to_cstring(interp, value);
  43     char   *t                 = str;
  44     char   *first_num_offset  = str;
  45     char   *second_num_offset = NULL;
  46     STRING *S;
  47
  48     INTVAL  i                 = 0;
  49     INTVAL  first_num_minus   = 0;
  50     INTVAL  second_num_minus  = 0;
  51     UINTVAL first_num_length, second_num_length;
  52
  53     /* walk the string and identify the real and imaginary parts */
  54
  55     if (*t == '-') {
  56         /* first number is negative */
  57         t++;
  58         first_num_minus = 1;
  59
  60         /* allow for an optional space */
  61         if (*t == ' ')
  62             t++;
  63         first_num_offset = t;
  64     }
  65
  66     /* skip digits */
  67     while (*t >= '0' && *t <= '9')
  68         t++;
  69
  70     if (*t == '.') {
  71         /* this number has a decimal point */
  72         t++;
  73
  74         /* skip digits */
  75         while (*t >= '0' && *t <= '9')
  76             t++;
  77     }
  78
  79     /* save the length of the real part */
  80     first_num_length = t - first_num_offset;
  81
  82     /* end of string; we only have a real part */
  83     if (*t == 0) {
  84         second_num_length = 0;
  85     }
  86     else if ((*t == 'i' || *t == 'j') && *(t+1) == 0) {
  87         /* there is only an imaginary part, so the first number was
  88             actually the imaginary part */
  89         second_num_length = first_num_length;
  90         first_num_length  = 0;
  91         second_num_offset = first_num_offset;
  92         second_num_minus  = first_num_minus;
  93         first_num_minus   = 0;
  94
  95         /* this is useful if there is no number for
  96             the imaginary part, like in "-i" */
  97         i = 1;
  98     }
  99     else {
 100         /* skip an optional space */
 101         if (*t == ' ')
 102             t++;
 103
 104         /* expect "+" or "-" and the imaginary part */
 105         if (*t == '+' || *t == '-') {
 106             /* save the sign */
 107             second_num_minus = (*t == '-');
 108             t++;
 109
 110             /* skip another optional space */
 111             if (*t == ' ')
 112                 t++;
 113
 114             /* save the beginning of the imaginary part */
 115             second_num_offset = t;
 116
 117             /* skip digits */
 118             while (*t >= '0' && *t <= '9')
 119                 t++;
 120
 121             if (*t == '.') {
 122                 /* this number has a decimal point */
 123                 t++;
 124
 125                 /* skip digits */
 126                 while (*t >= '0' && *t <= '9')
 127                     t++;
 128             }
 129
 130             /* save the length of the imaginary part */
 131             second_num_length = t - second_num_offset;
 132
 133             /* allow for one more optional space */
 134             if (*t == ' ')
 135                 t++;
 136
 137             /* verify that the string ends properly */
 138             if ((*t != 'i' && *t != 'j') || (*(t+1) != 0)) {
 139                 /* imaginary part does not end in 'i' or 'j' */
 140                 Parrot_ex_throw_from_c_args(interp, NULL,
 141                     EXCEPTION_INVALID_STRING_REPRESENTATION,
 142                     "Complex: malformed string");
 143             }
 144
 145             /* this is useful if there is no number for the
 146                 imaginary part, like in "2+i" */
 147             i = 1;
 148
 149             /* all is OK, save the number */
 150         }
 151         else {
 152             /* "+" or "-" not found: error */
 153             Parrot_str_free_cstring(str);
 154
 155             Parrot_ex_throw_from_c_args(interp, NULL,
 156                 EXCEPTION_INVALID_STRING_REPRESENTATION,
 157                 "Complex: malformed string");
 158         }
 159     }
 160
 161     /* now we have the offsets and the lengths we turn them into float values */
 162
 163     if (first_num_length) {
 164         /* there is a real part, interpret it */
 165         S   = Parrot_str_new(interp, first_num_offset, first_num_length);
 166         *re = Parrot_str_to_num(interp, S);
 167     }
 168     else {
 169         /* consider the real part 0.0 */
 170         *re = 0.0;
 171     }
 172
 173     if (second_num_length) {
 174         /* there is an imaginary part, interpret it */
 175         S   = Parrot_str_new(interp, second_num_offset, second_num_length);
 176         *im = Parrot_str_to_num(interp, S);
 177     }
 178     else {
 179         /* consider the imaginary part 0.0 */
 180         if (i) /* the string was something like "1+i" */
 181             *im = 1.0;
 182         else
 183             *im = 0.0;
 184     }
 185
 186     if (first_num_minus)
 187         *re = -*re;
 188
 189     if (second_num_minus)
 190         *im = -*im;
 191
 192     Parrot_str_free_cstring(str);
 193 }
 194
 195 static void
 196 int_check_divide_zero(PARROT_INTERP, INTVAL value) {
 197     if (value == 0)
 198         Parrot_ex_throw_from_c_args(interp, NULL, EXCEPTION_DIV_BY_ZERO,
 199             "Divide by zero");
 200 }
 201
 202 static void
 203 float_check_divide_zero(PARROT_INTERP, FLOATVAL value) {
 204     if (FLOAT_IS_ZERO(value))
 205         Parrot_ex_throw_from_c_args(interp, NULL, EXCEPTION_DIV_BY_ZERO,
 206             "Divide by zero");
 207 }
 208
 209 static void
 210 complex_check_divide_zero(PARROT_INTERP, PMC *value) {
 211     /* Throw an exception if we are dividing by zero. Check both the real part
 212      * and the imaginary part.*/
 213     if (FLOAT_IS_ZERO(VTABLE_get_number_keyed_int(interp, value, 0))
 214             && FLOAT_IS_ZERO(VTABLE_get_number_keyed_int(interp, value, 1)))
 215         Parrot_ex_throw_from_c_args(interp, NULL, EXCEPTION_DIV_BY_ZERO,
 216             "Divide by zero");
 217 }
 218
 219
 220 pmclass Complex auto_attrs {
 221
 222     ATTR FLOATVAL re; /* real part */
 223     ATTR FLOATVAL im; /* imaginary part */
 224
 225 /*
 226
 227 =item C<PMC *instantiate(PMC *sig)>
 228
 229 Create a new complex PMC with passed arguments according to pdd03.
 230
 231 =cut
 232
 233 */
 234
 235     VTABLE PMC *instantiate(PMC *sig) {
 236         return PMCNULL;
 237
 238         /* TODO -- really build this thing */
 239 #if 0
 240         PMC *res = pmc_new(INTERP, enum_class_Complex);
 241         FLOATVAL re    = 0.0;
 242         FLOATVAL im    = 0.0;
 243         int      argcI = REG_INT(interp, 1);
 244         int      argcS = REG_INT(interp, 2);
 245         int      argcP = REG_INT(interp, 3);
 246         int      argcN = REG_INT(interp, 4);
 247
 248         /*
 249          * we can only allow 0..2 arguments of one kind. For
 250          * mixed (e.g. N, P) args the order of arguments isn't fixed
 251          */
 252         if (argcP) {
 253             re = VTABLE_get_number(INTERP, REG_PMC(interp, 5));
 254
 255             if (argcP == 2)
 256                 im = VTABLE_get_number(INTERP, REG_PMC(interp, 6));
 257             /* TODO throw exception if argument mismatch */
 258         }
 259         else if (argcN) {
 260             re = REG_NUM(interp, 5);
 261
 262             if (argcN == 2)
 263                 im = REG_NUM(interp, 6);
 264         }
 265         else if (argcI) {
 266             re = REG_INT(interp, 5);
 267
 268             if (argcI == 2)
 269                 im = REG_INT(interp, 6);
 270         }
 271         else if (argcS == 1) {
 272             complex_parse_string(INTERP, &RE(res), &IM(res), REG_STR(interp, 5));
 273             return res;
 274         }
 275
 276         RE(res) = re;
 277         IM(res) = im;
 278         return res;
 279 #endif
 280     }
 281
 282 /*
 283
 284 =item C<opcode_t *invoke(void *next)>
 285
 286 Pythonic object constructor. SELF is a Complex Class object. Return a new
 287 C<complex> object according to 2.1. Built-in Functions.
 288
 289 =cut
 290
 291 */
 292     VTABLE opcode_t *invoke(void *next) {
 293         const int   argcP = REG_INT(interp, 3);
 294         PMC * const res   = pmc_new(INTERP, enum_class_Complex);
 295
 296         if (argcP == 1) {
 297             PMC * const arg = REG_PMC(interp, 5);
 298
 299             if (arg->vtable->base_type == enum_class_String)
 300                 VTABLE_set_string_native(INTERP, res, VTABLE_get_string(interp, arg));
 301             else
 302                 SET_ATTR_re(INTERP, res, VTABLE_get_number(INTERP, arg));
 303         }
 304         else if (argcP == 2) {
 305             SET_ATTR_re(INTERP, res, VTABLE_get_number(INTERP, REG_PMC(interp, 5)));
 306             SET_ATTR_im(INTERP, res, VTABLE_get_number(INTERP, REG_PMC(interp, 6)));
 307         }
 308
 309         REG_PMC(interp, 5) = res;
 310         return (opcode_t *)next;
 311     }
 312
 313 /*
 314
 315 =back
 316
 317 =head2 Methods
 318
 319 =over 4
 320
 321 =item C<void init()>
 322
 323 Initializes the complex number with the value 0+0i.
 324
 325 =item C<void init_pmc(PMC *initializer)>
 326
 327 Initializes the complex number with the specified initializer.
 328 The initializer can be a string PMC or a numeric array with (real, imag)
 329
 330 =item C<PMC *clone()>
 331
 332 Creates an identical copy of the complex number.
 333
 334 =cut
 335
 336 */
 337
 338     VTABLE void init() {
 339         SET_ATTR_re(INTERP, SELF, 0.0);
 340         SET_ATTR_im(INTERP, SELF, 0.0);
 341     }
 342
 343     VTABLE void init_pmc(PMC *initializer) {
 344         const INTVAL arg_type = VTABLE_type(interp, initializer);
 345         SELF.init();
 346         switch (arg_type) {
 347             case enum_class_String:
 348                 SELF.set_string_native(VTABLE_get_string(interp, initializer));
 349                 break;
 350             case enum_class_FixedFloatArray:
 351             case enum_class_ResizableFloatArray:
 352             case enum_class_FixedIntegerArray:
 353             case enum_class_ResizableIntegerArray:
 354                 if (VTABLE_get_integer(interp, initializer) == 2) {
 355                     FLOATVAL re = VTABLE_get_number_keyed_int(interp, initializer, 0);
 356                     FLOATVAL im = VTABLE_get_number_keyed_int(interp, initializer, 1);
 357                     SET_ATTR_re(INTERP, SELF, re);
 358                     SET_ATTR_im(INTERP, SELF, im);
 359                     break;
 360                 }
 361                 /* else let it fall to default */
 362             default:
 363                 if (VTABLE_isa(interp, initializer, CONST_STRING(interp, "String"))) {
 364                     STRING * s = VTABLE_get_string(interp, initializer);
 365                     SELF.set_string_native(s);
 366                 }
 367                 else {
 368                     Parrot_ex_throw_from_c_args(interp, NULL,
 369                         EXCEPTION_INVALID_OPERATION,
 370                         "Invalid Complex initializer");
 371                 }
 372         }
 373     }
 374
 375     VTABLE PMC *clone() {
 376         PMC * const dest = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
 377         FLOATVAL re, im;
 378
 379         GET_ATTR_re(INTERP, SELF, re);
 380         GET_ATTR_im(INTERP, SELF, im);
 381         SET_ATTR_re(INTERP, dest, re);
 382         SET_ATTR_im(INTERP, dest, im);
 383
 384         return dest;
 385     }
 386
 387 /*
 388
 389 =item C<INTVAL get_integer()>
 390
 391 Returns the modulus of the complex number as an integer.
 392
 393 =item C<FLOATVAL get_number()>
 394
 395 Returns the modulus of the complex number.
 396
 397 =item C<STRING *get_string()>
 398
 399 Returns the complex number as a string in the form C<a+bi>.
 400
 401 =item C<INTVAL get_bool()>
 402
 403 Returns true if the complex number is non-zero.
 404
 405 =cut
 406
 407 */
 408
 409     VTABLE INTVAL get_integer() {
 410         const FLOATVAL f = SELF.get_number();
 411         return (INTVAL)f;
 412     }
 413
 414     VTABLE FLOATVAL get_number() {
 415         FLOATVAL re, im;
 416         GET_ATTR_re(INTERP, SELF, re);
 417         GET_ATTR_im(INTERP, SELF, im);
 418         return sqrt(re * re + im * im);
 419     }
 420
 421     VTABLE STRING *get_string() {
 422         FLOATVAL re, im;
 423         GET_ATTR_re(INTERP, SELF, re);
 424         GET_ATTR_im(INTERP, SELF, im);
 425         return Parrot_sprintf_c(INTERP, "%vg%+vgi", re, im);
 426     }
 427
 428     VTABLE INTVAL get_bool() {
 429         FLOATVAL re, im;
 430         GET_ATTR_re(INTERP, SELF, re);
 431         GET_ATTR_im(INTERP, SELF, im);
 432         return !(FLOAT_IS_ZERO(re) && FLOAT_IS_ZERO(im));
 433     }
 434
 435 /*
 436
 437 =item C<INTVAL get_integer_keyed(PMC *key)>
 438
 439 =item C<INTVAL get_integer_keyed_str(STRING *key)>
 440
 441 =item C<FLOATVAL get_number_keyed(PMC *key)>
 442
 443 =item C<FLOATVAL get_number_keyed_str(STRING *key)>
 444
 445 =item C<PMC *get_pmc_keyed(PMC *key)>
 446
 447 =item C<PMC *get_pmc_keyed_str(STRING *key)>
 448
 449 Returns the requested number (real part for C<real> and imaginary for C<imag>).
 450
 451 =cut
 452
 453 */
 454
 455     VTABLE INTVAL get_integer_keyed(PMC *key) {
 456         STRING * const s = VTABLE_get_string(INTERP, key);
 457         return SELF.get_integer_keyed_str(s);
 458     }
 459
 460     VTABLE INTVAL get_integer_keyed_str(STRING *key) {
 461         const FLOATVAL f = SELF.get_number_keyed_str(key);
 462         return (INTVAL)f;
 463     }
 464
 465     VTABLE FLOATVAL get_number_keyed(PMC *key) {
 466         STRING * const s = VTABLE_get_string(INTERP, key);
 467         return SELF.get_number_keyed_str(s);
 468     }
 469
 470     VTABLE FLOATVAL get_number_keyed_str(STRING *key) {
 471         FLOATVAL value;
 472         if (Parrot_str_equal(INTERP, key, CONST_STRING(INTERP, "real"))) {
 473             GET_ATTR_re(INTERP, SELF, value);
 474         }
 475         else if (Parrot_str_equal(INTERP, key, CONST_STRING(INTERP, "imag"))) {
 476             GET_ATTR_im(INTERP, SELF, value);
 477         }
 478         else
 479             Parrot_ex_throw_from_c_args(interp, NULL,
 480                 EXCEPTION_INVALID_OPERATION, "Complex: key is neither 'real' or 'imag'");
 481         return value;
 482     }
 483
 484     VTABLE PMC *get_pmc_keyed(PMC *key) {
 485         if (VTABLE_isa(interp, key, CONST_STRING(interp, "Integer"))) {
 486             const INTVAL i = VTABLE_get_integer(interp, key);
 487             return SELF.get_pmc_keyed_int(i);
 488         }
 489         else {
 490             STRING * const s = VTABLE_get_string(INTERP, key);
 491             return SELF.get_pmc_keyed_str(s);
 492         }
 493     }
 494
 495     VTABLE PMC *get_pmc_keyed_str(STRING *key) {
 496         PMC * const    ret = pmc_new(INTERP, enum_class_Float);
 497         const FLOATVAL val = SELF.get_number_keyed_str(key);
 498         VTABLE_set_number_native(INTERP, ret, val);
 499         return ret;
 500     }
 501
 502 /*
 503
 504 =item C<PMC *get_pmc_keyed_int(INTVAL key)>
 505
 506 Returns the requested number (real part for C<0> and imaginary for C<1>).
 507
 508 =cut
 509
 510 */
 511
 512     VTABLE PMC *get_pmc_keyed_int(INTVAL key) {
 513         PMC * const    ret = pmc_new(INTERP, enum_class_Float);
 514         const FLOATVAL val = SELF.get_number_keyed_int(key);
 515         VTABLE_set_number_native(INTERP, ret, val);
 516         return ret;
 517     }
 518
 519 /*
 520
 521 =item C<FLOATVAL get_number_keyed_int(INTVAL key)>
 522
 523 Quick hack to emulate get_real() and get_imag():
 524
 525   key = 0 ... get real part
 526   key = 1 ... get imag part
 527
 528 =item C<void set_number_keyed_int(INTVAL key, FLOATVAL v)>
 529
 530 Set real or imag depending on key
 531
 532 =cut
 533
 534 */
 535
 536     VTABLE FLOATVAL get_number_keyed_int(INTVAL key) {
 537         FLOATVAL f;
 538         switch (key) {
 539             case 0:
 540                 GET_ATTR_re(INTERP, SELF, f);
 541                 break;
 542             case 1:
 543                 GET_ATTR_im(INTERP, SELF, f);
 544                 break;
 545             default:
 546                 Parrot_ex_throw_from_c_args(interp, NULL,
 547                     EXCEPTION_INVALID_OPERATION, "Complex: key must be 0 or 1");
 548         }
 549         return f;
 550     }
 551
 552     VTABLE void set_number_keyed_int(INTVAL key, FLOATVAL v) {
 553         switch (key) {
 554             case 0:
 555                 SET_ATTR_re(INTERP, SELF, v);
 556                 break;
 557             case 1:
 558                 SET_ATTR_im(INTERP, SELF, v);
 559                 break;
 560             default:
 561                 Parrot_ex_throw_from_c_args(interp, NULL,
 562                     EXCEPTION_INVALID_OPERATION, "Complex: key must be 0 or 1");
 563         }
 564     }
 565 /*
 566
 567 =item C<void set_string_native(STRING *value)>
 568
 569 Parses the string C<value> into a complex number; raises an exception
 570 on failure.
 571
 572 =item C<void set_pmc(PMC *value)>
 573
 574 if C<value> is a Complex PMC then the complex number is set to its
 575 value; otherwise C<value>'s string representation is parsed with
 576 C<set_string_native()>.
 577
 578 =item C<void set_integer_native(INTVAL value)>
 579
 580 =item C<void set_number_native(FLOATVAL value)>
 581
 582 Sets the real part of the complex number to C<value> and the imaginary
 583 part to C<0.0>
 584
 585 =cut
 586
 587 */
 588
 589     VTABLE void set_string_native(STRING *value) {
 590         FLOATVAL re, im;
 591         complex_parse_string(INTERP, &re, &im, value);
 592         SET_ATTR_re(INTERP, SELF, re);
 593         SET_ATTR_im(INTERP, SELF, im);
 594     }
 595
 596     VTABLE void set_pmc(PMC *value) {
 597         if (VTABLE_isa(INTERP, value, CONST_STRING(INTERP, "Complex"))) {
 598             FLOATVAL re, im;
 599             GET_ATTR_re(INTERP, value, re);
 600             GET_ATTR_im(INTERP, value, im);
 601             SET_ATTR_re(INTERP, SELF, re);
 602             SET_ATTR_im(INTERP, SELF, im);
 603         }
 604         else
 605             VTABLE_set_string_native(INTERP, SELF, VTABLE_get_string(INTERP, value));
 606     }
 607
 608     VTABLE void set_integer_native(INTVAL value) {
 609         SELF.set_number_native((FLOATVAL)value);
 610     }
 611
 612     VTABLE void set_number_native(FLOATVAL value) {
 613         SET_ATTR_re(INTERP, SELF, value);
 614         SET_ATTR_im(INTERP, SELF, 0.0);
 615     }
 616
 617 /*
 618
 619 =item C<void set_integer_keyed(PMC *key, INTVAL value)>
 620
 621 =item C<void set_integer_keyed_str(STRING *key, INTVAL value)>
 622
 623 =item C<void set_number_keyed(PMC *key, FLOATVAL value)>
 624
 625 =item C<void set_number_keyed_str(STRING *key, FLOATVAL value)>
 626
 627 =item C<void set_pmc_keyed(PMC *key, PMC *value)>
 628
 629 =item C<void set_pmc_keyed_str(STRING *key, PMC *value)>
 630
 631 Sets the requested number (real part for C<real> and imaginary for C<imag>)
 632 to C<value>.
 633
 634 =cut
 635
 636 */
 637
 638     VTABLE void set_integer_keyed(PMC *key, INTVAL value) {
 639         SELF.set_number_keyed(key, (FLOATVAL)value);
 640     }
 641
 642     VTABLE void set_integer_keyed_str(STRING *key, INTVAL value) {
 643         SELF.set_number_keyed_str(key, (FLOATVAL)value);
 644     }
 645
 646     VTABLE void set_number_keyed(PMC *key, FLOATVAL value) {
 647         if (VTABLE_isa(interp, key, CONST_STRING(interp, "Integer"))) {
 648             const INTVAL i = VTABLE_get_integer(interp, key);
 649             SELF.set_number_keyed_int(i, value);
 650         }
 651         else {
 652             STRING *s = VTABLE_get_string(INTERP, key);
 653             SELF.set_number_keyed_str(s, value);
 654         }
 655     }
 656
 657     VTABLE void set_number_keyed_str(STRING *key, FLOATVAL value) {
 658         if (Parrot_str_equal(INTERP, key, CONST_STRING(INTERP, "real"))) {
 659             SET_ATTR_re(INTERP, SELF, value);
 660         }
 661         else if (Parrot_str_equal(INTERP, key, CONST_STRING(INTERP, "imag"))) {
 662             SET_ATTR_im(INTERP, SELF, value);
 663         }
 664         else
 665             Parrot_ex_throw_from_c_args(interp, NULL,
 666                 EXCEPTION_INVALID_OPERATION, "Complex: key is neither 'real' or 'imag'");
 667     }
 668
 669     VTABLE void set_pmc_keyed(PMC *key, PMC *value) {
 670         const FLOATVAL f = VTABLE_get_number(INTERP, value);
 671         SELF.set_number_keyed(key, f);
 672     }
 673
 674     VTABLE void set_pmc_keyed_str(STRING *key, PMC *value) {
 675         const FLOATVAL f = VTABLE_get_number(INTERP, value);
 676         SELF.set_number_keyed_str(key, f);
 677     }
 678
 679 /*
 680
 681 =item C<PMC *add(PMC *value, PMC *dest)>
 682
 683 =item C<PMC *add_int(INTVAL value, PMC *dest)>
 684
 685 =item C<PMC *add_float(FLOATVAL value, PMC *dest)>
 686
 687 Adds C<value> to the complex number, placing the result in C<dest>.
 688
 689 =cut
 690
 691 */
 692
 693     MULTI PMC *add(Complex value, PMC *dest) {
 694         FLOATVAL self_re, self_im, val_re, val_im;
 695         dest = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
 696
 697         GET_ATTR_re(INTERP, SELF, self_re);
 698         GET_ATTR_im(INTERP, SELF, self_im);
 699
 700         GET_ATTR_re(INTERP, value, val_re);
 701         GET_ATTR_im(INTERP, value, val_im);
 702
 703         SET_ATTR_re(INTERP, dest, self_re + val_re);
 704         SET_ATTR_im(INTERP, dest, self_im + val_im);
 705
 706         return dest;
 707     }
 708
 709     MULTI PMC *add(DEFAULT value, PMC *dest) {
 710         FLOATVAL re, im;
 711         dest = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
 712
 713         GET_ATTR_re(INTERP, SELF, re);
 714         GET_ATTR_im(INTERP, SELF, im);
 715         SET_ATTR_re(INTERP, dest, re + VTABLE_get_number(INTERP, value));
 716         SET_ATTR_im(INTERP, dest, im);
 717
 718         return dest;
 719     }
 720
 721     VTABLE PMC *add_int(INTVAL value, PMC *dest) {
 722         return SELF.add_float((FLOATVAL)value, dest);
 723     }
 724
 725     VTABLE PMC *add_float(FLOATVAL value, PMC *dest) {
 726         FLOATVAL re, im;
 727         dest = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
 728
 729         GET_ATTR_re(INTERP, SELF, re);
 730         GET_ATTR_im(INTERP, SELF, im);
 731         SET_ATTR_re(INTERP, dest, re + value);
 732         SET_ATTR_im(INTERP, dest, im);
 733
 734         return dest;
 735     }
 736
 737     MULTI void i_add(Complex value) {
 738         FLOATVAL self_re, self_im, val_re, val_im;
 739
 740         GET_ATTR_re(INTERP, SELF, self_re);
 741         GET_ATTR_im(INTERP, SELF, self_im);
 742
 743         GET_ATTR_re(INTERP, value, val_re);
 744         GET_ATTR_im(INTERP, value, val_im);
 745
 746         SET_ATTR_re(INTERP, SELF, self_re + val_re);
 747         SET_ATTR_im(INTERP, SELF, self_im + val_im);
 748     }
 749
 750     MULTI void i_add(DEFAULT value) {
 751         FLOATVAL re;
 752
 753         GET_ATTR_re(INTERP, SELF, re);
 754         SET_ATTR_re(INTERP, SELF, re + VTABLE_get_number(INTERP, value));
 755     }
 756
 757     VTABLE void i_add_int(INTVAL value) {
 758         SELF.i_add_float((FLOATVAL)value);
 759     }
 760
 761     VTABLE void i_add_float(FLOATVAL value) {
 762         FLOATVAL re;
 763
 764         GET_ATTR_re(INTERP, SELF, re);
 765         SET_ATTR_re(INTERP, SELF, re + value);
 766     }
 767
 768 /*
 769
 770 =item C<PMC *subtract(PMC *value, PMC *dest)>
 771
 772 =item C<PMC *subtract_int(INTVAL value, PMC *dest)>
 773
 774 =item C<PMC *subtract_float(FLOATVAL value, PMC *dest)>
 775
 776 Subtracts C<value> from the complex number, placing the result in C<dest>.
 777
 778 =cut
 779
 780 */
 781
 782     MULTI PMC *subtract(Complex value, PMC *dest) {
 783         FLOATVAL self_re, self_im, val_re, val_im;
 784         dest = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
 785
 786         GET_ATTR_re(INTERP, SELF, self_re);
 787         GET_ATTR_im(INTERP, SELF, self_im);
 788
 789         GET_ATTR_re(INTERP, value, val_re);
 790         GET_ATTR_im(INTERP, value, val_im);
 791
 792         SET_ATTR_re(INTERP, dest, self_re - val_re);
 793         SET_ATTR_im(INTERP, dest, self_im - val_im);
 794
 795         return dest;
 796     }
 797
 798     MULTI PMC *subtract(DEFAULT value, PMC *dest) {
 799         FLOATVAL re, im;
 800         dest = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
 801
 802         GET_ATTR_re(INTERP, SELF, re);
 803         GET_ATTR_im(INTERP, SELF, im);
 804         SET_ATTR_re(INTERP, dest, re - VTABLE_get_number(INTERP, value));
 805         SET_ATTR_im(INTERP, dest, im);
 806
 807         return dest;
 808     }
 809
 810     VTABLE PMC *subtract_int(INTVAL value, PMC *dest) {
 811         return SELF.subtract_float((FLOATVAL)value, dest);
 812     }
 813
 814     VTABLE PMC *subtract_float(FLOATVAL value, PMC *dest) {
 815         FLOATVAL re, im;
 816         dest = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
 817
 818         GET_ATTR_re(INTERP, SELF, re);
 819         GET_ATTR_im(INTERP, SELF, im);
 820         SET_ATTR_re(INTERP, dest, re - value);
 821         SET_ATTR_im(INTERP, dest, im);
 822
 823         return dest;
 824     }
 825
 826     MULTI void i_subtract(Complex value) {
 827         FLOATVAL self_re, self_im, val_re, val_im;
 828
 829         GET_ATTR_re(INTERP, SELF, self_re);
 830         GET_ATTR_im(INTERP, SELF, self_im);
 831
 832         GET_ATTR_re(INTERP, value, val_re);
 833         GET_ATTR_im(INTERP, value, val_im);
 834
 835         SET_ATTR_re(INTERP, SELF, self_re - val_re);
 836         SET_ATTR_im(INTERP, SELF, self_im - val_im);
 837     }
 838
 839     MULTI void i_subtract(DEFAULT value) {
 840         FLOATVAL re;
 841
 842         GET_ATTR_re(INTERP, SELF, re);
 843         SET_ATTR_re(INTERP, SELF, re - VTABLE_get_number(INTERP, value));
 844     }
 845
 846     VTABLE void i_subtract_int(INTVAL value) {
 847         SELF.i_subtract_float((FLOATVAL) value);
 848     }
 849
 850     VTABLE void i_subtract_float(FLOATVAL value) {
 851         FLOATVAL re;
 852
 853         GET_ATTR_re(INTERP, SELF, re);
 854         SET_ATTR_re(INTERP, SELF, re - value);
 855     }
 856 /*
 857
 858 =item C<PMC *multiply(PMC *value, PMC *dest)>
 859
 860 =item C<PMC *multiply_int(INTVAL value, PMC *dest)>
 861
 862 =item C<PMC *multiply_float(FLOATVAL value, PMC *dest)>
 863
 864 Multiplies the complex number with C<value>, placing the result in C<dest>.
 865
 866 =item C<void i_multiply(PMC *value)>
 867
 868 =item C<void i_multiply_int(INTVAL value)>
 869
 870 =item C<void i_multiply_float(FLOATVAL value)>
 871
 872 Multiplies the complex number SELF inplace with C<value>.
 873
 874 =cut
 875
 876 */
 877
 878 /*
 879
 880   (a+ib)(c+id)=(ac-bd)+i((a+b)(c+d)-ac-bd).
 881   (a+bi)(c+di)=(ac-bd)+i(ad+bc)
 882
 883 */
 884     MULTI PMC *multiply(Complex value, PMC *dest) {
 885         FLOATVAL a, b, c, d;
 886
 887         dest = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
 888
 889         GET_ATTR_re(INTERP, SELF, a);
 890         GET_ATTR_im(INTERP, SELF, b);
 891         GET_ATTR_re(INTERP, value, c);
 892         GET_ATTR_im(INTERP, value, d);
 893         SET_ATTR_re(INTERP, dest, a * c - b * d);
 894         SET_ATTR_im(INTERP, dest, a * d + b * c);
 895
 896         return dest;
 897     }
 898
 899     MULTI PMC *multiply(DEFAULT value, PMC *dest) {
 900         FLOATVAL re, im;
 901         dest = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
 902
 903         GET_ATTR_re(INTERP, SELF, re);
 904         GET_ATTR_im(INTERP, SELF, im);
 905         SET_ATTR_re(INTERP, dest, re * VTABLE_get_number(INTERP, value));
 906         SET_ATTR_im(INTERP, dest, im * VTABLE_get_number(INTERP, value));
 907
 908         return dest;
 909     }
 910
 911     VTABLE PMC *multiply_int(INTVAL value, PMC *dest) {
 912         return SELF.multiply_float((FLOATVAL) value, dest);
 913     }
 914
 915     VTABLE PMC *multiply_float(FLOATVAL value, PMC *dest) {
 916         FLOATVAL re, im;
 917         dest = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
 918
 919         GET_ATTR_re(INTERP, SELF, re);
 920         GET_ATTR_im(INTERP, SELF, im);
 921         SET_ATTR_re(INTERP, dest, re * value);
 922         SET_ATTR_im(INTERP, dest, im * value);
 923
 924         return dest;
 925     }
 926
 927     MULTI void i_multiply(Complex value) {
 928         FLOATVAL a, b, c, d;
 929
 930         GET_ATTR_re(INTERP, SELF, a);
 931         GET_ATTR_im(INTERP, SELF, b);
 932         GET_ATTR_re(INTERP, value, c);
 933         GET_ATTR_im(INTERP, value, d);
 934         SET_ATTR_re(INTERP, SELF, a * c - b * d);
 935         SET_ATTR_im(INTERP, SELF, a * d + b * c);
 936     }
 937
 938     MULTI void i_multiply(DEFAULT value) {
 939         FLOATVAL re, im;
 940
 941         GET_ATTR_re(INTERP, SELF, re);
 942         GET_ATTR_im(INTERP, SELF, im);
 943         SET_ATTR_re(INTERP, SELF, re * VTABLE_get_number(INTERP, value));
 944         SET_ATTR_im(INTERP, SELF, im * VTABLE_get_number(INTERP, value));
 945     }
 946
 947     VTABLE void i_multiply_int(INTVAL value) {
 948         FLOATVAL re, im;
 949
 950         GET_ATTR_re(INTERP, SELF, re);
 951         GET_ATTR_im(INTERP, SELF, im);
 952         SET_ATTR_re(INTERP, SELF, re * value);
 953         SET_ATTR_im(INTERP, SELF, im * value);
 954     }
 955
 956     VTABLE void i_multiply_float(FLOATVAL value) {
 957         FLOATVAL re, im;
 958
 959         GET_ATTR_re(INTERP, SELF, re);
 960         GET_ATTR_im(INTERP, SELF, im);
 961         SET_ATTR_re(INTERP, SELF, re * value);
 962         SET_ATTR_im(INTERP, SELF, im * value);
 963     }
 964
 965 /*
 966
 967 =item C<PMC *divide(PMC *value, PMC *dest)>
 968
 969 =item C<PMC *divide_int(INTVAL value, PMC *dest)>
 970
 971 =item C<PMC *divide_float(FLOATVAL value, PMC *dest)>
 972
 973 Divide the complex number by C<value>, placing the result in C<dest>.
 974
 975 =item C<void i_divide(PMC *value, PMC *dest)>
 976
 977 =item C<void i_divide_int(INTVAL value, PMC *dest)>
 978
 979 =item C<void i_divide_float(FLOATVAL value, PMC *dest)>
 980
 981 Divide the complex number C<SELF> by C<value> inplace.
 982
 983 Throws divide by zero exception if divisor is zero.
 984
 985 =cut
 986
 987 TODO: for better fp precision
 988 http://docs.sun.com/source/806-3568/ncg_goldberg.html
 989 (a+ib)/(c+id) =
 990     (a + b(d/c)) / (c + d(d/c)) + i(b - a(d/c)) / (c + d(d/c)) if |d|<|c|
 991     (b + a(c/d)) / (d + c(c/d)) + i(-a + b(c/d)) / (d + c(c/d)) if |d|>=|c|
 992
 993 */
 994
 995     MULTI PMC *divide(Complex value, PMC *dest) {
 996         FLOATVAL mod, re, im;
 997         FLOATVAL self_re, self_im, val_re, val_im;
 998
 999         complex_check_divide_zero(INTERP, value);
1000         dest = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1001
1002         GET_ATTR_re(INTERP, SELF, self_re);
1003         GET_ATTR_im(INTERP, SELF, self_im);
1004         GET_ATTR_re(INTERP, value, val_re);
1005         GET_ATTR_im(INTERP, value, val_im);
1006
1007         /* a little speed optimisation: cache an intermediate number;
1008             I'm not sure the compiler does this */
1009
1010         if (self_im == 0.0 && val_im == 0.0) {
1011             re = self_re / val_re;
1012             im = 0.0;
1013         }
1014         else {
1015             mod = (val_re * val_re + val_im * val_im);
1016             re  = (self_re * val_re + self_im * val_im) / mod;
1017             im  = (self_im * val_re - self_re * val_im) / mod;
1018         }
1019
1020         SET_ATTR_re(INTERP, dest, re);
1021         SET_ATTR_im(INTERP, dest, im);
1022
1023         return dest;
1024     }
1025
1026     MULTI PMC *divide(DEFAULT value, PMC *dest) {
1027         FLOATVAL re, im;
1028         const FLOATVAL d = VTABLE_get_number(INTERP, value);
1029         float_check_divide_zero(INTERP, d);
1030         dest = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1031
1032         GET_ATTR_re(INTERP, SELF, re);
1033         GET_ATTR_im(INTERP, SELF, im);
1034
1035         SET_ATTR_re(INTERP, dest, re / d);
1036         SET_ATTR_im(INTERP, dest, im / d);
1037
1038         return dest;
1039     }
1040
1041     VTABLE PMC *divide_int(INTVAL value, PMC *dest) {
1042         FLOATVAL re, im;
1043         int_check_divide_zero(INTERP, value);
1044         dest = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1045
1046         GET_ATTR_re(INTERP, SELF, re);
1047         GET_ATTR_im(INTERP, SELF, im);
1048
1049         SET_ATTR_re(INTERP, dest, re / value);
1050         SET_ATTR_im(INTERP, dest, im / value);
1051
1052         return dest;
1053     }
1054
1055     VTABLE PMC *divide_float(FLOATVAL value, PMC *dest) {
1056         FLOATVAL re, im;
1057         float_check_divide_zero(INTERP, value);
1058         dest = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1059
1060         GET_ATTR_re(INTERP, SELF, re);
1061         GET_ATTR_im(INTERP, SELF, im);
1062
1063         SET_ATTR_re(INTERP, dest, re / value);
1064         SET_ATTR_im(INTERP, dest, im / value);
1065
1066         return dest;
1067     }
1068
1069     MULTI void i_divide(Complex value) {
1070         FLOATVAL re, im;
1071         FLOATVAL self_re, self_im, val_re, val_im;
1072
1073         complex_check_divide_zero(INTERP, value);
1074
1075         GET_ATTR_re(INTERP, SELF, self_re);
1076         GET_ATTR_im(INTERP, SELF, self_im);
1077         GET_ATTR_re(INTERP, value, val_re);
1078         GET_ATTR_im(INTERP, value, val_im);
1079
1080         if (self_im == 0.0 && val_im == 0.0) {
1081             re = self_re / val_re;
1082             im = 0.0;
1083         }
1084         else {
1085             /* a little speed optimisation: cache an intermediate number;
1086                I'm not sure the compiler does this */
1087             const FLOATVAL mod = (val_re * val_re + val_im * val_im);
1088             re  = (self_re * val_re + self_im * val_im) / mod;
1089             im  = (self_im * val_re - self_re * val_im) / mod;
1090         }
1091
1092         SET_ATTR_re(INTERP, SELF, re);
1093         SET_ATTR_im(INTERP, SELF, im);
1094
1095     }
1096
1097     MULTI void i_divide(DEFAULT value) {
1098         FLOATVAL re, im;
1099         const FLOATVAL d = VTABLE_get_number(INTERP, value);
1100         float_check_divide_zero(INTERP, d);
1101
1102         GET_ATTR_re(INTERP, SELF, re);
1103         GET_ATTR_im(INTERP, SELF, im);
1104
1105         SET_ATTR_re(INTERP, SELF, re / d);
1106         SET_ATTR_im(INTERP, SELF, im / d);
1107     }
1108
1109     VTABLE void i_divide_int(INTVAL value) {
1110         FLOATVAL re, im;
1111         int_check_divide_zero(INTERP, value);
1112
1113         GET_ATTR_re(INTERP, SELF, re);
1114         GET_ATTR_im(INTERP, SELF, im);
1115
1116         SET_ATTR_re(INTERP, SELF, re / value);
1117         SET_ATTR_im(INTERP, SELF, im / value);
1118     }
1119
1120     VTABLE void i_divide_float(FLOATVAL value) {
1121         FLOATVAL re, im;
1122         float_check_divide_zero(INTERP, value);
1123
1124         GET_ATTR_re(INTERP, SELF, re);
1125         GET_ATTR_im(INTERP, SELF, im);
1126
1127         SET_ATTR_re(INTERP, SELF, re / value);
1128         SET_ATTR_im(INTERP, SELF, im / value);
1129     }
1130
1131 /*
1132
1133 =item C<PMC *neg(PMC *dest)>
1134
1135 =item C<void neg()>
1136
1137 Set C<dest> to the negated value of C<SELF>.
1138
1139 =cut
1140
1141 */
1142
1143     VTABLE PMC *neg(PMC *dest) {
1144         FLOATVAL re, im;
1145         dest = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1146
1147         GET_ATTR_re(INTERP, SELF, re);
1148         GET_ATTR_im(INTERP, SELF, im);
1149
1150         SET_ATTR_re(INTERP, dest, -re);
1151         SET_ATTR_im(INTERP, dest, -im);
1152
1153         return dest;
1154     }
1155
1156     VTABLE void i_neg() {
1157         FLOATVAL re, im;
1158         GET_ATTR_re(INTERP, SELF, re);
1159         GET_ATTR_im(INTERP, SELF, im);
1160         SET_ATTR_re(INTERP, SELF, -re);
1161         SET_ATTR_im(INTERP, SELF, -im);
1162     }
1163
1164 /*
1165
1166 =item C<INTVAL is_equal(PMC *value)>
1167
1168 Compares the complex number with C<value> and returns true if they are equal.
1169
1170 =cut
1171
1172 */
1173
1174     MULTI INTVAL is_equal(Complex value) {
1175         FLOATVAL self_re, self_im, val_re, val_im;
1176         GET_ATTR_re(INTERP, SELF, self_re);
1177         GET_ATTR_im(INTERP, SELF, self_im);
1178         GET_ATTR_re(INTERP, value, val_re);
1179         GET_ATTR_im(INTERP, value, val_im);
1180         return (INTVAL)(self_re == val_re && self_im == val_im);
1181     }
1182
1183     MULTI INTVAL is_equal(DEFAULT value) {
1184         FLOATVAL re, im;
1185         GET_ATTR_re(INTERP, SELF, re);
1186         GET_ATTR_im(INTERP, SELF, im);
1187
1188         if (im != 0.0)
1189             return (INTVAL)0;
1190
1191         return (re == VTABLE_get_number(INTERP, value));
1192     }
1193
1194 /*
1195
1196 =item C<PMC *absolute(PMC *dest)>
1197
1198 =item C<void i_absolute()>
1199
1200 Sets C<dest> to the absolute value of SELF that is the distance from (0.0).
1201
1202 =cut
1203
1204 */
1205
1206 /*
1207
1208   TODO for better precision: hinted by vaxman according to "Numerical Recipes
1209   in Fortran 77", 2nd edition, Press, Vetterling, Teukolsky, Flannery,
1210   Cambridge University Press, 2001, pp. 171ff:
1211
1212
1213 |a+ib|=|a|*sqrt(1+(b/a)**2), if |a|>=|b|,
1214        |b|*sqrt(1+(a/b)**2)  else.
1215
1216 */
1217
1218     VTABLE PMC *absolute(PMC *dest) {
1219         FLOATVAL re, im, d;
1220         GET_ATTR_re(INTERP, SELF, re);
1221         GET_ATTR_im(INTERP, SELF, im);
1222         d = sqrt(re*re + im*im);
1223
1224         dest = pmc_new(INTERP,
1225             Parrot_get_ctx_HLL_type(INTERP, enum_class_Float));
1226
1227         VTABLE_set_number_native(INTERP, dest, d);
1228         return dest;
1229     }
1230
1231     VTABLE void i_absolute() {
1232         FLOATVAL re, im, d;
1233         GET_ATTR_re(INTERP, SELF, re);
1234         GET_ATTR_im(INTERP, SELF, im);
1235         d = sqrt(re*re + im*im);
1236         pmc_reuse(INTERP, SELF, enum_class_Float, 0);
1237         VTABLE_set_number_native(INTERP, SELF, d);
1238     }
1239
1240 /*
1241
1242 =item C<METHOD ln()>
1243
1244 Returns the natural logarithm of SELF as a PMC.
1245
1246 =cut
1247
1248 ln z = ln |z| + i arg(z)
1249 |x + iy| = sqrt(x^2 + y^2)
1250 arg(x + iy) = atan2(y, x)
1251
1252 Some special cases
1253 ln(-1) = pi i
1254 ln(0) = -inf
1255 ln(1) = 0
1256 ln(e) = 1
1257 ln(+-i) = +- (pi i)/2
1258
1259 */
1260
1261     METHOD ln() {
1262         PMC * const d  = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1263         FLOATVAL re, im, result_re, result_im;
1264         GET_ATTR_re(INTERP, SELF, re);
1265         GET_ATTR_im(INTERP, SELF, im);
1266
1267         /* This is necessary for atan2 to behave */
1268         if (im == -0.0)
1269             im = 0.0;
1270
1271         result_re = log(sqrt(re*re + im*im));
1272         if (re == 0.0 && im == 0.0) /* atan2(0, 0) not portable */
1273             result_im = 0.0;
1274         else
1275             result_im = atan2(im, re);
1276
1277         SET_ATTR_re(INTERP, d, result_re);
1278         SET_ATTR_im(INTERP, d, result_im);
1279
1280         RETURN(PMC *d);
1281     }
1282
1283 /*
1284
1285 =item C<METHOD exp()>
1286
1287 Returns e ^ SELF as a PMC.
1288
1289 =cut
1290
1291 exp(a + bi) = exp(a) * (cos(b) + i * sin(b))
1292
1293 */
1294
1295     METHOD exp() {
1296         PMC * const d  = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1297         FLOATVAL re, im, f;
1298         GET_ATTR_re(INTERP, SELF, re);
1299         GET_ATTR_im(INTERP, SELF, im);
1300
1301         f  = exp(re);
1302
1303         SET_ATTR_re(INTERP, d, f * cos(im));
1304
1305         /* If only sin(pi) worked. */
1306         if (im == 4.0 * atan(1.0)) {
1307             SET_ATTR_im(INTERP, d, 0.0);
1308         }
1309         else {
1310             SET_ATTR_im(INTERP, d, f * sin(im));
1311         }
1312
1313         RETURN(PMC *d);
1314     }
1315
1316 /*
1317
1318 =item C<METHOD PMC *sin()>
1319
1320 =item C<METHOD PMC *cos()>
1321
1322 =item C<METHOD PMC *tan()>
1323
1324 =item C<METHOD PMC *csc()>
1325
1326 =item C<METHOD PMC *sec()>
1327
1328 =item C<METHOD PMC *cot()>
1329
1330 Returns C<FUNC>(SELF).
1331
1332 =cut
1333
1334  => sin(a + bi) = sin(a)cosh(b)+i*cos(a)sinh(b)
1335     sin(z) = ((e ^ zi) - (e ^ -zi)) / (2i)
1336  => cos(a + bi) = cos(a) * cosh(b) - i * sin(a) * sinh(b)
1337     cos(z) = ((e ^ zi) + (e ^ -zi)) / 2
1338
1339     sin(iz) = i sinh(z)
1340     cos(iz) = cosh(z)
1341
1342     sinh(iz) = i sin(z)
1343     cosh(iz) = cos z
1344     sinh(a + bi) = sinh(a) * cos(b) + i * cosh(a) * sin(b)
1345     cosh(a + bi) = cosh(a) * cos(b) + i * sinh(a) * sin(b)
1346
1347 */
1348
1349     METHOD sin() {
1350         PMC * const d = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1351         FLOATVAL re, im, result_re, result_im;
1352         GET_ATTR_re(INTERP, SELF, re);
1353         GET_ATTR_im(INTERP, SELF, im);
1354
1355         if (FLOAT_IS_ZERO(im)) {
1356             result_re = sin(re);
1357             result_im = 0.0;
1358         }
1359         else if (FLOAT_IS_ZERO(re)) {
1360             result_re = 0.0;
1361             result_im = sinh(im);
1362         }
1363         else {
1364             result_re = sin(re) * cosh(im);
1365
1366             if (im == -0.0)
1367                 result_im = 0.0;
1368             else
1369                 result_im = cos(re) * sinh(im);
1370
1371         }
1372
1373         SET_ATTR_re(INTERP, d, result_re);
1374         SET_ATTR_im(INTERP, d, result_im);
1375
1376         RETURN(PMC *d);
1377     }
1378
1379     METHOD cos() {
1380         PMC * const d  = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1381         FLOATVAL re, im, result_re, result_im;
1382         GET_ATTR_re(INTERP, SELF, re);
1383         GET_ATTR_im(INTERP, SELF, im);
1384
1385         if (FLOAT_IS_ZERO(re)) {
1386             result_re = cosh(im);
1387             result_im = 0.0;
1388         }
1389         else if (FLOAT_IS_ZERO(im)) {
1390             result_re = cos(re);
1391             result_im = 0.0;
1392         }
1393         else {
1394             result_re = cos(re) * cosh(im);
1395             result_im = -1.0 * sin(re) * sinh(im);
1396         }
1397
1398         SET_ATTR_re(INTERP, d, result_re);
1399         SET_ATTR_im(INTERP, d, result_im);
1400
1401         RETURN(PMC *d);
1402     }
1403
1404     METHOD tan() {
1405         PMC * const d = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1406         PMC * const e = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1407
1408         (PMC *d) = PCCINVOKE(INTERP, SELF, "sin");
1409         (PMC *e) = PCCINVOKE(INTERP, SELF, "cos");
1410
1411         Parrot_Complex_multi_i_divide_Complex(INTERP, d, e);
1412
1413         RETURN(PMC *d);
1414     }
1415
1416     METHOD cot() {
1417         PMC * const d = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1418         PMC * const e = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1419         (PMC *d) = PCCINVOKE(INTERP, SELF, "cos");
1420         (PMC *e) = PCCINVOKE(INTERP, SELF, "sin");
1421
1422         Parrot_Complex_multi_i_divide_Complex(INTERP, d, e);
1423
1424         RETURN(PMC *d);
1425     }
1426
1427     METHOD sec() {
1428         PMC * const d = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1429         PMC * const e = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1430         SET_ATTR_re(INTERP, d, 1.0);
1431         SET_ATTR_im(INTERP, d, 0.0);
1432         (PMC *e) = PCCINVOKE(INTERP, SELF, "cos");
1433
1434         Parrot_Complex_multi_i_divide_Complex(INTERP, d, e);
1435
1436         RETURN(PMC *d);
1437     }
1438
1439     METHOD csc() {
1440         PMC * const d = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1441         PMC * const e = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1442
1443         SET_ATTR_re(INTERP, d, 1.0);
1444         SET_ATTR_im(INTERP, d, 0.0);
1445
1446         (PMC *e) = PCCINVOKE(INTERP, SELF, "sin");
1447
1448         Parrot_Complex_multi_i_divide_Complex(INTERP, d, e);
1449
1450         RETURN(PMC *d);
1451     }
1452
1453 /*
1454
1455 =item C<METHOD PMC *asin()>
1456
1457 =item C<METHOD PMC *acos()>
1458
1459 =item C<METHOD PMC *atan()>
1460
1461 =item C<METHOD PMC *acsc()>
1462
1463 =item C<METHOD PMC *asec()>
1464
1465 =item C<METHOD PMC *acot()>
1466
1467 Returns the inverse function of SELF.
1468
1469 =cut
1470
1471  => arcsin z = -i ln(iz + sqrt(1-z*z))
1472  => arccos z = pi/2 + i * ln(iz + sqrt(1 - z*z))
1473     arccos z = -i ln(z + sqrt(z*z-1))
1474  => arctan z = i/2 ln((i+z) / (i-z))
1475     arctan z = 1/2 i (ln(1-iz) - ln(1 + iz))
1476
1477  => acot(z) = atan(1 / z)
1478     acot(z) = i/2 (ln((z - i) / z) - ln((z + i) / z))
1479  => asec(z) = acos(1 / z)
1480     asec(z) = 1/2 pi + i ln(sqrt(1 - 1/zz) + i/z)
1481  => acsc(z) = asin(1 / z)
1482     acsc(z) = -i ln(sqrt(1 - 1/zz + i/z))
1483
1484 */
1485
1486     METHOD asin() {
1487         FLOATVAL d_re, d_im, e_re, e_im, self_re, self_im;
1488         PMC * const d = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1489         PMC *       e = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1490         GET_ATTR_re(INTERP, SELF, self_re);
1491         GET_ATTR_im(INTERP, SELF, self_im);
1492
1493         e      = Parrot_Complex_multi_multiply_Complex_PMC(INTERP, SELF, SELF, e);
1494         GET_ATTR_re(INTERP, e, e_re);
1495         GET_ATTR_im(INTERP, e, e_im);
1496         SET_ATTR_re(INTERP, e, 1.0 - e_re);
1497         SET_ATTR_im(INTERP, e, -e_im);
1498
1499         (PMC *d) = PCCINVOKE(INTERP, e, "sqrt");
1500         GET_ATTR_re(INTERP, d, d_re);
1501         GET_ATTR_im(INTERP, d, d_im);
1502         SET_ATTR_re(INTERP, d, d_re - self_im);
1503         SET_ATTR_im(INTERP, d, d_im + self_re);
1504
1505         (PMC *d) = PCCINVOKE(INTERP, d, "ln");
1506         GET_ATTR_re(INTERP, d, d_re);
1507         GET_ATTR_im(INTERP, d, d_im);
1508         SET_ATTR_re(INTERP, e, d_im);
1509         SET_ATTR_im(INTERP, e, d_re ? -d_re : 0.0);
1510
1511         RETURN(PMC *e);
1512     }
1513
1514     METHOD acos() {
1515         FLOATVAL d_re, d_im, e_re, e_im, self_re, self_im;
1516         PMC * const d = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1517         PMC *       e = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1518         GET_ATTR_re(INTERP, SELF, self_re);
1519         GET_ATTR_im(INTERP, SELF, self_im);
1520
1521         e = Parrot_Complex_multi_multiply_Complex_PMC(INTERP, SELF, SELF, e);
1522         GET_ATTR_re(INTERP, e, e_re);
1523         GET_ATTR_im(INTERP, e, e_im);
1524         SET_ATTR_re(INTERP, e, 1.0 - e_re);
1525         SET_ATTR_im(INTERP, e, -e_im);
1526
1527         (PMC *d) = PCCINVOKE(INTERP, e, "sqrt");
1528         GET_ATTR_re(INTERP, d, d_re);
1529         GET_ATTR_im(INTERP, d, d_im);
1530         SET_ATTR_re(INTERP, d, d_re + self_im);
1531         SET_ATTR_im(INTERP, d, d_im - self_re);
1532
1533         (PMC *e) = PCCINVOKE(INTERP, d, "ln");
1534         GET_ATTR_re(INTERP, e, e_re);
1535         GET_ATTR_im(INTERP, e, e_im);
1536         SET_ATTR_re(INTERP, d, e_im + 2.0 * atan(1.0));
1537         SET_ATTR_im(INTERP, d, e_re ? -e_re : 0.0);
1538
1539         RETURN(PMC *d);
1540     }
1541
1542     METHOD atan() {
1543         PMC * const d  = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1544         PMC * const e  = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1545         FLOATVAL re, im, d_re, d_im;
1546         GET_ATTR_re(INTERP, SELF, re);
1547         GET_ATTR_im(INTERP, SELF, im);
1548
1549         SET_ATTR_re(INTERP, d, re);
1550         SET_ATTR_im(INTERP, d, 1 + im);
1551         SET_ATTR_re(INTERP, e, -re);
1552         SET_ATTR_im(INTERP, e, 1 - im);
1553
1554         Parrot_Complex_multi_i_divide_Complex(INTERP, d, e);
1555
1556         (PMC *d) = PCCINVOKE(INTERP, d, "ln");
1557         GET_ATTR_re(INTERP, d, d_re);
1558         GET_ATTR_im(INTERP, d, d_im);
1559
1560         SET_ATTR_re(INTERP, e, (d_im ? d_im : -0.0) / -2.0);
1561         SET_ATTR_im(INTERP, e, d_re / 2.0);
1562
1563         RETURN(PMC *e);
1564     }
1565
1566     METHOD acot() {
1567         PMC * const d = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1568         PMC     *e;
1569         FLOATVAL re, im;
1570         GET_ATTR_re(INTERP, SELF, re);
1571         GET_ATTR_im(INTERP, SELF, im);
1572
1573         SET_ATTR_re(INTERP, d,  re / (re * re + im * im));
1574         SET_ATTR_im(INTERP, d, -im / (re * re + im * im));
1575
1576         (PMC *e) = PCCINVOKE(INTERP, d, "atan");
1577         RETURN(PMC *e);
1578     }
1579
1580     METHOD acsc() {
1581         PMC * const d = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1582         PMC     *e;
1583         FLOATVAL re, im;
1584         GET_ATTR_re(INTERP, SELF, re);
1585         GET_ATTR_im(INTERP, SELF, im);
1586
1587         SET_ATTR_re(INTERP, d,  re / (re * re + im * im));
1588         SET_ATTR_im(INTERP, d, -im / (re * re + im * im));
1589
1590         (PMC *e) = PCCINVOKE(INTERP, d, "asin");
1591         RETURN(PMC *e);
1592     }
1593
1594     METHOD asec() {
1595         PMC * const d = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1596         PMC     *e;
1597         FLOATVAL re, im;
1598         GET_ATTR_re(INTERP, SELF, re);
1599         GET_ATTR_im(INTERP, SELF, im);
1600
1601         SET_ATTR_re(INTERP, d,  re / (re * re + im * im));
1602         SET_ATTR_im(INTERP, d, -im / (re * re + im * im));
1603
1604         (PMC *e) = PCCINVOKE(INTERP, d, "acos");
1605         RETURN(PMC *e);
1606     }
1607
1608 /*
1609
1610 =item C<METHOD PMC *sinh()>
1611
1612 Returns the arctangent of SELF.
1613
1614 =item C<METHOD PMC *cosh()>
1615
1616 Returns the arcsine of SELF.
1617
1618 =item C<METHOD PMC *tanh()>
1619
1620 Returns the arccosine of SELF.
1621
1622 =cut
1623
1624 tanh(z) = sinh(z) / cosh(z)
1625
1626 */
1627
1628     METHOD sinh() {
1629         PMC * const d = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1630         FLOATVAL re, im;
1631         GET_ATTR_re(INTERP, SELF, re);
1632         GET_ATTR_im(INTERP, SELF, im);
1633
1634         SET_ATTR_re(INTERP, d, sinh(re) * cos(im));
1635         SET_ATTR_im(INTERP, d, im ? cosh(re) * sin(im) : 0.0);
1636
1637         RETURN(PMC *d);
1638     }
1639
1640     METHOD cosh() {
1641         PMC * const d = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1642         FLOATVAL re, im;
1643         GET_ATTR_re(INTERP, SELF, re);
1644         GET_ATTR_im(INTERP, SELF, im);
1645
1646         SET_ATTR_re(INTERP, d, cosh(re) * cos(im));
1647         if (re == 0.0 || im == 0.0) {
1648             SET_ATTR_im(INTERP, d, 0.0);
1649         }
1650         else {
1651             SET_ATTR_im(INTERP, d, sinh(re) * sin(im));
1652         }
1653
1654         RETURN(PMC *d);
1655     }
1656
1657     METHOD tanh() {
1658         PMC * const d = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1659         PMC * const e = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1660
1661         (PMC *d) = PCCINVOKE(INTERP, SELF, "sinh");
1662         (PMC *e) = PCCINVOKE(INTERP, SELF, "cosh");
1663
1664         Parrot_Complex_multi_i_divide_Complex(INTERP, d, e);
1665
1666         RETURN(PMC *d);
1667     }
1668
1669     METHOD coth() {
1670         FLOATVAL re, im;
1671         PMC * const d = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1672
1673         (PMC *d) = PCCINVOKE(INTERP, SELF, "tanh");
1674         GET_ATTR_re(INTERP, d, re);
1675         GET_ATTR_im(INTERP, d, im);
1676
1677         SET_ATTR_re(INTERP, d, re ?  re / (re * re + im * im) : 0.0);
1678         SET_ATTR_im(INTERP, d, im ? -im / (re * re + im * im) : 0.0);
1679
1680         RETURN(PMC *d);
1681     }
1682
1683     METHOD csch() {
1684         FLOATVAL re, im;
1685         PMC * const d = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1686
1687         (PMC *d) = PCCINVOKE(INTERP, SELF, "sinh");
1688         GET_ATTR_re(INTERP, d, re);
1689         GET_ATTR_im(INTERP, d, im);
1690
1691         SET_ATTR_re(INTERP, d, re ?  re / (re * re + im * im) : 0.0);
1692         SET_ATTR_im(INTERP, d, im ? -im / (re * re + im * im) : 0.0);
1693
1694         RETURN(PMC *d);
1695     }
1696
1697     METHOD sech() {
1698         FLOATVAL re, im;
1699         PMC * const d = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1700
1701         (PMC *d) = PCCINVOKE(INTERP, SELF, "cosh");
1702         GET_ATTR_re(INTERP, d, re);
1703         GET_ATTR_im(INTERP, d, im);
1704
1705         SET_ATTR_re(INTERP, d, re ?  re / (re * re + im * im) : 0.0);
1706         SET_ATTR_im(INTERP, d, im ? -im / (re * re + im * im) : 0.0);
1707
1708         RETURN(PMC *d);
1709     }
1710
1711 /*
1712
1713 =item C<METHOD PMC *asinh()>
1714
1715 =item C<METHOD PMC *acosh()>
1716
1717 =item C<METHOD PMC *atanh()>
1718
1719 =item C<METHOD PMC *acsch()>
1720
1721 =item C<METHOD PMC *asech()>
1722
1723 =item C<METHOD PMC *acoth()>
1724
1725 The inverse hyperbolic functions.  Currently all broken, but for
1726 C<func(a+bi) = c+di>, C<|c|> and C<|d|> will be correct, confusingly enough.
1727
1728 =cut
1729
1730 asinh z = -ln(sqrt(1+zz) - z)
1731 asinh z = ln(sqrt(zz + 1) + z)
1732
1733 asinh = i asin(-ix)
1734 acosh = i acos(x)
1735 atanh = i atan(-ix)
1736
1737 */
1738
1739     METHOD asinh() {
1740         FLOATVAL re, im;
1741         PMC * const d = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1742         PMC * const e = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1743
1744         GET_ATTR_re(INTERP, SELF, re);
1745         GET_ATTR_im(INTERP, SELF, im);
1746         SET_ATTR_re(INTERP, d,  im);
1747         SET_ATTR_im(INTERP, d, -re);
1748
1749         (PMC *d) = PCCINVOKE(INTERP, d, "asin");
1750         GET_ATTR_re(INTERP, d, re);
1751         GET_ATTR_im(INTERP, d, im);
1752         SET_ATTR_re(INTERP, e, -im);
1753         SET_ATTR_im(INTERP, e,  re);
1754
1755         RETURN(PMC *e);
1756     }
1757
1758     METHOD acosh() {
1759         FLOATVAL re, im;
1760         PMC * const d = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1761         PMC * const e = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1762
1763         (PMC *d) = PCCINVOKE(INTERP, SELF, "acos");
1764         GET_ATTR_re(INTERP, d, re);
1765         GET_ATTR_im(INTERP, d, im);
1766         SET_ATTR_re(INTERP, e, -im);
1767         SET_ATTR_im(INTERP, e,  re);
1768
1769         RETURN(PMC *e);
1770     }
1771
1772     METHOD atanh() {
1773         FLOATVAL re, im;
1774         PMC * const d = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1775         PMC * const e = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1776
1777         GET_ATTR_re(INTERP, SELF, re);
1778         GET_ATTR_im(INTERP, SELF, im);
1779         SET_ATTR_re(INTERP, d,  im);
1780         SET_ATTR_im(INTERP, d, -re);
1781
1782         (PMC *d) = PCCINVOKE(INTERP, d, "atan");
1783         GET_ATTR_re(INTERP, d, re);
1784         GET_ATTR_im(INTERP, d, im);
1785         SET_ATTR_re(INTERP, e, -im);
1786         SET_ATTR_im(INTERP, e,  re);
1787
1788         RETURN(PMC *e);
1789     }
1790
1791     METHOD acoth() {
1792         PMC * const d = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1793         PMC     *e;
1794         FLOATVAL re, im;
1795         GET_ATTR_re(INTERP, SELF, re);
1796         GET_ATTR_im(INTERP, SELF, im);
1797
1798         SET_ATTR_re(INTERP, d,  re / (re * re + im * im));
1799         SET_ATTR_im(INTERP, d, -im / (re * re + im * im));
1800
1801         (PMC *e) = PCCINVOKE(INTERP, d, "atanh");
1802         RETURN(PMC *e);
1803     }
1804
1805     METHOD acsch() {
1806         PMC * const d = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1807         PMC     *e;
1808         FLOATVAL re, im;
1809         GET_ATTR_re(INTERP, SELF, re);
1810         GET_ATTR_im(INTERP, SELF, im);
1811
1812         SET_ATTR_re(INTERP, d,  re / (re * re + im * im));
1813         SET_ATTR_im(INTERP, d, -im / (re * re + im * im));
1814
1815         (PMC *e) = PCCINVOKE(INTERP, d, "asinh");
1816         RETURN(PMC *e);
1817     }
1818
1819     METHOD asech() {
1820         PMC * const d = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1821         PMC     *e;
1822         FLOATVAL re, im;
1823         GET_ATTR_re(INTERP, SELF, re);
1824         GET_ATTR_im(INTERP, SELF, im);
1825
1826         SET_ATTR_re(INTERP, d,  re / (re * re + im * im));
1827         SET_ATTR_im(INTERP, d, -im / (re * re + im * im));
1828
1829         (PMC *e) = PCCINVOKE(INTERP, d, "acosh");
1830         RETURN(PMC *e);
1831     }
1832
1833 /*
1834
1835 =item C<PMC *pow(PMC *value, PMC *dest)>
1836
1837 Return SELF to the C<value>th power and return result in C<dest>.
1838
1839 =item C<METHOD PMC *sqrt()>
1840
1841 Return the square root of SELF.
1842
1843 =cut
1844
1845 TODO: mmd in other pmc's to allow .Integer ^ .Complex, etc.
1846 and i_pow, and pow_(float|int), etc
1847
1848 x ^ y = exp(y * ln x))
1849
1850 */
1851
1852     MULTI PMC *pow(Complex value, PMC *dest) {
1853         PMC *l = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1854         PMC *log;
1855
1856         dest = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1857
1858
1859         Parrot_PCCINVOKE(interp, SELF, CONST_STRING(interp, "ln"),
1860                     "->P", &log);
1861
1862         l = Parrot_Complex_multi_multiply_Complex_PMC(INTERP,
1863                     log, value, l);
1864
1865         Parrot_PCCINVOKE(interp, l, CONST_STRING(interp, "exp"),
1866                     "->P", &dest);
1867
1868         return dest;
1869     }
1870
1871     MULTI PMC *pow(DEFAULT value, PMC *dest) {
1872         PMC *l = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1873         PMC *log;
1874
1875         dest = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1876
1877         Parrot_PCCINVOKE(interp, SELF, CONST_STRING(interp, "ln"),
1878                      "->P", &log);
1879         l = Parrot_Complex_multi_multiply_DEFAULT_PMC(INTERP, log, value, l);
1880         Parrot_PCCINVOKE(interp, l, CONST_STRING(interp, "exp"),
1881                      "->P", &dest);
1882
1883         return dest;
1884     }
1885
1886     METHOD sqrt() {
1887         PMC * const result = pmc_new(INTERP, VTABLE_type(INTERP, SELF));
1888         const FLOATVAL absval = SELF.get_number();
1889         FLOATVAL sx, sy, rx, ry;
1890         GET_ATTR_re(INTERP, SELF, sx);
1891         GET_ATTR_im(INTERP, SELF, sy);
1892
1893         rx = sqrt((absval + sx) / 2);
1894         ry = sqrt((absval - sx) / 2);
1895         if (sy < 0)
1896             ry = -ry;
1897         SET_ATTR_re(INTERP, result, rx);
1898         SET_ATTR_im(INTERP, result, ry);
1899         RETURN(PMC *result);
1900     }
1901
1902 }
1903
1904 /*
1905
1906 =back
1907
1908 =cut
1909
1910 */
1911
1912 /*
1913  * Local variables:
1914  *   c-file-style: "parrot"
1915  * End:
1916  * vim: expandtab shiftwidth=4:
1917  */