path: root/lisp/mp/mpr.c

/*
 * Copyright (c) 2002 by The XFree86 Project, Inc.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a
 * copy of this software and associated documentation files (the "Software"),
 * to deal in the Software without restriction, including without limitation
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
 * and/or sell copies of the Software, and to permit persons to whom the
 * Software is furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *  
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL
 * THE XFREE86 PROJECT BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
 * WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF
 * OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 *
 * Except as contained in this notice, the name of the XFree86 Project shall
 * not be used in advertising or otherwise to promote the sale, use or other
 * dealings in this Software without prior written authorization from the
 * XFree86 Project.
 *
 * Author: Paulo César Pereira de Andrade
 */

/* $XFree86$ */

#include "mp.h"

/*
 * TODO:
 *  Implement a fast gcd and divexact for integers, so that this code
 * could be changed to do intermediary calculations faster using smaller
 * numbers.
 */

/*
 * Prototypes
 */
	/* do the hard work of mpr_add and mpr_sub */
static void mpr_addsub(mpr *rop, mpr *op1, mpr *op2, int sub);

	/* do the hard work of mpr_cmp and mpr_cmpabs */
static int mpr_docmp(mpr *op1, mpr *op2, int sign);

/*
 * Implementation
 */
void
mpr_init(mpr *op)
{
    op->num.digs = mp_malloc(sizeof(BNS));
    op->num.sign = 0;
    op->num.size = op->num.alloc = 1;
    op->num.digs[0] = 0;

    op->den.digs = mp_malloc(sizeof(BNS));
    op->den.sign = 0;
    op->den.size = op->den.alloc = 1;
    op->den.digs[0] = 1;
}

void
mpr_clear(mpr *op)
{
    op->num.sign = 0;
    op->num.size = op->num.alloc = 0;
    mp_free(op->num.digs);

    op->den.sign = 0;
    op->den.size = op->den.alloc = 0;
    mp_free(op->den.digs);
}

void
mpr_set(mpr *rop, mpr *op)
{
    if (rop != op) {
	mpi_set(mpr_num(rop), mpr_num(op));
	mpi_set(mpr_den(rop), mpr_den(op));
    }
}

void
mpr_seti(mpr *rop, long num, long den)
{
    mpi_seti(mpr_num(rop), num);
    mpi_seti(mpr_den(rop), den);
}

void
mpr_setd(mpr *rop, double d)
{
    double val, num;
    int e, sign;

    /* initialize */
    if (d < 0) {
	sign = 1;
	val = -d;
    }
    else {
	sign = 0;
	val = d;
    }

    val = frexp(val, &e);
    while (modf(val, &num) != 0.0 && val <= DBL_MAX / 2.0) {
	--e;
	val *= 2.0;
    }
    if (e >= 0) {
	mpi_setd(mpr_num(rop), d);
	mpi_seti(mpr_den(rop), 1);
    }
    else {
	mpi_setd(mpr_num(rop), sign ? -num : num);
	mpi_setd(mpr_den(rop), ldexp(1.0, -e));
    }
}

void
mpr_setstr(mpr *rop, char *str, int base)
{
    char *slash = strchr(str, '/');

    mpi_setstr(mpr_num(rop), str, base);
    if (slash != NULL)
	mpi_setstr(mpr_den(rop), slash + 1, base);
    else
	mpi_seti(mpr_den(rop), 1);
}

void
mpr_canonicalize(mpr *op)
{
    mpi gcd;

    memset(&gcd, '\0', sizeof(mpi));

    mpi_gcd(&gcd, mpr_num(op), mpr_den(op));
    if (mpi_cmpabsi(&gcd, 1)) {
	mpi_div(mpr_num(op), mpr_num(op), &gcd);
	mpi_div(mpr_den(op), mpr_den(op), &gcd);
    }

    if (op->den.sign) {
	op->num.sign = !op->num.sign;
	op->den.sign = 0;
    }

    mpi_clear(&gcd);
}

void
mpr_add(mpr *rop, mpr *op1, mpr *op2)
{
    mpr_addsub(rop, op1, op2, 0);
}

void
mpr_addi(mpr *rop, mpr *op1, long op2)
{
    mpi prod;

    memset(&prod, '\0', sizeof(mpi));

    mpi_muli(&prod, mpr_den(op1), op2);
    mpi_add(mpr_num(rop), mpr_num(op1), &prod);
    mpi_clear(&prod);
}

void
mpr_sub(mpr *rop, mpr *op1, mpr *op2)
{
    mpr_addsub(rop, op1, op2, 1);
}

void
mpr_subi(mpr *rop, mpr *op1, long op2)
{
    mpi prod;

    memset(&prod, '\0', sizeof(mpi));

    mpi_muli(&prod, mpr_den(op1), op2);
    mpi_sub(mpr_num(rop), mpr_num(op1), &prod);
    mpi_clear(&prod);
}

static void
mpr_addsub(mpr *rop, mpr *op1, mpr *op2, int sub)
{
    mpi prod1, prod2;

    memset(&prod1, '\0', sizeof(mpi));
    memset(&prod2, '\0', sizeof(mpi));

    mpi_mul(&prod1, mpr_num(op1), mpr_den(op2));
    mpi_mul(&prod2, mpr_num(op2), mpr_den(op1));

    if (sub)
	mpi_sub(mpr_num(rop), &prod1, &prod2);
    else
	mpi_add(mpr_num(rop), &prod1, &prod2);

    mpi_clear(&prod1);
    mpi_clear(&prod2);

    mpi_mul(mpr_den(rop), mpr_den(op1), mpr_den(op2));
}

void
mpr_mul(mpr *rop, mpr *op1, mpr *op2)
{
    /* check if temporary storage is required */
    if (op1 == op2 && rop == op1) {
	mpi prod;

	memset(&prod, '\0', sizeof(mpi));

	mpi_mul(&prod, mpr_num(op1), mpr_num(op2));
	mpi_mul(mpr_den(rop), mpr_den(op1), mpr_den(op2));
	mpi_set(mpr_num(rop), &prod);

	mpi_clear(&prod);
    }
    else {
	mpi_mul(mpr_num(rop), mpr_num(op1), mpr_num(op2));
	mpi_mul(mpr_den(rop), mpr_den(op1), mpr_den(op2));
    }
}

void
mpr_muli(mpr *rop, mpr *op1, long op2)
{
    mpi_muli(mpr_num(rop), mpr_num(op1), op2);
}

void
mpr_div(mpr *rop, mpr *op1, mpr *op2)
{
    /* check if temporary storage is required */
    if (op1 == op2 && rop == op1) {
	mpi prod;

	memset(&prod, '\0', sizeof(mpi));

	mpi_mul(&prod, mpr_num(op1), mpr_den(op2));
	mpi_mul(mpr_den(rop), mpr_num(op2), mpr_den(op1));
	mpi_set(mpr_num(rop), &prod);

	mpi_clear(&prod);
    }
    else {
	mpi_mul(mpr_num(rop), mpr_num(op1), mpr_den(op2));
	mpi_mul(mpr_den(rop), mpr_num(op2), mpr_den(op1));
    }
}

void
mpr_divi(mpr *rop, mpr *op1, long op2)
{
    mpi_muli(mpr_den(rop), mpr_den(op1), op2);
}

void
mpr_inv(mpr *rop, mpr *op)
{
    if (rop == op)
	mpi_swap(mpr_num(op), mpr_den(op));
    else {
	mpi_set(mpr_num(rop), mpr_den(op));
	mpi_set(mpr_den(rop), mpr_num(op));
    }
}

void
mpr_neg(mpr *rop, mpr *op)
{
    mpi_neg(mpr_num(rop), mpr_num(op));
    mpi_set(mpr_den(rop), mpr_den(op));
}

void
mpr_abs(mpr *rop, mpr *op)
{
    if (mpr_num(op)->sign)
	mpi_neg(mpr_num(rop), mpr_num(op));
    else
	mpi_set(mpr_num(rop), mpr_num(op));

    /* op may not be canonicalized */
    if (mpr_den(op)->sign)
	mpi_neg(mpr_den(rop), mpr_den(op));
    else
	mpi_set(mpr_den(rop), mpr_den(op));
}

int
mpr_cmp(mpr *op1, mpr *op2)
{
    return (mpr_docmp(op1, op2, 1));
}

int
mpr_cmpi(mpr *op1, long op2)
{
    int cmp;
    mpr rat;

    mpr_init(&rat);
    mpi_seti(mpr_num(&rat), op2);
    cmp = mpr_docmp(op1, &rat, 1);
    mpr_clear(&rat);

    return (cmp);
}

int
mpr_cmpabs(mpr *op1, mpr *op2)
{
    return (mpr_docmp(op1, op2, 0));
}

int
mpr_cmpabsi(mpr *op1, long op2)
{
    int cmp;
    mpr rat;

    mpr_init(&rat);
    mpi_seti(mpr_num(&rat), op2);
    cmp = mpr_docmp(op1, &rat, 0);
    mpr_clear(&rat);

    return (cmp);
}

static int
mpr_docmp(mpr *op1, mpr *op2, int sign)
{
    int cmp, neg;
    mpi prod1, prod2;

    neg = 0;
    if (sign) {
	/* if op1 is negative */
	if (mpr_num(op1)->sign ^ mpr_den(op1)->sign) {
	    /* if op2 is positive */
	    if (!(mpr_num(op2)->sign ^ mpr_den(op2)->sign))
		return (-1);
	    else
		neg = 1;
	}
	/* if op2 is negative */
	else if (mpr_num(op2)->sign ^ mpr_den(op2)->sign)
	    return (1);
	/* else same sign */
    }

    /* if denominators are equal, compare numerators */
    if (mpi_cmpabs(mpr_den(op1), mpr_den(op2)) == 0) {
	cmp = mpi_cmpabs(mpr_num(op1), mpr_num(op2));
	if (cmp == 0)
	    return (0);
	if (sign && neg)
	    return (cmp < 0 ? 1 : -1);
	return (cmp);
    }

    memset(&prod1, '\0', sizeof(mpi));
    memset(&prod2, '\0', sizeof(mpi));

    /* "divide" op1 by op2
     * if result is smaller than 1, op1 is smaller than op2 */
    mpi_mul(&prod1, mpr_num(op1), mpr_den(op2));
    mpi_mul(&prod2, mpr_num(op2), mpr_den(op1));

    cmp = mpi_cmpabs(&prod1, &prod2);

    mpi_clear(&prod1);
    mpi_clear(&prod2);

    if (sign && neg)
	return (cmp < 0 ? 1 : -1);
    return (cmp);
}

void
mpr_swap(mpr *op1, mpr *op2)
{
    if (op1 != op2) {
	mpr swap;

	memcpy(&swap, op1, sizeof(mpr));
	memcpy(op1, op2, sizeof(mpr));
	memcpy(op2, &swap, sizeof(mpr));
    }
}

int
mpr_fiti(mpr *op)
{
    return (mpi_fiti(mpr_num(op)) && mpi_fiti(mpr_den(op)));
}

double
mpr_getd(mpr *op)
{
    return (mpi_getd(mpr_num(op)) / mpi_getd(mpr_den(op)));
}

char *
mpr_getstr(char *str, mpr *op, int base)
{
    int len;

    if (str == NULL) {
	len = mpi_getsize(mpr_num(op), base) + mpr_num(op)->sign + 1 +
	      mpi_getsize(mpr_den(op), base) + mpr_den(op)->sign + 1;

	str = mp_malloc(len);
    }

    (void)mpi_getstr(str, mpr_num(op), base);
    len = strlen(str);
    str[len] = '/';
    (void)mpi_getstr(str + len + 1, mpr_den(op), base);

    return (str);
}