/*  arith12.c                         Copyright (C) 1990-2007 Codemist Ltd */

/*
 * Arithmetic functions... specials for Reduce (esp. factoriser)
 *
 */


/*
 * This code may be used and modified, and redistributed in binary
 * or source form, subject to the "CCL Public License", which should
 * accompany it. This license is a variant on the BSD license, and thus
 * permits use of code derived from this in either open and commercial
 * projects: but it does require that updates to this code be made
 * available back to the originators of the package.
 * Before merging other code in with this or linking this code
 * with other packages or libraries please check that the license terms
 * of the other material are compatible with those of this.
 */

/* Signature: 1114527e 18-Jan-2007 */


#include "headers.h"


#define FP_EVALUATE   1


Lisp_Object Lfrexp(Lisp_Object nil, Lisp_Object a)
{
    double d;
    int x;
    d = float_of_number(a);
    d = frexp(d, &x);
    if (d == 1.0) d = 0.5, x++;
    a = make_boxfloat(d, TYPE_DOUBLE_FLOAT);
    errexit();
    return Lcons(nil, fixnum_of_int((int32_t)x), a);
}

Lisp_Object Lmodular_difference(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
{
    int32_t r;
    CSL_IGNORE(nil);
    r = int_of_fixnum(a) - int_of_fixnum(b);
    if (r < 0) r += current_modulus;
    return onevalue(fixnum_of_int(r));
}

Lisp_Object Lmodular_minus(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    if (a != fixnum_of_int(0))
    {   int32_t r = current_modulus - int_of_fixnum(a);
        a = fixnum_of_int(r);
    }
    return onevalue(a);
}

Lisp_Object Lmodular_number(Lisp_Object nil, Lisp_Object a)
{
    int32_t r;
    a = Cremainder(a, fixnum_of_int(current_modulus));
    errexit();
    r = int_of_fixnum(a);
    if (r < 0) r += current_modulus;
    return onevalue(fixnum_of_int(r));
}

Lisp_Object Lmodular_plus(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
{
    int32_t r;
    CSL_IGNORE(nil);
    r = int_of_fixnum(a) + int_of_fixnum(b);
    if (r >= current_modulus) r -= current_modulus;
    return onevalue(fixnum_of_int(r));
}

Lisp_Object Lmodular_reciprocal(Lisp_Object nil, Lisp_Object n)
{
    int32_t a, b, x, y;
    CSL_IGNORE(nil);
    a = current_modulus;
    b = int_of_fixnum(n);
    x = 0;
    y = 1;
    if (b == 0) return aerror1("modular-reciprocal", n);
    if (b < 0) b = current_modulus - ((-b)%current_modulus);
    while (b != 1)
    {   int32_t w, t;
        if (b == 0)
            return aerror1("non-prime modulus in modular-reciprocal",
                           fixnum_of_int(current_modulus));
        w = a / b;
        t = b;
        b = a - b*w;
        a = t;
        t = y;
        y = x - y*w;
        x = t;
    }
    if (y < 0) y += current_modulus;
    return onevalue(fixnum_of_int(y));
}

Lisp_Object Lmodular_times(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
{
#ifndef HAVE_UINT64_T
    uint32_t h, l;
#endif
    uint32_t r, cm;
    int32_t aa, bb;
    CSL_IGNORE(nil);
    cm = (uint32_t)current_modulus;
    aa = int_of_fixnum(a);
    bb = int_of_fixnum(b);
/*
 * The constant 46341 is sqrt(2^31) suitable rounded - if my modulus
 * is no bigger than that then a and b will both be strictly smaller,
 * and hence a*b will be < 2^31 and hence in range for 32-bit signed
 * arithmetic.  I make this test because I expect Imultiply and Idivide
 * to be pretty painful, while regular C multiplication and division are
 * (probably!) much better.
 */
    if (cm <= 46341U) r = (aa * bb) % cm;
    else
    {
#ifdef HAVE_UINT64_T
        r = (uint32_t)(((uint64_t)aa * (uint64_t)bb) % (uint32_t)cm);
#else
        Dmultiply(h, l, aa, bb, 0);
        Ddivide(r, l, h, l, cm);
#endif
    }
    return onevalue(fixnum_of_int(r));
}

Lisp_Object Lmodular_quotient(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
{
    CSL_IGNORE(nil);
    push(a);
    b = Lmodular_reciprocal(nil, b);
    pop(a);
    errexit();
    return Lmodular_times(nil, a, b);
}

Lisp_Object Lmodular_expt(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
{
    int32_t x, r, p;
    uint32_t h, l;
    CSL_IGNORE(nil);
    x = int_of_fixnum(b);
    if (x == 0) return onevalue(fixnum_of_int(1));
    p = int_of_fixnum(a);
/*
 * I could play games here on half-length current_modulus and use
 * native C arithmetic, but I judge this case not to be quite that
 * critically important. Also on 64-bit machines I could do more
 * work in-line.
 */
    p = p % current_modulus; /* In case somebody is being silly! */
    while ((x & 1) == 0)
    {   Dmultiply(h, l, p, p, 0);
        Ddivide(p, l, h, l, current_modulus);
        x = x/2;
    }
    r = p;
    while (x != 1)
    {   Dmultiply(h, l, p, p, 0);
        Ddivide(p, l, h, l, current_modulus);
        x = x/2;
        if ((x & 1) != 0)
        {   Dmultiply(h, l, r, p, 0);
            Ddivide(r, l, h, l, current_modulus);
        }
    }
    return onevalue(fixnum_of_int(r));
}

Lisp_Object Lset_small_modulus(Lisp_Object nil, Lisp_Object a)
{
    int32_t r, old = current_modulus;
    CSL_IGNORE(nil);
    if (!is_fixnum(a)) return aerror1("set-small-modulus", a);
    r = int_of_fixnum(a);
/*
 * I COULD allow a small modulus of up to 2^27, but for compatibility
 * with Cambridge Lisp I will limit myself to 24 bits.
 */
    if (r > 0x00ffffff) return aerror1("set-small-modulus", a);
    current_modulus = r;
    return onevalue(fixnum_of_int(old));
}

Lisp_Object Liadd1(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    if (!is_fixnum(a)) return aerror1("iadd1", a);
    return onevalue((Lisp_Object)((int32_t)a + 0x10));
}

Lisp_Object Lidifference(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
{
    CSL_IGNORE(nil);
    if (!is_fixnum(a) || !is_fixnum(b)) return aerror2("idifference", a, b);
    return onevalue((Lisp_Object)((int32_t)a - (int32_t)b + TAG_FIXNUM));
}

/*
 * xdifference is provided just for the support of the CASE operator. It
 * subtracts its arguments but returns NIL if either argument is not
 * an small integer or if the result overflows.
 */

Lisp_Object Lxdifference(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
{
    int32_t r;
    if (!is_fixnum(a) || !is_fixnum(b)) return onevalue(nil);
    r = int_of_fixnum(a) - int_of_fixnum(b);
    if (r < -0x08000000 || r > 0x07ffffff) return onevalue(nil);
    return onevalue(fixnum_of_int(r));
}

Lisp_Object Ligreaterp(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
{
    CSL_IGNORE(nil);
    if (!is_fixnum(a) || !is_fixnum(b)) return aerror2("igreaterp", a, b);
    return onevalue(Lispify_predicate(a > b));
}

Lisp_Object Lilessp(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
{
    CSL_IGNORE(nil);
    if (!is_fixnum(a) || !is_fixnum(b)) return aerror2("ilessp", a, b);
    return onevalue(Lispify_predicate(a < b));
}

Lisp_Object Ligeq(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
{
    CSL_IGNORE(nil);
    if (!is_fixnum(a) || !is_fixnum(b)) return aerror2("igeq", a, b);
    return onevalue(Lispify_predicate(a >= b));
}

Lisp_Object Lileq(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
{
    CSL_IGNORE(nil);
    if (!is_fixnum(a) || !is_fixnum(b)) return aerror2("ileq", a, b);
    return onevalue(Lispify_predicate(a <= b));
}

static Lisp_Object MS_CDECL Lilogand(Lisp_Object nil, int nargs, ...)
{
    va_list a;
    Lisp_Object r;
    if (nargs == 0) return onevalue(fixnum_of_int(-1));
    if (nargs > ARG_CUT_OFF) return aerror("too many args for ilogand");
    CSL_IGNORE(nil);
    va_start(a, nargs);
    r = va_arg(a, Lisp_Object);
    if (!is_fixnum(r)) return aerror1("ilogand", r);
    while (--nargs != 0)
    {   Lisp_Object w = va_arg(a, Lisp_Object);
        if (!is_fixnum(w))
        {   va_end(a);
            return aerror1("ilogand", w);
        }
        r = (Lisp_Object)((int32_t)r & (int32_t)w);
    }
    va_end(a);
    return onevalue(r);
}

static Lisp_Object MS_CDECL Lilogor(Lisp_Object nil, int nargs, ...)
{
    va_list a;
    Lisp_Object r;
    if (nargs == 0) return onevalue(fixnum_of_int(0));
    if (nargs > ARG_CUT_OFF) return aerror("too many args for ilogor");
    CSL_IGNORE(nil);
    va_start(a, nargs);
    r = va_arg(a, Lisp_Object);
    if (!is_fixnum(r)) return aerror1("ilogor", r);
    while (--nargs != 0)
    {   Lisp_Object w = va_arg(a, Lisp_Object);
        if (!is_fixnum(w))
        {   va_end(a);
            return aerror1("ilogor", w);
        }
        r = (Lisp_Object)((int32_t)r | (int32_t)w);
    }
    va_end(a);
    return onevalue(r);
}

static Lisp_Object MS_CDECL Lilogxor(Lisp_Object nil, int nargs, ...)
{
    va_list a;
    Lisp_Object r;
    if (nargs == 0) return onevalue(fixnum_of_int(0));
    if (nargs > ARG_CUT_OFF) return aerror("too many args for ilogxor");
    CSL_IGNORE(nil);
    va_start(a, nargs);
    r = va_arg(a, Lisp_Object);
    if (!is_fixnum(r)) return aerror1("ilogxor", r);
    while (--nargs != 0)
    {   Lisp_Object w = va_arg(a, Lisp_Object);
        if (!is_fixnum(w))
        {   va_end(a);
            return aerror1("ilogxor", w);
        }
        r = (Lisp_Object)(((int32_t)r ^ (int32_t)w) + TAG_FIXNUM);
    }
    va_end(a);
    return onevalue(r);
}

static Lisp_Object Lilogand2(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
{
    CSL_IGNORE(nil);
    if (!is_fixnum(a) || !is_fixnum(b)) return aerror2("ilogand", a, b);
    return onevalue(a & b);    
}

static Lisp_Object Lilogor2(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
{
    CSL_IGNORE(nil);
    if (!is_fixnum(a) || !is_fixnum(b)) return aerror2("ilogor", a, b);
    return onevalue(a | b);    
}

static Lisp_Object Lilogxor2(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
{
    CSL_IGNORE(nil);
    if (!is_fixnum(a) || !is_fixnum(b)) return aerror2("ilogxor", a, b);
    return onevalue((a ^ b) + TAG_FIXNUM);    
}

Lisp_Object Limin(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
{
    CSL_IGNORE(nil);
    if (!is_fixnum(a) || !is_fixnum(b)) return aerror2("imin", a, b);
    return onevalue(a < b ? a : b);
}

Lisp_Object Limax(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
{
    CSL_IGNORE(nil);
    if (!is_fixnum(a) || !is_fixnum(b)) return aerror2("imax", a, b);
    return onevalue(a > b ? a : b);
}

Lisp_Object Liminus(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    if (!is_fixnum(a)) return aerror1("iminus", a);
    return onevalue((Lisp_Object)(2*TAG_FIXNUM - (int32_t)a));
}

Lisp_Object Liminusp(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    return onevalue(Lispify_predicate((int32_t)a < (int32_t)fixnum_of_int(0)));
}

static Lisp_Object MS_CDECL Liplus(Lisp_Object nil, int nargs, ...)
{
    va_list a;
    Lisp_Object r;
    if (nargs == 0) return onevalue(fixnum_of_int(0));
    if (nargs > ARG_CUT_OFF) return aerror("too many args for iplus");
    CSL_IGNORE(nil);
    va_start(a, nargs);
    r = va_arg(a, Lisp_Object);
    if (!is_fixnum(r)) return aerror1("iplus", r);
    while (--nargs != 0)
    {   Lisp_Object w = va_arg(a, Lisp_Object);
        if (!is_fixnum(w))
        {   va_end(a);
            return aerror1("iplus", w);
        }
        r = (Lisp_Object)((int32_t)r + (int32_t)w - TAG_FIXNUM);
    }
    va_end(a);
    return onevalue(r);
}

Lisp_Object Liplus2(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
{
    CSL_IGNORE(nil);
    if (!is_fixnum(a) || !is_fixnum(b)) return aerror2("iplus2", a, b);
    return onevalue((Lisp_Object)((int32_t)a + (int32_t)b - TAG_FIXNUM));
}

Lisp_Object Liquotient(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
{
    int32_t aa, bb, c;
    CSL_IGNORE(nil);
    if (!is_fixnum(a) || !is_fixnum(b) ||
        b == fixnum_of_int(0)) return aerror2("iquotient", a, b);
/* C does not define the exact behaviour of /, % on -ve args */
    aa = int_of_fixnum(a);
    bb = int_of_fixnum(b);
    c = aa % bb;
    if (aa < 0)
    {   if (c > 0) c -= bb;
    }
    else if (c < 0) c += bb;
    return onevalue(fixnum_of_int((aa-c)/bb));
}

Lisp_Object Liremainder(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
{
    int32_t aa, bb, c;
    CSL_IGNORE(nil);
    if (!is_fixnum(a) || !is_fixnum(b) ||
        b == fixnum_of_int(0)) return aerror2("iremainder", a, b);
/* C does not define the exact behaviour of /, % on -ve args */
    aa = int_of_fixnum(a);
    bb = int_of_fixnum(b);
    c = aa % bb;
    if (aa < 0)
    {   if (c > 0) c -= bb;
    }
    else if (c < 0) c += bb;
    return onevalue(fixnum_of_int(c));
}

Lisp_Object Lirightshift(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
{
    CSL_IGNORE(nil);
    if (!is_fixnum(a) || !is_fixnum(b)) return aerror2("irightshift", a, b);
    return onevalue(fixnum_of_int(int_of_fixnum(a) >> int_of_fixnum(b)));
}

Lisp_Object Lisub1(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    if (!is_fixnum(a)) return aerror1("isub1", a);
    return onevalue((Lisp_Object)((int32_t)a - 0x10));
}

static Lisp_Object MS_CDECL Litimes(Lisp_Object nil, int nargs, ...)
{
    va_list a;
    Lisp_Object rr;
    int32_t r;
    if (nargs == 0) return onevalue(fixnum_of_int(1));
    if (nargs > ARG_CUT_OFF) return aerror("too many args for itimes");
    CSL_IGNORE(nil);
    va_start(a, nargs);
    rr = va_arg(a, Lisp_Object);
    if (!is_fixnum(rr)) return aerror1("itimes", rr);
    r = int_of_fixnum(rr);
    while (nargs-- != 0)
    {   Lisp_Object w = va_arg(a, Lisp_Object);
        if (!is_fixnum(w))
        {   va_end(a);
            return aerror1("itimes", w);
        }
        r = r * int_of_fixnum(w);
    }
    va_end(a);
    return onevalue(fixnum_of_int(r));
}

Lisp_Object Litimes2(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
{
    CSL_IGNORE(nil);
    if (!is_fixnum(a) || !is_fixnum(b)) return aerror2("itimes2", a, b);
    return onevalue(fixnum_of_int(int_of_fixnum(a) * int_of_fixnum(b)));
}

Lisp_Object Lionep(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    return onevalue(Lispify_predicate((int32_t)a == (int32_t)fixnum_of_int(1)));
}

Lisp_Object Lizerop(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    return onevalue(Lispify_predicate((int32_t)a == (int32_t)fixnum_of_int(0)));
}

#ifdef FP_EVALUATE

static double fp_args[32];
static double fp_stack[16];

/* codes 0 to 31 just load up arguments */
#define FP_RETURN        32
#define FP_PLUS          33
#define FP_DIFFERENCE    34
#define FP_TIMES         35
#define FP_QUOTIENT      36
#define FP_MINUS         37
#define FP_SQUARE        38
#define FP_CUBE          39
#define FP_SIN           40
#define FP_COS           41
#define FP_TAN           42
#define FP_EXP           43
#define FP_LOG           44
#define FP_SQRT          45


static Lisp_Object Lfp_eval(Lisp_Object nil, Lisp_Object code,
                                             Lisp_Object args)
/*
 * The object of this code is to support fast evaluation of numeric
 * expressions.  The first argument is a vector of byte opcodes, while
 * the second arg is a list of floating point values whose value will (or
 * at least may) be used.  There are at most 32 values in this list.
 */
{
    int n = 0;
    double w;
    unsigned char *p;
    if (!is_vector(code)) return aerror("fp-evaluate");
    while (consp(args))
    {   fp_args[n++] = float_of_number(qcar(args));
        args = qcdr(args);
    }
    n = 0;
    p = &ucelt(code, 0);
    for (;;)
    {   int op = *p++;
/*
 * Opcodes 0 to 31 just load up the corresponding input value.
 */
        if (op < 32)
        {   fp_stack[n++] = fp_args[op];
            continue;
        }
        switch (op)
        {
    default:
            return aerror("Bad op in fp-evaluate");
    case FP_RETURN:
            args = make_boxfloat(fp_stack[0], TYPE_DOUBLE_FLOAT);
            errexit();
            return onevalue(args);
    case FP_PLUS:
            n--;
            fp_stack[n] += fp_stack[n-1];
            continue;
    case FP_DIFFERENCE:
            n--;
            fp_stack[n] -= fp_stack[n-1];
            continue;
    case FP_TIMES:
            n--;
            fp_stack[n] *= fp_stack[n-1];
            continue;
    case FP_QUOTIENT:
            n--;
            fp_stack[n] /= fp_stack[n-1];
            continue;
    case FP_MINUS:
            fp_stack[n] = -fp_stack[n];
            continue;
    case FP_SQUARE:
            fp_stack[n] *= fp_stack[n];
            continue;
    case FP_CUBE:
            w = fp_stack[n];
            w *= w;
            fp_stack[n] *= w;
            continue;
    case FP_SIN:
            fp_stack[n] = sin(fp_stack[n]);
            continue;
    case FP_COS:
            fp_stack[n] = cos(fp_stack[n]);
            continue;
    case FP_TAN:
            fp_stack[n] = tan(fp_stack[n]);
            continue;
    case FP_EXP:
            fp_stack[n] = exp(fp_stack[n]);
            continue;
    case FP_LOG:
            fp_stack[n] = log(fp_stack[n]);
            continue;
    case FP_SQRT:
            fp_stack[n] = sqrt(fp_stack[n]);
            continue;
        }
    }
}

#endif

setup_type const arith12_setup[] =
{
    {"frexp",                   Lfrexp, too_many_1, wrong_no_1},
    {"modular-difference",      too_few_2, Lmodular_difference, wrong_no_2},
    {"modular-minus",           Lmodular_minus, too_many_1, wrong_no_1},
    {"modular-number",          Lmodular_number, too_many_1, wrong_no_1},
    {"modular-plus",            too_few_2, Lmodular_plus, wrong_no_2},
    {"modular-quotient",        too_few_2, Lmodular_quotient, wrong_no_2},
    {"modular-reciprocal",      Lmodular_reciprocal, too_many_1, wrong_no_1},
    {"modular-times",           too_few_2, Lmodular_times, wrong_no_2},
    {"modular-expt",            too_few_2, Lmodular_expt, wrong_no_2},
    {"set-small-modulus",       Lset_small_modulus, too_many_1, wrong_no_1},
    {"iadd1",                   Liadd1, too_many_1, wrong_no_1},
    {"idifference",             too_few_2, Lidifference, wrong_no_2},
    {"xdifference",             too_few_2, Lxdifference, wrong_no_2},
    {"igeq",                    too_few_2, Ligeq, wrong_no_2},
    {"igreaterp",               too_few_2, Ligreaterp, wrong_no_2},
    {"ileq",                    too_few_2, Lileq, wrong_no_2},
    {"ilessp",                  too_few_2, Lilessp, wrong_no_2},
    {"ilogand",                 Lidentity, Lilogand2, Lilogand},
    {"ilogor",                  Lidentity, Lilogor2, Lilogor},
    {"ilogxor",                 Lidentity, Lilogxor2, Lilogxor},
    {"imax",                    too_few_2, Limax, wrong_no_2},
    {"imin",                    too_few_2, Limin, wrong_no_2},
    {"iminus",                  Liminus, too_many_1, wrong_no_1},
    {"iminusp",                 Liminusp, too_many_1, wrong_no_1},
    {"iplus",                   Lidentity, Liplus2, Liplus},
    {"iplus2",                  too_few_2, Liplus2, wrong_no_2},
    {"iquotient",               too_few_2, Liquotient, wrong_no_2},
    {"iremainder",              too_few_2, Liremainder, wrong_no_2},
    {"irightshift",             too_few_2, Lirightshift, wrong_no_2},
    {"isub1",                   Lisub1, too_many_1, wrong_no_1},
    {"itimes",                  Lidentity, Litimes2, Litimes},
    {"itimes2",                 too_few_2, Litimes2, wrong_no_2},
    {"ionep",                   Lionep, too_many_1, wrong_no_1},
    {"izerop",                  Lizerop, too_many_1, wrong_no_1},
#ifdef FP_EVALUATE
    {"fp-evaluate",             too_few_2, Lfp_eval, wrong_no_2},
#endif
    {NULL,                      0, 0, 0}
};

/* end of arith12.c */