glisp/src/lisp_math.c

#include "lisp_math.h"

#include "list.h"
#include "stack.h"

#include <errno.h>
#include <stdckdint.h>

LispVal *copy_lisp_gmp(LispVal *gmp) {
    assert(LISP_GMP_P(gmp));
    LispGmp *n = lisp_alloc_object(sizeof(LispGmp), TYPE_GMP);
    mpz_init_set(n->val, ((LispGmp *) gmp)->val);
    return n;
}

LispVal *make_gmp(fixnum_t value) {
    LispGmp *n = lisp_alloc_object(sizeof(LispGmp), TYPE_GMP);
    mpz_init_set_si(n->val, value);
    return n;
}

LispVal *make_gmp_unsigned(unsigned long value) {
    LispGmp *n = lisp_alloc_object(sizeof(LispGmp), TYPE_GMP);
    mpz_init_set_ui(n->val, value);
    return n;
}

LispVal *parse_gmp(const char *str, int base) {
    LispGmp *n = lisp_alloc_object(sizeof(LispGmp), TYPE_GMP);
    mpz_init(n->val);
    int res = mpz_set_str(n->val, str, base);
    if (res != 0) {
        free(n);
        return NULL;
    }
    return n;
}

LispVal *parse_number(const char *str, int base) {
    char *endptr;
    errno = 0;
    intmax_t n = strtoimax(str, &endptr, base);
    if ((n == INTMAX_MIN || n == INTMAX_MAX) && errno == ERANGE) {
        return parse_gmp(str, base);
    } else if (n < MOST_NEGATIVE_FIXNUM || n > MOST_POSITIVE_FIXNUM) {
        return parse_gmp(str, base);
    } else if (*endptr) { // trailing garbage
        return NULL;
    }
    return MAKE_FIXNUM(n);
}

DEFUN(fixnump, "fixnump", (LispVal * val), "(val)", "") {
    return FIXNUMP(val) ? Qt : Qnil;
}

DEFUN(integerp, "integerp", (LispVal * val), "(val)", "") {
    return LISP_GMP_P(val) ? Qt : Qnil;
}

DEFUN(floatp, "floatp", (LispVal * val), "(val)", "") {
    return LISP_FLOAT_P(val) ? Qt : Qnil;
}

DEFUN(numberp, "numberp", (LispVal * val), "(val)", "") {
    return FIXNUMP(val) || LISP_FLOAT_P(val) || LISP_GMP_P(val) ? Qt : Qnil;
}

enum Operator {
    OPERATOR_ADD,
    OPERATOR_SUB,
    OPERATOR_MUL,
    OPERATOR_DIV,
    OPERATOR_AND,
    OPERATOR_IOR,
    OPERATOR_XOR,
};

static inline LispVal *type_for_operator(enum Operator op) {
    switch (op) {
    case OPERATOR_ADD:
    case OPERATOR_SUB:
    case OPERATOR_MUL:
    case OPERATOR_DIV:
        return Qnumber;
    case OPERATOR_AND:
    case OPERATOR_IOR:
    case OPERATOR_XOR:
        return Qinteger;
    default:
        abort();
    }
}

static LispVal *float_math_driver(enum Operator op, bool first,
                                  lisp_float_t acc, LispVal *nums) {
    assert(!first || LISP_FLOAT_P(XCAR(nums)));
    switch (op) {
    case OPERATOR_AND:
    case OPERATOR_IOR:
    case OPERATOR_XOR:
        signal_type_error(XCAR(nums), Qinteger);
        break;
    default:
        break;
    }
    if (first) {
        acc = XLISP_FLOAT(XCAR(nums));
        nums = XCDR(nums);
    }
    DOLIST(cur, nums) {
        lisp_float_t n;
        if (FIXNUMP(cur)) {
            n = XFIXNUM(cur);
        } else if (LISP_FLOAT_P(cur)) {
            n = XLISP_FLOAT(cur);
        } else if (LISP_GMP_P(cur)) {
            n = mpz_get_d(((LispGmp *) cur)->val);
        } else {
            signal_type_error(cur, Qnumber);
        }
        switch (op) {
        case OPERATOR_ADD:
            acc += n;
            break;
        case OPERATOR_SUB:
            acc -= n;
            break;
        case OPERATOR_MUL:
            acc *= n;
            break;
        case OPERATOR_DIV:
            acc /= n;
            break;
        default:
            abort();
        }
    }
    return MAKE_LISP_FLOAT(acc);
}

static LispVal *gmp_math_driver(enum Operator op, bool first,
                                fixnum_t fixnum_acc, LispVal *nums) {
    assert(!first || LISP_GMP_P(XCAR(nums)));
    LispGmp *acc = make_gmp(fixnum_acc);
    if (first) {
        mpz_set(acc->val, ((LispGmp *) XCAR(nums))->val);
        nums = XCDR(nums);
        first = false;
    }
    DOTAILS(rest, nums) {
        LispVal *cur = XCAR(rest);
        if (FIXNUMP(cur)) {
            fixnum_t f = XFIXNUM(cur);
            switch (op) {
            case OPERATOR_ADD:
                if (f < 0) {
                    mpz_sub_ui(acc->val, acc->val, -f);
                } else {
                    mpz_add_ui(acc->val, acc->val, f);
                }
                break;
            case OPERATOR_SUB:
                if (f < 0) {
                    mpz_add_ui(acc->val, acc->val, -f);
                } else {
                    mpz_sub_ui(acc->val, acc->val, f);
                }
                break;
            case OPERATOR_MUL:
                mpz_mul_si(acc->val, acc->val, f);
                break;
            case OPERATOR_DIV:
                if (f == 0) {
                    lisp_signal(Qdivision_by_zero, Qnil);
                }
                if (f > 0) {
                    mpz_mul_si(acc->val, acc->val, -1);
                }
                mpz_tdiv_ui(acc->val, f);
                break;
#define BITWISE(o)                      \
    {                                   \
        mpz_t n;                        \
        mpz_init_set_si(n, f);          \
        mpz_##o(acc->val, acc->val, n); \
        mpz_clear(n);                   \
    }
            case OPERATOR_AND:
                BITWISE(and);
                break;
            case OPERATOR_IOR:
                BITWISE(ior);
                break;
            case OPERATOR_XOR:
                BITWISE(xor);
                break;
#undef BITWISE
            default:
                abort();
            }
        } else if (LISP_GMP_P(cur)) {
            mpz_srcptr cv = ((LispGmp *) cur)->val;
            switch (op) {
            case OPERATOR_ADD:
                mpz_add(acc->val, acc->val, cv);
                break;
            case OPERATOR_SUB:
                mpz_sub(acc->val, acc->val, cv);
                break;
            case OPERATOR_MUL:
                mpz_mul(acc->val, acc->val, cv);
                break;
            case OPERATOR_DIV:
                mpz_tdiv_q(acc->val, acc->val, cv);
                break;
            case OPERATOR_AND:
                mpz_and(acc->val, acc->val, cv);
                break;
            case OPERATOR_IOR:
                mpz_ior(acc->val, acc->val, cv);
                break;
            case OPERATOR_XOR:
                mpz_xor(acc->val, acc->val, cv);
                break;
            default:
                abort();
            }
        } else if (LISP_FLOAT_P(cur)) {
            return float_math_driver(op, first, mpz_get_d(acc->val), rest);
        } else {
            signal_type_error(cur, type_for_operator(op));
        }
    }
    if (mpz_cmp_si(acc->val, MOST_POSITIVE_FIXNUM) <= 0
        && mpz_cmp_si(acc->val, MOST_NEGATIVE_FIXNUM) >= 0) {
        return MAKE_FIXNUM(mpz_get_si(acc->val));
    }
    return acc;
}

static LispVal *math_driver(enum Operator op, LispVal *nums) {
    assert(!NILP(nums));
    bool first = true;
    fixnum_t acc;
    LispVal *rest;
    if (FIXNUMP(XCAR(nums))) {
        rest = XCDR(nums);
        acc = XFIXNUM(XCAR(nums));
        first = false;
    } else {
        rest = nums;
        acc = 0;
        goto change_type;
    }
    while (FIXNUMP(XCAR(rest))) {
        bool did_overflow = false;
        fixnum_t cur = XFIXNUM(XCAR(rest));
        fixnum_t out;
        switch (op) {
        case OPERATOR_ADD:
            did_overflow = ckd_add(&out, acc, cur);
            break;
        case OPERATOR_SUB:
            did_overflow = ckd_sub(&out, acc, cur);
            break;
        case OPERATOR_MUL:
            did_overflow = ckd_mul(&out, acc, cur);
            break;
        case OPERATOR_DIV:
            if (cur == 0) {
                lisp_signal(Qdivision_by_zero, Qnil);
            }
            out = acc / cur;
            break;
        case OPERATOR_AND:
            out = acc & cur;
            break;
        case OPERATOR_IOR:
            out = acc | cur;
            break;
        case OPERATOR_XOR:
            out = acc ^ cur;
            break;
        default:
            abort();
        }
        if (did_overflow || out > MOST_POSITIVE_FIXNUM
            || out < MOST_NEGATIVE_FIXNUM) {
            break;
        }
        acc = out;
        rest = XCDR(rest);
        if (NILP(rest)) {
            return MAKE_FIXNUM(acc);
        }
    }
change_type:
    return LISP_FLOAT_P(XCAR(rest)) ? float_math_driver(op, first, acc, rest)
                                    : gmp_math_driver(op, first, acc, rest);
}

DEFUN(plus, "+", (LispVal * nums), "(&rest nums)", "") {
    if (NILP(nums)) {
        return MAKE_FIXNUM(0);
    }
    return math_driver(OPERATOR_ADD, nums);
}

DEFUN(minus, "-", (LispVal * nums), "(&rest nums)", "") {
    if (NILP(nums)) {
        return MAKE_FIXNUM(0);
    } else if (NILP(XCDR(nums))) {
        LispVal *num = XCAR(nums);
        if (FIXNUMP(num)) {
            return MAKE_FIXNUM(-XFIXNUM(num));
        } else if (LISP_FLOAT_P(num)) {
            return MAKE_LISP_FLOAT(-XLISP_FLOAT(num));
        } else if (LISP_GMP_P(num)) {
            LispGmp *o = make_gmp(0);
            LispGmp *a = num;
            mpz_mul_si(o->val, a->val, -1);
            return o;
        }
        signal_type_error(num, Qnumber);
    }
    return math_driver(OPERATOR_SUB, nums);
}

DEFUN(times, "*", (LispVal * nums), "(&rest nums)", "") {
    if (NILP(nums)) {
        return MAKE_FIXNUM(1);
    }
    return math_driver(OPERATOR_MUL, nums);
}

DEFUN(divide, "/", (LispVal * num, LispVal *nums), "(num &rest nums)", "") {
    if (NILP(nums)) {
        if (FIXNUMP(num)) {
            fixnum_t fv = XFIXNUM(num);
            if (fv == 0) {
                lisp_signal(Qdivision_by_zero, Qnil);
            } else if (fv == 1) {
                return MAKE_FIXNUM(1);
            } else if (fv == -1) {
                return MAKE_FIXNUM(-1);
            }
            return MAKE_FIXNUM(0);
        } else if (LISP_GMP_P(num)) {
            mpz_srcptr n = ((LispGmp *) num)->val;
            if (mpz_cmp_ui(n, 0) == 0) {
                lisp_signal(Qdivision_by_zero, Qnil);
            } else if (mpz_cmp_ui(n, 1) == 0) {
                return MAKE_FIXNUM(1);
            } else if (mpz_cmp_ui(n, -1) == 0) {
                return MAKE_FIXNUM(-1);
            }
            return MAKE_FIXNUM(0);
        }
    }
    return math_driver(OPERATOR_DIV, CONS(num, nums));
}

DEFUN(logand, "logand", (LispVal * nums), "(&rest nums)", "") {
    if (NILP(nums)) {
        return MAKE_FIXNUM(-1);
    }
    return math_driver(OPERATOR_AND, nums);
}

DEFUN(logior, "logior", (LispVal * nums), "(&rest nums)", "") {
    if (NILP(nums)) {
        return MAKE_FIXNUM(0);
    }
    return math_driver(OPERATOR_IOR, nums);
}

DEFUN(logxor, "logxor", (LispVal * nums), "(&rest nums)", "") {
    if (NILP(nums)) {
        return MAKE_FIXNUM(0);
    }
    return math_driver(OPERATOR_XOR, nums);
}

enum CompResult {
    CR_UNORDERED = 0,
    CR_LT = 1,
    CR_EQ = 2,
    CR_GT = 4,
};

enum Comparison {
    COMP_EQ = CR_EQ,
    COMP_GT = CR_GT,
    COMP_LT = CR_LT,
    COMP_GE = CR_GT | CR_EQ,
    COMP_LE = CR_LT | CR_EQ,
};

static inline enum CompResult do_float_comparison(lisp_float_t f1,
                                                  LispVal *cur) {
    if (isnan(f1)) {
        return CR_UNORDERED;
    }
    if (LISP_FLOAT_P(cur)) {
        lisp_float_t f2 = XLISP_FLOAT(cur);
        if (isnan(f2)) {
            return CR_UNORDERED;
        }
        return f1 < f2 ? CR_LT : (f1 > f2 ? CR_GT : CR_EQ);
    } else if (FIXNUMP(cur)) {
        fixnum_t fn = XFIXNUM(cur);
        return f1 < fn ? CR_LT : (f1 > fn ? CR_GT : CR_EQ);
    } else if (LISP_GMP_P(cur)) {
        mpz_srcptr g = ((LispGmp *) cur)->val;
        switch (mpz_cmp_d(g, f1)) {
        case -1:
            return CR_GT;
        case 0:
            return CR_EQ;
        case 1:
            return CR_LT;
        default:
            abort();
        }
    }
    signal_type_error(cur, Qnumber);
}

static inline enum CompResult do_gmp_comparison(mpz_srcptr n1, LispVal *cur) {
    int res;
    if (LISP_FLOAT_P(cur)) {
        lisp_float_t f = XLISP_FLOAT(cur);
        if (isnan(f)) {
            return CR_UNORDERED;
        }
        res = mpz_cmp_d(n1, f);
    } else if (FIXNUMP(cur)) {
        res = mpz_cmp_si(n1, XFIXNUM(cur));
    } else if (LISP_GMP_P(cur)) {
        mpz_srcptr n2 = ((LispGmp *) cur)->val;
        res = mpz_cmp(n1, n2);
    } else {
        signal_type_error(cur, Qnumber);
    }
    switch (res) {
    case -1:
        return CR_LT;
    case 0:
        return CR_EQ;
    case 1:
        return CR_GT;
    default:
        abort();
    }
}

static inline enum CompResult do_fixnum_comparison(fixnum_t n1, LispVal *cur) {
    if (LISP_FLOAT_P(cur)) {
        lisp_float_t n2 = XLISP_FLOAT(cur);
        if (isnan(n2)) {
            return CR_UNORDERED;
        }
        return n1 < n2 ? CR_LT : (n1 > n2 ? CR_GT : CR_EQ);
    } else if (FIXNUMP(cur)) {
        fixnum_t n2 = XFIXNUM(cur);
        return n1 < n2 ? CR_LT : (n1 > n2 ? CR_GT : CR_EQ);
    } else if (LISP_GMP_P(cur)) {
        mpz_srcptr g = ((LispGmp *) cur)->val;
        switch (mpz_cmp_si(g, n1)) {
        case -1:
            return CR_GT;
        case 0:
            return CR_EQ;
        case 1:
            return CR_LT;
        default:
            abort();
        }
    }
    signal_type_error(cur, Qnumber);
}

static bool comp_driver(enum Comparison c, LispVal *first, LispVal *rest) {
    if (NILP(rest)) {
        return true;
    }
    DOLIST(cur, rest) {
        enum CompResult res;
        if (LISP_FLOAT_P(first)) {
            res = do_float_comparison(XLISP_FLOAT(first), cur);
        } else if (LISP_GMP_P(first)) {
            res = do_gmp_comparison(((LispGmp *) first)->val, cur);
        } else if (FIXNUMP(first)) {
            res = do_fixnum_comparison(XFIXNUM(first), cur);
        } else {
            signal_type_error(first, Qnumber);
        }
        if ((c & res) == 0) {
            return false;
        }
        first = cur;
    }
    return true;
}

DEFUN(math_eq, "=", (LispVal * num, LispVal *nums), "(num &rest nums)", "") {
    return comp_driver(COMP_EQ, num, nums) ? Qt : Qnil;
}

DEFUN(gt, ">", (LispVal * num, LispVal *nums), "(num &rest nums)", "") {
    return comp_driver(COMP_GT, num, nums) ? Qt : Qnil;
}

DEFUN(lt, "<", (LispVal * num, LispVal *nums), "(num &rest nums)", "") {
    return comp_driver(COMP_LT, num, nums) ? Qt : Qnil;
}

DEFUN(ge, ">=", (LispVal * num, LispVal *nums), "(num &rest nums)", "") {
    return comp_driver(COMP_GE, num, nums) ? Qt : Qnil;
}

DEFUN(le, "<=", (LispVal * num, LispVal *nums), "(num &rest nums)", "") {
    return comp_driver(COMP_LE, num, nums) ? Qt : Qnil;
}

DEFINE_SYMBOL(integer, "integer");
DEFINE_SYMBOL(number, "number");

DEFINE_SYMBOL(division_by_zero, "division-by-zero");
DEFINE_CONDITION_CLASS(division_by_zero, error);