
/* this file interfaces to arbitrary precision math routines in libmpfr

   Copyright (C) 2006 Dennis Furey

   This program is free software; you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation; either version 2, or (at your option)
   any later version.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program; if not, write to the Free Software Foundation,
   Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.

*/

#include <avm/common.h>
#include <avm/error.h>
#include <avm/lists.h>
#include <avm/compare.h>
#include <avm/listfuns.h>
#include <avm/chrcodes.h>
#include <avm/mpfr.h>
#include <time.h>

#ifndef HAVE_MEMMOVE
extern void 
*memmove(char *dest, const char *source, unsigned length)
#endif

#define DEFAULT_PREC 160   /* used by the dbl2mp function */

/* non-zero means static variables are initialized */
static int initialized = 0;

/* error messages as lists of lists of character representations */
static list shared_cell = NULL;
static list mpfr_error = NULL;
static list bad_mpfr_spec = NULL;
static list mpfr_overflow = NULL;
static list memory_overflow = NULL;
static list inf_string = NULL;
static list ninf_string = NULL;
static list nan_string = NULL;
static list unrecognized_function_name = NULL;

/* function names as lists of lists of character representations */
static list wild = NULL;
static list funs = NULL;

/* kinds of numbers */
static list unsigned_nan = NULL;
static list negative_inf = NULL;
static list negative_zero = NULL;
static list negative_regular = NULL;
static list positive_inf = NULL;
static list positive_zero = NULL;
static list positive_regular = NULL;

#if HAVE_MPFR
#if HAVE_GMP

static gmp_randstate_t state;

typedef int (*mpfr_mutator)(mpfr_t);
typedef int (*mpfr_unary_operator)(mpfr_t,mpfr_t,mp_rnd_t);
typedef int (*mpfr_binary_operator)(mpfr_t,mpfr_t,mpfr_t,mp_rnd_t);
typedef int (*mpfr_unary_predicate)(mpfr_t);
typedef int (*mpfr_binary_predicate)(mpfr_t,mpfr_t);
typedef int (*mpfr_binary_ui)(mpfr_t,mpfr_t,unsigned long,mp_rnd_t);
typedef int (*mpfr_constant)(mpfr_t,mp_rnd_t);

#endif
#endif


avm_mpfr_ptr
avm_mpfr_of_list(operand, message, fault)
     list operand;
     list *message;
     int *fault;

     /* This function is analogous to the avm_value_of_list function
        in listfuns.c but is specific to lists representing mpfr
        numbers. The operand is required to be in the form
        ((prec,(sign,kind)),((esign,exponent),mantissa)), where sign
        and esign are boolean, mantissa is a list of characters, and
        the others are naturals as lists of booleans lsb first. The
        value field of the mantissa caches the array of limbs used in
        the mpfr representation and the value field of the whole
        operand caches the record referenced by the avm_mpfr_ptr,
        whose mantissa field must point to the array referenced by the
        value field of the operand's mantissa. The caches need not be
        initialized on entry (i.e., they can be NULL pointers) but
        they will be initialized to the results generated by this
        function as a side effect. */
{
#if HAVE_MPFR
  int kind;
  mp_prec_t prec;
  mp_exp_t exponent;
  void *mantissa;
  avm_mpfr_ptr x;

  *message = NULL;
  exponent = prec = kind = 0;
  *fault = ! (operand?(operand->head?(operand->head->tail?(operand->tail?operand->tail->head:NULL):NULL):NULL):NULL);
  if (*fault ? NULL : operand->value)
    {
      *fault = !(operand->tail->tail->value);
      if (!*fault)
	*fault = (((void *) operand->tail->tail->value) != mpfr_custom_get_mantissa(*((avm_mpfr_ptr) operand->value)));
      if (!*fault)
	return ((avm_mpfr_ptr) operand->value);
    }
  if (!*fault)
    {
      prec = (mp_prec_t) avm_counter(operand->head->head);
      kind = (int) avm_counter(operand->head->tail->tail);
      if (operand->head->tail->head)
	{
	  kind = -kind;
	  *fault = (operand->head->tail->head->head ? 1 : !!(operand->head->tail->head->tail));
	}
      exponent = (mp_exp_t) avm_counter(operand->tail->head->tail);
      if (*fault ? NULL : operand->tail->head->head)
	{
	  exponent = -exponent;
	  *fault = (operand->tail->head->head->head ? 1 : !!(operand->tail->head->head->tail));
	}
      *fault = (*fault ? 1 : (((prec < MPFR_PREC_MIN) ? 1 : (prec > MPFR_PREC_MAX))));
    }
  if (!*fault)
    mantissa = (void *) avm_value_of_list(operand->tail->tail,message,fault);
  if (*fault)
    {
      *message = (*message ? *message : avm_copied(bad_mpfr_spec));
      return NULL;
    }
  x = (avm_mpfr_ptr) malloc(sizeof(mpfr_t));
  if (*fault = !x)
    {
      *message = avm_copied(memory_overflow);
      return NULL;
    }
  mpfr_custom_init_set(*x, kind, exponent, prec, mantissa);
  operand->value = (void *) x;
  return x;
#else
  *message = avm_copied(mpfr_error);
  *fault = 1;
  return NULL;
#endif
}







list
avm_list_of_mpfr(x, fault)
     avm_mpfr_ptr x;
     int *fault;

     /* This function is analogous to the avm_list_of_value function
	in listfuns.c but is specific to lists representing mpfr
	numbers. The list it returns is in the form
	((prec,(sign,kind)),((esign,exponent),mantissa)) as specified
	above, with caches initialized. */
{
#if HAVE_MPFR
  mp_prec_t prec;
  mp_exp_t exponent;
  void *mantissa;
  list prec_sign_kind,esign_exponent,result;

  mantissa = mpfr_custom_get_mantissa(*x);
  if (*fault = !mantissa)
    return avm_copied(bad_mpfr_spec);
  prec = mpfr_get_prec(*x);
  prec_sign_kind = NULL;
  switch (mpfr_custom_get_kind(*x))
    {
    case MPFR_NAN_KIND:
      prec_sign_kind = avm_recoverable_join(avm_natural((counter) prec),avm_copied(unsigned_nan));
      break;
    case MPFR_INF_KIND:
      prec_sign_kind = avm_recoverable_join(avm_natural((counter) prec),avm_copied(positive_inf));
      break;
    case MPFR_ZERO_KIND:
      prec_sign_kind = avm_recoverable_join(avm_natural((counter) prec),avm_copied(positive_zero));
      break;
    case MPFR_REGULAR_KIND:
      prec_sign_kind = avm_recoverable_join(avm_natural((counter) prec),avm_copied(positive_regular));
      break;
    case -MPFR_INF_KIND:
      prec_sign_kind = avm_recoverable_join(avm_natural((counter) prec),avm_copied(negative_inf));
      break;
    case -MPFR_ZERO_KIND:
      prec_sign_kind = avm_recoverable_join(avm_natural((counter) prec),avm_copied(negative_zero));
      break;
    case -MPFR_REGULAR_KIND:
      prec_sign_kind = avm_recoverable_join(avm_natural((counter) prec),avm_copied(negative_regular));
      break;
    default:
      *fault = 1;
    }
  if (*fault ? 1 : (*fault = (prec < MPFR_PREC_MIN ? 1 : (prec > MPFR_PREC_MAX))))
    {
      avm_dispose(prec_sign_kind);
      return avm_copied(bad_mpfr_spec);
    }
  if (*fault = ! (prec_sign_kind ? prec_sign_kind->head : NULL))
    {
      avm_dispose(prec_sign_kind);
      free(mantissa);
      free(x);
      return avm_copied(memory_overflow);
    }
  exponent = mpfr_custom_get_exp(*x);
  if (exponent == 0)
    esign_exponent = avm_copied(shared_cell);
  else
    {
      if (exponent < 0)
	esign_exponent = avm_recoverable_join(avm_copied(shared_cell),avm_natural((counter) -exponent));
      else
	esign_exponent = avm_recoverable_join(NULL,avm_natural((counter) exponent));
      *fault = ! (esign_exponent ? esign_exponent->tail : NULL);
    }
  result = (*fault ? NULL : avm_list_of_value((void *) mantissa,mpfr_custom_get_size(prec),fault));
  if (*fault = (*fault ? 1 : !result))
    {
      avm_dispose(prec_sign_kind);
      avm_dispose(esign_exponent);
      free(mantissa);
      free(*x);
      return avm_copied(memory_overflow);
    }
  result->discontiguous = 1;
  result = avm_recoverable_join(prec_sign_kind,avm_recoverable_join(esign_exponent,result));
  if (*fault = ! (result ? result->tail : NULL))
    {
      avm_dispose(result);
      free(mantissa);
      free(*x);
      return avm_copied(memory_overflow);
    }
  result->tail->discontiguous = 1;
  result->discontiguous = 1;
  mpfr_custom_move(*x,(void *) result->tail->tail->value);
  result->value = (void *) x;
  free(mantissa);
  return result;
#else
  *fault = 1;
  return avm_copied(mpfr_error);
#endif
}



#if HAVE_MPFR


static list
unary_evaluation(operator, operand, fault)
     mpfr_unary_operator operator;
     list operand;
     int *fault;

     /* performs a unary operation on mpfr numbers */
{
  avm_mpfr_ptr x;
  avm_mpfr_ptr y;
  list message;
  mp_prec_t prec;
  void *mantissa;

  message = NULL;
  x = avm_mpfr_of_list(operand,&message,fault);
  if (*fault = (*fault ? 1 : !!message))
    return message;
  if (!x)
    avm_internal_error(61);
  prec = mpfr_get_prec(*x);
  if (*fault = ! (mantissa = (void *) malloc(mpfr_custom_get_size(prec))))
    return avm_copied(memory_overflow);
  if (*fault = !(y = (avm_mpfr_ptr) malloc(sizeof(mpfr_t))))
    {
      free(mantissa);
      return avm_copied(memory_overflow);
    }
  mpfr_custom_init(mantissa,prec);
  mpfr_custom_init_set(*y,MPFR_NAN_KIND,0,prec,mantissa);
  (*operator)(*y,*x,GMP_RNDN);
  return avm_list_of_mpfr(y,fault);
}







static list
binary_ui_evaluation(operator, operand, fault)
     mpfr_binary_ui operator;
     list operand;
     int *fault;

     /* performs a binary operation on mpfr numbers with integers */
{
  avm_mpfr_ptr x;
  unsigned long y;
  avm_mpfr_ptr z;
  list message;
  mp_prec_t prec;
  void *mantissa;

  message = NULL;
  if (*fault = !operand)
    return avm_copied(bad_mpfr_spec);
  x = avm_mpfr_of_list(operand->head,&message,fault);
  if (*fault = (*fault ? 1 : !!message))
    return message;
  if (!x)
    avm_internal_error(63);
  y = (unsigned long) avm_counter(operand->tail);
  prec = mpfr_get_prec(*x);
  if (*fault = ! (mantissa = (void *) malloc(mpfr_custom_get_size(prec))))
    return avm_copied(memory_overflow);
  if (*fault = !(z = (avm_mpfr_ptr) malloc(sizeof(mpfr_t))))
    {
      free(mantissa);
      return avm_copied(memory_overflow);
    }
  mpfr_custom_init(mantissa,prec);
  mpfr_custom_init_set(*z,MPFR_NAN_KIND,0,prec,mantissa);
  (*operator)(*z,*x,y,GMP_RNDN);
  return avm_list_of_mpfr(z,fault);
}






static list
binary_evaluation(operator, operand, fault)
     mpfr_binary_operator operator;
     list operand;
     int *fault;

     /* This performs a binary operation on mpfr numbers. The
	precision of the result is the higher of the two input
	precisions. Although mathematically it makes more sense for it
	to be the lower precision, this convention is more convenient
	for functions with embedded constants (e.g., a function that
	adds 1 to everything) because they can be used on arguments of
	any precision without needing different versions. */
{
  avm_mpfr_ptr x,y,z;
  list message;
  mp_prec_t x_prec,y_prec,prec;
  void *mantissa;

  message = NULL;
  if (*fault = !operand)
    return avm_copied(bad_mpfr_spec);
  x = avm_mpfr_of_list(operand->head,&message,fault);
  if (*fault = (*fault ? 1 : !!message))
    return message;
  y = avm_mpfr_of_list(operand->tail,&message,fault);
  if (*fault = (*fault ? 1 : !!message))
    return message;
  if (x ? !y : 1)
    avm_internal_error(62);
  x_prec = mpfr_get_prec(*x);
  y_prec = mpfr_get_prec(*y);
  prec = ((x_prec < y_prec) ? y_prec : x_prec);
  if (*fault = ! (mantissa = (void *) malloc(mpfr_custom_get_size(prec))))
    return avm_copied(memory_overflow);
  if (*fault = !(z = (avm_mpfr_ptr) malloc(sizeof(*z))))
    {
      free(mantissa);
      return avm_copied(memory_overflow);
    }
  mpfr_custom_init(mantissa,prec);
  mpfr_custom_init_set(*z,MPFR_NAN_KIND,0,prec,mantissa);
  (*operator)(*z,*x,*y,GMP_RNDN);
  return avm_list_of_mpfr(z,fault);
}












static list
vid_evaluation(operand, fault)
     list operand;
     int *fault;

     /* special case of binary evaluation for the inverse quotient,
	which isn't defined in the mpfr library */
{
  avm_mpfr_ptr x,y,z;
  list message;
  mp_prec_t x_prec,y_prec,prec;
  void *mantissa;

  message = NULL;
  if (*fault = !operand)
    return avm_copied(bad_mpfr_spec);
  x = avm_mpfr_of_list(operand->head,&message,fault);
  if (*fault = (*fault ? 1 : !!message))
    return message;
  y = avm_mpfr_of_list(operand->tail,&message,fault);
  if (*fault = (*fault ? 1 : !!message))
    return message;
  if (x ? !y : 1)
    avm_internal_error(106);
  x_prec = mpfr_get_prec(*x);
  y_prec = mpfr_get_prec(*y);
  prec = ((x_prec < y_prec) ? y_prec : x_prec);
  if (*fault = ! (mantissa = (void *) malloc(mpfr_custom_get_size(prec))))
    return avm_copied(memory_overflow);
  if (*fault = !(z = (avm_mpfr_ptr) malloc(sizeof(*z))))
    {
      free(mantissa);
      return avm_copied(memory_overflow);
    }
  mpfr_custom_init(mantissa,prec);
  mpfr_custom_init_set(*z,MPFR_NAN_KIND,0,prec,mantissa);
  mpfr_div(*z,*y,*x,GMP_RNDN);
  return avm_list_of_mpfr(z,fault);
}












static list
bus_evaluation(operand, fault)
     list operand;
     int *fault;

     /* special case of binary evaluation for the inverse difference,
	which isn't defined in the mpfr library */
{
  avm_mpfr_ptr x,y,z;
  list message;
  mp_prec_t x_prec,y_prec,prec;
  void *mantissa;

  message = NULL;
  if (*fault = !operand)
    return avm_copied(bad_mpfr_spec);
  x = avm_mpfr_of_list(operand->head,&message,fault);
  if (*fault = (*fault ? 1 : !!message))
    return message;
  y = avm_mpfr_of_list(operand->tail,&message,fault);
  if (*fault = (*fault ? 1 : !!message))
    return message;
  if (x ? !y : 1)
    avm_internal_error(108);
  x_prec = mpfr_get_prec(*x);
  y_prec = mpfr_get_prec(*y);
  prec = ((x_prec < y_prec) ? y_prec : x_prec);
  if (*fault = ! (mantissa = (void *) malloc(mpfr_custom_get_size(prec))))
    return avm_copied(memory_overflow);
  if (*fault = !(z = (avm_mpfr_ptr) malloc(sizeof(*z))))
    {
      free(mantissa);
      return avm_copied(memory_overflow);
    }
  mpfr_custom_init(mantissa,prec);
  mpfr_custom_init_set(*z,MPFR_NAN_KIND,0,prec,mantissa);
  mpfr_sub(*z,*y,*x,GMP_RNDN);
  return avm_list_of_mpfr(z,fault);
}









static list
unary_predicate_evaluation(operator, operand, fault)
     mpfr_unary_predicate operator;
     list operand;
     int *fault;

     /* evaluates a unary predicate on mpfr numbers */
{
  avm_mpfr_ptr x;
  list message;

  message = NULL;
  x = avm_mpfr_of_list(operand,&message,fault);
  if (*fault = (*fault ? 1 : !!message))
    return message;
  return ((*operator)(*x) ? avm_copied(shared_cell) : NULL);
}








static list
binary_predicate_evaluation(operator, operand, fault)
     mpfr_binary_predicate operator;
     list operand;
     int *fault;

     /* evaluates a binary predicate on mpfr numbers */
{
  avm_mpfr_ptr x,y;
  list message;

  message = NULL;
  if (*fault = !operand)
    return avm_copied(bad_mpfr_spec);
  x = avm_mpfr_of_list(operand->head,&message,fault);
  if (*fault = (*fault ? 1 : !!message))
    return message;
  y = avm_mpfr_of_list(operand->tail,&message,fault);
  if (*fault = (*fault ? 1 : !!message))
    return message;
  return ((*operator)(*x,*y) ? avm_copied(shared_cell) : NULL);
}







static list
constant_evaluation(operator, operand, fault)
     mpfr_constant operator;
     list operand;
     int *fault;

     /* returns a list representing an mpfr constant evaluated to a
	specified precision */
{
  avm_mpfr_ptr y;
  list message;
  mp_prec_t prec;
  void *mantissa;

  prec = (mp_prec_t) avm_counter(operand);
  if (*fault = (prec < MPFR_PREC_MIN ? 1 : (prec > MPFR_PREC_MAX)))
    return avm_copied(bad_mpfr_spec);
  if (*fault = ! (mantissa = (void *) malloc(mpfr_custom_get_size(prec))))
    return avm_copied(memory_overflow);
  if (*fault = !(y = (avm_mpfr_ptr) malloc(sizeof(*y))))
    {
      free(mantissa);
      return avm_copied(memory_overflow);
    }
  mpfr_custom_init(mantissa,prec);
  mpfr_custom_init_set(*y,MPFR_NAN_KIND,0,prec,mantissa);
  (*operator)(*y,GMP_RNDN);
  return avm_list_of_mpfr(y,fault);
}





static list
infinity(direction, operand, fault)
     int direction;
     list operand;
     int *fault;
{
  avm_mpfr_ptr y;
  list message;
  mp_prec_t prec;
  void *mantissa;

  prec = (mp_prec_t) avm_counter(operand);
  if (*fault = (prec < MPFR_PREC_MIN ? 1 : (prec > MPFR_PREC_MAX)))
    return avm_copied(bad_mpfr_spec);
  if (*fault = ! (mantissa = (void *) malloc(mpfr_custom_get_size(prec))))
    return avm_copied(memory_overflow);
  if (*fault = !(y = (avm_mpfr_ptr) malloc(sizeof(*y))))
    {
      free(mantissa);
      return avm_copied(memory_overflow);
    }
  mpfr_custom_init(mantissa,prec);
  mpfr_custom_init_set(*y,MPFR_NAN_KIND,0,prec,mantissa);
  mpfr_set_inf(*y,direction);
  return avm_list_of_mpfr(y,fault);
}








static list
nan(operand, fault)
     list operand;
     int *fault;
{
  avm_mpfr_ptr y;
  list message;
  mp_prec_t prec;
  void *mantissa;

  prec = (mp_prec_t) avm_counter(operand);
  if (*fault = (prec < MPFR_PREC_MIN ? 1 : (prec > MPFR_PREC_MAX)))
    return avm_copied(bad_mpfr_spec);
  if (*fault = ! (mantissa = (void *) malloc(mpfr_custom_get_size(prec))))
    return avm_copied(memory_overflow);
  if (*fault = !(y = (avm_mpfr_ptr) malloc(sizeof(*y))))
    {
      free(mantissa);
      return avm_copied(memory_overflow);
    }
  mpfr_custom_init(mantissa,prec);
  mpfr_custom_init_set(*y,MPFR_NAN_KIND,0,prec,mantissa);
  mpfr_set_nan(*y);
  return avm_list_of_mpfr(y,fault);
}










static list
precision(operand, fault)
     list operand;
     int *fault;

     /* returns the precision of an mpfr number as a natural */
{
  if (*fault = !(operand ? operand->head : NULL))
    return avm_copied(bad_mpfr_spec);
  return avm_copied (operand->head->head);
}






static list
sin_cos(operand, fault)
     list operand;
     int *fault;

     /* returns the mpfr sin_cos function of an argument, currently
	the only mpfr function that returns a pair of values */
{
  avm_mpfr_ptr x,y,z;
  list message;
  mp_prec_t prec;
  void *x_mantissa;
  void *y_mantissa;
  list x_list,y_list,result;

  message = NULL;
  z = avm_mpfr_of_list(operand,&message,fault);
  if (*fault = (*fault ? 1 : !!message))
    return message;
  if (!z)
    avm_internal_error(64);
  prec = mpfr_get_prec(*z);
  if (*fault = ! (x_mantissa = (void *) malloc(mpfr_custom_get_size(prec))))
    return avm_copied(memory_overflow);
  if (*fault = !(x = (avm_mpfr_ptr) malloc(sizeof(*x))))
    {
      free(x_mantissa);
      return avm_copied(memory_overflow);
    }
  if (*fault = ! (y_mantissa = (void *) malloc(mpfr_custom_get_size(prec))))
    return avm_copied(memory_overflow);
  if (*fault = !(y = (avm_mpfr_ptr) malloc(sizeof(*y))))
    {
      free(x);
      free(x_mantissa);
      free(y_mantissa);
      return avm_copied(memory_overflow);
    }
  mpfr_custom_init(x_mantissa,prec);
  mpfr_custom_init(y_mantissa,prec);
  mpfr_custom_init_set(*x,MPFR_NAN_KIND,0,prec,x_mantissa);
  mpfr_custom_init_set(*y,MPFR_NAN_KIND,0,prec,y_mantissa);
  mpfr_sin_cos(*y,*x,*z,GMP_RNDN);
  x_list = avm_list_of_mpfr(x,fault);
  if (*fault)
    {
      free(y);
      free(y_mantissa);
      return x_list;
    }
  y_list = avm_list_of_mpfr(y,fault);
  if (*fault)
    {
      avm_dispose(x_list);
      return y_list;
    }
  result = avm_recoverable_join(y_list,x_list);
  return ((*fault = !result) ? avm_copied(memory_overflow) : result);
}







static list
mutator_evaluation(operator, operand, fault)
     mpfr_mutator operator;
     list operand;
     int *fault;

     /* implements a non-destructive interface to mpfr functions with
	no rounding mode parameter that modify their argument in place
	(e.g., nextabove and nextbelow) */
{
  avm_mpfr_ptr x;
  avm_mpfr_ptr y;
  list message;
  mp_prec_t prec;
  void *mantissa;

  message = NULL;
  x = avm_mpfr_of_list(operand,&message,fault);
  if (*fault = (*fault ? 1 : !!message))
    return message;
  if (!x)
    avm_internal_error(65);
  prec = mpfr_get_prec(*x);
  if (*fault = ! (mantissa = (void *) malloc(mpfr_custom_get_size(prec))))
    return avm_copied(memory_overflow);
  if (*fault = !(y = (avm_mpfr_ptr) malloc(sizeof(mpfr_t))))
    {
      free(mantissa);
      return avm_copied(memory_overflow);
    }
  mpfr_custom_init(mantissa,prec);
  mpfr_custom_init_set(*y,MPFR_NAN_KIND,0,prec,mantissa);
  mpfr_set(*y,*x,GMP_RNDN);
  (*operator)(*y);
  return avm_list_of_mpfr(y,fault);
}









static list
grow(direction, operand, fault)
     int direction;
     list operand;
     int *fault;

     /* returns a copy of the operand with different precision */
{
  avm_mpfr_ptr x,y;
  list message;
  mp_prec_t prec,x_prec;
  void *mantissa;
  int growth;

  message = NULL;
  if (*fault = !operand)
    return avm_copied(bad_mpfr_spec);
  x = avm_mpfr_of_list(operand->head,&message,fault);
  if (*fault = (*fault ? 1 : !!message))
    return message;
  if (!x)
    avm_internal_error(66);
  growth = (int) avm_counter(operand->tail);
  if (*fault = (operand->tail ? !growth : 0))
    return avm_copied(memory_overflow);
  x_prec = mpfr_get_prec(*x);
  if (direction > 0)
    {
      prec = x_prec + growth;
      prec = ((prec > MPFR_PREC_MAX) ? MPFR_PREC_MAX : prec);
      prec = ((prec < x_prec) ? MPFR_PREC_MAX : prec);
    }
  else
    {
      prec = x_prec - growth;
      prec = ((prec < MPFR_PREC_MIN) ? MPFR_PREC_MIN : prec);
      prec = ((prec > x_prec) ? MPFR_PREC_MIN : prec);
    }
  if (*fault = ! (mantissa = (void *) malloc(mpfr_custom_get_size(prec))))
    return avm_copied(memory_overflow);
  memset (mantissa, 0, mpfr_custom_get_size (prec));
  if (*fault = !(y = (avm_mpfr_ptr) malloc(sizeof(*y))))
    {
      free(mantissa);
      return avm_copied(memory_overflow);
    }
  mpfr_custom_init(mantissa,prec);
  mpfr_custom_init_set(*y,MPFR_ZERO_KIND,0,prec,mantissa);
  mpfr_set(*y,*x,GMP_RNDN);
  return avm_list_of_mpfr(y,fault);
}







static list
dbl2mp(operand, fault)
     list operand;
     int *fault;

     /* conversion from IEEE double precision to mpfr format with a
	specified precision; the operand should represent a double precision number;
	the precision of the result is a hard coded constant DEFAULT_PREC */
{
  avm_mpfr_ptr y;
  double *x;
  list message;
  void *mantissa;

  if (*fault = !operand)
    return avm_copied(bad_mpfr_spec);
  message = NULL;
  x = (double *) avm_value_of_list(operand,&message,fault);
  if (*fault = (*fault ? 1 : !!message))
    return avm_copied(message);
  if (*fault = ! (mantissa = (void *) malloc(mpfr_custom_get_size(DEFAULT_PREC))))
    return avm_copied(memory_overflow);
  if (*fault = !(y = (avm_mpfr_ptr) malloc(sizeof(*y))))
    {
      free(mantissa);
      return avm_copied(memory_overflow);
    }
  mpfr_custom_init(mantissa,DEFAULT_PREC);
  mpfr_custom_init_set(*y,MPFR_NAN_KIND,0,DEFAULT_PREC,mantissa);
  mpfr_set_d(*y,*x,GMP_RNDN);
  return avm_list_of_mpfr(y,fault);
}






static list
mp2dbl(operand, fault)
     list operand;
     int *fault;

     /* converts an mpfr number to a double */
{
  avm_mpfr_ptr x;
  double y;
  list message;

  message = NULL;
  x = avm_mpfr_of_list(operand, &message, fault);
  if (*fault = (*fault ? 1 : !!message))
    return message;
  if (!x)
    avm_internal_error(68);
  y = mpfr_get_d(*x,GMP_RNDN);
  return avm_list_of_value((void *) &y,sizeof(double),fault);
}







static list
str2mp(operand, fault)
     list operand;
     int *fault;

     /* operand should be of the form (prec,string) */
{
  char *string;
  list message;
  avm_mpfr_ptr y;
  mpfr_prec_t prec;
  void *mantissa;

  if (*fault = !operand)
    return avm_copied(bad_mpfr_spec);
  message = NULL;
  string = avm_standard_unstrung (operand->tail, &message, fault);
  if (*fault = (*fault ? 1 : !!message))
    return message;
  avm_dispose(message);
  prec = (mp_prec_t) avm_counter(operand->head);
  if (prec < MPFR_PREC_MIN)
    prec = MPFR_PREC_MIN;
  if (*fault = (prec > MPFR_PREC_MAX))
    {
      free(string);
      return avm_copied(bad_mpfr_spec);
    }
  if (*fault = ! (mantissa = (void *) malloc(mpfr_custom_get_size(prec))))
    {
      free(string);
      return avm_copied(memory_overflow);
    }
  if (*fault = !(y = (avm_mpfr_ptr) malloc(sizeof(*y))))
    {
      free(mantissa);
      free(string);
      return avm_copied(memory_overflow);
    }
  mpfr_custom_init(mantissa,prec);
  mpfr_custom_init_set(*y,MPFR_NAN_KIND,0,prec,mantissa);
  if (*fault = (mpfr_set_str(*y,string,0,GMP_RNDN) != 0))
    {
      free(string);
      free(mantissa);
      free(y);
      return avm_copied(bad_mpfr_spec);
    }
  free(string);
  return avm_list_of_mpfr(y,fault);
}







static list
mp2str(operand, fault)
     list operand;
     int *fault;

     /* converts an mpfr number to a string in standard exponential
        notation, with sufficiently many digits to express the
        precision */
{
  avm_mpfr_ptr x;
  list message;
  size_t n;
  char *str;
  mp_exp_t exponent;
  char exp_str[81];
  list result;

  message = NULL;
  x = avm_mpfr_of_list(operand,&message,fault);
  if (*fault = (*fault ? 1 : !!message))
    return message;
  if (!x)
    avm_internal_error(69);
  if (mpfr_nan_p(*x))
    return avm_copied(nan_string);
  else if (mpfr_inf_p(*x))
    {
      if (mpfr_sgn(*x) > 0)
	return avm_copied(inf_string);
      else
	return avm_copied(ninf_string);
    }
  n = 2 + ((mpfr_get_prec(*x) * 3) / 10);
  n = (n > 2 ? n : 3);
  if (*fault = ! (str = (char *) malloc(n + 84)))
    return avm_copied(memory_overflow);
  if (*fault = (mpfr_get_str(str,&exponent,10,n,*x,GMP_RNDN) != str))
    {
      free(str);
      return avm_copied(mpfr_error);
    }
  if (str[0] == '-')
    {
      memmove(&str[3],&str[2],strlen(str) - 1);
      str[2] = '.';
    }
  else
    {
      memmove(&str[2],&str[1],strlen(str));
      str[1] = '.';
    }
  exp_str[0] = 'E';
  exp_str[1] = '+';
  if ((exponent - 1 > 9) ? 1 : (exponent - 1 < -9))
    snprintf(&exp_str[(exponent < 0) ? 1 : 2],80,"%ld",exponent - 1);
  else if (exponent > 0)
    snprintf(&exp_str[2],80,"0%ld",exponent - 1);
  else if (mpfr_zero_p(*x))
    {
      exp_str[2] = '0';
      exp_str[3] = '0';
      exp_str[4] = '\0';
    }
  else 
    {
      exp_str[1] = '-';
      exp_str[2] = '0';
      snprintf(&exp_str[3],80,"%ld",1 - exponent);
    }
  strncat(str,exp_str,81);
  result = avm_recoverable_strung(str,fault);
  free(str);
  if (*fault = (*fault ? 1 : !result))
    {
      avm_dispose(result);
      return avm_copied(memory_overflow);
    }
  return result;
}







static list
nat2mp(operand, fault)
     list operand;
     int *fault;

     /* The operand represents a natural number as a list of bits lsb
	first. The result is a list representing an mpfr number of the
	implied precision. */
{
  avm_mpfr_ptr y;
  list temporary,reversal;
  mp_prec_t prec;
  void *mantissa;

  prec = (mp_prec_t) avm_recoverable_length (operand);
  if (*fault = (operand ? !prec : 0))
    return avm_copied (mpfr_overflow);
  if (prec < MPFR_PREC_MIN)
    prec = MPFR_PREC_MIN;
  if (*fault = (prec > MPFR_PREC_MAX))
    return avm_copied(mpfr_overflow);
  if (*fault = ! (mantissa = (void *) malloc(mpfr_custom_get_size(prec))))
    return avm_copied(memory_overflow);
  if (*fault = !(y = (avm_mpfr_ptr) malloc(sizeof(*y))))
    {
      free(mantissa);
      return avm_copied(memory_overflow);
    }
  mpfr_custom_init(mantissa,prec);
  mpfr_custom_init_set(*y,MPFR_ZERO_KIND,0,prec,mantissa);
  reversal = NULL;
  while (operand)
    {
      operand = (temporary = operand)->tail;
      temporary->tail = reversal;
      reversal = temporary;
    }
  while (reversal)
    {
      mpfr_mul_ui (*y, *y, 2, GMP_RNDN);
      if (reversal->head)
	mpfr_add_ui (*y, *y, 1, GMP_RNDN);
      reversal = (temporary = reversal)->tail;
      temporary->tail = operand;
      operand = temporary;
    }
  return avm_list_of_mpfr(y,fault);
}







#if HAVE_GMP


static list
approximately_equal(operand, fault)
     list operand;
     int* fault;

     /* operand should be of the form (prec,x,y), where prec is a
	natural. The result is a true list if x and y have the same
	exponent and agree in the most significant prec bits. */
{
  unsigned long prec;
  avm_mpfr_ptr x,y;
  list message;

  if (*fault = ! (operand ? operand->tail : NULL))
    return avm_copied(bad_mpfr_spec);
  prec = (unsigned long) avm_counter(operand->head);
  message = NULL;
  x = avm_mpfr_of_list(operand->tail->head,&message,fault);
  if (*fault = (*fault ? 1 : !!message))
    return message;
  y = avm_mpfr_of_list(operand->tail->tail,&message,fault);
  if (*fault = (*fault ? 1 : !!message))
    return message;
  return (mpfr_eq(*x,*y,prec) ? avm_copied(shared_cell) : NULL);
}






static list
urandomb(operand, fault)
     list operand;
     int *fault;

     /* the operand is a natural specifying the precision */
{
  avm_mpfr_ptr y;
  list message;
  mp_prec_t prec;
  void *mantissa;

  prec = (mp_prec_t) avm_counter(operand);
  if (*fault = (prec < MPFR_PREC_MIN ? 1 : (prec > MPFR_PREC_MAX)))
    return avm_copied(bad_mpfr_spec);
  if (*fault = ! (mantissa = (void *) malloc(mpfr_custom_get_size(prec))))
    return avm_copied(memory_overflow);
  if (*fault = !(y = (avm_mpfr_ptr) malloc(sizeof(*y))))
    {
      free(mantissa);
      return avm_copied(memory_overflow);
    }
  mpfr_custom_init(mantissa,prec);
  mpfr_custom_init_set(*y,MPFR_NAN_KIND,0,prec,mantissa);
  mpfr_urandomb(*y,state);
  return avm_list_of_mpfr(y,fault);
}






#endif /* HAVE_GMP */
#endif /* HAVE_MPFR */






list
avm_have_mpfr_call (function_name,fault)
     list function_name;
     int *fault;

/* this reports the availability of a function */
{
#if HAVE_MPFR
  list membership;
  list comparison;
  list result;

  if (!initialized)
    avm_initialize_mpfr ();
  if (*fault)
    return NULL;
  comparison = avm_binary_comparison (function_name, wild, fault);
  if (*fault)
    return comparison;
  if (comparison)
    {
      avm_dispose(comparison);
      return avm_copied(funs);
    }
  if (!(membership = avm_binary_membership (function_name, funs, fault)) ? 1 : *fault)
    return membership;
  avm_dispose(membership);
  return ((*fault = !(result = avm_recoverable_join(avm_copied(function_name),NULL))) ? avm_copied(memory_overflow) : result);
#endif
  return NULL;
}







list
avm_mpfr_call (function_name, argument, fault)
     list function_name;
     list argument;
     int *fault;
{
#if HAVE_MPFR
  list message;
  int function_number;

  if (*fault)
    return NULL;
  if (! initialized)
    avm_initialize_mpfr ();
  if (!(function_number = 0xff & (function_name ? function_name->characterization : 0)))
    {
      message = avm_position (function_name, funs, fault);
      if (*fault = (*fault ? 1 : !message))
	return (message ? message : avm_copied (unrecognized_function_name));
      function_number = message->characterization;
      function_name->characterization = function_number;
      avm_dispose (message);
    }
  switch (function_number)
    {
    case  1: return binary_evaluation((mpfr_binary_operator) &mpfr_add, argument, fault);
    case  2: return binary_evaluation((mpfr_binary_operator) &mpfr_sub, argument, fault);
    case  3: return binary_evaluation((mpfr_binary_operator) &mpfr_mul, argument, fault);
    case  4: return binary_evaluation((mpfr_binary_operator) &mpfr_div, argument, fault);
    case  5: return binary_evaluation((mpfr_binary_operator) &mpfr_pow, argument, fault);
    case  6: return binary_evaluation((mpfr_binary_operator) &mpfr_atan2, argument, fault);
    case  7: return binary_evaluation((mpfr_binary_operator) &mpfr_hypot, argument, fault);
    case  8: return binary_evaluation((mpfr_binary_operator) &mpfr_min, argument, fault);
    case  9: return binary_evaluation((mpfr_binary_operator) &mpfr_max, argument, fault);
    case 10: return unary_evaluation((mpfr_unary_operator) &mpfr_sqrt, argument, fault);
    case 11: return unary_evaluation((mpfr_unary_operator) &mpfr_sqr, argument, fault);
    case 12: return unary_evaluation((mpfr_unary_operator) &mpfr_cbrt, argument, fault);
    case 13: return unary_evaluation((mpfr_unary_operator) &mpfr_neg, argument, fault);
    case 14: return unary_evaluation((mpfr_unary_operator) &mpfr_abs, argument, fault);
    case 15: return unary_evaluation((mpfr_unary_operator) &mpfr_log, argument, fault);
    case 16: return unary_evaluation((mpfr_unary_operator) &mpfr_log2, argument, fault);
    case 17: return unary_evaluation((mpfr_unary_operator) &mpfr_log10, argument, fault);
    case 18: return unary_evaluation((mpfr_unary_operator) &mpfr_exp, argument, fault);
    case 19: return unary_evaluation((mpfr_unary_operator) &mpfr_exp2, argument, fault);
    case 20: return unary_evaluation((mpfr_unary_operator) &mpfr_exp10, argument, fault);
    case 21: return unary_evaluation((mpfr_unary_operator) &mpfr_cos, argument, fault);
    case 22: return unary_evaluation((mpfr_unary_operator) &mpfr_sin, argument, fault);
    case 23: return unary_evaluation((mpfr_unary_operator) &mpfr_tan, argument, fault);
    case 24: return unary_evaluation((mpfr_unary_operator) &mpfr_acos, argument, fault);
    case 25: return unary_evaluation((mpfr_unary_operator) &mpfr_asin, argument, fault);
    case 26: return unary_evaluation((mpfr_unary_operator) &mpfr_atan, argument, fault);
    case 27: return unary_evaluation((mpfr_unary_operator) &mpfr_cosh, argument, fault);
    case 28: return unary_evaluation((mpfr_unary_operator) &mpfr_sinh, argument, fault);
    case 29: return unary_evaluation((mpfr_unary_operator) &mpfr_tanh, argument, fault);
    case 30: return unary_evaluation((mpfr_unary_operator) &mpfr_acosh, argument, fault);
    case 31: return unary_evaluation((mpfr_unary_operator) &mpfr_asinh, argument, fault);
    case 32: return unary_evaluation((mpfr_unary_operator) &mpfr_atanh, argument, fault);
    case 33: return unary_evaluation((mpfr_unary_operator) &mpfr_log1p, argument, fault);
    case 34: return unary_evaluation((mpfr_unary_operator) &mpfr_expm1, argument, fault);
    case 35: return unary_evaluation((mpfr_unary_operator) &mpfr_eint, argument, fault);
    case 36: return unary_evaluation((mpfr_unary_operator) &mpfr_gamma, argument, fault);
    case 37: return unary_evaluation((mpfr_unary_operator) &mpfr_lngamma, argument, fault);
    case 38: return unary_evaluation((mpfr_unary_operator) &mpfr_erf, argument, fault);
    case 39: return unary_evaluation((mpfr_unary_operator) &mpfr_erfc, argument, fault);
    case 40: return unary_evaluation((mpfr_unary_operator) &mpfr_ceil, argument, fault);
    case 41: return unary_evaluation((mpfr_unary_operator) &mpfr_floor, argument, fault);
    case 42: return unary_evaluation((mpfr_unary_operator) &mpfr_round, argument, fault);
    case 43: return unary_evaluation((mpfr_unary_operator) &mpfr_trunc, argument, fault);
    case 44: return unary_evaluation((mpfr_unary_operator) &mpfr_frac, argument, fault);
    case 45: return binary_ui_evaluation((mpfr_binary_ui) &mpfr_root, argument, fault);
    case 46: return binary_ui_evaluation((mpfr_binary_ui) &mpfr_pow_ui, argument, fault);
    case 47: return binary_ui_evaluation((mpfr_binary_ui) &mpfr_mul_2ui, argument, fault);
    case 48: return binary_ui_evaluation((mpfr_binary_ui) &mpfr_div_2ui, argument, fault);
    case 49: return binary_predicate_evaluation((mpfr_binary_predicate) &mpfr_equal_p, argument, fault);
    case 50: return binary_predicate_evaluation((mpfr_binary_predicate) &mpfr_cmp_abs, argument, fault);
    case 51: return binary_predicate_evaluation((mpfr_binary_predicate) &mpfr_greater_p, argument, fault);
    case 52: return binary_predicate_evaluation((mpfr_binary_predicate) &mpfr_greaterequal_p, argument, fault);
    case 53: return binary_predicate_evaluation((mpfr_binary_predicate) &mpfr_less_p, argument, fault);
    case 54: return binary_predicate_evaluation((mpfr_binary_predicate) &mpfr_lessequal_p, argument, fault);
    case 55: return binary_predicate_evaluation((mpfr_binary_predicate) &mpfr_lessgreater_p, argument, fault);
    case 56: return unary_predicate_evaluation((mpfr_unary_predicate) &mpfr_nan_p, argument, fault);
    case 57: return unary_predicate_evaluation((mpfr_unary_predicate) &mpfr_inf_p, argument, fault);
    case 58: return unary_predicate_evaluation((mpfr_unary_predicate) &mpfr_number_p, argument, fault);
    case 59: return unary_predicate_evaluation((mpfr_unary_predicate) &mpfr_zero_p, argument, fault);
    case 60: return unary_predicate_evaluation((mpfr_unary_predicate) &mpfr_integer_p, argument, fault);
    case 61: return constant_evaluation((mpfr_constant) &mpfr_const_log2, argument, fault);
    case 62: return constant_evaluation((mpfr_constant) &mpfr_const_pi, argument, fault);
    case 63: return constant_evaluation((mpfr_constant) &mpfr_const_catalan, argument, fault);
    case 64: return precision(argument,fault);
    case 65: return sin_cos(argument,fault);
    case 66: return mutator_evaluation((mpfr_mutator) mpfr_nextabove, argument, fault);
    case 67: return mutator_evaluation((mpfr_mutator) mpfr_nextbelow, argument, fault);
    case 68: return grow(1, argument, fault);
    case 69: return grow(-1, argument, fault);
    case 70: return dbl2mp(argument, fault);
    case 71: return mp2dbl(argument, fault);
    case 72: return str2mp(argument, fault);
    case 73: return mp2str(argument, fault);
    case 74: return infinity(1, argument, fault);
    case 75: return infinity(-1, argument, fault);
    case 76: return nan(argument, fault);
    case 77: return vid_evaluation(argument, fault);
    case 78: return bus_evaluation(argument, fault);
    case 79: return nat2mp(argument, fault);
#if HAVE_GMP
    case 80: return approximately_equal(argument, fault);
    case 81: return urandomb(argument, fault);
#endif 
    }
#endif /* HAVE_MPFR */
  *fault = 1;
  return avm_copied (unrecognized_function_name);
}







void
avm_initialize_mpfr ()

     /* This initializes some static data structures. */

{
  char *funames[] = {
    "add",            /* binary */
    "sub",
    "mul",
    "div",
    "pow",
    "atan2",
    "hypot",
    "min",
    "max",
    "sqrt",            /* unary */
    "sqr",
    "cbrt",
    "neg",
    "abs",
    "log",
    "log2",
    "log10",
    "exp",
    "exp2",
    "exp10",
    "cos",
    "sin",
    "tan",
    "acos",
    "asin",
    "atan",
    "cosh",
    "sinh",
    "tanh",
    "acosh",
    "asinh",
    "atanh",
    "log1p",
    "expm1",
    "eint",
    "gamma",
    "lngamma",
    "erf",
    "erfc",
    "ceil",
    "floor",
    "round",
    "trunc",
    "frac",
    "root",           /* binary ui */
    "pow_ui",
    "mul_2ui",
    "div_2ui",
    "equal_p",            /* binary predicate */
    "unequal_abs",
    "greater_p",
    "greaterequal_p",
    "less_p",
    "lessequal_p",
    "lessgreater_p",
    "nan_p",          /* unary predicate */
    "inf_p",
    "number_p",
    "zero_p",
    "integer_p",
    "const_log2",     /* constants */
    "pi",
    "const_catalan",
    "prec",           /* other */
    "sin_cos",
    "nextabove",
    "nextbelow",
    "grow",
    "shrink",
    "dbl2mp",
    "mp2dbl",
    "str2mp",
    "mp2str",
    "inf",
    "ninf",
    "nan",
    "vid",
    "bus",
    "nat2mp",
#if HAVE_GMP
    "eq",
    "urandomb",
#endif
    NULL};
  list back;
  int string_number;
  unsigned long seed;

  if (initialized)
    return;
  initialized = 1;
  avm_initialize_lists ();
  avm_initialize_chrcodes ();
  wild = avm_strung("*");
  shared_cell = avm_join (NULL,NULL);
  memory_overflow = avm_join (avm_strung ("memory overflow"), NULL);
  mpfr_error = avm_join (avm_strung ("mpfr error"), NULL);
  bad_mpfr_spec = avm_join (avm_strung ("bad mpfr specification"), NULL);
  mpfr_overflow = avm_join (avm_strung ("mpfr overflow"), NULL);
  unrecognized_function_name = avm_join (avm_strung ("unrecognized mpfr function name"), NULL);
  unsigned_nan = avm_copied(shared_cell);
#if HAVE_MPFR
  negative_inf = avm_join(avm_copied(shared_cell),avm_natural(MPFR_INF_KIND));
  negative_zero = avm_join(avm_copied(shared_cell),avm_natural(MPFR_ZERO_KIND));
  negative_regular = avm_join(avm_copied(shared_cell),avm_natural(MPFR_REGULAR_KIND));
  positive_inf = avm_join(NULL,avm_natural(MPFR_INF_KIND));
  positive_zero = avm_join(NULL,avm_natural(MPFR_ZERO_KIND));
  positive_regular = avm_join(NULL,avm_natural(MPFR_REGULAR_KIND));
#endif
  inf_string = avm_strung("inf");
  ninf_string = avm_strung("-inf");
  nan_string = avm_strung("nan");
  string_number = 0;
  funs = back = NULL;
  while (funames[string_number])
    avm_enqueue (&funs, &back, avm_standard_strung (funames[string_number++]));
#if HAVE_GMP
#if HAVE_MPFR
  gmp_randinit_default(state);
  seed = time(0);
  gmp_randseed_ui(state,seed);
#endif
#endif
}





void
avm_count_mpfr ()

     /* This frees some static data structures as an aid to the
	detection of memory leaks. */
{

  if (!initialized)
    return;
  initialized = 0;
  avm_dispose (funs);
  avm_dispose (wild);
  avm_dispose (shared_cell);
  avm_dispose (mpfr_error);
  avm_dispose (bad_mpfr_spec);
  avm_dispose (mpfr_overflow);
  avm_dispose (memory_overflow);
  avm_dispose (unrecognized_function_name);
  avm_dispose (unsigned_nan);
#if HAVE_MPFR
  avm_dispose (negative_inf);
  avm_dispose (negative_zero);
  avm_dispose (negative_regular);
  avm_dispose (positive_inf);
  avm_dispose (positive_zero);
  avm_dispose (positive_regular);
#endif
  avm_dispose(inf_string);
  avm_dispose(ninf_string);
  avm_dispose(nan_string);
  inf_string = NULL;
  ninf_string = NULL;
  nan_string = NULL;
  unsigned_nan = NULL;
  negative_inf = NULL;
  negative_zero = NULL;
  negative_regular = NULL;
  positive_inf = NULL;
  positive_zero = NULL;
  positive_regular = NULL;
  funs = NULL;
  wild = NULL;
  shared_cell = NULL;
  mpfr_error = NULL;
  memory_overflow = NULL;
  unrecognized_function_name = NULL;
}
