437 lines
8.7 KiB
C
437 lines
8.7 KiB
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);
|
|
}
|