Blame isl-0.14/imath/imrat.c

Packit fb9d21
/*
Packit fb9d21
  Name:     imrat.c
Packit fb9d21
  Purpose:  Arbitrary precision rational arithmetic routines.
Packit fb9d21
  Author:   M. J. Fromberger <http://spinning-yarns.org/michael/>
Packit fb9d21
Packit fb9d21
  Copyright (C) 2002-2007 Michael J. Fromberger, All Rights Reserved.
Packit fb9d21
Packit fb9d21
  Permission is hereby granted, free of charge, to any person obtaining a copy
Packit fb9d21
  of this software and associated documentation files (the "Software"), to deal
Packit fb9d21
  in the Software without restriction, including without limitation the rights
Packit fb9d21
  to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
Packit fb9d21
  copies of the Software, and to permit persons to whom the Software is
Packit fb9d21
  furnished to do so, subject to the following conditions:
Packit fb9d21
Packit fb9d21
  The above copyright notice and this permission notice shall be included in
Packit fb9d21
  all copies or substantial portions of the Software.
Packit fb9d21
Packit fb9d21
  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
Packit fb9d21
  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
Packit fb9d21
  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL THE
Packit fb9d21
  AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
Packit fb9d21
  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
Packit fb9d21
  OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
Packit fb9d21
  SOFTWARE.
Packit fb9d21
 */
Packit fb9d21
Packit fb9d21
#include "imrat.h"
Packit fb9d21
#include <stdlib.h>
Packit fb9d21
#include <string.h>
Packit fb9d21
#include <ctype.h>
Packit fb9d21
#include <assert.h>
Packit fb9d21
Packit fb9d21
/* {{{ Useful macros */
Packit fb9d21
Packit fb9d21
#define TEMP(K) (temp + (K))
Packit fb9d21
#define SETUP(E, C) \
Packit fb9d21
do{if((res = (E)) != MP_OK) goto CLEANUP; ++(C);}while(0)
Packit fb9d21
Packit fb9d21
/* Argument checking:
Packit fb9d21
   Use CHECK() where a return value is required; NRCHECK() elsewhere */
Packit fb9d21
#define CHECK(TEST)   assert(TEST)
Packit fb9d21
#define NRCHECK(TEST) assert(TEST)
Packit fb9d21
Packit fb9d21
/* }}} */
Packit fb9d21
Packit fb9d21
/* Reduce the given rational, in place, to lowest terms and canonical form.
Packit fb9d21
   Zero is represented as 0/1, one as 1/1.  Signs are adjusted so that the sign
Packit fb9d21
   of the numerator is definitive. */
Packit fb9d21
static mp_result s_rat_reduce(mp_rat r);
Packit fb9d21
Packit fb9d21
/* Common code for addition and subtraction operations on rationals. */
Packit fb9d21
static mp_result s_rat_combine(mp_rat a, mp_rat b, mp_rat c,
Packit fb9d21
			       mp_result (*comb_f)(mp_int, mp_int, mp_int));
Packit fb9d21
Packit fb9d21
/* {{{ mp_rat_init(r) */
Packit fb9d21
Packit fb9d21
mp_result mp_rat_init(mp_rat r)
Packit fb9d21
{
Packit fb9d21
  mp_result res;
Packit fb9d21
Packit fb9d21
  if ((res = mp_int_init(MP_NUMER_P(r))) != MP_OK)
Packit fb9d21
    return res;
Packit fb9d21
  if ((res = mp_int_init(MP_DENOM_P(r))) != MP_OK) {
Packit fb9d21
    mp_int_clear(MP_NUMER_P(r));
Packit fb9d21
    return res;
Packit fb9d21
  }
Packit fb9d21
Packit fb9d21
  return mp_int_set_value(MP_DENOM_P(r), 1);
Packit fb9d21
}
Packit fb9d21
Packit fb9d21
/* }}} */
Packit fb9d21
Packit fb9d21
/* {{{ mp_rat_alloc() */
Packit fb9d21
Packit fb9d21
mp_rat mp_rat_alloc(void)
Packit fb9d21
{
Packit fb9d21
  mp_rat out = malloc(sizeof(*out));
Packit fb9d21
Packit fb9d21
  if (out != NULL) {
Packit fb9d21
    if (mp_rat_init(out) != MP_OK) {
Packit fb9d21
      free(out);
Packit fb9d21
      return NULL;
Packit fb9d21
    }
Packit fb9d21
  }
Packit fb9d21
Packit fb9d21
  return out;
Packit fb9d21
}
Packit fb9d21
Packit fb9d21
/* }}} */
Packit fb9d21
Packit fb9d21
/* {{{ mp_rat_reduce(r) */
Packit fb9d21
Packit fb9d21
mp_result mp_rat_reduce(mp_rat r) {
Packit fb9d21
  return s_rat_reduce(r);
Packit fb9d21
}
Packit fb9d21
Packit fb9d21
/* }}} */
Packit fb9d21
Packit fb9d21
/* {{{ mp_rat_init_size(r, n_prec, d_prec) */
Packit fb9d21
Packit fb9d21
mp_result mp_rat_init_size(mp_rat r, mp_size n_prec, mp_size d_prec)
Packit fb9d21
{
Packit fb9d21
  mp_result res;
Packit fb9d21
Packit fb9d21
  if ((res = mp_int_init_size(MP_NUMER_P(r), n_prec)) != MP_OK)
Packit fb9d21
    return res;
Packit fb9d21
  if ((res = mp_int_init_size(MP_DENOM_P(r), d_prec)) != MP_OK) {
Packit fb9d21
    mp_int_clear(MP_NUMER_P(r));
Packit fb9d21
    return res;
Packit fb9d21
  }
Packit fb9d21
  
Packit fb9d21
  return mp_int_set_value(MP_DENOM_P(r), 1);
Packit fb9d21
}
Packit fb9d21
Packit fb9d21
/* }}} */
Packit fb9d21
Packit fb9d21
/* {{{ mp_rat_init_copy(r, old) */
Packit fb9d21
Packit fb9d21
mp_result mp_rat_init_copy(mp_rat r, mp_rat old)
Packit fb9d21
{
Packit fb9d21
  mp_result res;
Packit fb9d21
Packit fb9d21
  if ((res = mp_int_init_copy(MP_NUMER_P(r), MP_NUMER_P(old))) != MP_OK)
Packit fb9d21
    return res;
Packit fb9d21
  if ((res = mp_int_init_copy(MP_DENOM_P(r), MP_DENOM_P(old))) != MP_OK) 
Packit fb9d21
    mp_int_clear(MP_NUMER_P(r));
Packit fb9d21
  
Packit fb9d21
  return res;
Packit fb9d21
}
Packit fb9d21
Packit fb9d21
/* }}} */
Packit fb9d21
Packit fb9d21
/* {{{ mp_rat_set_value(r, numer, denom) */
Packit fb9d21
Packit fb9d21
mp_result mp_rat_set_value(mp_rat r, mp_small numer, mp_small denom)
Packit fb9d21
{
Packit fb9d21
  mp_result res;
Packit fb9d21
Packit fb9d21
  if (denom == 0)
Packit fb9d21
    return MP_UNDEF;
Packit fb9d21
Packit fb9d21
  if ((res = mp_int_set_value(MP_NUMER_P(r), numer)) != MP_OK)
Packit fb9d21
    return res;
Packit fb9d21
  if ((res = mp_int_set_value(MP_DENOM_P(r), denom)) != MP_OK)
Packit fb9d21
    return res;
Packit fb9d21
Packit fb9d21
  return s_rat_reduce(r);
Packit fb9d21
}
Packit fb9d21
Packit fb9d21
/* }}} */
Packit fb9d21
Packit fb9d21
/* {{{ mp_rat_set_uvalue(r, numer, denom) */
Packit fb9d21
Packit fb9d21
mp_result mp_rat_set_uvalue(mp_rat r, mp_usmall numer, mp_usmall denom)
Packit fb9d21
{
Packit fb9d21
  mp_result res;
Packit fb9d21
Packit fb9d21
  if (denom == 0)
Packit fb9d21
    return MP_UNDEF;
Packit fb9d21
Packit fb9d21
  if ((res = mp_int_set_uvalue(MP_NUMER_P(r), numer)) != MP_OK)
Packit fb9d21
    return res;
Packit fb9d21
  if ((res = mp_int_set_uvalue(MP_DENOM_P(r), denom)) != MP_OK)
Packit fb9d21
    return res;
Packit fb9d21
Packit fb9d21
  return s_rat_reduce(r);
Packit fb9d21
}
Packit fb9d21
Packit fb9d21
/* }}} */
Packit fb9d21
Packit fb9d21
/* {{{ mp_rat_clear(r) */
Packit fb9d21
Packit fb9d21
void      mp_rat_clear(mp_rat r)
Packit fb9d21
{
Packit fb9d21
  mp_int_clear(MP_NUMER_P(r));
Packit fb9d21
  mp_int_clear(MP_DENOM_P(r));
Packit fb9d21
Packit fb9d21
}
Packit fb9d21
Packit fb9d21
/* }}} */
Packit fb9d21
Packit fb9d21
/* {{{ mp_rat_free(r) */
Packit fb9d21
Packit fb9d21
void      mp_rat_free(mp_rat r)
Packit fb9d21
{
Packit fb9d21
  NRCHECK(r != NULL);
Packit fb9d21
  
Packit fb9d21
  if (r->num.digits != NULL)
Packit fb9d21
    mp_rat_clear(r);
Packit fb9d21
Packit fb9d21
  free(r);
Packit fb9d21
}
Packit fb9d21
Packit fb9d21
/* }}} */
Packit fb9d21
Packit fb9d21
/* {{{ mp_rat_numer(r, z) */
Packit fb9d21
Packit fb9d21
mp_result mp_rat_numer(mp_rat r, mp_int z)
Packit fb9d21
{
Packit fb9d21
  return mp_int_copy(MP_NUMER_P(r), z);
Packit fb9d21
}
Packit fb9d21
Packit fb9d21
/* }}} */
Packit fb9d21
Packit fb9d21
/* {{{ mp_rat_numer_ref(r) */
Packit fb9d21
Packit fb9d21
mp_int mp_rat_numer_ref(mp_rat r)
Packit fb9d21
{
Packit fb9d21
  return MP_NUMER_P(r);
Packit fb9d21
}
Packit fb9d21
Packit fb9d21
/* }}} */
Packit fb9d21
Packit fb9d21
Packit fb9d21
/* {{{ mp_rat_denom(r, z) */
Packit fb9d21
Packit fb9d21
mp_result mp_rat_denom(mp_rat r, mp_int z)
Packit fb9d21
{
Packit fb9d21
  return mp_int_copy(MP_DENOM_P(r), z);
Packit fb9d21
}
Packit fb9d21
Packit fb9d21
/* }}} */
Packit fb9d21
Packit fb9d21
/* {{{ mp_rat_denom_ref(r) */
Packit fb9d21
Packit fb9d21
mp_int    mp_rat_denom_ref(mp_rat r)
Packit fb9d21
{
Packit fb9d21
  return MP_DENOM_P(r);
Packit fb9d21
}
Packit fb9d21
Packit fb9d21
/* }}} */
Packit fb9d21
Packit fb9d21
/* {{{ mp_rat_sign(r) */
Packit fb9d21
Packit fb9d21
mp_sign   mp_rat_sign(mp_rat r)
Packit fb9d21
{
Packit fb9d21
  return MP_SIGN(MP_NUMER_P(r));
Packit fb9d21
}
Packit fb9d21
Packit fb9d21
/* }}} */
Packit fb9d21
Packit fb9d21
/* {{{ mp_rat_copy(a, c) */
Packit fb9d21
Packit fb9d21
mp_result mp_rat_copy(mp_rat a, mp_rat c)
Packit fb9d21
{
Packit fb9d21
  mp_result res;
Packit fb9d21
Packit fb9d21
  if ((res = mp_int_copy(MP_NUMER_P(a), MP_NUMER_P(c))) != MP_OK)
Packit fb9d21
    return res;
Packit fb9d21
  
Packit fb9d21
  res = mp_int_copy(MP_DENOM_P(a), MP_DENOM_P(c));
Packit fb9d21
  return res;
Packit fb9d21
}
Packit fb9d21
Packit fb9d21
/* }}} */
Packit fb9d21
Packit fb9d21
/* {{{ mp_rat_zero(r) */
Packit fb9d21
Packit fb9d21
void      mp_rat_zero(mp_rat r)
Packit fb9d21
{
Packit fb9d21
  mp_int_zero(MP_NUMER_P(r));
Packit fb9d21
  mp_int_set_value(MP_DENOM_P(r), 1);
Packit fb9d21
  
Packit fb9d21
}
Packit fb9d21
Packit fb9d21
/* }}} */
Packit fb9d21
Packit fb9d21
/* {{{ mp_rat_abs(a, c) */
Packit fb9d21
Packit fb9d21
mp_result mp_rat_abs(mp_rat a, mp_rat c)
Packit fb9d21
{
Packit fb9d21
  mp_result res;
Packit fb9d21
Packit fb9d21
  if ((res = mp_int_abs(MP_NUMER_P(a), MP_NUMER_P(c))) != MP_OK)
Packit fb9d21
    return res;
Packit fb9d21
  
Packit fb9d21
  res = mp_int_abs(MP_DENOM_P(a), MP_DENOM_P(c));
Packit fb9d21
  return res;
Packit fb9d21
}
Packit fb9d21
Packit fb9d21
/* }}} */
Packit fb9d21
Packit fb9d21
/* {{{ mp_rat_neg(a, c) */
Packit fb9d21
Packit fb9d21
mp_result mp_rat_neg(mp_rat a, mp_rat c)
Packit fb9d21
{
Packit fb9d21
  mp_result res;
Packit fb9d21
Packit fb9d21
  if ((res = mp_int_neg(MP_NUMER_P(a), MP_NUMER_P(c))) != MP_OK)
Packit fb9d21
    return res;
Packit fb9d21
Packit fb9d21
  res = mp_int_copy(MP_DENOM_P(a), MP_DENOM_P(c));
Packit fb9d21
  return res;
Packit fb9d21
}
Packit fb9d21
Packit fb9d21
/* }}} */
Packit fb9d21
Packit fb9d21
/* {{{ mp_rat_recip(a, c) */
Packit fb9d21
Packit fb9d21
mp_result mp_rat_recip(mp_rat a, mp_rat c)
Packit fb9d21
{
Packit fb9d21
  mp_result res;
Packit fb9d21
Packit fb9d21
  if (mp_rat_compare_zero(a) == 0)
Packit fb9d21
    return MP_UNDEF;
Packit fb9d21
Packit fb9d21
  if ((res = mp_rat_copy(a, c)) != MP_OK)
Packit fb9d21
    return res;
Packit fb9d21
Packit fb9d21
  mp_int_swap(MP_NUMER_P(c), MP_DENOM_P(c));
Packit fb9d21
Packit fb9d21
  /* Restore the signs of the swapped elements */
Packit fb9d21
  {
Packit fb9d21
    mp_sign tmp = MP_SIGN(MP_NUMER_P(c));
Packit fb9d21
Packit fb9d21
    MP_SIGN(MP_NUMER_P(c)) = MP_SIGN(MP_DENOM_P(c));
Packit fb9d21
    MP_SIGN(MP_DENOM_P(c)) = tmp;
Packit fb9d21
  }
Packit fb9d21
Packit fb9d21
  return MP_OK;
Packit fb9d21
}
Packit fb9d21
Packit fb9d21
/* }}} */
Packit fb9d21
Packit fb9d21
/* {{{ mp_rat_add(a, b, c) */
Packit fb9d21
Packit fb9d21
mp_result mp_rat_add(mp_rat a, mp_rat b, mp_rat c)
Packit fb9d21
{
Packit fb9d21
  return s_rat_combine(a, b, c, mp_int_add);
Packit fb9d21
Packit fb9d21
}
Packit fb9d21
Packit fb9d21
/* }}} */
Packit fb9d21
Packit fb9d21
/* {{{ mp_rat_sub(a, b, c) */
Packit fb9d21
Packit fb9d21
mp_result mp_rat_sub(mp_rat a, mp_rat b, mp_rat c)
Packit fb9d21
{
Packit fb9d21
  return s_rat_combine(a, b, c, mp_int_sub);
Packit fb9d21
Packit fb9d21
}
Packit fb9d21
Packit fb9d21
/* }}} */
Packit fb9d21
Packit fb9d21
/* {{{ mp_rat_mul(a, b, c) */
Packit fb9d21
Packit fb9d21
mp_result mp_rat_mul(mp_rat a, mp_rat b, mp_rat c)
Packit fb9d21
{
Packit fb9d21
  mp_result res;
Packit fb9d21
Packit fb9d21
  if ((res = mp_int_mul(MP_NUMER_P(a), MP_NUMER_P(b), MP_NUMER_P(c))) != MP_OK)
Packit fb9d21
    return res;
Packit fb9d21
Packit fb9d21
  if (mp_int_compare_zero(MP_NUMER_P(c)) != 0) {
Packit fb9d21
    if ((res = mp_int_mul(MP_DENOM_P(a), MP_DENOM_P(b), MP_DENOM_P(c))) != MP_OK)
Packit fb9d21
      return res;
Packit fb9d21
  }
Packit fb9d21
Packit fb9d21
  return s_rat_reduce(c);
Packit fb9d21
}
Packit fb9d21
Packit fb9d21
/* }}} */
Packit fb9d21
Packit fb9d21
/* {{{ mp_int_div(a, b, c) */
Packit fb9d21
Packit fb9d21
mp_result mp_rat_div(mp_rat a, mp_rat b, mp_rat c)
Packit fb9d21
{
Packit fb9d21
  mp_result res = MP_OK;
Packit fb9d21
Packit fb9d21
  if (mp_rat_compare_zero(b) == 0)
Packit fb9d21
    return MP_UNDEF;
Packit fb9d21
Packit fb9d21
  if (c == a || c == b) {
Packit fb9d21
    mpz_t tmp;
Packit fb9d21
Packit fb9d21
    if ((res = mp_int_init(&tmp)) != MP_OK) return res;
Packit fb9d21
    if ((res = mp_int_mul(MP_NUMER_P(a), MP_DENOM_P(b), &tmp)) != MP_OK) 
Packit fb9d21
      goto CLEANUP;
Packit fb9d21
    if ((res = mp_int_mul(MP_DENOM_P(a), MP_NUMER_P(b), MP_DENOM_P(c))) != MP_OK)
Packit fb9d21
      goto CLEANUP;
Packit fb9d21
    res = mp_int_copy(&tmp, MP_NUMER_P(c));
Packit fb9d21
Packit fb9d21
  CLEANUP:
Packit fb9d21
    mp_int_clear(&tmp);
Packit fb9d21
  }
Packit fb9d21
  else {
Packit fb9d21
    if ((res = mp_int_mul(MP_NUMER_P(a), MP_DENOM_P(b), MP_NUMER_P(c))) != MP_OK)
Packit fb9d21
      return res;
Packit fb9d21
    if ((res = mp_int_mul(MP_DENOM_P(a), MP_NUMER_P(b), MP_DENOM_P(c))) != MP_OK)
Packit fb9d21
      return res;
Packit fb9d21
  }
Packit fb9d21
Packit fb9d21
  if (res != MP_OK)
Packit fb9d21
    return res;
Packit fb9d21
  else
Packit fb9d21
    return s_rat_reduce(c);
Packit fb9d21
}
Packit fb9d21
Packit fb9d21
/* }}} */
Packit fb9d21
Packit fb9d21
/* {{{ mp_rat_add_int(a, b, c) */
Packit fb9d21
Packit fb9d21
mp_result mp_rat_add_int(mp_rat a, mp_int b, mp_rat c)
Packit fb9d21
{
Packit fb9d21
  mpz_t tmp;
Packit fb9d21
  mp_result res;
Packit fb9d21
Packit fb9d21
  if ((res = mp_int_init_copy(&tmp, b)) != MP_OK)
Packit fb9d21
    return res;
Packit fb9d21
Packit fb9d21
  if ((res = mp_int_mul(&tmp, MP_DENOM_P(a), &tmp)) != MP_OK)
Packit fb9d21
    goto CLEANUP;
Packit fb9d21
Packit fb9d21
  if ((res = mp_rat_copy(a, c)) != MP_OK)
Packit fb9d21
    goto CLEANUP;
Packit fb9d21
Packit fb9d21
  if ((res = mp_int_add(MP_NUMER_P(c), &tmp, MP_NUMER_P(c))) != MP_OK)
Packit fb9d21
    goto CLEANUP;
Packit fb9d21
Packit fb9d21
  res = s_rat_reduce(c);
Packit fb9d21
Packit fb9d21
 CLEANUP:
Packit fb9d21
  mp_int_clear(&tmp);
Packit fb9d21
  return res;
Packit fb9d21
}
Packit fb9d21
Packit fb9d21
/* }}} */
Packit fb9d21
Packit fb9d21
/* {{{ mp_rat_sub_int(a, b, c) */
Packit fb9d21
Packit fb9d21
mp_result mp_rat_sub_int(mp_rat a, mp_int b, mp_rat c)
Packit fb9d21
{
Packit fb9d21
  mpz_t tmp;
Packit fb9d21
  mp_result res;
Packit fb9d21
Packit fb9d21
  if ((res = mp_int_init_copy(&tmp, b)) != MP_OK)
Packit fb9d21
    return res;
Packit fb9d21
Packit fb9d21
  if ((res = mp_int_mul(&tmp, MP_DENOM_P(a), &tmp)) != MP_OK)
Packit fb9d21
    goto CLEANUP;
Packit fb9d21
Packit fb9d21
  if ((res = mp_rat_copy(a, c)) != MP_OK)
Packit fb9d21
    goto CLEANUP;
Packit fb9d21
Packit fb9d21
  if ((res = mp_int_sub(MP_NUMER_P(c), &tmp, MP_NUMER_P(c))) != MP_OK)
Packit fb9d21
    goto CLEANUP;
Packit fb9d21
Packit fb9d21
  res = s_rat_reduce(c);
Packit fb9d21
Packit fb9d21
 CLEANUP:
Packit fb9d21
  mp_int_clear(&tmp);
Packit fb9d21
  return res;
Packit fb9d21
}
Packit fb9d21
Packit fb9d21
/* }}} */
Packit fb9d21
Packit fb9d21
/* {{{ mp_rat_mul_int(a, b, c) */
Packit fb9d21
Packit fb9d21
mp_result mp_rat_mul_int(mp_rat a, mp_int b, mp_rat c)
Packit fb9d21
{
Packit fb9d21
  mp_result res;
Packit fb9d21
Packit fb9d21
  if ((res = mp_rat_copy(a, c)) != MP_OK)
Packit fb9d21
    return res;
Packit fb9d21
Packit fb9d21
  if ((res = mp_int_mul(MP_NUMER_P(c), b, MP_NUMER_P(c))) != MP_OK)
Packit fb9d21
    return res;
Packit fb9d21
Packit fb9d21
  return s_rat_reduce(c);
Packit fb9d21
}
Packit fb9d21
Packit fb9d21
/* }}} */
Packit fb9d21
Packit fb9d21
/* {{{ mp_rat_div_int(a, b, c) */
Packit fb9d21
Packit fb9d21
mp_result mp_rat_div_int(mp_rat a, mp_int b, mp_rat c)
Packit fb9d21
{
Packit fb9d21
  mp_result res;
Packit fb9d21
Packit fb9d21
  if (mp_int_compare_zero(b) == 0)
Packit fb9d21
    return MP_UNDEF;
Packit fb9d21
Packit fb9d21
  if ((res = mp_rat_copy(a, c)) != MP_OK)
Packit fb9d21
    return res;
Packit fb9d21
Packit fb9d21
  if ((res = mp_int_mul(MP_DENOM_P(c), b, MP_DENOM_P(c))) != MP_OK)
Packit fb9d21
    return res;
Packit fb9d21
Packit fb9d21
  return s_rat_reduce(c);
Packit fb9d21
}
Packit fb9d21
Packit fb9d21
/* }}} */
Packit fb9d21
Packit fb9d21
/* {{{ mp_rat_expt(a, b, c) */
Packit fb9d21
Packit fb9d21
mp_result mp_rat_expt(mp_rat a, mp_small b, mp_rat c)
Packit fb9d21
{
Packit fb9d21
  mp_result  res;
Packit fb9d21
Packit fb9d21
  /* Special cases for easy powers. */
Packit fb9d21
  if (b == 0)
Packit fb9d21
    return mp_rat_set_value(c, 1, 1);
Packit fb9d21
  else if(b == 1)
Packit fb9d21
    return mp_rat_copy(a, c);
Packit fb9d21
Packit fb9d21
  /* Since rationals are always stored in lowest terms, it is not necessary to
Packit fb9d21
     reduce again when raising to an integer power. */
Packit fb9d21
  if ((res = mp_int_expt(MP_NUMER_P(a), b, MP_NUMER_P(c))) != MP_OK)
Packit fb9d21
    return res;
Packit fb9d21
Packit fb9d21
  return mp_int_expt(MP_DENOM_P(a), b, MP_DENOM_P(c));
Packit fb9d21
}
Packit fb9d21
Packit fb9d21
/* }}} */
Packit fb9d21
Packit fb9d21
/* {{{ mp_rat_compare(a, b) */
Packit fb9d21
Packit fb9d21
int       mp_rat_compare(mp_rat a, mp_rat b)
Packit fb9d21
{
Packit fb9d21
  /* Quick check for opposite signs.  Works because the sign of the numerator
Packit fb9d21
     is always definitive. */
Packit fb9d21
  if (MP_SIGN(MP_NUMER_P(a)) != MP_SIGN(MP_NUMER_P(b))) {
Packit fb9d21
    if (MP_SIGN(MP_NUMER_P(a)) == MP_ZPOS)
Packit fb9d21
      return 1;
Packit fb9d21
    else
Packit fb9d21
      return -1;
Packit fb9d21
  }
Packit fb9d21
  else {
Packit fb9d21
    /* Compare absolute magnitudes; if both are positive, the answer stands,
Packit fb9d21
       otherwise it needs to be reflected about zero. */
Packit fb9d21
    int cmp = mp_rat_compare_unsigned(a, b);
Packit fb9d21
Packit fb9d21
    if (MP_SIGN(MP_NUMER_P(a)) == MP_ZPOS)
Packit fb9d21
      return cmp;
Packit fb9d21
    else
Packit fb9d21
      return -cmp;
Packit fb9d21
  }
Packit fb9d21
}
Packit fb9d21
Packit fb9d21
/* }}} */
Packit fb9d21
Packit fb9d21
/* {{{ mp_rat_compare_unsigned(a, b) */
Packit fb9d21
Packit fb9d21
int       mp_rat_compare_unsigned(mp_rat a, mp_rat b)
Packit fb9d21
{
Packit fb9d21
  /* If the denominators are equal, we can quickly compare numerators without
Packit fb9d21
     multiplying.  Otherwise, we actually have to do some work. */
Packit fb9d21
  if (mp_int_compare_unsigned(MP_DENOM_P(a), MP_DENOM_P(b)) == 0)
Packit fb9d21
    return mp_int_compare_unsigned(MP_NUMER_P(a), MP_NUMER_P(b));
Packit fb9d21
Packit fb9d21
  else {
Packit fb9d21
    mpz_t  temp[2];
Packit fb9d21
    mp_result res;
Packit fb9d21
    int  cmp = INT_MAX, last = 0;
Packit fb9d21
Packit fb9d21
    /* t0 = num(a) * den(b), t1 = num(b) * den(a) */
Packit fb9d21
    SETUP(mp_int_init_copy(TEMP(last), MP_NUMER_P(a)), last);
Packit fb9d21
    SETUP(mp_int_init_copy(TEMP(last), MP_NUMER_P(b)), last);
Packit fb9d21
Packit fb9d21
    if ((res = mp_int_mul(TEMP(0), MP_DENOM_P(b), TEMP(0))) != MP_OK ||
Packit fb9d21
	(res = mp_int_mul(TEMP(1), MP_DENOM_P(a), TEMP(1))) != MP_OK)
Packit fb9d21
      goto CLEANUP;
Packit fb9d21
    
Packit fb9d21
    cmp = mp_int_compare_unsigned(TEMP(0), TEMP(1));
Packit fb9d21
    
Packit fb9d21
  CLEANUP:
Packit fb9d21
    while (--last >= 0)
Packit fb9d21
      mp_int_clear(TEMP(last));
Packit fb9d21
Packit fb9d21
    return cmp;
Packit fb9d21
  }
Packit fb9d21
}
Packit fb9d21
Packit fb9d21
/* }}} */
Packit fb9d21
Packit fb9d21
/* {{{ mp_rat_compare_zero(r) */
Packit fb9d21
Packit fb9d21
int       mp_rat_compare_zero(mp_rat r)
Packit fb9d21
{
Packit fb9d21
  return mp_int_compare_zero(MP_NUMER_P(r));
Packit fb9d21
}
Packit fb9d21
Packit fb9d21
/* }}} */
Packit fb9d21
Packit fb9d21
/* {{{ mp_rat_compare_value(r, n, d) */
Packit fb9d21
Packit fb9d21
int       mp_rat_compare_value(mp_rat r, mp_small n, mp_small d)
Packit fb9d21
{
Packit fb9d21
  mpq_t tmp;
Packit fb9d21
  mp_result res;
Packit fb9d21
  int  out = INT_MAX;
Packit fb9d21
Packit fb9d21
  if ((res = mp_rat_init(&tmp)) != MP_OK)
Packit fb9d21
    return out;
Packit fb9d21
  if ((res = mp_rat_set_value(&tmp, n, d)) != MP_OK)
Packit fb9d21
    goto CLEANUP;
Packit fb9d21
  
Packit fb9d21
  out = mp_rat_compare(r, &tmp);
Packit fb9d21
  
Packit fb9d21
 CLEANUP:
Packit fb9d21
  mp_rat_clear(&tmp);
Packit fb9d21
  return out;
Packit fb9d21
}
Packit fb9d21
Packit fb9d21
/* }}} */
Packit fb9d21
Packit fb9d21
/* {{{ mp_rat_is_integer(r) */
Packit fb9d21
Packit fb9d21
int       mp_rat_is_integer(mp_rat r)
Packit fb9d21
{
Packit fb9d21
  return (mp_int_compare_value(MP_DENOM_P(r), 1) == 0);
Packit fb9d21
}
Packit fb9d21
Packit fb9d21
/* }}} */
Packit fb9d21
Packit fb9d21
/* {{{ mp_rat_to_ints(r, *num, *den) */
Packit fb9d21
Packit fb9d21
mp_result mp_rat_to_ints(mp_rat r, mp_small *num, mp_small *den)
Packit fb9d21
{
Packit fb9d21
  mp_result res;
Packit fb9d21
Packit fb9d21
  if ((res = mp_int_to_int(MP_NUMER_P(r), num)) != MP_OK)
Packit fb9d21
    return res;
Packit fb9d21
Packit fb9d21
  res = mp_int_to_int(MP_DENOM_P(r), den);
Packit fb9d21
  return res;
Packit fb9d21
}
Packit fb9d21
Packit fb9d21
/* }}} */
Packit fb9d21
Packit fb9d21
/* {{{ mp_rat_to_string(r, radix, *str, limit) */
Packit fb9d21
Packit fb9d21
mp_result mp_rat_to_string(mp_rat r, mp_size radix, char *str, int limit)
Packit fb9d21
{
Packit fb9d21
  char *start;
Packit fb9d21
  int   len;
Packit fb9d21
  mp_result res;
Packit fb9d21
Packit fb9d21
  /* Write the numerator.  The sign of the rational number is written by the
Packit fb9d21
     underlying integer implementation. */
Packit fb9d21
  if ((res = mp_int_to_string(MP_NUMER_P(r), radix, str, limit)) != MP_OK)
Packit fb9d21
    return res;
Packit fb9d21
Packit fb9d21
  /* If the value is zero, don't bother writing any denominator */
Packit fb9d21
  if (mp_int_compare_zero(MP_NUMER_P(r)) == 0)
Packit fb9d21
    return MP_OK;
Packit fb9d21
  
Packit fb9d21
  /* Locate the end of the numerator, and make sure we are not going to exceed
Packit fb9d21
     the limit by writing a slash. */
Packit fb9d21
  len = strlen(str);
Packit fb9d21
  start = str + len;
Packit fb9d21
  limit -= len;
Packit fb9d21
  if(limit == 0)
Packit fb9d21
    return MP_TRUNC;
Packit fb9d21
Packit fb9d21
  *start++ = '/';
Packit fb9d21
  limit -= 1;
Packit fb9d21
  
Packit fb9d21
  res = mp_int_to_string(MP_DENOM_P(r), radix, start, limit);
Packit fb9d21
  return res;
Packit fb9d21
}
Packit fb9d21
Packit fb9d21
/* }}} */
Packit fb9d21
Packit fb9d21
/* {{{ mp_rat_to_decimal(r, radix, prec, *str, limit) */
Packit fb9d21
mp_result mp_rat_to_decimal(mp_rat r, mp_size radix, mp_size prec,
Packit fb9d21
                            mp_round_mode round, char *str, int limit)
Packit fb9d21
{
Packit fb9d21
  mpz_t temp[3];
Packit fb9d21
  mp_result res;
Packit fb9d21
  char *start = str;
Packit fb9d21
  int len, lead_0, left = limit, last = 0;
Packit fb9d21
    
Packit fb9d21
  SETUP(mp_int_init_copy(TEMP(last), MP_NUMER_P(r)), last);
Packit fb9d21
  SETUP(mp_int_init(TEMP(last)), last);
Packit fb9d21
  SETUP(mp_int_init(TEMP(last)), last);
Packit fb9d21
Packit fb9d21
  /* Get the unsigned integer part by dividing denominator into the absolute
Packit fb9d21
     value of the numerator. */
Packit fb9d21
  mp_int_abs(TEMP(0), TEMP(0));
Packit fb9d21
  if ((res = mp_int_div(TEMP(0), MP_DENOM_P(r), TEMP(0), TEMP(1))) != MP_OK)
Packit fb9d21
    goto CLEANUP;
Packit fb9d21
Packit fb9d21
  /* Now:  T0 = integer portion, unsigned;
Packit fb9d21
           T1 = remainder, from which fractional part is computed. */
Packit fb9d21
Packit fb9d21
  /* Count up leading zeroes after the radix point. */
Packit fb9d21
  for (lead_0 = 0; lead_0 < prec && mp_int_compare(TEMP(1), MP_DENOM_P(r)) < 0; 
Packit fb9d21
      ++lead_0) {
Packit fb9d21
    if ((res = mp_int_mul_value(TEMP(1), radix, TEMP(1))) != MP_OK)
Packit fb9d21
      goto CLEANUP;
Packit fb9d21
  }
Packit fb9d21
Packit fb9d21
  /* Multiply remainder by a power of the radix sufficient to get the right
Packit fb9d21
     number of significant figures. */
Packit fb9d21
  if (prec > lead_0) {
Packit fb9d21
    if ((res = mp_int_expt_value(radix, prec - lead_0, TEMP(2))) != MP_OK)
Packit fb9d21
      goto CLEANUP;
Packit fb9d21
    if ((res = mp_int_mul(TEMP(1), TEMP(2), TEMP(1))) != MP_OK)
Packit fb9d21
      goto CLEANUP;
Packit fb9d21
  }
Packit fb9d21
  if ((res = mp_int_div(TEMP(1), MP_DENOM_P(r), TEMP(1), TEMP(2))) != MP_OK)
Packit fb9d21
    goto CLEANUP;
Packit fb9d21
Packit fb9d21
  /* Now:  T1 = significant digits of fractional part;
Packit fb9d21
           T2 = leftovers, to use for rounding. 
Packit fb9d21
Packit fb9d21
     At this point, what we do depends on the rounding mode.  The default is
Packit fb9d21
     MP_ROUND_DOWN, for which everything is as it should be already.
Packit fb9d21
  */
Packit fb9d21
  switch (round) {
Packit fb9d21
    int cmp;
Packit fb9d21
Packit fb9d21
  case MP_ROUND_UP:
Packit fb9d21
    if (mp_int_compare_zero(TEMP(2)) != 0) {
Packit fb9d21
      if (prec == 0)
Packit fb9d21
	res = mp_int_add_value(TEMP(0), 1, TEMP(0));
Packit fb9d21
      else
Packit fb9d21
	res = mp_int_add_value(TEMP(1), 1, TEMP(1));
Packit fb9d21
    }
Packit fb9d21
    break;
Packit fb9d21
Packit fb9d21
  case MP_ROUND_HALF_UP:
Packit fb9d21
  case MP_ROUND_HALF_DOWN:
Packit fb9d21
    if ((res = mp_int_mul_pow2(TEMP(2), 1, TEMP(2))) != MP_OK)
Packit fb9d21
      goto CLEANUP;
Packit fb9d21
Packit fb9d21
    cmp = mp_int_compare(TEMP(2), MP_DENOM_P(r));    
Packit fb9d21
Packit fb9d21
    if (round == MP_ROUND_HALF_UP)
Packit fb9d21
      cmp += 1;
Packit fb9d21
Packit fb9d21
    if (cmp > 0) {
Packit fb9d21
      if (prec == 0)
Packit fb9d21
	res = mp_int_add_value(TEMP(0), 1, TEMP(0));
Packit fb9d21
      else
Packit fb9d21
	res = mp_int_add_value(TEMP(1), 1, TEMP(1));
Packit fb9d21
    }
Packit fb9d21
    break;
Packit fb9d21
    
Packit fb9d21
  case MP_ROUND_DOWN:
Packit fb9d21
    break;  /* No action required */
Packit fb9d21
Packit fb9d21
  default: 
Packit fb9d21
    return MP_BADARG; /* Invalid rounding specifier */
Packit fb9d21
  }
Packit fb9d21
Packit fb9d21
  /* The sign of the output should be the sign of the numerator, but if all the
Packit fb9d21
     displayed digits will be zero due to the precision, a negative shouldn't
Packit fb9d21
     be shown. */
Packit fb9d21
  if (MP_SIGN(MP_NUMER_P(r)) == MP_NEG &&
Packit fb9d21
      (mp_int_compare_zero(TEMP(0)) != 0 ||
Packit fb9d21
       mp_int_compare_zero(TEMP(1)) != 0)) {
Packit fb9d21
    *start++ = '-';
Packit fb9d21
    left -= 1;
Packit fb9d21
  }
Packit fb9d21
Packit fb9d21
  if ((res = mp_int_to_string(TEMP(0), radix, start, left)) != MP_OK)
Packit fb9d21
    goto CLEANUP;
Packit fb9d21
  
Packit fb9d21
  len = strlen(start);
Packit fb9d21
  start += len;
Packit fb9d21
  left -= len;
Packit fb9d21
  
Packit fb9d21
  if (prec == 0) 
Packit fb9d21
    goto CLEANUP;
Packit fb9d21
  
Packit fb9d21
  *start++ = '.';
Packit fb9d21
  left -= 1;
Packit fb9d21
  
Packit fb9d21
  if (left < prec + 1) {
Packit fb9d21
    res = MP_TRUNC;
Packit fb9d21
    goto CLEANUP;
Packit fb9d21
  }
Packit fb9d21
Packit fb9d21
  memset(start, '0', lead_0 - 1);
Packit fb9d21
  left -= lead_0;
Packit fb9d21
  start += lead_0 - 1;
Packit fb9d21
Packit fb9d21
  res = mp_int_to_string(TEMP(1), radix, start, left);
Packit fb9d21
Packit fb9d21
 CLEANUP:
Packit fb9d21
  while (--last >= 0)
Packit fb9d21
    mp_int_clear(TEMP(last));
Packit fb9d21
  
Packit fb9d21
  return res;
Packit fb9d21
}
Packit fb9d21
Packit fb9d21
/* }}} */
Packit fb9d21
Packit fb9d21
/* {{{ mp_rat_string_len(r, radix) */
Packit fb9d21
Packit fb9d21
mp_result mp_rat_string_len(mp_rat r, mp_size radix)
Packit fb9d21
{
Packit fb9d21
  mp_result n_len, d_len = 0;
Packit fb9d21
Packit fb9d21
  n_len = mp_int_string_len(MP_NUMER_P(r), radix);
Packit fb9d21
Packit fb9d21
  if (mp_int_compare_zero(MP_NUMER_P(r)) != 0)
Packit fb9d21
    d_len = mp_int_string_len(MP_DENOM_P(r), radix);
Packit fb9d21
Packit fb9d21
  /* Though simplistic, this formula is correct.  Space for the sign flag is
Packit fb9d21
     included in n_len, and the space for the NUL that is counted in n_len
Packit fb9d21
     counts for the separator here.  The space for the NUL counted in d_len
Packit fb9d21
     counts for the final terminator here. */
Packit fb9d21
Packit fb9d21
  return n_len + d_len;
Packit fb9d21
Packit fb9d21
}
Packit fb9d21
Packit fb9d21
/* }}} */
Packit fb9d21
Packit fb9d21
/* {{{ mp_rat_decimal_len(r, radix, prec) */
Packit fb9d21
Packit fb9d21
mp_result mp_rat_decimal_len(mp_rat r, mp_size radix, mp_size prec)
Packit fb9d21
{
Packit fb9d21
  int  z_len, f_len;
Packit fb9d21
Packit fb9d21
  z_len = mp_int_string_len(MP_NUMER_P(r), radix);
Packit fb9d21
  
Packit fb9d21
  if (prec == 0)
Packit fb9d21
    f_len = 1; /* terminator only */
Packit fb9d21
  else
Packit fb9d21
    f_len = 1 + prec + 1; /* decimal point, digits, terminator */
Packit fb9d21
  
Packit fb9d21
  return z_len + f_len;
Packit fb9d21
}
Packit fb9d21
Packit fb9d21
/* }}} */
Packit fb9d21
Packit fb9d21
/* {{{ mp_rat_read_string(r, radix, *str) */
Packit fb9d21
Packit fb9d21
mp_result mp_rat_read_string(mp_rat r, mp_size radix, const char *str)
Packit fb9d21
{
Packit fb9d21
  return mp_rat_read_cstring(r, radix, str, NULL);
Packit fb9d21
}
Packit fb9d21
Packit fb9d21
/* }}} */
Packit fb9d21
Packit fb9d21
/* {{{ mp_rat_read_cstring(r, radix, *str, **end) */
Packit fb9d21
Packit fb9d21
mp_result mp_rat_read_cstring(mp_rat r, mp_size radix, const char *str, 
Packit fb9d21
			      char **end)
Packit fb9d21
{
Packit fb9d21
  mp_result res;
Packit fb9d21
  char *endp;
Packit fb9d21
Packit fb9d21
  if ((res = mp_int_read_cstring(MP_NUMER_P(r), radix, str, &endp)) != MP_OK &&
Packit fb9d21
      (res != MP_TRUNC))
Packit fb9d21
    return res;
Packit fb9d21
Packit fb9d21
  /* Skip whitespace between numerator and (possible) separator */
Packit fb9d21
  while (isspace((unsigned char) *endp))
Packit fb9d21
    ++endp;
Packit fb9d21
  
Packit fb9d21
  /* If there is no separator, we will stop reading at this point. */
Packit fb9d21
  if (*endp != '/') {
Packit fb9d21
    mp_int_set_value(MP_DENOM_P(r), 1);
Packit fb9d21
    if (end != NULL)
Packit fb9d21
      *end = endp;
Packit fb9d21
    return res;
Packit fb9d21
  }
Packit fb9d21
  
Packit fb9d21
  ++endp; /* skip separator */
Packit fb9d21
  if ((res = mp_int_read_cstring(MP_DENOM_P(r), radix, endp, end)) != MP_OK)
Packit fb9d21
    return res;
Packit fb9d21
  
Packit fb9d21
  /* Make sure the value is well-defined */
Packit fb9d21
  if (mp_int_compare_zero(MP_DENOM_P(r)) == 0)
Packit fb9d21
    return MP_UNDEF;
Packit fb9d21
Packit fb9d21
  /* Reduce to lowest terms */
Packit fb9d21
  return s_rat_reduce(r);
Packit fb9d21
}
Packit fb9d21
Packit fb9d21
/* }}} */
Packit fb9d21
Packit fb9d21
/* {{{ mp_rat_read_ustring(r, radix, *str, **end) */
Packit fb9d21
Packit fb9d21
/* Read a string and figure out what format it's in.  The radix may be supplied
Packit fb9d21
   as zero to use "default" behaviour.
Packit fb9d21
Packit fb9d21
   This function will accept either a/b notation or decimal notation.
Packit fb9d21
 */
Packit fb9d21
mp_result mp_rat_read_ustring(mp_rat r, mp_size radix, const char *str, 
Packit fb9d21
			      char **end)
Packit fb9d21
{
Packit fb9d21
  char      *endp;
Packit fb9d21
  mp_result  res;
Packit fb9d21
Packit fb9d21
  if (radix == 0)
Packit fb9d21
    radix = 10;  /* default to decimal input */
Packit fb9d21
Packit fb9d21
  if ((res = mp_rat_read_cstring(r, radix, str, &endp)) != MP_OK) {
Packit fb9d21
    if (res == MP_TRUNC) {
Packit fb9d21
      if (*endp == '.')
Packit fb9d21
	res = mp_rat_read_cdecimal(r, radix, str, &endp);
Packit fb9d21
    }
Packit fb9d21
    else
Packit fb9d21
      return res;
Packit fb9d21
  }
Packit fb9d21
Packit fb9d21
  if (end != NULL)
Packit fb9d21
    *end = endp;
Packit fb9d21
Packit fb9d21
  return res;
Packit fb9d21
}
Packit fb9d21
Packit fb9d21
/* }}} */
Packit fb9d21
Packit fb9d21
/* {{{ mp_rat_read_decimal(r, radix, *str) */
Packit fb9d21
Packit fb9d21
mp_result mp_rat_read_decimal(mp_rat r, mp_size radix, const char *str)
Packit fb9d21
{
Packit fb9d21
  return mp_rat_read_cdecimal(r, radix, str, NULL);
Packit fb9d21
}
Packit fb9d21
Packit fb9d21
/* }}} */
Packit fb9d21
Packit fb9d21
/* {{{ mp_rat_read_cdecimal(r, radix, *str, **end) */
Packit fb9d21
Packit fb9d21
mp_result mp_rat_read_cdecimal(mp_rat r, mp_size radix, const char *str, 
Packit fb9d21
			       char **end)
Packit fb9d21
{
Packit fb9d21
  mp_result res;
Packit fb9d21
  mp_sign   osign;
Packit fb9d21
  char *endp;
Packit fb9d21
Packit fb9d21
  while (isspace((unsigned char) *str))
Packit fb9d21
    ++str;
Packit fb9d21
  
Packit fb9d21
  switch (*str) {
Packit fb9d21
  case '-':
Packit fb9d21
    osign = MP_NEG;
Packit fb9d21
    break;
Packit fb9d21
  default:
Packit fb9d21
    osign = MP_ZPOS;
Packit fb9d21
  }
Packit fb9d21
  
Packit fb9d21
  if ((res = mp_int_read_cstring(MP_NUMER_P(r), radix, str, &endp)) != MP_OK &&
Packit fb9d21
     (res != MP_TRUNC))
Packit fb9d21
    return res;
Packit fb9d21
Packit fb9d21
  /* This needs to be here. */
Packit fb9d21
  (void) mp_int_set_value(MP_DENOM_P(r), 1);
Packit fb9d21
Packit fb9d21
  if (*endp != '.') {
Packit fb9d21
    if (end != NULL)
Packit fb9d21
      *end = endp;
Packit fb9d21
    return res;
Packit fb9d21
  }
Packit fb9d21
Packit fb9d21
  /* If the character following the decimal point is whitespace or a sign flag,
Packit fb9d21
     we will consider this a truncated value.  This special case is because
Packit fb9d21
     mp_int_read_string() will consider whitespace or sign flags to be valid
Packit fb9d21
     starting characters for a value, and we do not want them following the
Packit fb9d21
     decimal point.
Packit fb9d21
Packit fb9d21
     Once we have done this check, it is safe to read in the value of the
Packit fb9d21
     fractional piece as a regular old integer.
Packit fb9d21
  */
Packit fb9d21
  ++endp;
Packit fb9d21
  if (*endp == '\0') {
Packit fb9d21
    if (end != NULL)
Packit fb9d21
      *end = endp;
Packit fb9d21
    return MP_OK;
Packit fb9d21
  }
Packit fb9d21
  else if(isspace((unsigned char) *endp) || *endp == '-' || *endp == '+') {
Packit fb9d21
    return MP_TRUNC;
Packit fb9d21
  }
Packit fb9d21
  else {
Packit fb9d21
    mpz_t  frac;
Packit fb9d21
    mp_result save_res;
Packit fb9d21
    char  *save = endp;
Packit fb9d21
    int    num_lz = 0;
Packit fb9d21
Packit fb9d21
    /* Make a temporary to hold the part after the decimal point. */
Packit fb9d21
    if ((res = mp_int_init(&frac)) != MP_OK)
Packit fb9d21
      return res;
Packit fb9d21
    
Packit fb9d21
    if ((res = mp_int_read_cstring(&frac, radix, endp, &endp)) != MP_OK &&
Packit fb9d21
       (res != MP_TRUNC))
Packit fb9d21
      goto CLEANUP;
Packit fb9d21
Packit fb9d21
    /* Save this response for later. */
Packit fb9d21
    save_res = res;
Packit fb9d21
Packit fb9d21
    if (mp_int_compare_zero(&frac) == 0)
Packit fb9d21
      goto FINISHED;
Packit fb9d21
Packit fb9d21
    /* Discard trailing zeroes (somewhat inefficiently) */
Packit fb9d21
    while (mp_int_divisible_value(&frac, radix))
Packit fb9d21
      if ((res = mp_int_div_value(&frac, radix, &frac, NULL)) != MP_OK)
Packit fb9d21
	goto CLEANUP;
Packit fb9d21
    
Packit fb9d21
    /* Count leading zeros after the decimal point */
Packit fb9d21
    while (save[num_lz] == '0')
Packit fb9d21
      ++num_lz;
Packit fb9d21
Packit fb9d21
    /* Find the least power of the radix that is at least as large as the
Packit fb9d21
       significant value of the fractional part, ignoring leading zeroes.  */
Packit fb9d21
    (void) mp_int_set_value(MP_DENOM_P(r), radix); 
Packit fb9d21
    
Packit fb9d21
    while (mp_int_compare(MP_DENOM_P(r), &frac) < 0) {
Packit fb9d21
      if ((res = mp_int_mul_value(MP_DENOM_P(r), radix, MP_DENOM_P(r))) != MP_OK)
Packit fb9d21
	goto CLEANUP;
Packit fb9d21
    }
Packit fb9d21
    
Packit fb9d21
    /* Also shift by enough to account for leading zeroes */
Packit fb9d21
    while (num_lz > 0) {
Packit fb9d21
      if ((res = mp_int_mul_value(MP_DENOM_P(r), radix, MP_DENOM_P(r))) != MP_OK)
Packit fb9d21
	goto CLEANUP;
Packit fb9d21
Packit fb9d21
      --num_lz;
Packit fb9d21
    }
Packit fb9d21
Packit fb9d21
    /* Having found this power, shift the numerator leftward that many, digits,
Packit fb9d21
       and add the nonzero significant digits of the fractional part to get the
Packit fb9d21
       result. */
Packit fb9d21
    if ((res = mp_int_mul(MP_NUMER_P(r), MP_DENOM_P(r), MP_NUMER_P(r))) != MP_OK)
Packit fb9d21
      goto CLEANUP;
Packit fb9d21
    
Packit fb9d21
    { /* This addition needs to be unsigned. */
Packit fb9d21
      MP_SIGN(MP_NUMER_P(r)) = MP_ZPOS;
Packit fb9d21
      if ((res = mp_int_add(MP_NUMER_P(r), &frac, MP_NUMER_P(r))) != MP_OK)
Packit fb9d21
	goto CLEANUP;
Packit fb9d21
Packit fb9d21
      MP_SIGN(MP_NUMER_P(r)) = osign;
Packit fb9d21
    }
Packit fb9d21
    if ((res = s_rat_reduce(r)) != MP_OK)
Packit fb9d21
      goto CLEANUP;
Packit fb9d21
Packit fb9d21
    /* At this point, what we return depends on whether reading the fractional
Packit fb9d21
       part was truncated or not.  That information is saved from when we
Packit fb9d21
       called mp_int_read_string() above. */
Packit fb9d21
  FINISHED:
Packit fb9d21
    res = save_res;
Packit fb9d21
    if (end != NULL)
Packit fb9d21
      *end = endp;
Packit fb9d21
Packit fb9d21
  CLEANUP:
Packit fb9d21
    mp_int_clear(&frac);
Packit fb9d21
Packit fb9d21
    return res;
Packit fb9d21
  }
Packit fb9d21
}
Packit fb9d21
Packit fb9d21
/* }}} */
Packit fb9d21
Packit fb9d21
/* Private functions for internal use.  Make unchecked assumptions about format
Packit fb9d21
   and validity of inputs. */
Packit fb9d21
Packit fb9d21
/* {{{ s_rat_reduce(r) */
Packit fb9d21
Packit fb9d21
static mp_result s_rat_reduce(mp_rat r)
Packit fb9d21
{
Packit fb9d21
  mpz_t gcd;
Packit fb9d21
  mp_result res = MP_OK;
Packit fb9d21
Packit fb9d21
  if (mp_int_compare_zero(MP_NUMER_P(r)) == 0) {
Packit fb9d21
    mp_int_set_value(MP_DENOM_P(r), 1);
Packit fb9d21
    return MP_OK;
Packit fb9d21
  }
Packit fb9d21
Packit fb9d21
  /* If the greatest common divisor of the numerator and denominator is greater
Packit fb9d21
     than 1, divide it out. */
Packit fb9d21
  if ((res = mp_int_init(&gcd)) != MP_OK)
Packit fb9d21
    return res;
Packit fb9d21
Packit fb9d21
  if ((res = mp_int_gcd(MP_NUMER_P(r), MP_DENOM_P(r), &gcd)) != MP_OK)
Packit fb9d21
    goto CLEANUP;
Packit fb9d21
Packit fb9d21
  if (mp_int_compare_value(&gcd, 1) != 0) {
Packit fb9d21
    if ((res = mp_int_div(MP_NUMER_P(r), &gcd, MP_NUMER_P(r), NULL)) != MP_OK)
Packit fb9d21
      goto CLEANUP;
Packit fb9d21
    if ((res = mp_int_div(MP_DENOM_P(r), &gcd, MP_DENOM_P(r), NULL)) != MP_OK)
Packit fb9d21
      goto CLEANUP;
Packit fb9d21
  }
Packit fb9d21
Packit fb9d21
  /* Fix up the signs of numerator and denominator */
Packit fb9d21
  if (MP_SIGN(MP_NUMER_P(r)) == MP_SIGN(MP_DENOM_P(r)))
Packit fb9d21
    MP_SIGN(MP_NUMER_P(r)) = MP_SIGN(MP_DENOM_P(r)) = MP_ZPOS;
Packit fb9d21
  else {
Packit fb9d21
    MP_SIGN(MP_NUMER_P(r)) = MP_NEG;
Packit fb9d21
    MP_SIGN(MP_DENOM_P(r)) = MP_ZPOS;
Packit fb9d21
  }
Packit fb9d21
Packit fb9d21
 CLEANUP:
Packit fb9d21
  mp_int_clear(&gcd;;
Packit fb9d21
Packit fb9d21
  return res;
Packit fb9d21
}
Packit fb9d21
Packit fb9d21
/* }}} */
Packit fb9d21
Packit fb9d21
/* {{{ s_rat_combine(a, b, c, comb_f) */
Packit fb9d21
Packit fb9d21
static mp_result s_rat_combine(mp_rat a, mp_rat b, mp_rat c, 
Packit fb9d21
			       mp_result (*comb_f)(mp_int, mp_int, mp_int))
Packit fb9d21
{
Packit fb9d21
  mp_result res;
Packit fb9d21
Packit fb9d21
  /* Shortcut when denominators are already common */
Packit fb9d21
  if (mp_int_compare(MP_DENOM_P(a), MP_DENOM_P(b)) == 0) {
Packit fb9d21
    if ((res = (comb_f)(MP_NUMER_P(a), MP_NUMER_P(b), MP_NUMER_P(c))) != MP_OK)
Packit fb9d21
      return res;
Packit fb9d21
    if ((res = mp_int_copy(MP_DENOM_P(a), MP_DENOM_P(c))) != MP_OK)
Packit fb9d21
      return res;
Packit fb9d21
    
Packit fb9d21
    return s_rat_reduce(c);
Packit fb9d21
  }
Packit fb9d21
  else {
Packit fb9d21
    mpz_t  temp[2];
Packit fb9d21
    int    last = 0;
Packit fb9d21
Packit fb9d21
    SETUP(mp_int_init_copy(TEMP(last), MP_NUMER_P(a)), last);
Packit fb9d21
    SETUP(mp_int_init_copy(TEMP(last), MP_NUMER_P(b)), last);
Packit fb9d21
    
Packit fb9d21
    if ((res = mp_int_mul(TEMP(0), MP_DENOM_P(b), TEMP(0))) != MP_OK)
Packit fb9d21
      goto CLEANUP;
Packit fb9d21
    if ((res = mp_int_mul(TEMP(1), MP_DENOM_P(a), TEMP(1))) != MP_OK)
Packit fb9d21
      goto CLEANUP;
Packit fb9d21
    if ((res = (comb_f)(TEMP(0), TEMP(1), MP_NUMER_P(c))) != MP_OK)
Packit fb9d21
      goto CLEANUP;
Packit fb9d21
Packit fb9d21
    res = mp_int_mul(MP_DENOM_P(a), MP_DENOM_P(b), MP_DENOM_P(c));
Packit fb9d21
Packit fb9d21
  CLEANUP:
Packit fb9d21
    while (--last >= 0) 
Packit fb9d21
      mp_int_clear(TEMP(last));
Packit fb9d21
Packit fb9d21
    if (res == MP_OK)
Packit fb9d21
      return s_rat_reduce(c);
Packit fb9d21
    else
Packit fb9d21
      return res;
Packit fb9d21
  }
Packit fb9d21
}
Packit fb9d21
Packit fb9d21
/* }}} */
Packit fb9d21
Packit fb9d21
/* Here there be dragons */