complex.pmc

/*
Copyright (C) 2004-2007, The Perl Foundation.
$Id$

=head1 NAME

src/pmc/complex.pmc - Complex Numbers PMC Class

=head1 DESCRIPTION

C<Complex> provides a representation of complex numbers. It handles
string parsing/generating and basic mathematical operations.

=head2 Functions

=over 4

=cut

Equations used are sometimes listed.  At times, multiple equations are given,
but those starting with => are the ones used

*/

#include "parrot/parrot.h"


#define RE(obj) (((FLOATVAL*)PMC_struct_val(obj))[0])
#define IM(obj) (((FLOATVAL*)PMC_struct_val(obj))[1])


/*

=item C<static void
complex_parse_string(Interp* interp,
        FLOATVAL *re, FLOATVAL *im, STRING *value)>

Parses the string in C<value> to produce a complex number, represented
by the real (C<*re>) and imaginary (C<*im>) parts. Raises an exception
if it cannot understand the string.  The string should be of the form
C<a+bi> with optional spaces around C<+> and before C<i>. You can also
use C<j> instead of C<i>.

=cut

*/

static void
complex_parse_string(Interp* interp,
        FLOATVAL *re, FLOATVAL *im, STRING *value) {
    char *str;
    INTVAL first_num_length, second_num_length;
    INTVAL first_num_minus, second_num_minus, i;
    char *first_num_offset, *second_num_offset = NULL, *t;
    STRING *S;

    t = str = string_to_cstring(interp, value);
    i = 0;
    first_num_offset = str;
    first_num_minus = second_num_minus = 0;

    /* walk the string and identify the real and imaginary parts */

    if (*t == '-') {
        /* first number is negative */
        t++;
        first_num_minus = 1;
        if (*t == ' ') t++; /* allow for an optional space */
        first_num_offset = t;
    }
    while (*t >= '0' && *t <= '9') t++; /* skip digits */
    if (*t == '.') {
        /* this number has a decimal point */
        t++;
        while (*t >= '0' && *t <= '9') t++; /* skip digits */
    }
    /* save the length of the real part */
    first_num_length = t - first_num_offset;

    if (*t == 0) {
        /* end of string; we only have a real part */
        second_num_length = 0;
    }
    else if ((*t == 'i' || *t == 'j') && *(t+1) == 0) {
        /* there is only an imaginary part, so the first number was
            actually the imaginary part */
        second_num_length = first_num_length;
        first_num_length = 0;
        second_num_offset = first_num_offset;
        second_num_minus = first_num_minus;
        first_num_minus = 0;
        /* this is useful if there is no number for
            the imaginary part, like in "-i" */
        i = 1;
    }
    else {
        if (*t == ' ') t++; /* skip an optional space */
        /* expect "+" or "-" and the imaginary part */
        if (*t == '+' || *t == '-') {
            second_num_minus = (*t == '-'); /* save the sign */
            t++;
            if (*t == ' ') t++; /* skip another optional space */

            /* save the beginning of the imaginary part */
            second_num_offset = t;
            while (*t >= '0' && *t <= '9') t++; /* skip digits */
            if (*t == '.') {
                /* this number has a decimal point */
                t++;
                while (*t >= '0' && *t <= '9') t++; /* skip digits */
            }
            /* save the length of the imaginary part */
            second_num_length = t - second_num_offset;

            if (*t == ' ') t++; /* allow for one more optional space */

            /* verify that the string ends properly */
            if ((*t != 'i' && *t != 'j') || (*(t+1) != 0)) {
                /* imaginary part does not end in 'i' or 'j' */
                real_exception(interp, NULL, E_ValueError,
                    "Complex: malformed string");
                first_num_length = second_num_length = 0;
            }
            /* this is useful if there is no number for the
                imaginary part, like in "2+i" */
            i = 1;

            /* all is OK, save the number */
        }
        else {
            /* "+" or "-" not found: error */
            real_exception(interp, NULL, E_ValueError,
                "Complex: malformed string");
            first_num_length = second_num_length = 0;
        }
    }

    /* now we have the offsets and the lengths;
        we turn them into float values */

    if (first_num_length) {
        /* there is a real part, interpret it */
        S = string_from_cstring(interp,
                first_num_offset, first_num_length);
        *re = string_to_num(interp, S);
    }
    else {
        /* consider the real part 0.0 */
        *re = 0.0;
    }

    if (second_num_length) {
        /* there is an imaginary part, interpret it */
        S = string_from_cstring(interp,
                second_num_offset, second_num_length);
        *im = string_to_num(interp, S);
    }
    else {
        /* consider the imaginary part 0.0 */
        if (i) /* the string was something like "1+i" */
            *im = 1.0;
        else
            *im = 0.0;
    }
    if (first_num_minus) *re = - *re;
    if (second_num_minus) *im = - *im;

    string_cstring_free(str);
}

/*

=item C<static FLOATVAL*
complex_locate_keyed_num(Interp* interp, PMC* self, STRING *key)>

Interpret the string C<key>; valid keys are C<real> and C<imag>,
representing the real and imaginary parts of the complex number.

=cut

*/

static FLOATVAL*
complex_locate_keyed_num(Interp* interp, PMC* self, STRING *key) {
    /* do imag first since real can be read much faster anyway */
    STRING *imag = string_from_cstring(interp, "imag", 4);
    STRING *real;

    if (0 == string_equal(interp, key, imag))
        return &IM(self);
    real = string_from_cstring(interp, "real", 4);
    if (0 == string_equal(interp, key, real))
        return &RE(self);
    real_exception(interp, NULL, E_KeyError,
        "Complex: key is neither 'real' or 'imag'");
    return NULL;
}

static void
int_check_divide_zero(PARROT_INTERP, INTVAL value)
{
    if (value == 0) {
        real_exception(interp, NULL, E_ZeroDivisionError,
                       "Divide by zero");
    }
}

static void
float_check_divide_zero(PARROT_INTERP, FLOATVAL value)
{
    if (value == 0.0) {
        real_exception(interp, NULL, E_ZeroDivisionError,
                       "Divide by zero");
    }
}

static void
complex_check_divide_zero(PARROT_INTERP, PMC *value)
{
    /* Throw an exception if we are dividing by zero. */
    if (RE(value) == 0.0 && IM(value) == 0.0) {
        real_exception(interp, NULL, E_ZeroDivisionError,
                       "Divide by zero");
    }
}


pmclass Complex {

/*

=item C<PMC* instantiate(PMC* sig)>

Create a new complex PMC with passed arguments according to pdd03.

=cut

*/

    PMC* instantiate(PMC* sig) {
        PMC *res = pmc_new(INTERP, enum_class_Complex);
        FLOATVAL re = 0.0, im = 0.0;
        int argcI = REG_INT(interp, 1);
        int argcS = REG_INT(interp, 2);
        int argcP = REG_INT(interp, 3);
        int argcN = REG_INT(interp, 4);

        return PMCNULL;  /* TODO */
        /*
         * we can only allow 0..2 arguments of one kind. For
         * mixed (e.g. N,P) args the order of arguments isn't fixed
         */
        if (argcP) {
            re = VTABLE_get_number(INTERP, REG_PMC(interp, 5));
            if (argcP == 2)
                im = VTABLE_get_number(INTERP, REG_PMC(interp, 6));
            /*
             * TODO throw exception if argument mismatch
             */
        }
        else if (argcN) {
            re = REG_NUM(interp, 5);
            if (argcN == 2)
                im = REG_NUM(interp, 6);
        }
        else if (argcI) {
            re = REG_INT(interp, 5);
            if (argcI == 2)
                im = REG_INT(interp, 6);
        }
        else if (argcS == 1) {
            complex_parse_string(INTERP, &RE(res), &IM(res), REG_STR(interp, 5));
            return res;
        }
        RE(res) = re;
        IM(res) = im;
        return res;
    }

/*

=item C<opcode_t* invoke(void* next)>

Pythonic object constructor. SELF is a Complex Class object. Return a new
C<complex> object according to 2.1. Built-in Functions.

=cut

*/
    opcode_t* invoke(void *next) {
        const int argcP = REG_INT(interp, 3);
        PMC * const res  = pmc_new(INTERP, enum_class_Complex);

        if (argcP == 1) {
            PMC * const arg = REG_PMC(interp, 5);

            if (arg->vtable->base_type == enum_class_String)
                VTABLE_set_string_native(INTERP, res, PMC_str_val(arg));
            else
                RE(res) = VTABLE_get_number(INTERP, arg);
        }
        else if (argcP == 2) {
            RE(res) = VTABLE_get_number(INTERP, REG_PMC(interp, 5));
            IM(res) = VTABLE_get_number(INTERP, REG_PMC(interp, 6));
        }

        REG_PMC(interp, 5) = res;
        return (opcode_t *)next;
    }

/*

=back

=head2 Methods

=over 4

=item C<void init()>

Initializes the complex number with the value 0+0i.

=item C<void init_pmc(PMC* initializer)>

Initializes the complex number with the specified values.
(not implemented)

=item C<void destroy()>

Cleans up.

=item C<PMC* clone()>

Creates an identical copy of the complex number.

=cut

*/

    void init() {
        /* XXX should check if mem_sys_allocate failed */
        FLOATVAL* data = (FLOATVAL*)mem_sys_allocate(2 * sizeof (FLOATVAL));
        PMC_struct_val(SELF) = data;
        PObj_active_destroy_SET(SELF);
        RE(SELF) = IM(SELF) = 0.0;
    }

    void init_pmc(PMC* initializer) {
        /* XXX not implemented */
        DYNSELF.init();
    }

    void destroy() {
        mem_sys_free(PMC_struct_val(SELF));
        PMC_struct_val(SELF) = NULL;
    }

    void morph(INTVAL type) {
        if (SELF->vtable->base_type == type)
            return;
        SUPER(type);
    }

    PMC* clone() {
        PMC* dest = pmc_new_noinit(INTERP, SELF->vtable->base_type);
        FLOATVAL* data = (FLOATVAL*)mem_sys_allocate(2 * sizeof (FLOATVAL));
        PMC_struct_val(dest) = data;
        PObj_active_destroy_SET(dest);
        RE(dest) = RE(SELF);
        IM(dest) = IM(SELF);
        return dest;
    }

/*

=item C<INTVAL get_integer()>

Returns the modulus of the complex number as an integer.

=item C<FLOATVAL get_number()>

Returns the modulus of the complex number.

=item C<STRING* get_string()>

Returns the complex number as a string in the form C<a+bi>.

=item C<INTVAL get_bool()>

Returns true if the complex number is non-zero.

=cut

*/

    INTVAL get_integer() {
        return (INTVAL)(DYNSELF.get_number());
    }

    FLOATVAL get_number() {
        return sqrt(RE(SELF)*RE(SELF) + IM(SELF)*IM(SELF));
    }

    STRING* get_string() {
        STRING *s;
            s = Parrot_sprintf_c(INTERP,
                    "%vg%+vgi", RE(SELF), IM(SELF));
        return s;
    }

    INTVAL get_bool() {
        return (INTVAL)(RE(SELF) != 0.0 || IM(SELF) != 0.0);
    }

/*

=item C<INTVAL get_integer_keyed(PMC* key)>

=item C<INTVAL get_integer_keyed_str(STRING* key)>

=item C<FLOATVAL get_number_keyed(PMC* key)>

=item C<FLOATVAL get_number_keyed_str(STRING* key)>

=item C<PMC* get_pmc_keyed(PMC* key)>

=item C<PMC* get_pmc_keyed_str(STRING* key)>

Returns the requested number (real part for C<real> and imaginary for C<imag>).

=cut

*/

    INTVAL get_integer_keyed(PMC* key) {
        STRING * const s = VTABLE_get_string(INTERP, key);
        return DYNSELF.get_integer_keyed_str(s);
    }

    INTVAL get_integer_keyed_str(STRING* key) {
        return (INTVAL)(DYNSELF.get_number_keyed_str(key));
    }

    FLOATVAL get_number_keyed(PMC* key) {
        STRING * const s = VTABLE_get_string(INTERP, key);
        return DYNSELF.get_number_keyed_str(s);
    }

    FLOATVAL get_number_keyed_str(STRING* key) {
        const FLOATVAL * const num = complex_locate_keyed_num(INTERP, SELF, key);
        if (num)
            return *num;
        return 0.0;
    }

    PMC* get_pmc_keyed(PMC* key) {
        STRING * const s = VTABLE_get_string(INTERP, key);
        return DYNSELF.get_pmc_keyed_str(s);
    }

    PMC* get_pmc_keyed_str(STRING* key) {
        PMC * const ret = pmc_new(INTERP, enum_class_Float);
        const FLOATVAL val = DYNSELF.get_number_keyed_str(key);
        VTABLE_set_number_native(INTERP, ret, val);
        return ret;
    }

/*

=item C<PMC* get_pmc_keyed_int(INTVAL key)>

Returns the requested number (real part for C<0> and imaginary for C<1>).

=cut

*/

    PMC* get_pmc_keyed_int(INTVAL key) {
        PMC * const ret = pmc_new(INTERP, enum_class_Float);
        const FLOATVAL val = DYNSELF.get_number_keyed_int(key);
        VTABLE_set_number_native(INTERP, ret, val);
        return ret;
    }

/*

=item C<FLOATVAL get_number_keyed_int(INTVAL key)>

Quick hack to emulate get_real() and get_imag():

  key = 0 ... get real part
  key = 1 ... get imag part

=item C<void set_number_keyed_int(INTVAL key, FLOATVAL v)>

Set real or imag depending on key

=cut

*/

    FLOATVAL get_number_keyed_int(INTVAL key) {
        switch (key) {
            case 0:
                return RE(SELF);
            case 1:
                return IM(SELF);
            default:
                real_exception(interp, NULL, E_KeyError,
                        "Complex: key must be 0 or 1");
        }
        return 0.0;
    }

    void set_number_keyed_int(INTVAL key, FLOATVAL v) {
        switch (key) {
            case 0:
                RE(SELF) = v;
                break;
            case 1:
                IM(SELF) = v;
                break;
            default:
                real_exception(interp, NULL, E_KeyError,
                        "Complex: key must be 0 or 1");
        }
    }
/*

=item C<void set_string_native(STRING* value)>

Parses the string C<value> into a complex number; raises an exception
on failure.

=item C<void set_pmc(PMC* value)>

if C<value> is a Complex PMC then the complex number is set to its
value; otherwise C<value>'s string representation is parsed with
C<set_string_native()>.

=item C<void set_integer_native(INTVAL value)>

=item C<void set_number_native(FLOATVAL value)>

Sets the real part of the complex number to C<value> and the imaginary
part to C<0.0>

=cut

*/

    void set_string_native(STRING* value) {
        complex_parse_string(INTERP, &RE(SELF), &IM(SELF), value);
    }

    void set_pmc(PMC* value) {
        if (value->vtable->base_type == enum_class_Complex) {
            RE(SELF) = RE(value);
            IM(SELF) = IM(value);
        }
        else {
            DYNSELF.set_string_native(VTABLE_get_string(INTERP, value));
        }
    }

    void set_integer_native(INTVAL value) {
        DYNSELF.set_number_native(value);
    }

    void set_number_native(FLOATVAL value) {
        RE(SELF) = value;
        IM(SELF) = 0.0;
    }

/*

=item C<void set_integer_keyed(PMC* key, INTVAL value)>

=item C<void set_integer_keyed_str(STRING* key, INTVAL value)>

=item C<void set_number_keyed(PMC* key, FLOATVAL value)>

=item C<void set_number_keyed_str(STRING* key, FLOATVAL value)>

=item C<void set_pmc_keyed(PMC* key, PMC* value)>

=item C<void set_pmc_keyed_str(STRING* key, PMC* value)>

Sets the requested number (real part for C<real> and imaginary for C<imag>)
to C<value>.

=cut

*/

    void set_integer_keyed(PMC* key, INTVAL value) {
        DYNSELF.set_number_keyed(key, value);
    }

    void set_integer_keyed_str(STRING* key, INTVAL value) {
        DYNSELF.set_number_keyed_str(key, value);
    }

    void set_number_keyed(PMC* key, FLOATVAL value) {
        STRING* s = VTABLE_get_string(INTERP, key);
        DYNSELF.set_number_keyed_str(s, value);
    }

    void set_number_keyed_str(STRING* key, FLOATVAL value) {
        FLOATVAL * const num = complex_locate_keyed_num(INTERP, SELF, key);
        if (num)
            *num = value;
    }

    void set_pmc_keyed(PMC* key, PMC* value) {
        const FLOATVAL f = VTABLE_get_number(INTERP, value);
        DYNSELF.set_number_keyed(key, f);
    }

    void set_pmc_keyed_str(STRING* key, PMC* value) {
        const FLOATVAL f = VTABLE_get_number(INTERP, value);
        DYNSELF.set_number_keyed_str(key, f);
    }

/*

=item C<PMC* add(PMC* value, PMC* dest)>

=item C<PMC* add_int(INTVAL value, PMC* dest)>

=item C<PMC* add_float(FLOATVAL value, PMC* dest)>

Adds C<value> to the complex number, placing the result in C<dest>.

=cut

*/

    PMC* add(PMC* value, PMC* dest) {
MMD_Complex: {
            if (dest)
                VTABLE_morph(INTERP, dest, SELF->vtable->base_type);
            else {
                dest = pmc_new(INTERP, SELF->vtable->base_type);
            }
            RE(dest) = RE(SELF) + RE(value);
            IM(dest) = IM(SELF) + IM(value);
            return dest;
        }
MMD_DEFAULT: {
            if (dest)
                VTABLE_morph(INTERP, dest, SELF->vtable->base_type);
            else {
                dest = pmc_new(INTERP, SELF->vtable->base_type);
            }
            RE(dest) = RE(SELF) + VTABLE_get_number(INTERP, value);
            IM(dest) = IM(SELF);
            return dest;
        }
    }

    PMC* add_int(INTVAL value, PMC* dest) {
        if (dest)
            VTABLE_morph(INTERP, dest, SELF->vtable->base_type);
        else {
            dest = pmc_new(INTERP, SELF->vtable->base_type);
        }
        RE(dest) = RE(SELF) + value;
        IM(dest) = IM(SELF);
        return dest;
    }

    PMC* add_float(FLOATVAL value, PMC* dest) {
        if (dest)
            VTABLE_morph(INTERP, dest, SELF->vtable->base_type);
        else {
            dest = pmc_new(INTERP, SELF->vtable->base_type);
        }
        VTABLE_morph(INTERP, dest, SELF->vtable->base_type);
        RE(dest) = RE(SELF) + value;
        IM(dest) = IM(SELF);
        return dest;
    }

    void i_add(PMC* value) {
MMD_Complex: {
            RE(SELF) += RE(value);
            IM(SELF) += IM(value);
        }
MMD_DEFAULT: {
            RE(SELF) += VTABLE_get_number(INTERP, value);
        }
    }

    void i_add_int(INTVAL value) {
        RE(SELF) += (FLOATVAL)value;
    }
    void i_add_float(FLOATVAL value) {
        RE(SELF) += value;
    }

/*

=item C<PMC* subtract(PMC* value, PMC* dest)>

=item C<PMC* subtract_int(INTVAL value, PMC* dest)>

=item C<PMC* subtract_float(FLOATVAL value, PMC* dest)>

Subtracts C<value> from the complex number, placing the result in C<dest>.

=cut

*/

    PMC* subtract(PMC* value, PMC* dest) {
MMD_Complex: {
            if (dest)
                VTABLE_morph(INTERP, dest, SELF->vtable->base_type);
            else {
                dest = pmc_new(INTERP, SELF->vtable->base_type);
            }
            RE(dest) = RE(SELF) - RE(value);
            IM(dest) = IM(SELF) - IM(value);
            return dest;
        }
MMD_DEFAULT: {
            if (dest)
                VTABLE_morph(INTERP, dest, SELF->vtable->base_type);
            else {
                dest = pmc_new(INTERP, SELF->vtable->base_type);
            }
            RE(dest) = RE(SELF) - VTABLE_get_number(INTERP, value);
            IM(dest) = IM(SELF);
            return dest;
        }
    }

    PMC* subtract_int(INTVAL value, PMC* dest) {
        if (dest)
            VTABLE_morph(INTERP, dest, SELF->vtable->base_type);
        else {
            dest = pmc_new(INTERP, SELF->vtable->base_type);
        }
        RE(dest) = RE(SELF) - value;
        IM(dest) = IM(SELF);
        return dest;
    }

    PMC* subtract_float(FLOATVAL value, PMC* dest) {
        if (dest)
            VTABLE_morph(INTERP, dest, SELF->vtable->base_type);
        else {
            dest = pmc_new(INTERP, SELF->vtable->base_type);
        }
        RE(dest) = RE(SELF) - value;
        IM(dest) = IM(SELF);
        return dest;
    }

    void i_subtract(PMC* value) {
MMD_Complex: {
            RE(SELF) -= RE(value);
            IM(SELF) -= IM(value);
        }
MMD_DEFAULT: {
            RE(SELF) -= VTABLE_get_number(INTERP, value);
        }
    }

    void i_subtract_int(INTVAL value) {
        RE(SELF) -= (FLOATVAL)value;
    }
    void i_subtract_float(FLOATVAL value) {
        RE(SELF) -= value;
    }
/*

=item C<PMC* multiply(PMC* value, PMC* dest)>

=item C<PMC* multiply_int(INTVAL value, PMC* dest)>

=item C<PMC* multiply_float(FLOATVAL value, PMC* dest)>

Multiplies the complex number with C<value>, placing the result in C<dest>.

=item C<void i_multiply(PMC* value)>

=item C<void i_multiply_int(INTVAL value)>

=item C<void i_multiply_float(FLOATVAL value)>

Multiplies the complex number SELF inplace with C<value>.

=cut

*/

/*

  TODO for better precision:

  (a+ib)(c+id)=(ac-bd)+i((a+b)(c+d)-ac-bd).
  (a+bi)(c+di)=(ac-bd)+i(ad+bc)

*/
    PMC* multiply(PMC* value, PMC* dest) {
MMD_Complex: {
            FLOATVAL re = RE(SELF) * RE(value) - IM(SELF) * IM(value);
            FLOATVAL im = IM(SELF) * RE(value) + RE(SELF) * IM(value);
            if (dest)
                VTABLE_morph(INTERP, dest, SELF->vtable->base_type);
            else {
                dest = pmc_new(INTERP, SELF->vtable->base_type);
            }
            RE(dest) = re;
            IM(dest) = im;
            return dest;
        }
MMD_DEFAULT: {
            FLOATVAL re = RE(SELF) * VTABLE_get_number(INTERP, value);
            FLOATVAL im = IM(SELF) * VTABLE_get_number(INTERP, value);
            if (dest)
                VTABLE_morph(INTERP, dest, SELF->vtable->base_type);
            else {
                dest = pmc_new(INTERP, SELF->vtable->base_type);
            }
            RE(dest) = re;
            IM(dest) = im;
            return dest;
        }
    }

    PMC* multiply_int(INTVAL value, PMC* dest) {
        if (dest)
            VTABLE_morph(INTERP, dest, SELF->vtable->base_type);
        else {
            dest = pmc_new(INTERP, SELF->vtable->base_type);
        }
        RE(dest) = RE(SELF) * value;
        IM(dest) = IM(SELF) * value;
        return dest;
    }

    PMC* multiply_float(FLOATVAL value, PMC* dest) {
        if (dest)
            VTABLE_morph(INTERP, dest, SELF->vtable->base_type);
        else {
            dest = pmc_new(INTERP, SELF->vtable->base_type);
        }
        RE(dest) = RE(SELF) * value;
        IM(dest) = IM(SELF) * value;
        return dest;
    }

    void i_multiply(PMC* value) {
MMD_Complex: {
            FLOATVAL re = RE(SELF) * RE(value) - IM(SELF) * IM(value);
            FLOATVAL im = IM(SELF) * RE(value) + RE(SELF) * IM(value);
            RE(SELF) = re;
            IM(SELF) = im;
        }
MMD_DEFAULT: {
            RE(SELF) *= VTABLE_get_number(INTERP, value);
            IM(SELF) *= VTABLE_get_number(INTERP, value);
        }
    }

    void i_multiply_int(INTVAL value) {
        RE(SELF) *= value;
        IM(SELF) *= value;
    }

    void i_multiply_float(FLOATVAL value) {
        RE(SELF) *= value;
        IM(SELF) *= value;
    }

/*

=item C<PMC* divide(PMC* value, PMC* dest)>

=item C<PMC* divide_int(INTVAL value, PMC* dest)>

=item C<PMC* divide_float(FLOATVAL value, PMC* dest)>

Divide the complex number by C<value>, placing the result in C<dest>.

=item C<void i_divide(PMC* value, PMC* dest)>

=item C<void i_divide_int(INTVAL value, PMC* dest)>

=item C<void i_divide_float(FLOATVAL value, PMC* dest)>

Divide the complex number C<SELF> by C<value> inplace.

Throws divide by zero exception if divisor is zero.

=cut

TODO: for better fp precision
http://docs.sun.com/source/806-3568/ncg_goldberg.html
(a+ib)/(c+id) =
    (a + b(d/c)) / (c + d(d/c)) + i(b - a(d/c)) / (c + d(d/c)) if |d|<|c|
    (b + a(c/d)) / (d + c(c/d)) + i(-a + b(c/d)) / (d + c(c/d)) if |d|>=|c|

*/

    PMC* divide(PMC* value, PMC* dest) {
MMD_Complex: {
            FLOATVAL mod, re, im;

            complex_check_divide_zero(INTERP, value);
            if (dest)
                VTABLE_morph(INTERP, dest, SELF->vtable->base_type);
            else {
                dest = pmc_new(INTERP, SELF->vtable->base_type);
            }
            /* a little speed optimisation: cache an intermediate number;
                I'm not sure the compiler does this */
            mod = (RE(value) * RE(value) + IM(value) * IM(value));
            re = (RE(SELF) * RE(value) + IM(SELF) * IM(value)) / mod;
            im = (IM(SELF) * RE(value) - RE(SELF) * IM(value)) / mod;
            RE(dest) = re;
            IM(dest) = im;
            return dest;
        }
MMD_DEFAULT: {
            FLOATVAL d = VTABLE_get_number(INTERP, value);
            float_check_divide_zero(INTERP, d);
            if (dest)
                VTABLE_morph(INTERP, dest, SELF->vtable->base_type);
            else {
                dest = pmc_new(INTERP, SELF->vtable->base_type);
            }
            RE(dest) = RE(SELF) / d;
            IM(dest) = IM(SELF) / d;
            return dest;
        }
    }

    PMC* divide_int(INTVAL value, PMC* dest) {
        int_check_divide_zero(INTERP, value);
        if (dest)
            VTABLE_morph(INTERP, dest, SELF->vtable->base_type);
        else {
            dest = pmc_new(INTERP, SELF->vtable->base_type);
        }
        RE(dest) = RE(SELF) / value;
        IM(dest) = IM(SELF) / value;
        return dest;
    }

    PMC* divide_float(FLOATVAL value, PMC* dest) {
        float_check_divide_zero(INTERP, value);
        if (dest)
            VTABLE_morph(INTERP, dest, SELF->vtable->base_type);
        else {
            dest = pmc_new(INTERP, SELF->vtable->base_type);
        }
        RE(dest) = RE(SELF) / value;
        IM(dest) = IM(SELF) / value;
        return dest;
    }

    void i_divide(PMC* value) {
MMD_Complex: {
            FLOATVAL mod, re, im;

            complex_check_divide_zero(INTERP, value);
            /* a little speed optimisation: cache an intermediate number;
                I'm not sure the compiler does this */
            mod = (RE(value) * RE(value) + IM(value) * IM(value));
            re = (RE(SELF) * RE(value) + IM(SELF) * IM(value)) / mod;
            im = (IM(SELF) * RE(value) - RE(SELF) * IM(value)) / mod;
            RE(SELF) = re;
            IM(SELF) = im;
        }
MMD_DEFAULT: {
            FLOATVAL d = VTABLE_get_number(INTERP, value);
            float_check_divide_zero(INTERP, d);
            RE(SELF) /= d;
            IM(SELF) /= d;
        }
    }

    void i_divide_int(INTVAL value) {
        int_check_divide_zero(INTERP, value);
        RE(SELF) /= value;
        IM(SELF) /= value;
    }

    void i_divide_float(FLOATVAL value) {
        float_check_divide_zero(INTERP, value);
        RE(SELF) /= value;
        IM(SELF) /= value;
    }


/*

=item C<PMC* neg(PMC *dest)>

=item C<void neg()>

Set C<dest> to the negated value of C<SELF>.

=cut

*/

    PMC* neg(PMC* dest) {
        if (!dest)
            dest = pmc_new(INTERP, SELF->vtable->base_type);
        else
            VTABLE_morph(INTERP, dest, SELF->vtable->base_type);
        RE(dest) = - RE(SELF);
        IM(dest) = - IM(SELF);
        return dest;
    }

    void i_neg() {
        RE(SELF) = - RE(SELF);
        IM(SELF) = - IM(SELF);
    }

/*

=item C<INTVAL is_equal(PMC* value)>

Compares the complex number with C<value> and returns true if they are equal.

=cut

*/

    INTVAL is_equal(PMC* value) {
MMD_Complex: {
            return (INTVAL)(
                    RE(SELF) == RE(value) &&
                    IM(SELF) == IM(value));
        }
MMD_DEFAULT: {
            if (IM(SELF) != 0.0)
                return (INTVAL)0;
            return (RE(SELF) == VTABLE_get_number(INTERP, value));
        }
    }

/*

=item C<PMC* absolute(PMC* dest)>

=item C<void i_absolute()>

Sets C<dest> to the absolute value of SELF that is the distance from (0.0).

=cut

*/

/*

  TODO for better precision: hinted by vaxman according to "Numerical Recipes
  in Fortran 77", 2nd edition, Press, Vetterling, Teukolsky, Flannery,
  Cambridge University Press, 2001, pp. 171ff:


|a+ib|=|a|*sqrt(1+(b/a)**2), if |a|>=|b|,
       |b|*sqrt(1+(a/b)**2)  else.

*/

    PMC* absolute(PMC *dest) {
        FLOATVAL d = sqrt(RE(SELF)*RE(SELF) + IM(SELF)*IM(SELF));
        if (!dest)
            dest = pmc_new(INTERP,
                Parrot_get_ctx_HLL_type(INTERP, enum_class_Float));
        VTABLE_set_number_native(INTERP, dest, d);
        return dest;
    }

    void i_absolute() {
        FLOATVAL d = sqrt(RE(SELF)*RE(SELF) + IM(SELF)*IM(SELF));
        VTABLE_morph(INTERP, SELF, enum_class_Float);   /* XXX */
        VTABLE_set_number_native(INTERP, SELF, d);
    }

/*

=item C<METHOD PMC* ln()>

Returns the natural logorithm of SELF.

=cut

ln z = ln |z| + i arg(z)
|x + iy| = sqrt(x^2 + y^2)
arg(x + iy) = atan2(y, x)

Some special cases
ln(-1) = pi i
ln(0) = -inf
ln(1) = 0
ln(e) = 1
ln(+-i) = +- (pi i)/2

*/

    METHOD PMC* ln() {
        FLOATVAL re, im;
        PMC *d = pmc_new(INTERP, SELF->vtable->base_type);
        re = RE(SELF);
        im = IM(SELF);
        /* This is necessary for atan2 to behave */
        if (im == -0.0) im = 0.0;
        RE(d) = log(sqrt(re*re + im*im));
        IM(d) = atan2(im, re);
        return d;
    }

/*

=item C<METHOD PMC* exp()>

Returns e ^ SELF.

=cut

exp(a + bi) = exp(a) * (cos(b) + i * sin(b))

*/

    METHOD PMC* exp() {
        FLOATVAL re, im, f;
        PMC *d = pmc_new(INTERP, SELF->vtable->base_type);
        re = RE(SELF);
        im = IM(SELF);
        f = exp(re);
        RE(d) = f * cos(im);
        /* If only sin(pi) worked. */
        if (im == 4.0 * atan(1.0))
            IM(d) = 0.0;
        else
            IM(d) = f * sin(im);
        return d;
    }

/*

=item C<METHOD PMC* sin()>

=item C<METHOD PMC* cos()>

=item C<METHOD PMC* tan()>

=item C<METHOD PMC* csc()>

=item C<METHOD PMC* sec()>

=item C<METHOD PMC* cot()>

Returns C<FUNC>(SELF).

=cut

 => sin(a + bi) = sin(a)cosh(b)+i*cos(a)sinh(b)
    sin(z) = ((e ^ zi) - (e ^ -zi)) / (2i)
 => cos(a + bi) = cos(a) * cosh(b) - i * sin(a) * sinh(b)
    cos(z) = ((e ^ zi) + (e ^ -zi)) / 2

    sin(iz) = i sinh(z)
    cos(iz) = cosh(z)

    sinh(iz) = i sin(z)
    cosh(iz) = cos z
    sinh(a + bi) = sinh(a) * cos(b) + i * cosh(a) * sin(b)
    cosh(a + bi) = cosh(a) * cos(b) + i * sinh(a) * sin(b)

*/

    METHOD PMC* sin() {
        FLOATVAL re, im;
        PMC *d = pmc_new(INTERP, SELF->vtable->base_type);
        re = RE(SELF);
        im = IM(SELF);
        if (im == 0.0) {
            RE(d) = sin(re);
            IM(d) = 0.0;
            return d;
        }
        if (re == 0.0) {
            RE(d) = 0.0;
            IM(d) = sinh(im);
            return d;
        }
        RE(d) = sin(re) * cosh(im);
        if (im == 0.0)
            IM(d) = im;
        else
            IM(d) = cos(re) * sinh(im);
        return d;
    }

    METHOD PMC* cos() {
        FLOATVAL re, im;
        PMC *d = pmc_new(INTERP, SELF->vtable->base_type);
        re = RE(SELF);
        im = IM(SELF);
        if (re == 0.0) {
            RE(d) = cosh(im);
            IM(d) = 0.0;
            return d;
        }
        if (im == 0.0) {
            RE(d) = cos(re);
            IM(d) = 0.0;
            return d;
        }
        RE(d) = cos(re) * cosh(im);
        IM(d) = -1.0 * sin(re) * sinh(im);
        return d;
    }

    METHOD PMC* tan() {
        PMC *d = pmc_new(INTERP, SELF->vtable->base_type);
        PMC *e = pmc_new(INTERP, SELF->vtable->base_type);
        d = Parrot_Complex_nci_sin(INTERP, SELF);
        e = Parrot_Complex_nci_cos(INTERP, SELF);
        Parrot_Complex_i_divide_Complex(INTERP, d, e);
        return d;
    }

    METHOD PMC* cot() {
        PMC *d = pmc_new(INTERP, SELF->vtable->base_type);
        PMC *e = pmc_new(INTERP, SELF->vtable->base_type);
        d = Parrot_Complex_nci_cos(INTERP, SELF);
        e = Parrot_Complex_nci_sin(INTERP, SELF);
        Parrot_Complex_i_divide_Complex(INTERP, d, e);
        return d;
    }

    METHOD PMC* sec() {
        PMC *d = pmc_new(INTERP, SELF->vtable->base_type);
        PMC *e = pmc_new(INTERP, SELF->vtable->base_type);
        RE(d) = 1.0;
        IM(d) = 0.0;
        e = Parrot_Complex_nci_cos(INTERP, SELF);
        Parrot_Complex_i_divide_Complex(INTERP, d, e);
        return d;
    }

    METHOD PMC* csc() {
        PMC *d = pmc_new(INTERP, SELF->vtable->base_type);
        PMC *e = pmc_new(INTERP, SELF->vtable->base_type);
        RE(d) = 1.0;
        IM(d) = 0.0;
        e = Parrot_Complex_nci_sin(INTERP, SELF);
        Parrot_Complex_i_divide_Complex(INTERP, d, e);
        return d;
    }

/*

=item C<METHOD PMC* asin()>

=item C<METHOD PMC* acos()>

=item C<METHOD PMC* atan()>

=item C<METHOD PMC* acsc()>

=item C<METHOD PMC* asec()>

=item C<METHOD PMC* acot()>

Returns the inverse function of SELF.

=cut

 => arcsin z = -i ln(iz + sqrt(1-z*z))
 => arccos z = pi/2 + i * ln(iz + sqrt(1 - z*z))
    arccos z = -i ln(z + sqrt(z*z-1))
 => arctan z = i/2 ln((i+z) / (i-z))
    arctan z = 1/2 i (ln(1-iz) - ln(1 + iz))

 => acot(z) = atan(1 / z)
    acot(z) = i/2 (ln((z - i) / z) - ln((z + i) / z))
 => asec(z) = acos(1 / z)
    asec(z) = 1/2 pi + i ln(sqrt(1 - 1/zz) + i/z)
 => acsc(z) = asin(1 / z)
    acsc(z) = -i ln(sqrt(1 - 1/zz + i/z))

*/

    METHOD PMC* asin() {
        PMC *d = pmc_new(INTERP, SELF->vtable->base_type);
        PMC *e = pmc_new(INTERP, SELF->vtable->base_type);
        e = Parrot_Complex_multiply_Complex(INTERP, SELF, SELF, e);
        RE(e) = 1.0 - RE(e);
        IM(e) = -IM(e);
        d = Parrot_Complex_nci_sqrt(INTERP, e);
        RE(d) -= IM(SELF);
        IM(d) += RE(SELF);
        d = Parrot_Complex_nci_ln(INTERP, d);
        RE(e) = IM(d);
        IM(e) = -RE(d);
        return e;
    }

    METHOD PMC* acos() {
        PMC *d = pmc_new(INTERP, SELF->vtable->base_type);
        PMC *e = pmc_new(INTERP, SELF->vtable->base_type);
        e = Parrot_Complex_multiply_Complex(INTERP, SELF, SELF, e);
        RE(e) = 1.0 - RE(e);
        IM(e) = -IM(e);
        d = Parrot_Complex_nci_sqrt(INTERP, e);
        RE(d) += IM(SELF);
        IM(d) -= RE(SELF);
        e = Parrot_Complex_nci_ln(INTERP, d);
        RE(d) = IM(e) + 2.0 * atan(1);
        IM(d) = -RE(e);
        return d;
    }

    METHOD PMC* atan() {
        FLOATVAL re, im;
        PMC *d = pmc_new(INTERP, SELF->vtable->base_type);
        PMC *e = pmc_new(INTERP, SELF->vtable->base_type);
        re = RE(SELF);
        im = IM(SELF);
        RE(d) = re;
        IM(d) = 1 + im;
        RE(e) = -re;
        IM(e) = 1 - im;
        Parrot_Complex_i_divide_Complex(INTERP, d, e);
        d = Parrot_Complex_nci_ln(INTERP, d);
        RE(e) = IM(d) / -2.0;
        IM(e) = RE(d) / 2.0;
        return e;
    }


    METHOD PMC* acot() {
        FLOATVAL a, b;
        PMC *d = pmc_new(INTERP, SELF->vtable->base_type);
        a = RE(SELF);
        b = IM(SELF);
        RE(d) = a/(a*a + b*b);
        IM(d) = -b/(a*a + b*b);
        d = Parrot_Complex_nci_atan(INTERP, d);
        return d;
    }

    METHOD PMC* acsc() {
        FLOATVAL a, b;
        PMC *d = pmc_new(INTERP, SELF->vtable->base_type);
        a = RE(SELF);
        b = IM(SELF);
        RE(d) = a/(a*a + b*b);
        IM(d) = -b/(a*a + b*b);
        d = Parrot_Complex_nci_asin(INTERP, d);
        return d;
    }

    METHOD PMC* asec() {
        FLOATVAL a, b;
        PMC *d = pmc_new(INTERP, SELF->vtable->base_type);
        a = RE(SELF);
        b = IM(SELF);
        RE(d) = a/(a*a + b*b);
        IM(d) = -b/(a*a + b*b);
        d = Parrot_Complex_nci_acos(INTERP, d);
        return d;
    }

/*

=item C<METHOD PMC* sinh()>

Returns the arctangent of SELF.

=item C<METHOD PMC* cosh()>

Returns the arcsine of SELF.

=item C<METHOD PMC* tanh()>

Returns the arccosine of SELF.

=cut

tanh(z) = sinh(z) / cosh(z)

*/

    METHOD PMC* sinh() {
        FLOATVAL a, b;
        PMC *d = pmc_new(INTERP, SELF->vtable->base_type);
        a = RE(SELF);
        b = IM(SELF);
        RE(d) = sinh(a) * cos(b);
        IM(d) = cosh(a) * sin(b);
        return d;
    }

    METHOD PMC* cosh() {
        FLOATVAL a, b;
        PMC *d = pmc_new(INTERP, SELF->vtable->base_type);
        a = RE(SELF);
        b = IM(SELF);
        RE(d) = cosh(a) * cos(b);
        IM(d) = sinh(a) * sin(b);
        return d;
    }

    METHOD PMC* tanh() {
        PMC *d = pmc_new(INTERP, SELF->vtable->base_type);
        PMC *e = pmc_new(INTERP, SELF->vtable->base_type);
        d = Parrot_Complex_nci_sinh(INTERP, SELF);
        e = Parrot_Complex_nci_cosh(INTERP, SELF);
        Parrot_Complex_i_divide_Complex(INTERP, d, e);
        return d;
    }

    METHOD PMC* coth() {
        FLOATVAL a, b;
        PMC *d = pmc_new(INTERP, SELF->vtable->base_type);
        d = Parrot_Complex_nci_tanh(INTERP, SELF);
        a = RE(d);
        b = IM(d);
        RE(d) = a/(a*a + b*b);
        IM(d) = -b/(a*a + b*b);
        return d;
    }

    METHOD PMC* csch() {
        FLOATVAL a, b;
        PMC *d = pmc_new(INTERP, SELF->vtable->base_type);
        d = Parrot_Complex_nci_sinh(INTERP, SELF);
        a = RE(d);
        b = IM(d);
        RE(d) = a/(a*a + b*b);
        IM(d) = -b/(a*a + b*b);
        return d;
    }

    METHOD PMC* sech() {
        FLOATVAL a, b;
        PMC *d = pmc_new(INTERP, SELF->vtable->base_type);
        d = Parrot_Complex_nci_cosh(INTERP, SELF);
        a = RE(d);
        b = IM(d);
        RE(d) = a/(a*a + b*b);
        IM(d) = -b/(a*a + b*b);
        return d;
    }

/*

=item C<METHOD PMC* asinh()>

=item C<METHOD PMC* acosh()>

=item C<METHOD PMC* atanh()>

=item C<METHOD PMC* acsch()>

=item C<METHOD PMC* asech()>

=item C<METHOD PMC* acoth()>

The inverse hyperbolic functions.  Currently all broken, but for
C<func(a+bi) = c+di>, C<|c|> and C<|d|> will be correct, confusingly enough.

=cut

asinh z = -ln(sqrt(1+zz) - z)
asinh z = ln(sqrt(zz + 1) + z)

asinh = i asin(-ix)
acosh = i acos(x)
atanh = i atan(-ix)

*/

    METHOD PMC* asinh() {
        PMC *d = pmc_new(INTERP, SELF->vtable->base_type);
        PMC *e = pmc_new(INTERP, SELF->vtable->base_type);
        RE(d) = IM(SELF);
        IM(d) = -RE(SELF);
        d = Parrot_Complex_nci_asin(INTERP, d);
        RE(e) = -IM(d);
        IM(e) = RE(d);
        return e;
    }

    METHOD PMC* acosh() {
        PMC *d = pmc_new(INTERP, SELF->vtable->base_type);
        PMC *e = pmc_new(INTERP, SELF->vtable->base_type);
        d = Parrot_Complex_nci_acos(INTERP, SELF);
        RE(e) = -IM(d);
        IM(e) = RE(d);
        return e;
    }

    METHOD PMC* atanh() {
        PMC *d = pmc_new(INTERP, SELF->vtable->base_type);
        PMC *e = pmc_new(INTERP, SELF->vtable->base_type);
        RE(d) = IM(SELF);
        IM(d) = -RE(SELF);
        d = Parrot_Complex_nci_atan(INTERP, d);
        RE(e) = -IM(d);
        IM(e) = RE(d);
        return e;
    }

    METHOD PMC* acoth() {
        FLOATVAL a, b;
        PMC *d = pmc_new(INTERP, SELF->vtable->base_type);
        a = RE(SELF);
        b = IM(SELF);
        RE(d) = a/(a*a + b*b);
        IM(d) = -b/(a*a + b*b);
        d = Parrot_Complex_nci_atanh(INTERP, d);
        return d;
    }

    METHOD PMC* acsch() {
        FLOATVAL a, b;
        PMC *d = pmc_new(INTERP, SELF->vtable->base_type);
        a = RE(SELF);
        b = IM(SELF);
        RE(d) = a/(a*a + b*b);
        IM(d) = -b/(a*a + b*b);
        d = Parrot_Complex_nci_asinh(INTERP, d);
        return d;
    }

    METHOD PMC* asech() {
        FLOATVAL a, b;
        PMC *d = pmc_new(INTERP, SELF->vtable->base_type);
        a = RE(SELF);
        b = IM(SELF);
        RE(d) = a/(a*a + b*b);
        IM(d) = -b/(a*a + b*b);
        d = Parrot_Complex_nci_acosh(INTERP, d);
        return d;
    }

/*

=item C<PMC* pow(PMC* value, PMC* dest)>

Return SELF to the C<value>th power and return result in C<dest>.

=item C<METHOD PMC* sqrt()>

Return the square root of SELF.

=cut

TODO: mmd in other pmc's to allow .Integer ^ .Complex, etc.
and i_pow, and pow_(float|int), etc

x ^ y = exp(y * ln x))

*/

    PMC* pow(PMC* value, PMC* dest) {
MMD_Complex: {
                 PMC *l = pmc_new(INTERP, SELF->vtable->base_type);
                 if (dest)
                     VTABLE_morph(INTERP, dest, SELF->vtable->base_type);
                 else
                     dest = pmc_new(INTERP, SELF->vtable->base_type);
                 l = Parrot_Complex_multiply_Complex(INTERP,
                         Parrot_Complex_nci_ln(INTERP, SELF), value, l);
                 dest = Parrot_Complex_nci_exp(INTERP, l);
                 return dest;
             }
MMD_DEFAULT: {
                 PMC *l = pmc_new(INTERP, SELF->vtable->base_type);
                 if (dest)
                     VTABLE_morph(INTERP, dest, SELF->vtable->base_type);
                 else
                     dest = pmc_new(INTERP, SELF->vtable->base_type);
                 l = Parrot_Complex_multiply(INTERP,
                         Parrot_Complex_nci_ln(INTERP, SELF), value, l);
                 dest = Parrot_Complex_nci_exp(INTERP, l);
                 return dest;
             }
    }

    /* sqrt(x) = exp(ln(x)/2) */
    METHOD PMC* sqrt() {
        PMC *d = pmc_new(INTERP, SELF->vtable->base_type);
        d = Parrot_Complex_nci_ln(INTERP, SELF);
        RE(d) *= 0.5;
        IM(d) *= 0.5;
        d = Parrot_Complex_nci_exp(INTERP, d);
        return d;
    }

}

/*

=back

=cut

*/

/*
 * Local variables:
 *   c-file-style: "parrot"
 * End:
 * vim: expandtab shiftwidth=4:
 */