nqp_bigint.ops

BEGIN_OPS_PREAMBLE

/* Parroty includes. */
#include "parrot/parrot.h"
#include "parrot/extend.h"
#include "parrot/dynext.h"
#include "../6model/sixmodelobject.h"
#include "../6model/reprs/P6bigint.h"

/* The ID of the bigint REPR. */
static INTVAL bigint_repr_id = 0;

static mp_int * get_bigint(PARROT_INTERP, PMC *obj) {
    struct SixModel_REPROps *r = REPR(obj);
    if (r->ID == bigint_repr_id)
        return &((P6bigintInstance *)PMC_data(obj))->body.i;
    else
        return &((P6bigintBody *)r->box_funcs->get_boxed_ref(interp, STABLE(obj), OBJECT_BODY(obj),
            bigint_repr_id))->i;
}

static FLOATVAL mp_get_double(mp_int *a) {
    FLOATVAL d    = 0.0;
    FLOATVAL sign = SIGN(a) == MP_NEG ? -1.0 : 1.0;
    int i;
    if (USED(a) == 0)
        return d;
    if (USED(a) == 1)
        return sign * (FLOATVAL) DIGIT(a, 0);

    mp_clamp(a);
    i = USED(a) - 1;
    d = (FLOATVAL) DIGIT(a, i);
    i--;
    if (i == -1) {
        return sign * d;
    }
    d *= pow(2.0, DIGIT_BIT);
    d += (FLOATVAL) DIGIT(a, i);

    d *= pow(2.0, DIGIT_BIT * i);
    return sign * d;
}


static void from_num(FLOATVAL d, mp_int *a) {
    FLOATVAL d_digit = pow(2, DIGIT_BIT);
    FLOATVAL da      = fabs(d);
    FLOATVAL upper;
    FLOATVAL lower;
    int      digits  = 0;

    mp_zero(a);

    while (da > d_digit * d_digit) {;
        da /= d_digit;
        digits++;
    }
    mp_grow(a, digits + 2);

    /* populate the top 2 digits */
    upper = da / d_digit;
    lower = fmod(da, d_digit);
    if (upper >= 1) {
        mp_set_long(a, (unsigned long) upper);
        mp_mul_2d(a, DIGIT_BIT , a);
    } else {
        a->used = 1;
    }
    DIGIT(a, 0) = (mp_digit) lower;

    /* shift the rest */
    mp_mul_2d(a, DIGIT_BIT * digits, a);
    if (d < 0)
        mp_neg(a, a);
    mp_clamp(a);
    mp_shrink(a);
}

static void grow_and_negate(mp_int *a, int size, mp_int *b) {
    int i;
    int actual_size = MAX(size, USED(a));
    mp_zero(b);
    mp_grow(b, actual_size);
    USED(b) = actual_size;
    for (i = 0; i < actual_size; i++) {
        DIGIT(b, i) = (~DIGIT(a, i)) & MP_MASK;
    }
    mp_add_d(b, 1, b);
}

static void two_complement_bitop(mp_int *a, mp_int *b, mp_int *c,
        int (*mp_bitop)(mp_int *, mp_int *, mp_int *)) {
    mp_int d;
    if (SIGN(a) ^ SIGN(b)) {
        /* exactly one of them is negative, so need to perform
         * some magic. tommath stores a sign bit, but Perl 6 expects
         * 2's complement */
        mp_init(&d);
        if (MP_NEG == SIGN(a)) {
            grow_and_negate(a, USED(b), &d);
            mp_bitop(&d, b, c);
        } else {
            grow_and_negate(b, USED(a), &d);
            mp_bitop(a, &d, c);
        }
        if (DIGIT(c, USED(c) - 1) & (MP_MASK > 1)) {
            grow_and_negate(c, c->used, &d);
            mp_copy(&d, c);
            mp_neg(c, c);
        }
        mp_clear(&d);
    } else {
        mp_bitop(a, b, c);
    }

}


END_OPS_PREAMBLE

/*
 * The ops in here mostly just delegate off to libtommath. nqp_bigint_setup must be
 * called before you create any types using the bigint representation.
 */

inline op nqp_bigint_setup() :base_core {
    /* Register the bigint representation. */
    if (!bigint_repr_id)
        bigint_repr_id = REGISTER_DYNAMIC_REPR(interp,
            Parrot_str_new_constant(interp, "P6bigint"),
            P6bigint_initialize);
}

inline op nqp_bigint_add(out PMC, in PMC, in PMC, in PMC) :base_core {
    mp_int *a = get_bigint(interp, $2);
    mp_int *b = get_bigint(interp, $3);
    $1 = REPR($4)->allocate(interp, STABLE($4));
    REPR($1)->initialize(interp, STABLE($1), OBJECT_BODY($1));
    mp_add(a, b, get_bigint(interp, $1));
}

inline op nqp_bigint_sub(out PMC, in PMC, in PMC, in PMC) :base_core {
    mp_int *a = get_bigint(interp, $2);
    mp_int *b = get_bigint(interp, $3);
    $1 = REPR($4)->allocate(interp, STABLE($4));
    REPR($1)->initialize(interp, STABLE($1), OBJECT_BODY($1));
    mp_sub(a, b, get_bigint(interp, $1));
}

inline op nqp_bigint_mul(out PMC, in PMC, in PMC, in PMC) :base_core {
    mp_int *a = get_bigint(interp, $2);
    mp_int *b = get_bigint(interp, $3);
    $1 = REPR($4)->allocate(interp, STABLE($4));
    REPR($1)->initialize(interp, STABLE($1), OBJECT_BODY($1));
    mp_mul(a, b, get_bigint(interp, $1));
}

inline op nqp_bigint_div(out PMC, in PMC, in PMC, in PMC) :base_core {
    mp_int *a = get_bigint(interp, $2);
    mp_int *b = get_bigint(interp, $3);
    int result;
    $1 = REPR($4)->allocate(interp, STABLE($4));
    REPR($1)->initialize(interp, STABLE($1), OBJECT_BODY($1));
    result = mp_div(a, b, get_bigint(interp, $1), NULL);
    if (result == MP_VAL)
        Parrot_ex_throw_from_c_args(interp, NULL, EXCEPTION_INVALID_OPERATION,
            "Divide by zero");
}

inline op nqp_bigint_mod(out PMC, in PMC, in PMC, in PMC) :base_core {
    mp_int *a = get_bigint(interp, $2);
    mp_int *b = get_bigint(interp, $3);
    $1 = REPR($4)->allocate(interp, STABLE($4));
    REPR($1)->initialize(interp, STABLE($1), OBJECT_BODY($1));
    mp_mod(a, b, get_bigint(interp, $1));
}

inline op nqp_bigint_exp_mod(out PMC, in PMC, in PMC, in PMC, in PMC) :base_core {
    mp_int *a = get_bigint(interp, $2);
    mp_int *b = get_bigint(interp, $3);
    mp_int *c = get_bigint(interp, $4);
    $1 = REPR($5)->allocate(interp, STABLE($5));
    REPR($1)->initialize(interp, STABLE($1), OBJECT_BODY($1));
    mp_exptmod(a, b, c, get_bigint(interp, $1));
}


inline op nqp_bigint_neg(out PMC, in PMC, in PMC) :base_core {
    mp_int *a = get_bigint(interp, $2);
    $1 = REPR($3)->allocate(interp, STABLE($3));
    REPR($1)->initialize(interp, STABLE($1), OBJECT_BODY($1));
    mp_neg(a, get_bigint(interp, $1));
}
inline op nqp_bigint_abs(out PMC, in PMC, in PMC) :base_core {
    mp_int *a = get_bigint(interp, $2);
    $1 = REPR($3)->allocate(interp, STABLE($3));
    REPR($1)->initialize(interp, STABLE($1), OBJECT_BODY($1));
    mp_abs(a, get_bigint(interp, $1));
}

inline op nqp_bigint_cmp(out INT, in PMC, in PMC) :base_core {
    mp_int *a = get_bigint(interp, $2);
    mp_int *b = get_bigint(interp, $3);
    $1 = mp_cmp(a, b);
}

inline op nqp_bigint_bool(out INT, in PMC) :base_core {
    $1 = !mp_iszero(get_bigint(interp, $2));
}

inline op nqp_bigint_eq(out INT, in PMC, in PMC) :base_core {
    mp_int *a = get_bigint(interp, $2);
    mp_int *b = get_bigint(interp, $3);
    $1 = MP_EQ == mp_cmp(a, b);
}

inline op nqp_bigint_ne(out INT, in PMC, in PMC) :base_core {
    mp_int *a = get_bigint(interp, $2);
    mp_int *b = get_bigint(interp, $3);
    $1 = MP_EQ != mp_cmp(a, b);
}

inline op nqp_bigint_gt(out INT, in PMC, in PMC) :base_cor {
    mp_int *a = get_bigint(interp, $2);
    mp_int *b = get_bigint(interp, $3);
    $1 = MP_GT == mp_cmp(a, b);
}

inline op nqp_bigint_ge(out INT, in PMC, in PMC) :base_core {
    mp_int *a = get_bigint(interp, $2);
    mp_int *b = get_bigint(interp, $3);
    $1 = MP_LT != mp_cmp(a, b);
}

inline op nqp_bigint_lt(out INT, in PMC, in PMC) :base_core {
    mp_int *a = get_bigint(interp, $2);
    mp_int *b = get_bigint(interp, $3);
    $1 = MP_LT == mp_cmp(a, b);
}

inline op nqp_bigint_le(out INT, in PMC, in PMC) :base_core {
    mp_int *a = get_bigint(interp, $2);
    mp_int *b = get_bigint(interp, $3);
    $1 = MP_GT != mp_cmp(a, b);
}

inline op nqp_bigint_gcd(out PMC, in PMC, in PMC, in PMC) :base_core {
    mp_int *a = get_bigint(interp, $2);
    mp_int *b = get_bigint(interp, $3);
    $1 = REPR($4)->allocate(interp, STABLE($4));
    REPR($1)->initialize(interp, STABLE($1), OBJECT_BODY($1));
    mp_gcd(a, b, get_bigint(interp, $1));
}

inline op nqp_bigint_lcm(out PMC, in PMC, in PMC, in PMC) :base_core {
    mp_int *a = get_bigint(interp, $2);
    mp_int *b = get_bigint(interp, $3);
    $1 = REPR($4)->allocate(interp, STABLE($4));
    REPR($1)->initialize(interp, STABLE($1), OBJECT_BODY($1));
    mp_lcm(a, b, get_bigint(interp, $1));
}

inline op nqp_bigint_from_str(out PMC, in STR, in PMC) :base_core {
    char *buf = Parrot_str_cstring(interp, $2);
    $1 = REPR($3)->allocate(interp, STABLE($3));
    REPR($1)->initialize(interp, STABLE($1), OBJECT_BODY($1));
    mp_read_radix(get_bigint(interp, $1), buf, 10);
}

inline op nqp_bigint_to_str(out STR, in PMC) :base_core {
    mp_int *i = get_bigint(interp, $2);
    int len;
    char *buf;
    mp_radix_size(i, 10, &len);
    buf = mem_sys_allocate(len);
    mp_toradix_n(i, buf, 10, len);
    /* len - 1 because buf is \0-terminated */
    $1 = Parrot_str_new(interp, buf, len - 1);
    mem_sys_free(buf);
}

inline op nqp_bigint_to_str_base(out STR, in PMC, in INT) :base_core {
    mp_int *i = get_bigint(interp, $2);
    int len;
    char *buf;
    mp_radix_size(i, $3, &len);
    buf = mem_sys_allocate(len);
    mp_toradix_n(i, buf, $3, len);
    /* len - 1 because buf is \0-terminated */
    $1 = Parrot_str_new(interp, buf, len - 1);
    mem_sys_free(buf);
}

inline op nqp_bigint_to_num(out NUM, in PMC) :base_core {
    mp_int  *a = get_bigint(interp, $2);
    $1         = mp_get_double(a);
}

inline op nqp_bigint_from_num(out PMC, in NUM, in PMC) :base_core {
    $1 = REPR($3)->allocate(interp, STABLE($3));
    REPR($1)->initialize(interp, STABLE($1), OBJECT_BODY($1));
    from_num($2, get_bigint(interp, $1));
}

inline op nqp_bigint_shr(out PMC, in PMC, in INT, in PMC) :base_core {
    mp_int *a = get_bigint(interp, $2);
    $1 = REPR($4)->allocate(interp, STABLE($4));
    REPR($1)->initialize(interp, STABLE($1), OBJECT_BODY($1));
    mp_div_2d(a, $3, get_bigint(interp, $1), NULL);
}

inline op nqp_bigint_shl(out PMC, in PMC, in INT, in PMC) :base_core {
    mp_int *b;
    mp_int *a = get_bigint(interp, $2);
    $1 = REPR($4)->allocate(interp, STABLE($4));
    REPR($1)->initialize(interp, STABLE($1), OBJECT_BODY($1));
    b = get_bigint(interp, $1);
    mp_mul_2d(a, $3, b);
}

inline op nqp_bigint_band(out PMC, in PMC, in PMC, in PMC) :base_core {
    mp_int *a = get_bigint(interp, $2);
    mp_int *b = get_bigint(interp, $3);
    $1 = REPR($4)->allocate(interp, STABLE($4));
    REPR($1)->initialize(interp, STABLE($1), OBJECT_BODY($1));
    two_complement_bitop(a, b, get_bigint(interp, $1), mp_and);
}
inline op nqp_bigint_bor(out PMC, in PMC, in PMC, in PMC) :base_core {
    mp_int *a = get_bigint(interp, $2);
    mp_int *b = get_bigint(interp, $3);
    $1 = REPR($4)->allocate(interp, STABLE($4));
    REPR($1)->initialize(interp, STABLE($1), OBJECT_BODY($1));
    two_complement_bitop(a, b, get_bigint(interp, $1), mp_or);
}
inline op nqp_bigint_bxor(out PMC, in PMC, in PMC, in PMC) :base_core {
    mp_int *a = get_bigint(interp, $2);
    mp_int *b = get_bigint(interp, $3);
    $1 = REPR($4)->allocate(interp, STABLE($4));
    REPR($1)->initialize(interp, STABLE($1), OBJECT_BODY($1));
    two_complement_bitop(a, b, get_bigint(interp, $1), mp_xor);
}
inline op nqp_bigint_bnot(out PMC, in PMC, in PMC) :base_core {
    mp_int *a = get_bigint(interp, $2);
    mp_int *b;
    $1 = REPR($3)->allocate(interp, STABLE($3));
    REPR($1)->initialize(interp, STABLE($1), OBJECT_BODY($1));
    b = get_bigint(interp, $1);

    /* 2s complement: add 1, negate */
    mp_add_d(a, 1, b);
    mp_neg(b, b);
}


/*

=item nqp_bigint_radix(out, radix, str, pos, flag, type)

Convert string $3 into a number starting at offset $4 and using radix $2.
The result of the conversion returns a FixedPMCArray of size 3 with objects
of type $6, which is either a bigint or some type that boxes one.

The contents of $6 are the return value, the base to divide the value by if
it is something after a decimal dot, and the new position.

The $5 flags is a bitmask that modifies the parse and/or result:
    0x01: negate the result (useful if you've already parsed a minus)
    0x02: parse a leading +/- and negate the result on -
    0x04: parse trailing zeroes but do not include in result
          (for parsing values after a decimal point)

=cut

*/

inline op nqp_bigint_radix(out PMC, in INT, in STR, in INT, in INT, in PMC) :base_core {
    PMC      *out;
    INTVAL   radix  = $2;
    STRING   *str   = $3;
    INTVAL   zpos   = $4;
    INTVAL   flags  = $5;
    INTVAL   chars  = Parrot_str_length(interp, str);
    int      neg    = 0;
    INTVAL   ch;
    mp_int   zvalue;
    mp_int   zbase;
    PMC      *value_obj;
    mp_int   *value;
    PMC      *base_obj;
    mp_int   *base;
    PMC      *pos_obj;
    INTVAL   pos    = -1;
    
    if (radix > 36) {
        Parrot_ex_throw_from_c_args(interp, NULL, EXCEPTION_INVALID_OPERATION,
            "Cannot convert radix of %d (max 36)", radix);
    }
    
    mp_init(&zvalue);
    mp_init(&zbase);
    mp_set_int(&zbase, 1);
    
    value_obj = REPR($6)->allocate(interp, STABLE($6));
    REPR(value_obj)->initialize(interp, STABLE(value_obj), OBJECT_BODY(value_obj));
    value = get_bigint(interp, value_obj);
    
    base_obj = REPR($6)->allocate(interp, STABLE($6));
    REPR(base_obj)->initialize(interp, STABLE(base_obj), OBJECT_BODY(base_obj));
    base = get_bigint(interp, base_obj);
    mp_set_int(base, 1);

    ch = (zpos < chars) ? STRING_ord(interp, str, zpos) : 0;
    if ((flags & 0x02) && (ch == '+' || ch == '-')) {
        neg = (ch == '-');
        zpos++;
        ch = (zpos < chars) ? STRING_ord(interp, str, zpos) : 0;
    }
    while (zpos < chars) {
        if (ch >= '0' && ch <= '9') ch = ch - '0';
        else if (ch >= 'a' && ch <= 'z') ch = ch - 'a' + 10;
        else if (ch >= 'A' && ch <= 'Z') ch = ch - 'A' + 10;
        else break;
        if (ch >= radix) break;
        mp_mul_d(&zvalue, radix, &zvalue);
        mp_add_d(&zvalue, ch, &zvalue);
        mp_mul_d(&zbase, radix, &zbase);
        zpos++; pos = zpos;
        if (ch != 0 || !(flags & 0x04)) { mp_copy(&zvalue, value); mp_copy(&zbase, base); }
        if (zpos >= chars) break;
        ch = STRING_ord(interp, str, zpos);
        if (ch != '_') continue;
        zpos++;
        if (zpos >= chars) break;
        ch = STRING_ord(interp, str, zpos);
    }
    
    mp_clear(&zvalue);
    mp_clear(&zbase);
    
    pos_obj = REPR($6)->allocate(interp, STABLE($6));
    REPR(pos_obj)->initialize(interp, STABLE(pos_obj), OBJECT_BODY(pos_obj));
    REPR(pos_obj)->box_funcs->set_int(interp, STABLE(pos_obj), OBJECT_BODY(pos_obj), pos);

    if (neg || flags & 0x01) { mp_neg(value, value); }
    out = Parrot_pmc_new(interp, enum_class_FixedPMCArray);
    VTABLE_set_integer_native(interp, out, 3);
    VTABLE_set_pmc_keyed_int(interp, out, 0, value_obj);
    VTABLE_set_pmc_keyed_int(interp, out, 1, base_obj);
    VTABLE_set_pmc_keyed_int(interp, out, 2, pos_obj);
    $1 = out;
}

/* calculates $1 = $2 ** $3
 * if it either overflows ($3 being too big), or $2 is negative,
 * a float is returned. $4 should contain the type object to box the
 * float into, $5 the type object to box the bigint into.
 */
inline op nqp_bigint_pow(out PMC, in PMC, in PMC, in PMC, in PMC) :base_core {
    mp_digit exponent_d = 0;
    mp_int *exponent = get_bigint(interp, $3);
    mp_int *base     = get_bigint(interp, $2);
    int cmp          = mp_cmp_d(exponent, 0);
    if (cmp == MP_EQ || MP_EQ == mp_cmp_d(base, 1)) {
        /* $x ** 0 or 1 ** $x */
        $1 = REPR($5)->allocate(interp, STABLE($5));
        REPR($1)->initialize(interp, STABLE($1), OBJECT_BODY($1));
        mp_set_int(get_bigint(interp, $1), 1);
    }
    else if (cmp == MP_GT) {
        exponent_d = mp_get_int(exponent);
        if (MP_GT == mp_cmp_d(exponent, exponent_d)) {
            /* the exponent is larger than what fits into an int register...
             * that's scary, and should be treated with care */


            /* XXX a bit ugly that it reuses cmp, but safe for now */
            cmp = mp_cmp_d(base, 0)
            if (MP_EQ == cmp || MP_EQ == mp_cmp_d(base, 1)) {
                /* 0 ** $big_number and 1 ** big_number are easy to do: */
                $1 = REPR($2)->allocate(interp, STABLE($2));
                REPR($1)->initialize(interp, STABLE($1), OBJECT_BODY($1));
                mp_copy(base, get_bigint(interp, $1));
            }
            else {
                FLOATVAL ZERO = 0.0;
                $1 = REPR($4)->allocate(interp, STABLE($4));
                REPR($1)->initialize(interp, STABLE($1), OBJECT_BODY($1));
                /* TODO: better ways to create +- Inf */
                if (MP_GT == cmp) {
                    REPR($1)->box_funcs->set_num(interp, STABLE($1), OBJECT_BODY($1), (FLOATVAL) 1.0/ZERO);
                }
                else {
                    REPR($1)->box_funcs->set_num(interp, STABLE($1), OBJECT_BODY($1), (FLOATVAL) -1.0/ZERO);
                }
            }
        }
        else {
            /* since the exponent fits into a digit, mp_expt_d is fine */
            $1 = REPR($5)->allocate(interp, STABLE($5));
            REPR($1)->initialize(interp, STABLE($1), OBJECT_BODY($1));
            mp_expt_d(get_bigint(interp, $2), exponent_d, get_bigint(interp, $1));
        }
    }
    else {
        FLOATVAL f_base = mp_get_double(base);
        FLOATVAL f_exp  = mp_get_double(exponent);
        $1 = REPR($4)->allocate(interp, STABLE($4));
        REPR($1)->initialize(interp, STABLE($1), OBJECT_BODY($1));
        REPR($1)->box_funcs->set_num(interp, STABLE($1), OBJECT_BODY($1), pow(f_base, f_exp));
    }
}

/* returns 1 if $2 is too large to fit into an INTVAL without loss of
   information */
inline op nqp_bigint_is_big(out INT, in PMC) :base_core {
    mp_int *a = get_bigint(interp, $2);
    $1 = a->used > 1;
    /* XXX somebody please check that on a 32 bit platform */
    if ( sizeof(INTVAL) * 8 < DIGIT_BIT && $1 == 0 && DIGIT(a, 0) & ~0x7FFFFFFFUL)
        $1 = 1;
}