From: Steven G. Kargl <kargls@comcast.net>
Date: Fri, 6 Aug 2004 20:36:05 +0000 (+0000)
Subject: arith.c: Add #define for model numbers.
X-Git-Url: https://git.libre-soc.org/?a=commitdiff_plain;h=f8e566e5253e4b0fe2dcd477fbc35ca5576cc7bc;p=gcc.git

arith.c: Add #define for model numbers.

2004-08-06  Steven G. Kargl  <kargls@comcast.net>

	* arith.c: Add #define for model numbers.  Remove global GMP variables.
	(natural_logarithm,common_logarithm,exponential,sine,
	cosine,arctangent,hypercos,hypersine ): Remove.
	(gfc_mpfr_to_mpz,gfc_set_model_kind,gfc_set_model): New functions.
	(arctangent2,gfc_arith_init_1,gfc_arith_done_1
	gfc_check_real_range, gfc_constant_result, gfc_range_check,
	gfc_arith_uminus,gfc_arith_plus, gfc_arith_minus, gfc_arith_times,
	gfc_arith_divide,complex_reciprocal,complex_pow_ui,
	gfc_arith_power,gfc_compare_expr,compare_complex,gfc_convert_real,
	gfc_convert_complex,gfc_int2real,gfc_int2complex,
	gfc_real2int,gfc_real2real,gfc_real2complex,
	gfc_complex2int,gfc_complex2real,gfc_complex2complex): Convert GMP
	to MPFR, use new functions.
	* arith.h: Remove extern global variables.
	(natural_logarithm,common_logarithm,exponential, sine, cosine,
	arctangent,hypercos,hypersine): Remove prototypes.
	(arctangent2): Update prototype from GMP to MPFR.
	(gfc_mpfr_to_mpz, gfc_set_model_kind,gfc_set_model): Add prototypes.
	* dump-parse-tree.c (gfc_show_expr): Convert GMP to MPFR.
	* expr.c (free_expr0,gfc_copy_expr): Convert GMP to MPFR.
	* gfortran.h (GFC_REAL_BITS): Remove.
	(arith): Add ARITH_NAN.
	Include mpfr.h.  Define GFC_RND_MODE.
	Rename GCC_GFORTRAN_H GFC_GFC_H.
	(gfc_expr): Convert GMP to MPFR.
	* module.c: Add arith.h, correct type in comment.
	(mio_gmp_real): Convert GMP to MPFR.
	(mio_expr):  Use gfc_set_model_kind().
	* primary.c:  Update copyright date with 2004.
	(match_real_constant,match_const_complex_part): Convert GMP to MPFR.
	* simplify.c: Remove global GMP variables
	(gfc_simplify_abs,gfc_simplify_acos,gfc_simplify_aimag,
	gfc_simplify_aint,gfc_simplify_dint,gfc_simplify_anint,
	gfc_simplify_dnint,gfc_simplify_asin,gfc_simplify_atan,
	gfc_simplify_atan2,gfc_simplify_ceiling,simplify_cmplx,
	gfc_simplify_conjg,gfc_simplify_cos,gfc_simplify_cosh,
	gfc_simplify_dim,gfc_simplify_dprod,gfc_simplify_epsilon,
	gfc_simplify_exp,gfc_simplify_exponent,gfc_simplify_floor,
	gfc_simplify_fraction,gfc_simplify_huge,gfc_simplify_int,
	gfc_simplify_ifix,gfc_simplify_idint,gfc_simplify_log,
	gfc_simplify_log10,simplify_min_max,gfc_simplify_mod,
	gfc_simplify_modulo,gfc_simplify_nearest,simplify_nint,
	gfc_simplify_rrspacing,gfc_simplify_scale,
	gfc_simplify_set_exponent,gfc_simplify_sign,gfc_simplify_sin,
	gfc_simplify_sinh,gfc_simplify_spacing,gfc_simplify_sqrt,
	gfc_simplify_tan,gfc_simplify_tanh,gfc_simplify_tiny,
	gfc_simplify_init_1,gfc_simplify_done_1):  Convert GMP to MPFR.
	Use new functions.
	* trans-const.c (gfc_conv_mpfr_to_tree): Rename from
	gfc_conv_mpf_to_tree.  Convert it to use MPFR
	(gfc_conv_constant_to_tree): Use it.
	* trans-const.h: Update prototype for gfc_conv_mpfr_to_tree().
	* trans-intrinsic.c: Add arith.h, remove gmp.h
	(gfc_conv_intrinsic_aint,gfc_conv_intrinsic_mod): Convert GMP to MPFR.

From-SVN: r85652
---

diff --git a/gcc/fortran/ChangeLog b/gcc/fortran/ChangeLog
index a3e1480a15f..4d0b3309ffc 100644
--- a/gcc/fortran/ChangeLog
+++ b/gcc/fortran/ChangeLog
@@ -1,3 +1,60 @@
+2004-08-06  Steven G. Kargl  <kargls@comcast.net>
+
+	* arith.c: Add #define for model numbers.  Remove global GMP variables.
+	(natural_logarithm,common_logarithm,exponential,sine,
+	cosine,arctangent,hypercos,hypersine ): Remove.
+	(gfc_mpfr_to_mpz,gfc_set_model_kind,gfc_set_model): New functions.
+	(arctangent2,gfc_arith_init_1,gfc_arith_done_1
+	gfc_check_real_range, gfc_constant_result, gfc_range_check,
+	gfc_arith_uminus,gfc_arith_plus, gfc_arith_minus, gfc_arith_times,
+	gfc_arith_divide,complex_reciprocal,complex_pow_ui,
+	gfc_arith_power,gfc_compare_expr,compare_complex,gfc_convert_real,
+	gfc_convert_complex,gfc_int2real,gfc_int2complex,
+	gfc_real2int,gfc_real2real,gfc_real2complex,
+	gfc_complex2int,gfc_complex2real,gfc_complex2complex): Convert GMP
+	to MPFR, use new functions.
+	* arith.h: Remove extern global variables.
+	(natural_logarithm,common_logarithm,exponential, sine, cosine,
+	arctangent,hypercos,hypersine): Remove prototypes.
+	(arctangent2): Update prototype from GMP to MPFR.
+	(gfc_mpfr_to_mpz, gfc_set_model_kind,gfc_set_model): Add prototypes.
+	* dump-parse-tree.c (gfc_show_expr): Convert GMP to MPFR.
+	* expr.c (free_expr0,gfc_copy_expr): Convert GMP to MPFR.
+	* gfortran.h (GFC_REAL_BITS): Remove.
+	(arith): Add ARITH_NAN.
+	Include mpfr.h.  Define GFC_RND_MODE.
+	Rename GCC_GFORTRAN_H GFC_GFC_H.
+	(gfc_expr): Convert GMP to MPFR.
+	* module.c: Add arith.h, correct type in comment.
+	(mio_gmp_real): Convert GMP to MPFR.
+	(mio_expr):  Use gfc_set_model_kind().
+	* primary.c:  Update copyright date with 2004.
+	(match_real_constant,match_const_complex_part): Convert GMP to MPFR.
+	* simplify.c: Remove global GMP variables
+	(gfc_simplify_abs,gfc_simplify_acos,gfc_simplify_aimag,
+	gfc_simplify_aint,gfc_simplify_dint,gfc_simplify_anint,
+	gfc_simplify_dnint,gfc_simplify_asin,gfc_simplify_atan,
+	gfc_simplify_atan2,gfc_simplify_ceiling,simplify_cmplx,
+	gfc_simplify_conjg,gfc_simplify_cos,gfc_simplify_cosh,
+	gfc_simplify_dim,gfc_simplify_dprod,gfc_simplify_epsilon,
+	gfc_simplify_exp,gfc_simplify_exponent,gfc_simplify_floor,
+	gfc_simplify_fraction,gfc_simplify_huge,gfc_simplify_int,
+	gfc_simplify_ifix,gfc_simplify_idint,gfc_simplify_log,
+	gfc_simplify_log10,simplify_min_max,gfc_simplify_mod,
+	gfc_simplify_modulo,gfc_simplify_nearest,simplify_nint,
+	gfc_simplify_rrspacing,gfc_simplify_scale,
+	gfc_simplify_set_exponent,gfc_simplify_sign,gfc_simplify_sin,
+	gfc_simplify_sinh,gfc_simplify_spacing,gfc_simplify_sqrt,
+	gfc_simplify_tan,gfc_simplify_tanh,gfc_simplify_tiny,
+	gfc_simplify_init_1,gfc_simplify_done_1):  Convert GMP to MPFR.
+	Use new functions.
+	* trans-const.c (gfc_conv_mpfr_to_tree): Rename from
+	gfc_conv_mpf_to_tree.  Convert it to use MPFR
+	(gfc_conv_constant_to_tree): Use it.
+	* trans-const.h: Update prototype for gfc_conv_mpfr_to_tree().
+	* trans-intrinsic.c: Add arith.h, remove gmp.h
+	(gfc_conv_intrinsic_aint,gfc_conv_intrinsic_mod): Convert GMP to MPFR.
+
 2004-08-06  Victor Leikehman  <lei@il.ibm.com>
 	Paul Brook  <paul@codesourcery.com>
 
diff --git a/gcc/fortran/arith.c b/gcc/fortran/arith.c
index b6aec5b951d..03ee14c0999 100644
--- a/gcc/fortran/arith.c
+++ b/gcc/fortran/arith.c
@@ -32,8 +32,6 @@ Software Foundation, 59 Temple Place - Suite 330, Boston, MA
 #include "gfortran.h"
 #include "arith.h"
 
-mpf_t pi, half_pi, two_pi, e;
-
 /* The gfc_(integer|real)_kinds[] structures have everything the front
    end needs to know about integers and real numbers on the target.
    Other entries of the structure are calculated from these values.
@@ -69,9 +67,31 @@ gfc_logical_info gfc_logical_kinds[] = {
   DEF_GFC_LOGICAL_KIND (0,  0)
 };
 
+
+/* IEEE-754 uses 1.xEe representation whereas the fortran standard
+   uses 0.xEe representation.  Hence the exponents below are biased
+   by one.  */
+
+#define GFC_SP_KIND      4
+#define GFC_SP_PREC     24   /* p    =   24, IEEE-754  */
+#define GFC_SP_EMIN   -125   /* emin = -126, IEEE-754  */
+#define GFC_SP_EMAX    128   /* emin =  127, IEEE-754  */
+
+/* Double precision model numbers.  */
+#define GFC_DP_KIND      8
+#define GFC_DP_PREC     53   /* p    =    53, IEEE-754  */
+#define GFC_DP_EMIN  -1021   /* emin = -1022, IEEE-754  */
+#define GFC_DP_EMAX   1024   /* emin =  1023, IEEE-754  */
+
+/* Quad precision model numbers.  Not used.  */
+#define GFC_QP_KIND     16
+#define GFC_QP_PREC    113   /* p    =    113, IEEE-754  */
+#define GFC_QP_EMIN -16381   /* emin = -16382, IEEE-754  */
+#define GFC_QP_EMAX  16384   /* emin =  16383, IEEE-754  */
+
 gfc_real_info gfc_real_kinds[] = {
-  DEF_GFC_REAL_KIND (4, 2, 24,  -125,  128),
-  DEF_GFC_REAL_KIND (8, 2, 53, -1021, 1024),
+  DEF_GFC_REAL_KIND (GFC_SP_KIND, 2, GFC_SP_PREC, GFC_SP_EMIN, GFC_SP_EMAX),
+  DEF_GFC_REAL_KIND (GFC_DP_KIND, 2, GFC_DP_PREC, GFC_DP_EMIN, GFC_DP_EMAX),
   DEF_GFC_REAL_KIND (0, 0,  0,     0,    0)
 };
 
@@ -82,440 +102,67 @@ gfc_real_info gfc_real_kinds[] = {
 int gfc_index_integer_kind;
 
 
-/* Compute the natural log of arg.
-
-   We first get the argument into the range 0.5 to 1.5 by successive
-   multiplications or divisions by e.  Then we use the series:
-
-     ln(x) = (x-1) - (x-1)^2/2 + (x-1)^3/3 - (x-1)^4/4 + ...
-
-   Because we are expanding in powers of (x-1), and 0.5 < x < 1.5, we
-   have -0.5 < (x-1) < 0.5.  Ignoring the harmonic term, this means
-   that each term is at most 1/(2^i), meaning one bit is gained per
-   iteration.
-
-   Not very efficient, but it doesn't have to be.  */
+/* MPFR does not have a direct replacement for mpz_set_f() from GMP.
+   It's easily implemented with a few calls though.  */
 
 void
-natural_logarithm (mpf_t * arg, mpf_t * result)
+gfc_mpfr_to_mpz(mpz_t z, mpfr_t x)
 {
-  mpf_t x, xp, t, log;
-  int i, p;
-
-  mpf_init_set (x, *arg);
-  mpf_init (t);
-
-  p = 0;
-
-  mpf_set_str (t, "0.5", 10);
-  while (mpf_cmp (x, t) < 0)
-    {
-      mpf_mul (x, x, e);
-      p--;
-    }
-
-  mpf_set_str (t, "1.5", 10);
-  while (mpf_cmp (x, t) > 0)
-    {
-      mpf_div (x, x, e);
-      p++;
-    }
-
-  mpf_sub_ui (x, x, 1);
-  mpf_init_set_ui (log, 0);
-  mpf_init_set_ui (xp, 1);
-
-  for (i = 1; i < GFC_REAL_BITS; i++)
-    {
-      mpf_mul (xp, xp, x);
-      mpf_div_ui (t, xp, i);
+  mp_exp_t e;
 
-      if (i % 2 == 0)
-	mpf_sub (log, log, t);
-      else
-	mpf_add (log, log, t);
-    }
-
-  /* Add in the log (e^p) = p */
-
-  if (p < 0)
-    mpf_sub_ui (log, log, -p);
+  e = mpfr_get_z_exp (z, x);
+  if (e > 0)
+    mpz_mul_2exp (z, z, e);
   else
-    mpf_add_ui (log, log, p);
-
-  mpf_clear (x);
-  mpf_clear (xp);
-  mpf_clear (t);
-
-  mpf_set (*result, log);
-  mpf_clear (log);
-}
-
-
-/* Calculate the common logarithm of arg.  We use the natural
-   logarithm of arg and of 10:
-
-   log10(arg) = log(arg)/log(10)  */
-
-void
-common_logarithm (mpf_t * arg, mpf_t * result)
-{
-  mpf_t i10, log10;
-
-  natural_logarithm (arg, result);
-
-  mpf_init_set_ui (i10, 10);
-  mpf_init (log10);
-  natural_logarithm (&i10, &log10);
-
-  mpf_div (*result, *result, log10);
-  mpf_clear (i10);
-  mpf_clear (log10);
+    mpz_tdiv_q_2exp (z, z, -e);
+  if (mpfr_sgn (x) < 0)
+    mpz_neg (z, z);
 }
 
-/* Calculate exp(arg).
-
-   We use a reduction of the form
-
-     x = Nln2 + r
 
-   Then we obtain exp(r) from the Maclaurin series.
-   exp(x) is then recovered from the identity
-
-     exp(x) = 2^N*exp(r).  */
+/* Set the model number precision by the requested KIND.  */
 
 void
-exponential (mpf_t * arg, mpf_t * result)
+gfc_set_model_kind (int kind)
 {
-  mpf_t two, ln2, power, q, r, num, denom, term, x, xp;
-  int i;
-  long n;
-  unsigned long p, mp;
-
-
-  mpf_init_set (x, *arg);
-
-  if (mpf_cmp_ui (x, 0) == 0)
-    {
-      mpf_set_ui (*result, 1);
-    }
-  else if (mpf_cmp_ui (x, 1) == 0)
-    {
-      mpf_set (*result, e);
-    }
-  else
-    {
-      mpf_init_set_ui (two, 2);
-      mpf_init (ln2);
-      mpf_init (q);
-      mpf_init (r);
-      mpf_init (power);
-      mpf_init (term);
-
-      natural_logarithm (&two, &ln2);
-
-      mpf_div (q, x, ln2);
-      mpf_floor (power, q);
-      mpf_mul (q, power, ln2);
-      mpf_sub (r, x, q);
-
-      mpf_init_set_ui (xp, 1);
-      mpf_init_set_ui (num, 1);
-      mpf_init_set_ui (denom, 1);
-
-      for (i = 1; i <= GFC_REAL_BITS + 10; i++)
+  switch (kind)
 	{
-	  mpf_mul (num, num, r);
-	  mpf_mul_ui (denom, denom, i);
-	  mpf_div (term, num, denom);
-	  mpf_add (xp, xp, term);
-	}
-
-      /* Reconstruction step */
-      n = (long) mpf_get_d (power);
-
-      if (n > 0)
-	{
-	  p = (unsigned int) n;
-	  mpf_mul_2exp (*result, xp, p);
-	}
-      else
-	{
-	  mp = (unsigned int) (-n);
-	  mpf_div_2exp (*result, xp, mp);
-	}
-
-      mpf_clear (two);
-      mpf_clear (ln2);
-      mpf_clear (q);
-      mpf_clear (r);
-      mpf_clear (power);
-      mpf_clear (num);
-      mpf_clear (denom);
-      mpf_clear (term);
-      mpf_clear (xp);
-    }
-
-  mpf_clear (x);
-}
-
-
-/* Calculate sin(arg).
-
-   We use a reduction of the form
-
-     x= N*2pi + r
-
-   Then we obtain sin(r) from the Maclaurin series.  */
-
-void
-sine (mpf_t * arg, mpf_t * result)
-{
-  mpf_t factor, q, r, num, denom, term, x, xp;
-  int i, sign;
-
-  mpf_init_set (x, *arg);
-
-  /* Special case (we do not treat multiples of pi due to roundoff issues) */
-  if (mpf_cmp_ui (x, 0) == 0)
-    {
-      mpf_set_ui (*result, 0);
-    }
-  else
-    {
-      mpf_init (q);
-      mpf_init (r);
-      mpf_init (factor);
-      mpf_init (term);
-
-      mpf_div (q, x, two_pi);
-      mpf_floor (factor, q);
-      mpf_mul (q, factor, two_pi);
-      mpf_sub (r, x, q);
-
-      mpf_init_set_ui (xp, 0);
-      mpf_init_set_ui (num, 1);
-      mpf_init_set_ui (denom, 1);
-
-      sign = -1;
-      for (i = 1; i < GFC_REAL_BITS + 10; i++)
-	{
-	  mpf_mul (num, num, r);
-	  mpf_mul_ui (denom, denom, i);
-	  if (i % 2 == 0)
-	    continue;
-
-	  sign = -sign;
-	  mpf_div (term, num, denom);
-	  if (sign > 0)
-	    mpf_add (xp, xp, term);
-	  else
-	    mpf_sub (xp, xp, term);
-	}
-
-      mpf_set (*result, xp);
-
-      mpf_clear (q);
-      mpf_clear (r);
-      mpf_clear (factor);
-      mpf_clear (num);
-      mpf_clear (denom);
-      mpf_clear (term);
-      mpf_clear (xp);
-    }
-
-  mpf_clear (x);
-}
-
-
-/* Calculate cos(arg).
-
-   Similar to sine.  */
-
-void
-cosine (mpf_t * arg, mpf_t * result)
-{
-  mpf_t factor, q, r, num, denom, term, x, xp;
-  int i, sign;
-
-  mpf_init_set (x, *arg);
-
-  /* Special case (we do not treat multiples of pi due to roundoff issues) */
-  if (mpf_cmp_ui (x, 0) == 0)
-    {
-      mpf_set_ui (*result, 1);
-    }
-  else
-    {
-      mpf_init (q);
-      mpf_init (r);
-      mpf_init (factor);
-      mpf_init (term);
-
-      mpf_div (q, x, two_pi);
-      mpf_floor (factor, q);
-      mpf_mul (q, factor, two_pi);
-      mpf_sub (r, x, q);
-
-      mpf_init_set_ui (xp, 1);
-      mpf_init_set_ui (num, 1);
-      mpf_init_set_ui (denom, 1);
-
-      sign = 1;
-      for (i = 1; i < GFC_REAL_BITS + 10; i++)
-	{
-	  mpf_mul (num, num, r);
-	  mpf_mul_ui (denom, denom, i);
-	  if (i % 2 != 0)
-	    continue;
-
-	  sign = -sign;
-	  mpf_div (term, num, denom);
-	  if (sign > 0)
-	    mpf_add (xp, xp, term);
-	  else
-	    mpf_sub (xp, xp, term);
-	}
-      mpf_set (*result, xp);
-
-      mpf_clear (q);
-      mpf_clear (r);
-      mpf_clear (factor);
-      mpf_clear (num);
-      mpf_clear (denom);
-      mpf_clear (term);
-      mpf_clear (xp);
+    case GFC_SP_KIND:
+      mpfr_set_default_prec (GFC_SP_PREC);
+      break;
+    case GFC_DP_KIND:
+      mpfr_set_default_prec (GFC_DP_PREC);
+      break;
+    case GFC_QP_KIND:
+      mpfr_set_default_prec (GFC_QP_PREC);
+      break;
+    default:
+      gfc_internal_error ("gfc_set_model_kind(): Bad model number");
     }
-
-  mpf_clear (x);
 }
 
 
-/* Calculate atan(arg).
-
-   Similar to sine but requires special handling for x near 1.  */
+/* Set the model number precision from mpfr_t x.  */
 
 void
-arctangent (mpf_t * arg, mpf_t * result)
+gfc_set_model (mpfr_t x)
 {
-  mpf_t absval, convgu, convgl, num, term, x, xp;
-  int i, sign;
-
-  mpf_init_set (x, *arg);
-
-  /* Special cases */
-  if (mpf_cmp_ui (x, 0) == 0)
+  switch (mpfr_get_prec (x))
     {
-      mpf_set_ui (*result, 0);
-    }
-  else if (mpf_cmp_ui (x, 1) == 0)
-    {
-      mpf_init (num);
-      mpf_div_ui (num, half_pi, 2);
-      mpf_set (*result, num);
-      mpf_clear (num);
-    }
-  else if (mpf_cmp_si (x, -1) == 0)
-    {
-      mpf_init (num);
-      mpf_div_ui (num, half_pi, 2);
-      mpf_neg (*result, num);
-      mpf_clear (num);
-    }
-  else
-    {				/* General cases */
-
-      mpf_init (absval);
-      mpf_abs (absval, x);
-
-      mpf_init_set_d (convgu, 1.5);
-      mpf_init_set_d (convgl, 0.5);
-      mpf_init_set_ui (num, 1);
-      mpf_init (term);
-
-      if (mpf_cmp (absval, convgl) < 0)
-	{
-	  mpf_init_set_ui (xp, 0);
-	  sign = -1;
-	  for (i = 1; i < GFC_REAL_BITS + 10; i++)
-	    {
-	      mpf_mul (num, num, absval);
-	      if (i % 2 == 0)
-		continue;
-
-	      sign = -sign;
-	      mpf_div_ui (term, num, i);
-	      if (sign > 0)
-		mpf_add (xp, xp, term);
-	      else
-		mpf_sub (xp, xp, term);
-	    }
-	}
-      else if (mpf_cmp (absval, convgu) >= 0)
-	{
-	  mpf_init_set (xp, half_pi);
-	  sign = 1;
-	  for (i = 1; i < GFC_REAL_BITS + 10; i++)
-	    {
-	      mpf_div (num, num, absval);
-	      if (i % 2 == 0)
-		continue;
-
-	      sign = -sign;
-	      mpf_div_ui (term, num, i);
-	      if (sign > 0)
-		mpf_add (xp, xp, term);
-	      else
-		mpf_sub (xp, xp, term);
-	    }
-	}
-      else
-	{
-	  mpf_init_set_ui (xp, 0);
-
-	  mpf_sub_ui (num, absval, 1);
-	  mpf_add_ui (term, absval, 1);
-	  mpf_div (absval, num, term);
-
-	  mpf_set_ui (num, 1);
-
-	  sign = -1;
-	  for (i = 1; i < GFC_REAL_BITS + 10; i++)
-	    {
-	      mpf_mul (num, num, absval);
-	      if (i % 2 == 0)
-		continue;
-	      sign = -sign;
-	      mpf_div_ui (term, num, i);
-	      if (sign > 0)
-		mpf_add (xp, xp, term);
-	      else
-		mpf_sub (xp, xp, term);
-	    }
-
-	  mpf_div_ui (term, half_pi, 2);
-	  mpf_add (xp, term, xp);
-	}
-
-      /* This makes sure to preserve the identity arctan(-x) = -arctan(x)
-         and improves accuracy to boot.  */
-
-      if (mpf_cmp_ui (x, 0) > 0)
-	mpf_set (*result, xp);
-      else
-	mpf_neg (*result, xp);
-
-      mpf_clear (absval);
-      mpf_clear (convgl);
-      mpf_clear (convgu);
-      mpf_clear (num);
-      mpf_clear (term);
-      mpf_clear (xp);
+    case GFC_SP_PREC:
+      mpfr_set_default_prec (GFC_SP_PREC);
+      break;
+    case GFC_DP_PREC:
+      mpfr_set_default_prec (GFC_DP_PREC);
+      break;
+    case GFC_QP_PREC:
+      mpfr_set_default_prec (GFC_QP_PREC);
+      break;
+    default:
+      gfc_internal_error ("gfc_set_model(): Bad model number");
     }
-  mpf_clear (x);
 }
 
-
 /* Calculate atan2 (y, x)
 
 atan2(y, x) = atan(y/x)				if x > 0,
@@ -525,97 +172,46 @@ atan2(y, x) = atan(y/x)				if x > 0,
 */
 
 void
-arctangent2 (mpf_t * y, mpf_t * x, mpf_t * result)
+arctangent2 (mpfr_t y, mpfr_t x, mpfr_t result)
 {
-  mpf_t t;
+  int i;
+  mpfr_t t;
 
-  mpf_init (t);
+  gfc_set_model (y);
+  mpfr_init (t);
 
-  switch (mpf_sgn (*x))
+  i = mpfr_sgn(x);
+
+  if (i > 0)
     {
-    case 1:
-      mpf_div (t, *y, *x);
-      arctangent (&t, result);
-      break;
-    case -1:
-      mpf_div (t, *y, *x);
-      mpf_abs (t, t);
-      arctangent (&t, &t);
-      mpf_sub (*result, pi, t);
-      if (mpf_sgn (*y) == -1)
-	mpf_neg (*result, *result);
-      break;
-    case 0:
-      if (mpf_sgn (*y) == 0)
-	mpf_set_ui (*result, 0);
+      mpfr_div (t, y, x, GFC_RND_MODE);
+      mpfr_atan (result, t, GFC_RND_MODE);
+    }
+  else if (i < 0)
+    {
+      mpfr_const_pi (result, GFC_RND_MODE);
+      mpfr_div (t, y, x, GFC_RND_MODE);
+      mpfr_abs (t, t, GFC_RND_MODE);
+      mpfr_atan (t, t, GFC_RND_MODE);
+      mpfr_sub (result, result, t, GFC_RND_MODE);
+      if (mpfr_sgn (y) < 0)
+	mpfr_neg (result, result, GFC_RND_MODE);
+    }
+      else
+	{
+      if (mpfr_sgn (y) == 0)
+	mpfr_set_ui (result, 0, GFC_RND_MODE);
       else
 	{
-	  mpf_set (*result, half_pi);
-	  if (mpf_sgn (*y) == -1)
-	    mpf_neg (*result, *result);
+          mpfr_const_pi (result, GFC_RND_MODE);
+          mpfr_div_ui (result, result, 2, GFC_RND_MODE);
+	  if (mpfr_sgn (y) < 0)
+	    mpfr_neg (result, result, GFC_RND_MODE);
 	}
-       break;
     }
-  mpf_clear (t);
-}
-
-/* Calculate cosh(arg). */
-
-void
-hypercos (mpf_t * arg, mpf_t * result)
-{
-  mpf_t neg, term1, term2, x, xp;
-
-  mpf_init_set (x, *arg);
 
-  mpf_init (neg);
-  mpf_init (term1);
-  mpf_init (term2);
-  mpf_init (xp);
+  mpfr_clear (t);
 
-  mpf_neg (neg, x);
-
-  exponential (&x, &term1);
-  exponential (&neg, &term2);
-
-  mpf_add (xp, term1, term2);
-  mpf_div_ui (*result, xp, 2);
-
-  mpf_clear (neg);
-  mpf_clear (term1);
-  mpf_clear (term2);
-  mpf_clear (x);
-  mpf_clear (xp);
-}
-
-
-/* Calculate sinh(arg). */
-
-void
-hypersine (mpf_t * arg, mpf_t * result)
-{
-  mpf_t neg, term1, term2, x, xp;
-
-  mpf_init_set (x, *arg);
-
-  mpf_init (neg);
-  mpf_init (term1);
-  mpf_init (term2);
-  mpf_init (xp);
-
-  mpf_neg (neg, x);
-
-  exponential (&x, &term1);
-  exponential (&neg, &term2);
-
-  mpf_sub (xp, term1, term2);
-  mpf_div_ui (*result, xp, 2);
-
-  mpf_clear (neg);
-  mpf_clear (term1);
-  mpf_clear (term2);
-  mpf_clear (x);
-  mpf_clear (xp);
 }
 
 
@@ -638,6 +234,9 @@ gfc_arith_error (arith code)
     case ARITH_UNDERFLOW:
       p = "Arithmetic underflow";
       break;
+    case ARITH_NAN:
+      p = "Arithmetic NaN";
+      break;
     case ARITH_DIV0:
       p = "Division by zero";
       break;
@@ -662,72 +261,17 @@ gfc_arith_init_1 (void)
 {
   gfc_integer_info *int_info;
   gfc_real_info *real_info;
-  mpf_t a, b;
+  mpfr_t a, b, c;
   mpz_t r;
-  int i, n, limit;
-
-  /* Set the default precision for GMP computations.  */
-  mpf_set_default_prec (GFC_REAL_BITS + 30);
-
-  /* Calculate e, needed by the natural_logarithm() subroutine.  */
-  mpf_init (b);
-  mpf_init_set_ui (e, 0);
-  mpf_init_set_ui (a, 1);
-
-  for (i = 1; i < 100; i++)
-    {
-      mpf_add (e, e, a);
-      mpf_div_ui (a, a, i);	/* 1/(i!) */
-    }
-
-  /* Calculate pi, 2pi, pi/2, and -pi/2, needed for trigonometric
-     functions.
-
-     We use the Bailey, Borwein and Plouffe formula:
-
-       pi = \sum{n=0}^\infty (1/16)^n [4/(8n+1) - 2/(8n+4) - 1/(8n+5) - 1/(8n+6)]
-
-     which gives about four bits per iteration.  */
-
-  mpf_init_set_ui (pi, 0);
-
-  mpf_init (two_pi);
-  mpf_init (half_pi);
-
-  limit = (GFC_REAL_BITS / 4) + 10;	/* (1/16)^n gives 4 bits per iteration */
-
-  for (n = 0; n < limit; n++)
-    {
-      mpf_set_ui (b, 4);
-      mpf_div_ui (b, b, 8 * n + 1);	/* 4/(8n+1) */
-
-      mpf_set_ui (a, 2);
-      mpf_div_ui (a, a, 8 * n + 4);	/* 2/(8n+4) */
-      mpf_sub (b, b, a);
-
-      mpf_set_ui (a, 1);
-      mpf_div_ui (a, a, 8 * n + 5);	/* 1/(8n+5) */
-      mpf_sub (b, b, a);
-
-      mpf_set_ui (a, 1);
-      mpf_div_ui (a, a, 8 * n + 6);	/* 1/(8n+6) */
-      mpf_sub (b, b, a);
-
-      mpf_set_ui (a, 16);
-      mpf_pow_ui (a, a, n);	/* 16^n */
-
-      mpf_div (b, b, a);
+  int i;
 
-      mpf_add (pi, pi, b);
-    }
+  gfc_set_model_kind (GFC_QP_KIND);
 
-  mpf_mul_ui (two_pi, pi, 2);
-  mpf_div_ui (half_pi, pi, 2);
+  mpfr_init (a);
+  mpz_init (r);
 
   /* Convert the minimum/maximum values for each kind into their
      GNU MP representation.  */
-  mpz_init (r);
-
   for (int_info = gfc_integer_kinds; int_info->kind != 0; int_info++)
     {
       /* Huge */
@@ -751,59 +295,76 @@ gfc_arith_init_1 (void)
       mpz_add_ui (int_info->max_int, int_info->max_int, 1);
 
       /* Range */
-      mpf_set_z (a, int_info->huge);
-      common_logarithm (&a, &a);
-      mpf_trunc (a, a);
-      mpz_set_f (r, a);
+      mpfr_set_z (a, int_info->huge, GFC_RND_MODE);
+      mpfr_log10 (a, a, GFC_RND_MODE);
+      mpfr_trunc (a, a);
+      gfc_mpfr_to_mpz (r, a);
       int_info->range = mpz_get_si (r);
     }
 
-  /*  mpf_set_default_prec(GFC_REAL_BITS); */
+  mpfr_clear (a);
+
   for (real_info = gfc_real_kinds; real_info->kind != 0; real_info++)
     {
-      /* Huge */
-      mpf_set_ui (a, real_info->radix);
-      mpf_set_ui (b, real_info->radix);
+      gfc_set_model_kind (real_info->kind);
 
-      mpf_pow_ui (a, a, real_info->max_exponent);
-      mpf_pow_ui (b, b, real_info->max_exponent - real_info->digits);
+      mpfr_init (a);
+      mpfr_init (b);
+      mpfr_init (c);
 
-      mpf_init (real_info->huge);
-      mpf_sub (real_info->huge, a, b);
+      /* huge(x) = (1 - b**(-p)) * b**(emax-1) * b  */
+      /* a = 1 - b**(-p) */
+      mpfr_set_ui (a, 1, GFC_RND_MODE);
+      mpfr_set_ui (b, real_info->radix, GFC_RND_MODE);
+      mpfr_pow_si (b, b, -real_info->digits, GFC_RND_MODE);
+      mpfr_sub (a, a, b, GFC_RND_MODE);
 
-      /* Tiny */
-      mpf_set_ui (b, real_info->radix);
-      mpf_pow_ui (b, b, 1 - real_info->min_exponent);
+      /* c = b**(emax-1) */
+      mpfr_set_ui (b, real_info->radix, GFC_RND_MODE);
+      mpfr_pow_ui (c, b, real_info->max_exponent - 1, GFC_RND_MODE);
 
-      mpf_init (real_info->tiny);
-      mpf_ui_div (real_info->tiny, 1, b);
+      /* a = a * c = (1 - b**(-p)) * b**(emax-1) */
+      mpfr_mul (a, a, c, GFC_RND_MODE);
 
-      /* Epsilon */
-      mpf_set_ui (b, real_info->radix);
-      mpf_pow_ui (b, b, real_info->digits - 1);
+      /* a = (1 - b**(-p)) * b**(emax-1) * b */
+      mpfr_mul_ui (a, a, real_info->radix, GFC_RND_MODE);
 
-      mpf_init (real_info->epsilon);
-      mpf_ui_div (real_info->epsilon, 1, b);
+      mpfr_init (real_info->huge);
+      mpfr_set (real_info->huge, a, GFC_RND_MODE);
 
-      /* Range */
-      common_logarithm (&real_info->huge, &a);
-      common_logarithm (&real_info->tiny, &b);
-      mpf_neg (b, b);
+      /* tiny(x) = b**(emin-1) */
+      mpfr_set_ui (b, real_info->radix, GFC_RND_MODE);
+      mpfr_pow_si (b, b, real_info->min_exponent - 1, GFC_RND_MODE);
+
+      mpfr_init (real_info->tiny);
+      mpfr_set (real_info->tiny, b, GFC_RND_MODE);
+
+      /* epsilon(x) = b**(1-p) */
+      mpfr_set_ui (b, real_info->radix, GFC_RND_MODE);
+      mpfr_pow_si (b, b, 1 - real_info->digits, GFC_RND_MODE);
+
+      mpfr_init (real_info->epsilon);
+      mpfr_set (real_info->epsilon, b, GFC_RND_MODE);
 
-      if (mpf_cmp (a, b) > 0)
-	mpf_set (a, b);		/* a = min(a, b) */
+      /* range(x) = int(min(log10(huge(x)), -log10(tiny)) */
+      mpfr_log10 (a, real_info->huge, GFC_RND_MODE);
+      mpfr_log10 (b, real_info->tiny, GFC_RND_MODE);
+      mpfr_neg (b, b, GFC_RND_MODE);
 
-      mpf_trunc (a, a);
-      mpz_set_f (r, a);
+      if (mpfr_cmp (a, b) > 0)
+	mpfr_set (a, b, GFC_RND_MODE);		/* a = min(a, b) */
+
+      mpfr_trunc (a, a);
+      gfc_mpfr_to_mpz (r, a);
       real_info->range = mpz_get_si (r);
 
-      /* Precision */
-      mpf_set_ui (a, real_info->radix);
-      common_logarithm (&a, &a);
+      /* precision(x) = int((p - 1) * log10(b)) + k */
+      mpfr_set_ui (a, real_info->radix, GFC_RND_MODE);
+      mpfr_log10 (a, a, GFC_RND_MODE);
 
-      mpf_mul_ui (a, a, real_info->digits - 1);
-      mpf_trunc (a, a);
-      mpz_set_f (r, a);
+      mpfr_mul_ui (a, a, real_info->digits - 1, GFC_RND_MODE);
+      mpfr_trunc (a, a);
+      gfc_mpfr_to_mpz (r, a);
       real_info->precision = mpz_get_si (r);
 
       /* If the radix is an integral power of 10, add one to the
@@ -811,11 +372,13 @@ gfc_arith_init_1 (void)
       for (i = 10; i <= real_info->radix; i *= 10)
 	if (i == real_info->radix)
 	  real_info->precision++;
+
+      mpfr_clear (a);
+      mpfr_clear (b);
+      mpfr_clear (c);
     }
 
   mpz_clear (r);
-  mpf_clear (a);
-  mpf_clear (b);
 }
 
 
@@ -827,12 +390,6 @@ gfc_arith_done_1 (void)
   gfc_integer_info *ip;
   gfc_real_info *rp;
 
-  mpf_clear (e);
-
-  mpf_clear (pi);
-  mpf_clear (half_pi);
-  mpf_clear (two_pi);
-
   for (ip = gfc_integer_kinds; ip->kind; ip++)
     {
       mpz_clear (ip->min_int);
@@ -842,9 +399,9 @@ gfc_arith_done_1 (void)
 
   for (rp = gfc_real_kinds; rp->kind; rp++)
     {
-      mpf_clear (rp->epsilon);
-      mpf_clear (rp->huge);
-      mpf_clear (rp->tiny);
+      mpfr_clear (rp->epsilon);
+      mpfr_clear (rp->huge);
+      mpfr_clear (rp->tiny);
     }
 }
 
@@ -1022,34 +579,35 @@ gfc_check_integer_range (mpz_t p, int kind)
    ARITH_UNDERFLOW.  */
 
 static arith
-gfc_check_real_range (mpf_t p, int kind)
+gfc_check_real_range (mpfr_t p, int kind)
 {
   arith retval;
-  mpf_t q;
+  mpfr_t q;
   int i;
 
-  mpf_init (q);
-  mpf_abs (q, p);
-
   i = validate_real (kind);
   if (i == -1)
     gfc_internal_error ("gfc_check_real_range(): Bad kind");
 
+  gfc_set_model (p);
+  mpfr_init (q);
+  mpfr_abs (q, p, GFC_RND_MODE);
+
   retval = ARITH_OK;
-  if (mpf_sgn (q) == 0)
+  if (mpfr_sgn (q) == 0)
     goto done;
 
-  if (mpf_cmp (q, gfc_real_kinds[i].huge) == 1)
+  if (mpfr_cmp (q, gfc_real_kinds[i].huge) > 0)
     {
       retval = ARITH_OVERFLOW;
       goto done;
     }
 
-  if (mpf_cmp (q, gfc_real_kinds[i].tiny) == -1)
+  if (mpfr_cmp (q, gfc_real_kinds[i].tiny) < 0)
     retval = ARITH_UNDERFLOW;
 
 done:
-  mpf_clear (q);
+  mpfr_clear (q);
 
   return retval;
 }
@@ -1081,12 +639,14 @@ gfc_constant_result (bt type, int kind, locus * where)
       break;
 
     case BT_REAL:
-      mpf_init (result->value.real);
+      gfc_set_model_kind (kind);
+      mpfr_init (result->value.real);
       break;
 
     case BT_COMPLEX:
-      mpf_init (result->value.complex.r);
-      mpf_init (result->value.complex.i);
+      gfc_set_model_kind (kind);
+      mpfr_init (result->value.complex.r);
+      mpfr_init (result->value.complex.i);
       break;
 
     default:
@@ -1173,9 +733,7 @@ gfc_arith_neqv (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
 
 
 /* Make sure a constant numeric expression is within the range for
-   its type and kind.  GMP is doing 130 bit arithmetic, so an UNDERFLOW
-   is numerically zero for REAL(4) and REAL(8) types.  Reset the value(s)
-   to exactly 0 for UNDERFLOW.  Note that there's also a gfc_check_range(),
+   its type and kind.  Note that there's also a gfc_check_range(),
    but that one deals with the intrinsic RANGE function.  */
 
 arith
@@ -1192,18 +750,18 @@ gfc_range_check (gfc_expr * e)
     case BT_REAL:
       rc = gfc_check_real_range (e->value.real, e->ts.kind);
       if (rc == ARITH_UNDERFLOW)
-        mpf_set_ui (e->value.real, 0);
+        mpfr_set_ui (e->value.real, 0, GFC_RND_MODE);
       break;
 
     case BT_COMPLEX:
       rc = gfc_check_real_range (e->value.complex.r, e->ts.kind);
       if (rc == ARITH_UNDERFLOW)
-        mpf_set_ui (e->value.complex.r, 0);
+        mpfr_set_ui (e->value.complex.r, 0, GFC_RND_MODE);
       if (rc == ARITH_OK || rc == ARITH_UNDERFLOW)
         {
           rc = gfc_check_real_range (e->value.complex.i, e->ts.kind);
           if (rc == ARITH_UNDERFLOW)
-            mpf_set_ui (e->value.complex.i, 0);
+            mpfr_set_ui (e->value.complex.i, 0, GFC_RND_MODE);
         }
 
       break;
@@ -1244,12 +802,12 @@ gfc_arith_uminus (gfc_expr * op1, gfc_expr ** resultp)
       break;
 
     case BT_REAL:
-      mpf_neg (result->value.real, op1->value.real);
+      mpfr_neg (result->value.real, op1->value.real, GFC_RND_MODE);
       break;
 
     case BT_COMPLEX:
-      mpf_neg (result->value.complex.r, op1->value.complex.r);
-      mpf_neg (result->value.complex.i, op1->value.complex.i);
+      mpfr_neg (result->value.complex.r, op1->value.complex.r, GFC_RND_MODE);
+      mpfr_neg (result->value.complex.i, op1->value.complex.i, GFC_RND_MODE);
       break;
 
     default:
@@ -1289,15 +847,16 @@ gfc_arith_plus (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
       break;
 
     case BT_REAL:
-      mpf_add (result->value.real, op1->value.real, op2->value.real);
+      mpfr_add (result->value.real, op1->value.real, op2->value.real,
+               GFC_RND_MODE);
       break;
 
     case BT_COMPLEX:
-      mpf_add (result->value.complex.r, op1->value.complex.r,
-	       op2->value.complex.r);
+      mpfr_add (result->value.complex.r, op1->value.complex.r,
+	       op2->value.complex.r, GFC_RND_MODE);
 
-      mpf_add (result->value.complex.i, op1->value.complex.i,
-	       op2->value.complex.i);
+      mpfr_add (result->value.complex.i, op1->value.complex.i,
+	       op2->value.complex.i, GFC_RND_MODE);
       break;
 
     default:
@@ -1337,16 +896,16 @@ gfc_arith_minus (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
       break;
 
     case BT_REAL:
-      mpf_sub (result->value.real, op1->value.real, op2->value.real);
+      mpfr_sub (result->value.real, op1->value.real, op2->value.real,
+                GFC_RND_MODE);
       break;
 
     case BT_COMPLEX:
-      mpf_sub (result->value.complex.r, op1->value.complex.r,
-	       op2->value.complex.r);
-
-      mpf_sub (result->value.complex.i, op1->value.complex.i,
-	       op2->value.complex.i);
+      mpfr_sub (result->value.complex.r, op1->value.complex.r,
+	       op2->value.complex.r, GFC_RND_MODE);
 
+      mpfr_sub (result->value.complex.i, op1->value.complex.i,
+	       op2->value.complex.i, GFC_RND_MODE);
       break;
 
     default:
@@ -1375,7 +934,7 @@ static arith
 gfc_arith_times (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
 {
   gfc_expr *result;
-  mpf_t x, y;
+  mpfr_t x, y;
   arith rc;
 
   result = gfc_constant_result (op1->ts.type, op1->ts.kind, &op1->where);
@@ -1387,23 +946,28 @@ gfc_arith_times (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
       break;
 
     case BT_REAL:
-      mpf_mul (result->value.real, op1->value.real, op2->value.real);
+      mpfr_mul (result->value.real, op1->value.real, op2->value.real,
+               GFC_RND_MODE);
       break;
 
     case BT_COMPLEX:
-      mpf_init (x);
-      mpf_init (y);
 
-      mpf_mul (x, op1->value.complex.r, op2->value.complex.r);
-      mpf_mul (y, op1->value.complex.i, op2->value.complex.i);
-      mpf_sub (result->value.complex.r, x, y);
+      /* FIXME:  possible numericals problem.  */
 
-      mpf_mul (x, op1->value.complex.r, op2->value.complex.i);
-      mpf_mul (y, op1->value.complex.i, op2->value.complex.r);
-      mpf_add (result->value.complex.i, x, y);
+      gfc_set_model (op1->value.complex.r);
+      mpfr_init (x);
+      mpfr_init (y);
 
-      mpf_clear (x);
-      mpf_clear (y);
+      mpfr_mul (x, op1->value.complex.r, op2->value.complex.r, GFC_RND_MODE);
+      mpfr_mul (y, op1->value.complex.i, op2->value.complex.i, GFC_RND_MODE);
+      mpfr_sub (result->value.complex.r, x, y, GFC_RND_MODE);
+
+      mpfr_mul (x, op1->value.complex.r, op2->value.complex.i, GFC_RND_MODE);
+      mpfr_mul (y, op1->value.complex.i, op2->value.complex.r, GFC_RND_MODE);
+      mpfr_add (result->value.complex.i, x, y, GFC_RND_MODE);
+
+      mpfr_clear (x);
+      mpfr_clear (y);
 
       break;
 
@@ -1433,7 +997,7 @@ static arith
 gfc_arith_divide (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
 {
   gfc_expr *result;
-  mpf_t x, y, div;
+  mpfr_t x, y, div;
   arith rc;
 
   rc = ARITH_OK;
@@ -1454,44 +1018,51 @@ gfc_arith_divide (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
       break;
 
     case BT_REAL:
-      if (mpf_sgn (op2->value.real) == 0)
+      /* FIXME: MPFR correctly generates NaN.  This may not be needed.  */
+      if (mpfr_sgn (op2->value.real) == 0)
 	{
 	  rc = ARITH_DIV0;
 	  break;
 	}
 
-      mpf_div (result->value.real, op1->value.real, op2->value.real);
+      mpfr_div (result->value.real, op1->value.real, op2->value.real,
+               GFC_RND_MODE);
       break;
 
     case BT_COMPLEX:
-      if (mpf_sgn (op2->value.complex.r) == 0
-	  && mpf_sgn (op2->value.complex.i) == 0)
+      /* FIXME: MPFR correctly generates NaN.  This may not be needed.  */
+      if (mpfr_sgn (op2->value.complex.r) == 0
+	  && mpfr_sgn (op2->value.complex.i) == 0)
 	{
 	  rc = ARITH_DIV0;
 	  break;
 	}
 
-      mpf_init (x);
-      mpf_init (y);
-      mpf_init (div);
+      gfc_set_model (op1->value.complex.r);
+      mpfr_init (x);
+      mpfr_init (y);
+      mpfr_init (div);
 
-      mpf_mul (x, op2->value.complex.r, op2->value.complex.r);
-      mpf_mul (y, op2->value.complex.i, op2->value.complex.i);
-      mpf_add (div, x, y);
+      /* FIXME: possible numerical problems.  */
+      mpfr_mul (x, op2->value.complex.r, op2->value.complex.r, GFC_RND_MODE);
+      mpfr_mul (y, op2->value.complex.i, op2->value.complex.i, GFC_RND_MODE);
+      mpfr_add (div, x, y, GFC_RND_MODE);
 
-      mpf_mul (x, op1->value.complex.r, op2->value.complex.r);
-      mpf_mul (y, op1->value.complex.i, op2->value.complex.i);
-      mpf_add (result->value.complex.r, x, y);
-      mpf_div (result->value.complex.r, result->value.complex.r, div);
+      mpfr_mul (x, op1->value.complex.r, op2->value.complex.r, GFC_RND_MODE);
+      mpfr_mul (y, op1->value.complex.i, op2->value.complex.i, GFC_RND_MODE);
+      mpfr_add (result->value.complex.r, x, y, GFC_RND_MODE);
+      mpfr_div (result->value.complex.r, result->value.complex.r, div,
+                GFC_RND_MODE);
 
-      mpf_mul (x, op1->value.complex.i, op2->value.complex.r);
-      mpf_mul (y, op1->value.complex.r, op2->value.complex.i);
-      mpf_sub (result->value.complex.i, x, y);
-      mpf_div (result->value.complex.i, result->value.complex.i, div);
+      mpfr_mul (x, op1->value.complex.i, op2->value.complex.r, GFC_RND_MODE);
+      mpfr_mul (y, op1->value.complex.r, op2->value.complex.i, GFC_RND_MODE);
+      mpfr_sub (result->value.complex.i, x, y, GFC_RND_MODE);
+      mpfr_div (result->value.complex.i, result->value.complex.i, div,
+                GFC_RND_MODE);
 
-      mpf_clear (x);
-      mpf_clear (y);
-      mpf_clear (div);
+      mpfr_clear (x);
+      mpfr_clear (y);
+      mpfr_clear (div);
 
       break;
 
@@ -1523,30 +1094,31 @@ gfc_arith_divide (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
 static void
 complex_reciprocal (gfc_expr * op)
 {
-  mpf_t mod, a, result_r, result_i;
-
-  mpf_init (mod);
-  mpf_init (a);
+  mpfr_t mod, a, re, im;
 
-  mpf_mul (mod, op->value.complex.r, op->value.complex.r);
-  mpf_mul (a, op->value.complex.i, op->value.complex.i);
-  mpf_add (mod, mod, a);
+  gfc_set_model (op->value.complex.r);
+  mpfr_init (mod);
+  mpfr_init (a);
+  mpfr_init (re);
+  mpfr_init (im);
 
-  mpf_init (result_r);
-  mpf_div (result_r, op->value.complex.r, mod);
+  /* FIXME:  another possible numerical problem.  */
+  mpfr_mul (mod, op->value.complex.r, op->value.complex.r, GFC_RND_MODE);
+  mpfr_mul (a, op->value.complex.i, op->value.complex.i, GFC_RND_MODE);
+  mpfr_add (mod, mod, a, GFC_RND_MODE);
 
-  mpf_init (result_i);
-  mpf_neg (result_i, op->value.complex.i);
-  mpf_div (result_i, result_i, mod);
+  mpfr_div (re, op->value.complex.r, mod, GFC_RND_MODE);
 
-  mpf_set (op->value.complex.r, result_r);
-  mpf_set (op->value.complex.i, result_i);
+  mpfr_neg (im, op->value.complex.i, GFC_RND_MODE);
+  mpfr_div (im, im, mod, GFC_RND_MODE);
 
-  mpf_clear (result_r);
-  mpf_clear (result_i);
+  mpfr_set (op->value.complex.r, re, GFC_RND_MODE);
+  mpfr_set (op->value.complex.i, im, GFC_RND_MODE);
 
-  mpf_clear (mod);
-  mpf_clear (a);
+  mpfr_clear (re);
+  mpfr_clear (im);
+  mpfr_clear (mod);
+  mpfr_clear (a);
 }
 
 
@@ -1555,32 +1127,37 @@ complex_reciprocal (gfc_expr * op)
 static void
 complex_pow_ui (gfc_expr * base, int power, gfc_expr * result)
 {
-  mpf_t temp_r, temp_i, a;
+  mpfr_t re, im, a;
 
-  mpf_set_ui (result->value.complex.r, 1);
-  mpf_set_ui (result->value.complex.i, 0);
+  gfc_set_model (base->value.complex.r);
+  mpfr_init (re);
+  mpfr_init (im);
+  mpfr_init (a);
 
-  mpf_init (temp_r);
-  mpf_init (temp_i);
-  mpf_init (a);
+  mpfr_set_ui (result->value.complex.r, 1, GFC_RND_MODE);
+  mpfr_set_ui (result->value.complex.i, 0, GFC_RND_MODE);
 
   for (; power > 0; power--)
     {
-      mpf_mul (temp_r, base->value.complex.r, result->value.complex.r);
-      mpf_mul (a, base->value.complex.i, result->value.complex.i);
-      mpf_sub (temp_r, temp_r, a);
+      mpfr_mul (re, base->value.complex.r, result->value.complex.r,
+                GFC_RND_MODE);
+      mpfr_mul (a, base->value.complex.i, result->value.complex.i,
+                GFC_RND_MODE);
+      mpfr_sub (re, re, a, GFC_RND_MODE);
 
-      mpf_mul (temp_i, base->value.complex.r, result->value.complex.i);
-      mpf_mul (a, base->value.complex.i, result->value.complex.r);
-      mpf_add (temp_i, temp_i, a);
+      mpfr_mul (im, base->value.complex.r, result->value.complex.i,
+                GFC_RND_MODE);
+      mpfr_mul (a, base->value.complex.i, result->value.complex.r,
+                GFC_RND_MODE);
+      mpfr_add (im, im, a, GFC_RND_MODE);
 
-      mpf_set (result->value.complex.r, temp_r);
-      mpf_set (result->value.complex.i, temp_i);
+      mpfr_set (result->value.complex.r, re, GFC_RND_MODE);
+      mpfr_set (result->value.complex.i, im, GFC_RND_MODE);
     }
 
-  mpf_clear (temp_r);
-  mpf_clear (temp_i);
-  mpf_clear (a);
+  mpfr_clear (re);
+  mpfr_clear (im);
+  mpfr_clear (a);
 }
 
 
@@ -1592,7 +1169,7 @@ gfc_arith_power (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
   int power, apower;
   gfc_expr *result;
   mpz_t unity_z;
-  mpf_t unity_f;
+  mpfr_t unity_f;
   arith rc;
 
   rc = ARITH_OK;
@@ -1611,25 +1188,23 @@ gfc_arith_power (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
 	    rc = ARITH_0TO0;
 	  else
 	    mpz_set_ui (result->value.integer, 1);
-
 	  break;
 
 	case BT_REAL:
-	  if (mpf_sgn (op1->value.real) == 0)
+	  if (mpfr_sgn (op1->value.real) == 0)
 	    rc = ARITH_0TO0;
 	  else
-	    mpf_set_ui (result->value.real, 1);
-
+	    mpfr_set_ui (result->value.real, 1, GFC_RND_MODE);
 	  break;
 
 	case BT_COMPLEX:
-	  if (mpf_sgn (op1->value.complex.r) == 0
-	      && mpf_sgn (op1->value.complex.i) == 0)
+	  if (mpfr_sgn (op1->value.complex.r) == 0
+	      && mpfr_sgn (op1->value.complex.i) == 0)
 	    rc = ARITH_0TO0;
 	  else
 	    {
-	      mpf_set_ui (result->value.complex.r, 1);
-	      mpf_set_ui (result->value.complex.i, 0);
+	      mpfr_set_ui (result->value.complex.r, 1, GFC_RND_MODE);
+	      mpfr_set_ui (result->value.complex.i, 0, GFC_RND_MODE);
 	    }
 
 	  break;
@@ -1638,8 +1213,7 @@ gfc_arith_power (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
 	  gfc_internal_error ("gfc_arith_power(): Bad base");
 	}
     }
-
-  if (power != 0)
+  else
     {
       apower = power;
       if (power < 0)
@@ -1661,22 +1235,24 @@ gfc_arith_power (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
 	  break;
 
 	case BT_REAL:
-	  mpf_pow_ui (result->value.real, op1->value.real, apower);
+	  mpfr_pow_ui (result->value.real, op1->value.real, apower,
+                       GFC_RND_MODE);
 
 	  if (power < 0)
 	    {
-	      mpf_init_set_ui (unity_f, 1);
-	      mpf_div (result->value.real, unity_f, result->value.real);
-	      mpf_clear (unity_f);
+              gfc_set_model (op1->value.real);
+	      mpfr_init (unity_f);
+	      mpfr_set_ui (unity_f, 1, GFC_RND_MODE);
+	      mpfr_div (result->value.real, unity_f, result->value.real,
+                        GFC_RND_MODE);
+	      mpfr_clear (unity_f);
 	    }
-
 	  break;
 
 	case BT_COMPLEX:
 	  complex_pow_ui (op1, apower, result);
 	  if (power < 0)
 	    complex_reciprocal (result);
-
 	  break;
 
 	default:
@@ -1748,7 +1324,7 @@ gfc_compare_expr (gfc_expr * op1, gfc_expr * op2)
       break;
 
     case BT_REAL:
-      rc = mpf_cmp (op1->value.real, op2->value.real);
+      rc = mpfr_cmp (op1->value.real, op2->value.real);
       break;
 
     case BT_CHARACTER:
@@ -1775,8 +1351,8 @@ static int
 compare_complex (gfc_expr * op1, gfc_expr * op2)
 {
 
-  return (mpf_cmp (op1->value.complex.r, op2->value.complex.r) == 0
-	  && mpf_cmp (op1->value.complex.i, op2->value.complex.i) == 0);
+  return (mpfr_cmp (op1->value.complex.r, op2->value.complex.r) == 0
+	  && mpfr_cmp (op1->value.complex.i, op2->value.complex.i) == 0);
 }
 
 
@@ -2544,12 +2120,12 @@ gfc_convert_real (const char *buffer, int kind, locus * where)
   const char *t;
 
   e = gfc_constant_result (BT_REAL, kind, where);
-  /* a leading plus is allowed, but not by mpf_set_str */
+  /* A leading plus is allowed in Fortran, but not by mpfr_set_str */
   if (buffer[0] == '+')
     t = buffer + 1;
   else
     t = buffer;
-  mpf_set_str (e->value.real, t, 10);
+  mpfr_set_str (e->value.real, t, 10, GFC_RND_MODE);
 
   return e;
 }
@@ -2564,8 +2140,8 @@ gfc_convert_complex (gfc_expr * real, gfc_expr * imag, int kind)
   gfc_expr *e;
 
   e = gfc_constant_result (BT_COMPLEX, kind, &real->where);
-  mpf_set (e->value.complex.r, real->value.real);
-  mpf_set (e->value.complex.i, imag->value.real);
+  mpfr_set (e->value.complex.r, real->value.real, GFC_RND_MODE);
+  mpfr_set (e->value.complex.i, imag->value.real, GFC_RND_MODE);
 
   return e;
 }
@@ -2621,7 +2197,7 @@ gfc_int2real (gfc_expr * src, int kind)
 
   result = gfc_constant_result (BT_REAL, kind, &src->where);
 
-  mpf_set_z (result->value.real, src->value.integer);
+  mpfr_set_z (result->value.real, src->value.integer, GFC_RND_MODE);
 
   if ((rc = gfc_check_real_range (result->value.real, kind)) != ARITH_OK)
     {
@@ -2644,8 +2220,8 @@ gfc_int2complex (gfc_expr * src, int kind)
 
   result = gfc_constant_result (BT_COMPLEX, kind, &src->where);
 
-  mpf_set_z (result->value.complex.r, src->value.integer);
-  mpf_set_ui (result->value.complex.i, 0);
+  mpfr_set_z (result->value.complex.r, src->value.integer, GFC_RND_MODE);
+  mpfr_set_ui (result->value.complex.i, 0, GFC_RND_MODE);
 
   if ((rc = gfc_check_real_range (result->value.complex.r, kind)) != ARITH_OK)
     {
@@ -2668,7 +2244,7 @@ gfc_real2int (gfc_expr * src, int kind)
 
   result = gfc_constant_result (BT_INTEGER, kind, &src->where);
 
-  mpz_set_f (result->value.integer, src->value.real);
+  gfc_mpfr_to_mpz (result->value.integer, src->value.real);
 
   if ((rc = gfc_check_integer_range (result->value.integer, kind))
       != ARITH_OK)
@@ -2692,7 +2268,7 @@ gfc_real2real (gfc_expr * src, int kind)
 
   result = gfc_constant_result (BT_REAL, kind, &src->where);
 
-  mpf_set (result->value.real, src->value.real);
+  mpfr_set (result->value.real, src->value.real, GFC_RND_MODE);
 
   rc = gfc_check_real_range (result->value.real, kind);
 
@@ -2700,7 +2276,7 @@ gfc_real2real (gfc_expr * src, int kind)
     {
       if (gfc_option.warn_underflow)
         gfc_warning ("%s at %L", gfc_arith_error (rc), &src->where);
-      mpf_set_ui(result->value.real, 0);
+      mpfr_set_ui(result->value.real, 0, GFC_RND_MODE);
     }
   else if (rc != ARITH_OK)
     {
@@ -2723,8 +2299,8 @@ gfc_real2complex (gfc_expr * src, int kind)
 
   result = gfc_constant_result (BT_COMPLEX, kind, &src->where);
 
-  mpf_set (result->value.complex.r, src->value.real);
-  mpf_set_ui (result->value.complex.i, 0);
+  mpfr_set (result->value.complex.r, src->value.real, GFC_RND_MODE);
+  mpfr_set_ui (result->value.complex.i, 0, GFC_RND_MODE);
 
   rc = gfc_check_real_range (result->value.complex.r, kind);
 
@@ -2732,7 +2308,7 @@ gfc_real2complex (gfc_expr * src, int kind)
     {
       if (gfc_option.warn_underflow)
         gfc_warning ("%s at %L", gfc_arith_error (rc), &src->where);
-      mpf_set_ui(result->value.complex.r, 0);
+      mpfr_set_ui(result->value.complex.r, 0, GFC_RND_MODE);
     }
   else if (rc != ARITH_OK)
     {
@@ -2755,7 +2331,7 @@ gfc_complex2int (gfc_expr * src, int kind)
 
   result = gfc_constant_result (BT_INTEGER, kind, &src->where);
 
-  mpz_set_f (result->value.integer, src->value.complex.r);
+  gfc_mpfr_to_mpz(result->value.integer, src->value.complex.r);
 
   if ((rc = gfc_check_integer_range (result->value.integer, kind))
       != ARITH_OK)
@@ -2779,7 +2355,7 @@ gfc_complex2real (gfc_expr * src, int kind)
 
   result = gfc_constant_result (BT_REAL, kind, &src->where);
 
-  mpf_set (result->value.real, src->value.complex.r);
+  mpfr_set (result->value.real, src->value.complex.r, GFC_RND_MODE);
 
   rc = gfc_check_real_range (result->value.real, kind);
 
@@ -2787,7 +2363,7 @@ gfc_complex2real (gfc_expr * src, int kind)
     {
       if (gfc_option.warn_underflow)
         gfc_warning ("%s at %L", gfc_arith_error (rc), &src->where);
-      mpf_set_ui(result->value.real, 0);
+      mpfr_set_ui(result->value.real, 0, GFC_RND_MODE);
     }
   if (rc != ARITH_OK)
     {
@@ -2810,8 +2386,8 @@ gfc_complex2complex (gfc_expr * src, int kind)
 
   result = gfc_constant_result (BT_COMPLEX, kind, &src->where);
 
-  mpf_set (result->value.complex.r, src->value.complex.r);
-  mpf_set (result->value.complex.i, src->value.complex.i);
+  mpfr_set (result->value.complex.r, src->value.complex.r, GFC_RND_MODE);
+  mpfr_set (result->value.complex.i, src->value.complex.i, GFC_RND_MODE);
 
   rc = gfc_check_real_range (result->value.complex.r, kind);
 
@@ -2819,7 +2395,7 @@ gfc_complex2complex (gfc_expr * src, int kind)
     {
       if (gfc_option.warn_underflow)
         gfc_warning ("%s at %L", gfc_arith_error (rc), &src->where);
-      mpf_set_ui(result->value.complex.r, 0);
+      mpfr_set_ui(result->value.complex.r, 0, GFC_RND_MODE);
     }
   else if (rc != ARITH_OK)
     {
@@ -2834,7 +2410,7 @@ gfc_complex2complex (gfc_expr * src, int kind)
     {
       if (gfc_option.warn_underflow)
         gfc_warning ("%s at %L", gfc_arith_error (rc), &src->where);
-      mpf_set_ui(result->value.complex.i, 0);
+      mpfr_set_ui(result->value.complex.i, 0, GFC_RND_MODE);
     }
   else if (rc != ARITH_OK)
     {
diff --git a/gcc/fortran/arith.h b/gcc/fortran/arith.h
index a1fdb1afd6b..1a718d4ea4c 100644
--- a/gcc/fortran/arith.h
+++ b/gcc/fortran/arith.h
@@ -24,19 +24,14 @@ Software Foundation, 59 Temple Place - Suite 330, Boston, MA
 
 #include "gfortran.h"
 
-/* Constants calculated during initialization.  */
-extern mpf_t pi, half_pi, two_pi, e;
-
-/* Calculate mathematically interesting functions.  */
-void natural_logarithm (mpf_t *, mpf_t *);
-void common_logarithm (mpf_t *, mpf_t *);
-void exponential (mpf_t *, mpf_t *);
-void sine (mpf_t *, mpf_t *);
-void cosine (mpf_t *, mpf_t *);
-void arctangent (mpf_t *, mpf_t *);
-void arctangent2 (mpf_t *, mpf_t *, mpf_t *);
-void hypercos (mpf_t *, mpf_t *);
-void hypersine (mpf_t *, mpf_t *);
+/* MPFR does not have mpfr_atan2(), which needs to return the principle
+   value of atan2().  MPFR also does not have the conversion of a mpfr_t
+   to a mpz_t, so declare a function for this as well.  */
+
+void arctangent2 (mpfr_t, mpfr_t, mpfr_t);
+void gfc_mpfr_to_mpz(mpz_t, mpfr_t);
+void gfc_set_model_kind (int);
+void gfc_set_model (mpfr_t);
 
 /* Return a constant result of a given type and kind, with locus.  */
 gfc_expr *gfc_constant_result (bt, int, locus *);
diff --git a/gcc/fortran/dump-parse-tree.c b/gcc/fortran/dump-parse-tree.c
index 8d23c908ed0..1c948d94253 100644
--- a/gcc/fortran/dump-parse-tree.c
+++ b/gcc/fortran/dump-parse-tree.c
@@ -363,7 +363,7 @@ gfc_show_expr (gfc_expr * p)
 	  break;
 
 	case BT_REAL:
-	  mpf_out_str (stdout, 10, 0, p->value.real);
+	  mpfr_out_str (stdout, 10, 0, p->value.real, GFC_RND_MODE);
 	  if (p->ts.kind != gfc_default_real_kind ())
 	    gfc_status ("_%d", p->ts.kind);
 	  break;
@@ -388,13 +388,13 @@ gfc_show_expr (gfc_expr * p)
 	case BT_COMPLEX:
 	  gfc_status ("(complex ");
 
-	  mpf_out_str (stdout, 10, 0, p->value.complex.r);
+	  mpfr_out_str (stdout, 10, 0, p->value.complex.r, GFC_RND_MODE);
 	  if (p->ts.kind != gfc_default_complex_kind ())
 	    gfc_status ("_%d", p->ts.kind);
 
 	  gfc_status (" ");
 
-	  mpf_out_str (stdout, 10, 0, p->value.complex.i);
+	  mpfr_out_str (stdout, 10, 0, p->value.complex.i, GFC_RND_MODE);
 	  if (p->ts.kind != gfc_default_complex_kind ())
 	    gfc_status ("_%d", p->ts.kind);
 
diff --git a/gcc/fortran/expr.c b/gcc/fortran/expr.c
index 74b785a5175..adff08e2070 100644
--- a/gcc/fortran/expr.c
+++ b/gcc/fortran/expr.c
@@ -154,7 +154,7 @@ free_expr0 (gfc_expr * e)
 	  break;
 
 	case BT_REAL:
-	  mpf_clear (e->value.real);
+	  mpfr_clear (e->value.real);
 	  break;
 
 	case BT_CHARACTER:
@@ -162,8 +162,8 @@ free_expr0 (gfc_expr * e)
 	  break;
 
 	case BT_COMPLEX:
-	  mpf_clear (e->value.complex.r);
-	  mpf_clear (e->value.complex.i);
+	  mpfr_clear (e->value.complex.r);
+	  mpfr_clear (e->value.complex.i);
 	  break;
 
 	default:
@@ -365,12 +365,17 @@ gfc_copy_expr (gfc_expr * p)
 	  break;
 
 	case BT_REAL:
-	  mpf_init_set (q->value.real, p->value.real);
+          gfc_set_model_kind (q->ts.kind);
+          mpfr_init (q->value.real);
+	  mpfr_set (q->value.real, p->value.real, GFC_RND_MODE);
 	  break;
 
 	case BT_COMPLEX:
-	  mpf_init_set (q->value.complex.r, p->value.complex.r);
-	  mpf_init_set (q->value.complex.i, p->value.complex.i);
+          gfc_set_model_kind (q->ts.kind);
+          mpfr_init (q->value.complex.r);
+          mpfr_init (q->value.complex.i);
+	  mpfr_set (q->value.complex.r, p->value.complex.r, GFC_RND_MODE);
+	  mpfr_set (q->value.complex.i, p->value.complex.i, GFC_RND_MODE);
 	  break;
 
 	case BT_CHARACTER:
diff --git a/gcc/fortran/gfortran.h b/gcc/fortran/gfortran.h
index 3ea8bb6431b..533479c63cd 100644
--- a/gcc/fortran/gfortran.h
+++ b/gcc/fortran/gfortran.h
@@ -59,7 +59,6 @@ char *alloca ();
 /* Major control parameters.  */
 
 #define GFC_MAX_SYMBOL_LEN 63
-#define GFC_REAL_BITS 100	/* Number of bits in g95's floating point numbers.  */
 #define GFC_MAX_LINE 132	/* Characters beyond this are not seen.  */
 #define GFC_MAX_DIMENSIONS 7	/* Maximum dimensions in an array.  */
 #define GFC_LETTERS 26		/* Number of letters in the alphabet.  */
@@ -184,7 +183,7 @@ extern mstring intrinsic_operators[];
 
 /* Arithmetic results.  */
 typedef enum
-{ ARITH_OK = 1, ARITH_OVERFLOW, ARITH_UNDERFLOW,
+{ ARITH_OK = 1, ARITH_OVERFLOW, ARITH_UNDERFLOW, ARITH_NAN,
   ARITH_DIV0, ARITH_0TO0, ARITH_INCOMMENSURATE
 }
 arith;
@@ -930,6 +929,8 @@ gfc_intrinsic_sym;
    EXPR_ARRAY      An array constructor.  */
 
 #include <gmp.h>
+#include <mpfr.h>
+#define GFC_RND_MODE GMP_RNDN
 
 typedef struct gfc_expr
 {
@@ -953,13 +954,14 @@ typedef struct gfc_expr
 
   union
   {
-    mpz_t integer;
-    mpf_t real;
     int logical;
+    mpz_t integer;
+
+    mpfr_t real;
 
     struct
     {
-      mpf_t r, i;
+      mpfr_t r, i;
     }
     complex;
 
@@ -1023,7 +1025,7 @@ typedef struct
   int kind, radix, digits, min_exponent, max_exponent;
 
   int range, precision;
-  mpf_t epsilon, huge, tiny;
+  mpfr_t epsilon, huge, tiny;
 }
 gfc_real_info;
 
@@ -1555,7 +1557,6 @@ match gfc_intrinsic_sub_interface (gfc_code *, int);
 
 /* simplify.c */
 void gfc_simplify_init_1 (void);
-void gfc_simplify_done_1 (void);
 
 /* match.c -- FIXME */
 void gfc_free_iterator (gfc_iterator *, int);
diff --git a/gcc/fortran/intrinsic.c b/gcc/fortran/intrinsic.c
index 022f1044e8e..659b507f6c5 100644
--- a/gcc/fortran/intrinsic.c
+++ b/gcc/fortran/intrinsic.c
@@ -1492,6 +1492,11 @@ add_functions (void)
              gfc_check_rand, NULL, NULL,
              i, BT_INTEGER, 4, 0);
 
+  /* Compatibility with HP FORTRAN 77/iX Reference.  Note, rand() and 
+     ran() use slightly different shoddy multiplicative congruential 
+     PRNG.  */
+  make_alias ("ran");
+
   make_generic ("rand", GFC_ISYM_RAND);
 
   add_sym_1 ("range", 0, 1, BT_INTEGER, di,
diff --git a/gcc/fortran/misc.c b/gcc/fortran/misc.c
index 3ef95db31e5..159a42150dd 100644
--- a/gcc/fortran/misc.c
+++ b/gcc/fortran/misc.c
@@ -309,7 +309,6 @@ gfc_done_1 (void)
 
   gfc_scanner_done_1 ();
   gfc_intrinsic_done_1 ();
-  gfc_simplify_done_1 ();
   gfc_iresolve_done_1 ();
   gfc_arith_done_1 ();
 }
diff --git a/gcc/fortran/module.c b/gcc/fortran/module.c
index 17254ce0f54..a9d0fa66c02 100644
--- a/gcc/fortran/module.c
+++ b/gcc/fortran/module.c
@@ -71,6 +71,7 @@ Software Foundation, 59 Temple Place - Suite 330, Boston, MA
 #include <time.h>
 
 #include "gfortran.h"
+#include "arith.h"
 #include "match.h"
 #include "parse.h" /* FIXME */
 
@@ -519,7 +520,7 @@ gfc_match_use (void)
 	tail->next = new;
       tail = new;
 
-      /* See what kind of interface we're dealing with.  Asusume it is
+      /* See what kind of interface we're dealing with.  Assume it is
          not an operator.  */
       new->operator = INTRINSIC_NONE;
       if (gfc_match_generic_spec (&type, name, &operator) == MATCH_ERROR)
@@ -2245,7 +2246,7 @@ mio_gmp_integer (mpz_t * integer)
 
 
 static void
-mio_gmp_real (mpf_t * real)
+mio_gmp_real (mpfr_t * real)
 {
   mp_exp_t exponent;
   char *p;
@@ -2255,14 +2256,14 @@ mio_gmp_real (mpf_t * real)
       if (parse_atom () != ATOM_STRING)
 	bad_module ("Expected real string");
 
-      mpf_init (*real);
-      mpf_set_str (*real, atom_string, -16);
+      mpfr_init (*real);
+      mpfr_set_str (*real, atom_string, 16, GFC_RND_MODE);
       gfc_free (atom_string);
 
     }
   else
     {
-      p = mpf_get_str (NULL, &exponent, 16, 0, *real);
+      p = mpfr_get_str (NULL, &exponent, 16, 0, *real, GFC_RND_MODE);
       atom_string = gfc_getmem (strlen (p) + 20);
 
       sprintf (atom_string, "0.%s@%ld", p, exponent);
@@ -2507,10 +2508,12 @@ mio_expr (gfc_expr ** ep)
 	  break;
 
 	case BT_REAL:
+          gfc_set_model_kind (e->ts.kind);
 	  mio_gmp_real (&e->value.real);
 	  break;
 
 	case BT_COMPLEX:
+          gfc_set_model_kind (e->ts.kind);
 	  mio_gmp_real (&e->value.complex.r);
 	  mio_gmp_real (&e->value.complex.i);
 	  break;
diff --git a/gcc/fortran/primary.c b/gcc/fortran/primary.c
index d054bfdb55b..eb5dc337f1d 100644
--- a/gcc/fortran/primary.c
+++ b/gcc/fortran/primary.c
@@ -1,5 +1,5 @@
 /* Primary expression subroutines
-   Copyright (C) 2000, 2001, 2002 Free Software Foundation, Inc.
+   Copyright (C) 2000, 2001, 2002, 2004 Free Software Foundation, Inc.
    Contributed by Andy Vaught
 
 This file is part of GNU G95.
@@ -436,7 +436,7 @@ done:
   buffer = alloca (count + 1);
   memset (buffer, '\0', count + 1);
 
-  /* Hack for mpf_init_set_str().  */
+  /* Hack for mpfr_set_str().  */
   p = buffer;
   while (count > 0)
     {
@@ -497,7 +497,7 @@ done:
     case ARITH_UNDERFLOW:
       if (gfc_option.warn_underflow)
         gfc_warning ("Real constant underflows its kind at %C");
-      mpf_set_ui(e->value.real, 0);
+      mpfr_set_ui (e->value.real, 0, GFC_RND_MODE);
       break;
 
     default:
@@ -1076,12 +1076,12 @@ done:
   buffer = alloca (count + 1);
   memset (buffer, '\0', count + 1);
 
-  /* Hack for mpf_init_set_str().  */
+  /* Hack for mpfr_set_str().  */
   p = buffer;
   while (count > 0)
     {
       c = gfc_next_char ();
-      if (c == 'd')
+      if (c == 'd' || c == 'q')
 	c = 'e';
       *p++ = c;
       count--;
diff --git a/gcc/fortran/simplify.c b/gcc/fortran/simplify.c
index d67b5c68ace..0a32d6f5cfc 100644
--- a/gcc/fortran/simplify.c
+++ b/gcc/fortran/simplify.c
@@ -30,9 +30,6 @@ Software Foundation, 59 Temple Place - Suite 330, Boston, MA
 #include "arith.h"
 #include "intrinsic.h"
 
-static mpf_t mpf_zero, mpf_half, mpf_one;
-static mpz_t mpz_zero;
-
 gfc_expr gfc_bad_expr;
 
 
@@ -148,7 +145,7 @@ gfc_expr *
 gfc_simplify_abs (gfc_expr * e)
 {
   gfc_expr *result;
-  mpf_t a, b;
+  mpfr_t a, b;
 
   if (e->expr_type != EXPR_CONSTANT)
     return NULL;
@@ -166,7 +163,7 @@ gfc_simplify_abs (gfc_expr * e)
     case BT_REAL:
       result = gfc_constant_result (BT_REAL, e->ts.kind, &e->where);
 
-      mpf_abs (result->value.real, e->value.real);
+      mpfr_abs (result->value.real, e->value.real, GFC_RND_MODE);
 
       result = range_check (result, "ABS");
       break;
@@ -174,17 +171,17 @@ gfc_simplify_abs (gfc_expr * e)
     case BT_COMPLEX:
       result = gfc_constant_result (BT_REAL, e->ts.kind, &e->where);
 
-      mpf_init (a);
-      mpf_mul (a, e->value.complex.r, e->value.complex.r);
-
-      mpf_init (b);
-      mpf_mul (b, e->value.complex.i, e->value.complex.i);
-
-      mpf_add (a, a, b);
-      mpf_sqrt (result->value.real, a);
+      gfc_set_model_kind (e->ts.kind);
+      mpfr_init (a);
+      mpfr_init (b);
+      /* FIXME:  Possible numerical problems.  */
+      mpfr_mul (a, e->value.complex.r, e->value.complex.r, GFC_RND_MODE);
+      mpfr_mul (b, e->value.complex.i, e->value.complex.i, GFC_RND_MODE);
+      mpfr_add (a, a, b, GFC_RND_MODE);
+      mpfr_sqrt (result->value.real, a, GFC_RND_MODE);
 
-      mpf_clear (a);
-      mpf_clear (b);
+      mpfr_clear (a);
+      mpfr_clear (b);
 
       result = range_check (result, "CABS");
       break;
@@ -231,12 +228,11 @@ gfc_expr *
 gfc_simplify_acos (gfc_expr * x)
 {
   gfc_expr *result;
-  mpf_t negative, square, term;
 
   if (x->expr_type != EXPR_CONSTANT)
     return NULL;
 
-  if (mpf_cmp_si (x->value.real, 1) > 0 || mpf_cmp_si (x->value.real, -1) < 0)
+  if (mpfr_cmp_si (x->value.real, 1) > 0 || mpfr_cmp_si (x->value.real, -1) < 0)
     {
       gfc_error ("Argument of ACOS at %L must be between -1 and 1",
 		 &x->where);
@@ -245,33 +241,7 @@ gfc_simplify_acos (gfc_expr * x)
 
   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
 
-  if (mpf_cmp_si (x->value.real, 1) == 0)
-    {
-      mpf_set_ui (result->value.real, 0);
-      return range_check (result, "ACOS");
-    }
-
-  if (mpf_cmp_si (x->value.real, -1) == 0)
-    {
-      mpf_set (result->value.real, pi);
-      return range_check (result, "ACOS");
-    }
-
-  mpf_init (negative);
-  mpf_init (square);
-  mpf_init (term);
-
-  mpf_pow_ui (square, x->value.real, 2);
-  mpf_ui_sub (term, 1, square);
-  mpf_sqrt (term, term);
-  mpf_div (term, x->value.real, term);
-  mpf_neg (term, term);
-  arctangent (&term, &negative);
-  mpf_add (result->value.real, half_pi, negative);
-
-  mpf_clear (negative);
-  mpf_clear (square);
-  mpf_clear (term);
+  mpfr_acos (result->value.real, x->value.real, GFC_RND_MODE);
 
   return range_check (result, "ACOS");
 }
@@ -370,7 +340,7 @@ gfc_simplify_aimag (gfc_expr * e)
     return NULL;
 
   result = gfc_constant_result (BT_REAL, e->ts.kind, &e->where);
-  mpf_set (result->value.real, e->value.complex.i);
+  mpfr_set (result->value.real, e->value.complex.i, GFC_RND_MODE);
 
   return range_check (result, "AIMAG");
 }
@@ -391,7 +361,7 @@ gfc_simplify_aint (gfc_expr * e, gfc_expr * k)
 
   rtrunc = gfc_copy_expr (e);
 
-  mpf_trunc (rtrunc->value.real, e->value.real);
+  mpfr_trunc (rtrunc->value.real, e->value.real);
 
   result = gfc_real2real (rtrunc, kind);
   gfc_free_expr (rtrunc);
@@ -410,7 +380,7 @@ gfc_simplify_dint (gfc_expr * e)
 
   rtrunc = gfc_copy_expr (e);
 
-  mpf_trunc (rtrunc->value.real, e->value.real);
+  mpfr_trunc (rtrunc->value.real, e->value.real);
 
   result = gfc_real2real (rtrunc, gfc_default_double_kind ());
   gfc_free_expr (rtrunc);
@@ -425,6 +395,7 @@ gfc_simplify_anint (gfc_expr * e, gfc_expr * k)
 {
   gfc_expr *rtrunc, *result;
   int kind, cmp;
+  mpfr_t half;
 
   kind = get_kind (BT_REAL, k, "ANINT", e->ts.kind);
   if (kind == -1)
@@ -437,22 +408,27 @@ gfc_simplify_anint (gfc_expr * e, gfc_expr * k)
 
   rtrunc = gfc_copy_expr (e);
 
-  cmp = mpf_cmp_ui (e->value.real, 0);
+  cmp = mpfr_cmp_ui (e->value.real, 0);
+
+  gfc_set_model_kind (kind);
+  mpfr_init (half);
+  mpfr_set_str (half, "0.5", 10, GFC_RND_MODE);
 
   if (cmp > 0)
     {
-      mpf_add (rtrunc->value.real, e->value.real, mpf_half);
-      mpf_trunc (result->value.real, rtrunc->value.real);
+      mpfr_add (rtrunc->value.real, e->value.real, half, GFC_RND_MODE);
+      mpfr_trunc (result->value.real, rtrunc->value.real);
     }
   else if (cmp < 0)
     {
-      mpf_sub (rtrunc->value.real, e->value.real, mpf_half);
-      mpf_trunc (result->value.real, rtrunc->value.real);
+      mpfr_sub (rtrunc->value.real, e->value.real, half, GFC_RND_MODE);
+      mpfr_trunc (result->value.real, rtrunc->value.real);
     }
   else
-    mpf_set_ui (result->value.real, 0);
+    mpfr_set_ui (result->value.real, 0, GFC_RND_MODE);
 
   gfc_free_expr (rtrunc);
+  mpfr_clear (half);
 
   return range_check (result, "ANINT");
 }
@@ -463,6 +439,7 @@ gfc_simplify_dnint (gfc_expr * e)
 {
   gfc_expr *rtrunc, *result;
   int cmp;
+  mpfr_t half;
 
   if (e->expr_type != EXPR_CONSTANT)
     return NULL;
@@ -472,22 +449,27 @@ gfc_simplify_dnint (gfc_expr * e)
 
   rtrunc = gfc_copy_expr (e);
 
-  cmp = mpf_cmp_ui (e->value.real, 0);
+  cmp = mpfr_cmp_ui (e->value.real, 0);
+
+  gfc_set_model_kind (gfc_default_double_kind ());
+  mpfr_init (half);
+  mpfr_set_str (half, "0.5", 10, GFC_RND_MODE);
 
   if (cmp > 0)
     {
-      mpf_add (rtrunc->value.real, e->value.real, mpf_half);
-      mpf_trunc (result->value.real, rtrunc->value.real);
+      mpfr_add (rtrunc->value.real, e->value.real, half, GFC_RND_MODE);
+      mpfr_trunc (result->value.real, rtrunc->value.real);
     }
   else if (cmp < 0)
     {
-      mpf_sub (rtrunc->value.real, e->value.real, mpf_half);
-      mpf_trunc (result->value.real, rtrunc->value.real);
+      mpfr_sub (rtrunc->value.real, e->value.real, half, GFC_RND_MODE);
+      mpfr_trunc (result->value.real, rtrunc->value.real);
     }
   else
-    mpf_set_ui (result->value.real, 0);
+    mpfr_set_ui (result->value.real, 0, GFC_RND_MODE);
 
   gfc_free_expr (rtrunc);
+  mpfr_clear (half);
 
   return range_check (result, "DNINT");
 }
@@ -497,12 +479,11 @@ gfc_expr *
 gfc_simplify_asin (gfc_expr * x)
 {
   gfc_expr *result;
-  mpf_t negative, square, term;
 
   if (x->expr_type != EXPR_CONSTANT)
     return NULL;
 
-  if (mpf_cmp_si (x->value.real, 1) > 0 || mpf_cmp_si (x->value.real, -1) < 0)
+  if (mpfr_cmp_si (x->value.real, 1) > 0 || mpfr_cmp_si (x->value.real, -1) < 0)
     {
       gfc_error ("Argument of ASIN at %L must be between -1 and 1",
 		 &x->where);
@@ -511,32 +492,7 @@ gfc_simplify_asin (gfc_expr * x)
 
   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
 
-  if (mpf_cmp_si (x->value.real, 1) == 0)
-    {
-      mpf_set (result->value.real, half_pi);
-      return range_check (result, "ASIN");
-    }
-
-  if (mpf_cmp_si (x->value.real, -1) == 0)
-    {
-      mpf_init (negative);
-      mpf_neg (negative, half_pi);
-      mpf_set (result->value.real, negative);
-      mpf_clear (negative);
-      return range_check (result, "ASIN");
-    }
-
-  mpf_init (square);
-  mpf_init (term);
-
-  mpf_pow_ui (square, x->value.real, 2);
-  mpf_ui_sub (term, 1, square);
-  mpf_sqrt (term, term);
-  mpf_div (term, x->value.real, term);
-  arctangent (&term, &result->value.real);
-
-  mpf_clear (square);
-  mpf_clear (term);
+  mpfr_asin(result->value.real, x->value.real, GFC_RND_MODE);
 
   return range_check (result, "ASIN");
 }
@@ -552,7 +508,7 @@ gfc_simplify_atan (gfc_expr * x)
 
   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
 
-  arctangent (&x->value.real, &result->value.real);
+  mpfr_atan(result->value.real, x->value.real, GFC_RND_MODE);
 
   return range_check (result, "ATAN");
 
@@ -569,17 +525,16 @@ gfc_simplify_atan2 (gfc_expr * y, gfc_expr * x)
 
   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
 
-
-  if (mpf_sgn (y->value.real) == 0 && mpf_sgn (x->value.real) == 0)
+  if (mpfr_sgn (y->value.real) == 0 && mpfr_sgn (x->value.real) == 0)
     {
       gfc_error
-	("If first argument of ATAN2 %L is zero, the second argument "
+	("If first argument of ATAN2 %L is zero, then the second argument "
 	  "must not be zero", &x->where);
       gfc_free_expr (result);
       return &gfc_bad_expr;
     }
 
-  arctangent2 (&y->value.real, &x->value.real, &result->value.real);
+  arctangent2 (y->value.real, x->value.real, result->value.real);
 
   return range_check (result, "ATAN2");
 
@@ -635,8 +590,8 @@ gfc_simplify_ceiling (gfc_expr * e, gfc_expr * k)
 
   ceil = gfc_copy_expr (e);
 
-  mpf_ceil (ceil->value.real, e->value.real);
-  mpz_set_f (result->value.integer, ceil->value.real);
+  mpfr_ceil (ceil->value.real, e->value.real);
+  gfc_mpfr_to_mpz(result->value.integer, ceil->value.real);
 
   gfc_free_expr (ceil);
 
@@ -684,21 +639,21 @@ simplify_cmplx (const char *name, gfc_expr * x, gfc_expr * y, int kind)
 
   result = gfc_constant_result (BT_COMPLEX, kind, &x->where);
 
-  mpf_set_ui (result->value.complex.i, 0);
+  mpfr_set_ui (result->value.complex.i, 0, GFC_RND_MODE);
 
   switch (x->ts.type)
     {
     case BT_INTEGER:
-      mpf_set_z (result->value.complex.r, x->value.integer);
+      mpfr_set_z (result->value.complex.r, x->value.integer, GFC_RND_MODE);
       break;
 
     case BT_REAL:
-      mpf_set (result->value.complex.r, x->value.real);
+      mpfr_set (result->value.complex.r, x->value.real, GFC_RND_MODE);
       break;
 
     case BT_COMPLEX:
-      mpf_set (result->value.complex.r, x->value.complex.r);
-      mpf_set (result->value.complex.i, x->value.complex.i);
+      mpfr_set (result->value.complex.r, x->value.complex.r, GFC_RND_MODE);
+      mpfr_set (result->value.complex.i, x->value.complex.i, GFC_RND_MODE);
       break;
 
     default:
@@ -710,11 +665,11 @@ simplify_cmplx (const char *name, gfc_expr * x, gfc_expr * y, int kind)
       switch (y->ts.type)
 	{
 	case BT_INTEGER:
-	  mpf_set_z (result->value.complex.i, y->value.integer);
+	  mpfr_set_z (result->value.complex.i, y->value.integer, GFC_RND_MODE);
 	  break;
 
 	case BT_REAL:
-	  mpf_set (result->value.complex.i, y->value.real);
+	  mpfr_set (result->value.complex.i, y->value.real, GFC_RND_MODE);
 	  break;
 
 	default:
@@ -752,7 +707,7 @@ gfc_simplify_conjg (gfc_expr * e)
     return NULL;
 
   result = gfc_copy_expr (e);
-  mpf_neg (result->value.complex.i, result->value.complex.i);
+  mpfr_neg (result->value.complex.i, result->value.complex.i, GFC_RND_MODE);
 
   return range_check (result, "CONJG");
 }
@@ -762,7 +717,7 @@ gfc_expr *
 gfc_simplify_cos (gfc_expr * x)
 {
   gfc_expr *result;
-  mpf_t xp, xq;
+  mpfr_t xp, xq;
 
   if (x->expr_type != EXPR_CONSTANT)
     return NULL;
@@ -772,23 +727,24 @@ gfc_simplify_cos (gfc_expr * x)
   switch (x->ts.type)
     {
     case BT_REAL:
-      cosine (&x->value.real, &result->value.real);
+      mpfr_cos (result->value.real, x->value.real, GFC_RND_MODE);
       break;
     case BT_COMPLEX:
-      mpf_init (xp);
-      mpf_init (xq);
+      gfc_set_model_kind (x->ts.kind);
+      mpfr_init (xp);
+      mpfr_init (xq);
 
-      cosine (&x->value.complex.r, &xp);
-      hypercos (&x->value.complex.i, &xq);
-      mpf_mul (result->value.complex.r, xp, xq);
+      mpfr_cos  (xp, x->value.complex.r, GFC_RND_MODE);
+      mpfr_cosh (xq, x->value.complex.i, GFC_RND_MODE);
+      mpfr_mul(result->value.complex.r, xp, xq, GFC_RND_MODE);
 
-      sine (&x->value.complex.r, &xp);
-      hypersine (&x->value.complex.i, &xq);
-      mpf_mul (xp, xp, xq);
-      mpf_neg (result->value.complex.i, xp);
+      mpfr_sin  (xp, x->value.complex.r, GFC_RND_MODE);
+      mpfr_sinh (xq, x->value.complex.i, GFC_RND_MODE);
+      mpfr_mul (xp, xp, xq, GFC_RND_MODE);
+      mpfr_neg (result->value.complex.i, xp, GFC_RND_MODE );
 
-      mpf_clear (xp);
-      mpf_clear (xq);
+      mpfr_clear (xp);
+      mpfr_clear (xq);
       break;
     default:
       gfc_internal_error ("in gfc_simplify_cos(): Bad type");
@@ -809,7 +765,7 @@ gfc_simplify_cosh (gfc_expr * x)
 
   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
 
-  hypercos (&x->value.real, &result->value.real);
+  mpfr_cosh (result->value.real, x->value.real, GFC_RND_MODE);
 
   return range_check (result, "COSH");
 }
@@ -902,15 +858,15 @@ gfc_simplify_dim (gfc_expr * x, gfc_expr * y)
       if (mpz_cmp (x->value.integer, y->value.integer) > 0)
 	mpz_sub (result->value.integer, x->value.integer, y->value.integer);
       else
-	mpz_set (result->value.integer, mpz_zero);
+	mpz_set_ui (result->value.integer, 0);
 
       break;
 
     case BT_REAL:
-      if (mpf_cmp (x->value.real, y->value.real) > 0)
-	mpf_sub (result->value.real, x->value.real, y->value.real);
+      if (mpfr_cmp (x->value.real, y->value.real) > 0)
+	mpfr_sub (result->value.real, x->value.real, y->value.real, GFC_RND_MODE);
       else
-	mpf_set (result->value.real, mpf_zero);
+	mpfr_set_ui (result->value.real, 0, GFC_RND_MODE);
 
       break;
 
@@ -925,7 +881,7 @@ gfc_simplify_dim (gfc_expr * x, gfc_expr * y)
 gfc_expr *
 gfc_simplify_dprod (gfc_expr * x, gfc_expr * y)
 {
-  gfc_expr *mult1, *mult2, *result;
+  gfc_expr *a1, *a2, *result;
 
   if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
     return NULL;
@@ -933,13 +889,13 @@ gfc_simplify_dprod (gfc_expr * x, gfc_expr * y)
   result =
     gfc_constant_result (BT_REAL, gfc_default_double_kind (), &x->where);
 
-  mult1 = gfc_real2real (x, gfc_default_double_kind ());
-  mult2 = gfc_real2real (y, gfc_default_double_kind ());
+  a1 = gfc_real2real (x, gfc_default_double_kind ());
+  a2 = gfc_real2real (y, gfc_default_double_kind ());
 
-  mpf_mul (result->value.real, mult1->value.real, mult2->value.real);
+  mpfr_mul (result->value.real, a1->value.real, a2->value.real, GFC_RND_MODE);
 
-  gfc_free_expr (mult1);
-  gfc_free_expr (mult2);
+  gfc_free_expr (a1);
+  gfc_free_expr (a2);
 
   return range_check (result, "DPROD");
 }
@@ -957,7 +913,7 @@ gfc_simplify_epsilon (gfc_expr * e)
 
   result = gfc_constant_result (BT_REAL, e->ts.kind, &e->where);
 
-  mpf_set (result->value.real, gfc_real_kinds[i].epsilon);
+  mpfr_set (result->value.real, gfc_real_kinds[i].epsilon, GFC_RND_MODE);
 
   return range_check (result, "EPSILON");
 }
@@ -967,86 +923,30 @@ gfc_expr *
 gfc_simplify_exp (gfc_expr * x)
 {
   gfc_expr *result;
-  mpf_t xp, xq;
-  double ln2, absval, rhuge;
+  mpfr_t xp, xq;
 
   if (x->expr_type != EXPR_CONSTANT)
     return NULL;
 
   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
 
-  /* Exactitude doesn't matter here */
-  ln2 = .6931472;
-  rhuge = ln2 * mpz_get_d (gfc_integer_kinds[0].huge);
-
   switch (x->ts.type)
     {
     case BT_REAL:
-      absval = mpf_get_d (x->value.real);
-      if (absval < 0)
-	absval = -absval;
-      if (absval > rhuge)
-	{
-	  /* Underflow (set arg to zero) if x is negative and its
-	     magnitude is greater than the maximum C long int times
-	     ln2, because the exponential method in arith.c will fail
-	     for such values.  */
-	  if (mpf_cmp_ui (x->value.real, 0) < 0)
-	    {
-	      if (pedantic == 1)
-		gfc_warning_now
-		  ("Argument of EXP at %L is negative and too large, "
-		   "setting result to zero", &x->where);
-	      mpf_set_ui (result->value.real, 0);
-	      return range_check (result, "EXP");
-	    }
-	  /* Overflow if magnitude of x is greater than C long int
-             huge times ln2.  */
-	  else
-	    {
-	      gfc_error ("Argument of EXP at %L too large", &x->where);
-	      gfc_free_expr (result);
-	      return &gfc_bad_expr;
-	    }
-	}
-      exponential (&x->value.real, &result->value.real);
+      mpfr_exp(result->value.real, x->value.real, GFC_RND_MODE);
       break;
 
     case BT_COMPLEX:
-      /* Using Euler's formula.  */
-      absval = mpf_get_d (x->value.complex.r);
-      if (absval < 0)
-	absval = -absval;
-      if (absval > rhuge)
-	{
-	  if (mpf_cmp_ui (x->value.complex.r, 0) < 0)
-	    {
-	      if (pedantic == 1)
-		gfc_warning_now
-		  ("Real part of argument of EXP at %L is negative "
-		   "and too large, setting result to zero", &x->where);
-
-	      mpf_set_ui (result->value.complex.r, 0);
-	      mpf_set_ui (result->value.complex.i, 0);
-	      return range_check (result, "EXP");
-	    }
-	  else
-	    {
-	      gfc_error ("Real part of argument of EXP at %L too large",
-			 &x->where);
-	      gfc_free_expr (result);
-	      return &gfc_bad_expr;
-	    }
-	}
-      mpf_init (xp);
-      mpf_init (xq);
-      exponential (&x->value.complex.r, &xq);
-      cosine (&x->value.complex.i, &xp);
-      mpf_mul (result->value.complex.r, xq, xp);
-      sine (&x->value.complex.i, &xp);
-      mpf_mul (result->value.complex.i, xq, xp);
-      mpf_clear (xp);
-      mpf_clear (xq);
+      gfc_set_model_kind (x->ts.kind);
+      mpfr_init (xp);
+      mpfr_init (xq);
+      mpfr_exp (xq, x->value.complex.r, GFC_RND_MODE);
+      mpfr_cos (xp, x->value.complex.i, GFC_RND_MODE);
+      mpfr_mul (result->value.complex.r, xq, xp, GFC_RND_MODE);
+      mpfr_sin (xp, x->value.complex.i, GFC_RND_MODE);
+      mpfr_mul (result->value.complex.i, xq, xp, GFC_RND_MODE);
+      mpfr_clear (xp);
+      mpfr_clear (xq);
       break;
 
     default:
@@ -1056,11 +956,11 @@ gfc_simplify_exp (gfc_expr * x)
   return range_check (result, "EXP");
 }
 
-
+/* FIXME:  MPFR should be able to do this better */
 gfc_expr *
 gfc_simplify_exponent (gfc_expr * x)
 {
-  mpf_t i2, absv, ln2, lnx;
+  mpfr_t i2, absv, ln2, lnx, zero;
   gfc_expr *result;
 
   if (x->expr_type != EXPR_CONSTANT)
@@ -1069,31 +969,39 @@ gfc_simplify_exponent (gfc_expr * x)
   result = gfc_constant_result (BT_INTEGER, gfc_default_integer_kind (),
 				&x->where);
 
-  if (mpf_cmp (x->value.real, mpf_zero) == 0)
+  gfc_set_model (x->value.real);
+  mpfr_init (zero);
+  mpfr_set_ui (zero, 0, GFC_RND_MODE);
+
+  if (mpfr_cmp (x->value.real, zero) == 0)
     {
       mpz_set_ui (result->value.integer, 0);
+      mpfr_clear (zero);
       return result;
     }
 
-  mpf_init_set_ui (i2, 2);
-  mpf_init (absv);
-  mpf_init (ln2);
-  mpf_init (lnx);
+  mpfr_init (i2);
+  mpfr_init (absv);
+  mpfr_init (ln2);
+  mpfr_init (lnx);
+
+  mpfr_set_ui (i2, 2, GFC_RND_MODE);
 
-  natural_logarithm (&i2, &ln2);
+  mpfr_log (ln2, i2, GFC_RND_MODE); 
+  mpfr_abs (absv, x->value.real, GFC_RND_MODE);
+  mpfr_log (lnx, absv, GFC_RND_MODE); 
 
-  mpf_abs (absv, x->value.real);
-  natural_logarithm (&absv, &lnx);
+  mpfr_div (lnx, lnx, ln2, GFC_RND_MODE);
+  mpfr_trunc (lnx, lnx);
+  mpfr_add_ui (lnx, lnx, 1, GFC_RND_MODE);
 
-  mpf_div (lnx, lnx, ln2);
-  mpf_trunc (lnx, lnx);
-  mpf_add_ui (lnx, lnx, 1);
-  mpz_set_f (result->value.integer, lnx);
+  gfc_mpfr_to_mpz (result->value.integer, lnx);
 
-  mpf_clear (i2);
-  mpf_clear (ln2);
-  mpf_clear (lnx);
-  mpf_clear (absv);
+  mpfr_clear (i2);
+  mpfr_clear (ln2);
+  mpfr_clear (lnx);
+  mpfr_clear (absv);
+  mpfr_clear (zero);
 
   return range_check (result, "EXPONENT");
 }
@@ -1116,7 +1024,7 @@ gfc_expr *
 gfc_simplify_floor (gfc_expr * e, gfc_expr * k)
 {
   gfc_expr *result;
-  mpf_t floor;
+  mpfr_t floor;
   int kind;
 
   kind = get_kind (BT_REAL, k, "FLOOR", gfc_default_real_kind ());
@@ -1128,10 +1036,13 @@ gfc_simplify_floor (gfc_expr * e, gfc_expr * k)
 
   result = gfc_constant_result (BT_INTEGER, kind, &e->where);
 
-  mpf_init (floor);
-  mpf_floor (floor, e->value.real);
-  mpz_set_f (result->value.integer, floor);
-  mpf_clear (floor);
+  gfc_set_model_kind (kind);
+  mpfr_init (floor);
+  mpfr_floor (floor, e->value.real);
+
+  gfc_mpfr_to_mpz (result->value.integer, floor);
+
+  mpfr_clear (floor);
 
   return range_check (result, "FLOOR");
 }
@@ -1141,7 +1052,7 @@ gfc_expr *
 gfc_simplify_fraction (gfc_expr * x)
 {
   gfc_expr *result;
-  mpf_t i2, absv, ln2, lnx, pow2;
+  mpfr_t i2, absv, ln2, lnx, pow2, zero;
   unsigned long exp2;
 
   if (x->expr_type != EXPR_CONSTANT)
@@ -1149,37 +1060,44 @@ gfc_simplify_fraction (gfc_expr * x)
 
   result = gfc_constant_result (BT_REAL, x->ts.kind, &x->where);
 
-  if (mpf_cmp (x->value.real, mpf_zero) == 0)
+  gfc_set_model_kind (x->ts.kind);
+  mpfr_init (zero);
+  mpfr_set_ui (zero, 0, GFC_RND_MODE);
+
+  if (mpfr_cmp (x->value.real, zero) == 0)
     {
-      mpf_set (result->value.real, mpf_zero);
+      mpfr_set (result->value.real, zero, GFC_RND_MODE);
+      mpfr_clear (zero);
       return result;
     }
 
-  mpf_init_set_ui (i2, 2);
-  mpf_init (absv);
-  mpf_init (ln2);
-  mpf_init (lnx);
-  mpf_init (pow2);
+  mpfr_init (i2);
+  mpfr_init (absv);
+  mpfr_init (ln2);
+  mpfr_init (lnx);
+  mpfr_init (pow2);
 
-  natural_logarithm (&i2, &ln2);
+  mpfr_set_ui (i2, 2, GFC_RND_MODE);
 
-  mpf_abs (absv, x->value.real);
-  natural_logarithm (&absv, &lnx);
+  mpfr_log (ln2, i2, GFC_RND_MODE);
+  mpfr_abs (absv, x->value.real, GFC_RND_MODE);
+  mpfr_log (lnx, absv, GFC_RND_MODE);
 
-  mpf_div (lnx, lnx, ln2);
-  mpf_trunc (lnx, lnx);
-  mpf_add_ui (lnx, lnx, 1);
+  mpfr_div (lnx, lnx, ln2, GFC_RND_MODE);
+  mpfr_trunc (lnx, lnx);
+  mpfr_add_ui (lnx, lnx, 1, GFC_RND_MODE);
 
-  exp2 = (unsigned long) mpf_get_d (lnx);
-  mpf_pow_ui (pow2, i2, exp2);
+  exp2 = (unsigned long) mpfr_get_d (lnx, GFC_RND_MODE);
+  mpfr_pow_ui (pow2, i2, exp2, GFC_RND_MODE);
 
-  mpf_div (result->value.real, absv, pow2);
+  mpfr_div (result->value.real, absv, pow2, GFC_RND_MODE);
 
-  mpf_clear (i2);
-  mpf_clear (ln2);
-  mpf_clear (absv);
-  mpf_clear (lnx);
-  mpf_clear (pow2);
+  mpfr_clear (i2);
+  mpfr_clear (ln2);
+  mpfr_clear (absv);
+  mpfr_clear (lnx);
+  mpfr_clear (pow2);
+  mpfr_clear (zero);
 
   return range_check (result, "FRACTION");
 }
@@ -1204,7 +1122,7 @@ gfc_simplify_huge (gfc_expr * e)
       break;
 
     case BT_REAL:
-      mpf_set (result->value.real, gfc_real_kinds[i].huge);
+      mpfr_set (result->value.real, gfc_real_kinds[i].huge, GFC_RND_MODE);
       break;
 
     bad_type:
@@ -1605,16 +1523,16 @@ gfc_simplify_int (gfc_expr * e, gfc_expr * k)
 
     case BT_REAL:
       rtrunc = gfc_copy_expr (e);
-      mpf_trunc (rtrunc->value.real, e->value.real);
-      mpz_set_f (result->value.integer, rtrunc->value.real);
+      mpfr_trunc (rtrunc->value.real, e->value.real);
+      gfc_mpfr_to_mpz (result->value.integer, rtrunc->value.real);
       gfc_free_expr (rtrunc);
       break;
 
     case BT_COMPLEX:
       rpart = gfc_complex2real (e, kind);
       rtrunc = gfc_copy_expr (rpart);
-      mpf_trunc (rtrunc->value.real, rpart->value.real);
-      mpz_set_f (result->value.integer, rtrunc->value.real);
+      mpfr_trunc (rtrunc->value.real, rpart->value.real);
+      gfc_mpfr_to_mpz (result->value.integer, rtrunc->value.real);
       gfc_free_expr (rpart);
       gfc_free_expr (rtrunc);
       break;
@@ -1642,8 +1560,8 @@ gfc_simplify_ifix (gfc_expr * e)
 
   rtrunc = gfc_copy_expr (e);
 
-  mpf_trunc (rtrunc->value.real, e->value.real);
-  mpz_set_f (result->value.integer, rtrunc->value.real);
+  mpfr_trunc (rtrunc->value.real, e->value.real);
+  gfc_mpfr_to_mpz (result->value.integer, rtrunc->value.real);
 
   gfc_free_expr (rtrunc);
   return range_check (result, "IFIX");
@@ -1663,8 +1581,8 @@ gfc_simplify_idint (gfc_expr * e)
 
   rtrunc = gfc_copy_expr (e);
 
-  mpf_trunc (rtrunc->value.real, e->value.real);
-  mpz_set_f (result->value.integer, rtrunc->value.real);
+  mpfr_trunc (rtrunc->value.real, e->value.real);
+  gfc_mpfr_to_mpz (result->value.integer, rtrunc->value.real);
 
   gfc_free_expr (rtrunc);
   return range_check (result, "IDINT");
@@ -2000,52 +1918,60 @@ gfc_expr *
 gfc_simplify_log (gfc_expr * x)
 {
   gfc_expr *result;
-  mpf_t xr, xi;
+  mpfr_t xr, xi, zero;
 
   if (x->expr_type != EXPR_CONSTANT)
     return NULL;
 
   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
 
+  gfc_set_model_kind (x->ts.kind);
+  mpfr_init (zero);
+  mpfr_set_ui (zero, 0, GFC_RND_MODE);
+
   switch (x->ts.type)
     {
     case BT_REAL:
-      if (mpf_cmp (x->value.real, mpf_zero) <= 0)
+      if (mpfr_cmp (x->value.real, zero) <= 0)
 	{
 	  gfc_error
 	    ("Argument of LOG at %L cannot be less than or equal to zero",
 	     &x->where);
 	  gfc_free_expr (result);
+          mpfr_clear (zero);
 	  return &gfc_bad_expr;
 	}
 
-      natural_logarithm (&x->value.real, &result->value.real);
+      mpfr_log(result->value.real, x->value.real, GFC_RND_MODE);
+      mpfr_clear (zero);
       break;
 
     case BT_COMPLEX:
-      if ((mpf_cmp (x->value.complex.r, mpf_zero) == 0)
-	  && (mpf_cmp (x->value.complex.i, mpf_zero) == 0))
+      if ((mpfr_cmp (x->value.complex.r, zero) == 0)
+	  && (mpfr_cmp (x->value.complex.i, zero) == 0))
 	{
 	  gfc_error ("Complex argument of LOG at %L cannot be zero",
 		     &x->where);
 	  gfc_free_expr (result);
+          mpfr_clear (zero);
 	  return &gfc_bad_expr;
 	}
 
-      mpf_init (xr);
-      mpf_init (xi);
+      mpfr_init (xr);
+      mpfr_init (xi);
 
-      arctangent2 (&x->value.complex.i, &x->value.complex.r,
-	&result->value.complex.i);
+      arctangent2 (x->value.complex.i, x->value.complex.r,
+	           result->value.complex.i);
 
-      mpf_mul (xr, x->value.complex.r, x->value.complex.r);
-      mpf_mul (xi, x->value.complex.i, x->value.complex.i);
-      mpf_add (xr, xr, xi);
-      mpf_sqrt (xr, xr);
-      natural_logarithm (&xr, &result->value.complex.r);
+      mpfr_mul (xr, x->value.complex.r, x->value.complex.r, GFC_RND_MODE);
+      mpfr_mul (xi, x->value.complex.i, x->value.complex.i, GFC_RND_MODE);
+      mpfr_add (xr, xr, xi, GFC_RND_MODE);
+      mpfr_sqrt (xr, xr, GFC_RND_MODE);
+      mpfr_log (result->value.complex.r, xr, GFC_RND_MODE);
 
-      mpf_clear (xr);
-      mpf_clear (xi);
+      mpfr_clear (xr);
+      mpfr_clear (xi);
+      mpfr_clear (zero);
 
       break;
 
@@ -2061,21 +1987,28 @@ gfc_expr *
 gfc_simplify_log10 (gfc_expr * x)
 {
   gfc_expr *result;
+  mpfr_t zero;
 
   if (x->expr_type != EXPR_CONSTANT)
     return NULL;
 
-  if (mpf_cmp (x->value.real, mpf_zero) <= 0)
+  gfc_set_model_kind (x->ts.kind);
+  mpfr_init (zero);
+  mpfr_set_ui (zero, 0, GFC_RND_MODE);
+
+  if (mpfr_cmp (x->value.real, zero) <= 0)
     {
       gfc_error
 	("Argument of LOG10 at %L cannot be less than or equal to zero",
 	 &x->where);
+      mpfr_clear (zero);
       return &gfc_bad_expr;
     }
 
   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
 
-  common_logarithm (&x->value.real, &result->value.real);
+  mpfr_log10 (result->value.real, x->value.real, GFC_RND_MODE);
+  mpfr_clear (zero);
 
   return range_check (result, "LOG10");
 }
@@ -2142,9 +2075,10 @@ simplify_min_max (gfc_expr * expr, int sign)
 	  break;
 
 	case BT_REAL:
-	  if (mpf_cmp (arg->expr->value.real, extremum->expr->value.real) *
+	  if (mpfr_cmp (arg->expr->value.real, extremum->expr->value.real) *
 	      sign > 0)
-	    mpf_set (extremum->expr->value.real, arg->expr->value.real);
+	    mpfr_set (extremum->expr->value.real, arg->expr->value.real,
+                      GFC_RND_MODE);
 
 	  break;
 
@@ -2235,7 +2169,7 @@ gfc_expr *
 gfc_simplify_mod (gfc_expr * a, gfc_expr * p)
 {
   gfc_expr *result;
-  mpf_t quot, iquot, term;
+  mpfr_t quot, iquot, term;
 
   if (a->expr_type != EXPR_CONSTANT || p->expr_type != EXPR_CONSTANT)
     return NULL;
@@ -2256,7 +2190,7 @@ gfc_simplify_mod (gfc_expr * a, gfc_expr * p)
       break;
 
     case BT_REAL:
-      if (mpf_cmp_ui (p->value.real, 0) == 0)
+      if (mpfr_cmp_ui (p->value.real, 0) == 0)
 	{
 	  /* Result is processor-dependent.  */
 	  gfc_error ("Second argument of MOD at %L is zero", &p->where);
@@ -2264,18 +2198,19 @@ gfc_simplify_mod (gfc_expr * a, gfc_expr * p)
 	  return &gfc_bad_expr;
 	}
 
-      mpf_init (quot);
-      mpf_init (iquot);
-      mpf_init (term);
+      gfc_set_model_kind (a->ts.kind);
+      mpfr_init (quot);
+      mpfr_init (iquot);
+      mpfr_init (term);
 
-      mpf_div (quot, a->value.real, p->value.real);
-      mpf_trunc (iquot, quot);
-      mpf_mul (term, iquot, p->value.real);
-      mpf_sub (result->value.real, a->value.real, term);
+      mpfr_div (quot, a->value.real, p->value.real, GFC_RND_MODE);
+      mpfr_trunc (iquot, quot);
+      mpfr_mul (term, iquot, p->value.real, GFC_RND_MODE);
+      mpfr_sub (result->value.real, a->value.real, term, GFC_RND_MODE);
 
-      mpf_clear (quot);
-      mpf_clear (iquot);
-      mpf_clear (term);
+      mpfr_clear (quot);
+      mpfr_clear (iquot);
+      mpfr_clear (term);
       break;
 
     default:
@@ -2290,7 +2225,7 @@ gfc_expr *
 gfc_simplify_modulo (gfc_expr * a, gfc_expr * p)
 {
   gfc_expr *result;
-  mpf_t quot, iquot, term;
+  mpfr_t quot, iquot, term;
 
   if (a->expr_type != EXPR_CONSTANT || p->expr_type != EXPR_CONSTANT)
     return NULL;
@@ -2313,7 +2248,7 @@ gfc_simplify_modulo (gfc_expr * a, gfc_expr * p)
       break;
 
     case BT_REAL:
-      if (mpf_cmp_ui (p->value.real, 0) == 0)
+      if (mpfr_cmp_ui (p->value.real, 0) == 0)
 	{
 	  /* Result is processor-dependent.  */
 	  gfc_error ("Second argument of MODULO at %L is zero", &p->where);
@@ -2321,19 +2256,20 @@ gfc_simplify_modulo (gfc_expr * a, gfc_expr * p)
 	  return &gfc_bad_expr;
 	}
 
-      mpf_init (quot);
-      mpf_init (iquot);
-      mpf_init (term);
+      gfc_set_model_kind (a->ts.kind);
+      mpfr_init (quot);
+      mpfr_init (iquot);
+      mpfr_init (term);
 
-      mpf_div (quot, a->value.real, p->value.real);
-      mpf_floor (iquot, quot);
-      mpf_mul (term, iquot, p->value.real);
+      mpfr_div (quot, a->value.real, p->value.real, GFC_RND_MODE);
+      mpfr_floor (iquot, quot);
+      mpfr_mul (term, iquot, p->value.real, GFC_RND_MODE);
 
-      mpf_clear (quot);
-      mpf_clear (iquot);
-      mpf_clear (term);
+      mpfr_clear (quot);
+      mpfr_clear (iquot);
+      mpfr_clear (term);
 
-      mpf_sub (result->value.real, a->value.real, term);
+      mpfr_sub (result->value.real, a->value.real, term, GFC_RND_MODE);
       break;
 
     default:
@@ -2376,7 +2312,7 @@ gfc_simplify_nearest (gfc_expr * x, gfc_expr * s)
 
   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
 
-  val = mpf_get_d (x->value.real);
+  val = mpfr_get_d (x->value.real, GFC_RND_MODE);
   p = gfc_real_kinds[k].digits;
 
   eps = 1.;
@@ -2387,32 +2323,32 @@ gfc_simplify_nearest (gfc_expr * x, gfc_expr * s)
 
   /* TODO we should make sure that 'float' matches kind 4 */
   match_float = gfc_real_kinds[k].kind == 4;
-  if (mpf_cmp_ui (s->value.real, 0) > 0)
+  if (mpfr_cmp_ui (s->value.real, 0) > 0)
     {
       if (match_float)
 	{
 	  rval = (float) val;
 	  rval = rval + eps;
-	  mpf_set_d (result->value.real, rval);
+	  mpfr_set_d (result->value.real, rval, GFC_RND_MODE);
 	}
       else
 	{
 	  val = val + eps;
-	  mpf_set_d (result->value.real, val);
+	  mpfr_set_d (result->value.real, val, GFC_RND_MODE);
 	}
     }
-  else if (mpf_cmp_ui (s->value.real, 0) < 0)
+  else if (mpfr_cmp_ui (s->value.real, 0) < 0)
     {
       if (match_float)
 	{
 	  rval = (float) val;
 	  rval = rval - eps;
-	  mpf_set_d (result->value.real, rval);
+	  mpfr_set_d (result->value.real, rval, GFC_RND_MODE);
 	}
       else
 	{
 	  val = val - eps;
-	  mpf_set_d (result->value.real, val);
+	  mpfr_set_d (result->value.real, val, GFC_RND_MODE);
 	}
     }
   else
@@ -2432,6 +2368,7 @@ simplify_nint (const char *name, gfc_expr * e, gfc_expr * k)
 {
   gfc_expr *rtrunc, *itrunc, *result;
   int kind, cmp;
+  mpfr_t half;
 
   kind = get_kind (BT_INTEGER, k, name, gfc_default_integer_kind ());
   if (kind == -1)
@@ -2445,25 +2382,30 @@ simplify_nint (const char *name, gfc_expr * e, gfc_expr * k)
   rtrunc = gfc_copy_expr (e);
   itrunc = gfc_copy_expr (e);
 
-  cmp = mpf_cmp_ui (e->value.real, 0);
+  cmp = mpfr_cmp_ui (e->value.real, 0);
+
+  gfc_set_model (e->value.real);
+  mpfr_init (half);
+  mpfr_set_str (half, "0.5", 10, GFC_RND_MODE);
 
   if (cmp > 0)
     {
-      mpf_add (rtrunc->value.real, e->value.real, mpf_half);
-      mpf_trunc (itrunc->value.real, rtrunc->value.real);
+      mpfr_add (rtrunc->value.real, e->value.real, half, GFC_RND_MODE);
+      mpfr_trunc (itrunc->value.real, rtrunc->value.real);
     }
   else if (cmp < 0)
     {
-      mpf_sub (rtrunc->value.real, e->value.real, mpf_half);
-      mpf_trunc (itrunc->value.real, rtrunc->value.real);
+      mpfr_sub (rtrunc->value.real, e->value.real, half, GFC_RND_MODE);
+      mpfr_trunc (itrunc->value.real, rtrunc->value.real);
     }
   else
-    mpf_set_ui (itrunc->value.real, 0);
+    mpfr_set_ui (itrunc->value.real, 0, GFC_RND_MODE);
 
-  mpz_set_f (result->value.integer, itrunc->value.real);
+  gfc_mpfr_to_mpz (result->value.integer, itrunc->value.real);
 
   gfc_free_expr (itrunc);
   gfc_free_expr (rtrunc);
+  mpfr_clear (half);
 
   return range_check (result, name);
 }
@@ -2937,7 +2879,7 @@ gfc_expr *
 gfc_simplify_rrspacing (gfc_expr * x)
 {
   gfc_expr *result;
-  mpf_t i2, absv, ln2, lnx, frac, pow2;
+  mpfr_t i2, absv, ln2, lnx, frac, pow2, zero;
   unsigned long exp2;
   int i, p;
 
@@ -2952,41 +2894,48 @@ gfc_simplify_rrspacing (gfc_expr * x)
 
   p = gfc_real_kinds[i].digits;
 
-  if (mpf_cmp (x->value.real, mpf_zero) == 0)
+  gfc_set_model_kind (x->ts.kind);
+  mpfr_init (zero);
+  mpfr_set_ui (zero, 0, GFC_RND_MODE);
+
+  if (mpfr_cmp (x->value.real, zero) == 0)
     {
-      mpf_ui_div (result->value.real, 1, gfc_real_kinds[i].tiny);
+      mpfr_ui_div (result->value.real, 1, gfc_real_kinds[i].tiny, GFC_RND_MODE);
+      mpfr_clear (zero);
       return result;
     }
 
-  mpf_init_set_ui (i2, 2);
-  mpf_init (ln2);
-  mpf_init (absv);
-  mpf_init (lnx);
-  mpf_init (frac);
-  mpf_init (pow2);
+  mpfr_init (i2);
+  mpfr_init (ln2);
+  mpfr_init (absv);
+  mpfr_init (lnx);
+  mpfr_init (frac);
+  mpfr_init (pow2);
 
-  natural_logarithm (&i2, &ln2);
+  mpfr_set_ui (i2, 2, GFC_RND_MODE);
 
-  mpf_abs (absv, x->value.real);
-  natural_logarithm (&absv, &lnx);
+  mpfr_log (ln2, i2, GFC_RND_MODE);
+  mpfr_abs (absv, x->value.real, GFC_RND_MODE);
+  mpfr_log (lnx, absv, GFC_RND_MODE);
 
-  mpf_div (lnx, lnx, ln2);
-  mpf_trunc (lnx, lnx);
-  mpf_add_ui (lnx, lnx, 1);
+  mpfr_div (lnx, lnx, ln2, GFC_RND_MODE);
+  mpfr_trunc (lnx, lnx);
+  mpfr_add_ui (lnx, lnx, 1, GFC_RND_MODE);
 
-  exp2 = (unsigned long) mpf_get_d (lnx);
-  mpf_pow_ui (pow2, i2, exp2);
-  mpf_div (frac, absv, pow2);
+  exp2 = (unsigned long) mpfr_get_d (lnx, GFC_RND_MODE);
+  mpfr_pow_ui (pow2, i2, exp2, GFC_RND_MODE);
+  mpfr_div (frac, absv, pow2, GFC_RND_MODE);
 
   exp2 = (unsigned long) p;
-  mpf_mul_2exp (result->value.real, frac, exp2);
+  mpfr_mul_2exp (result->value.real, frac, exp2, GFC_RND_MODE);
 
-  mpf_clear (i2);
-  mpf_clear (ln2);
-  mpf_clear (absv);
-  mpf_clear (lnx);
-  mpf_clear (frac);
-  mpf_clear (pow2);
+  mpfr_clear (i2);
+  mpfr_clear (ln2);
+  mpfr_clear (absv);
+  mpfr_clear (lnx);
+  mpfr_clear (frac);
+  mpfr_clear (pow2);
+  mpfr_clear (zero);
 
   return range_check (result, "RRSPACING");
 }
@@ -2996,7 +2945,7 @@ gfc_expr *
 gfc_simplify_scale (gfc_expr * x, gfc_expr * i)
 {
   int k, neg_flag, power, exp_range;
-  mpf_t scale, radix;
+  mpfr_t scale, radix;
   gfc_expr *result;
 
   if (x->expr_type != EXPR_CONSTANT || i->expr_type != EXPR_CONSTANT)
@@ -3004,9 +2953,9 @@ gfc_simplify_scale (gfc_expr * x, gfc_expr * i)
 
   result = gfc_constant_result (BT_REAL, x->ts.kind, &x->where);
 
-  if (mpf_sgn (x->value.real) == 0)
+  if (mpfr_sgn (x->value.real) == 0)
     {
-      mpf_set_ui (result->value.real, 0);
+      mpfr_set_ui (result->value.real, 0, GFC_RND_MODE);
       return result;
     }
 
@@ -3035,17 +2984,19 @@ gfc_simplify_scale (gfc_expr * x, gfc_expr * i)
       power = -power;
     }
 
-  mpf_init_set_ui (radix, gfc_real_kinds[k].radix);
-  mpf_init (scale);
-  mpf_pow_ui (scale, radix, power);
+  gfc_set_model_kind (x->ts.kind);
+  mpfr_init (scale);
+  mpfr_init (radix);
+  mpfr_set_ui (radix, gfc_real_kinds[k].radix, GFC_RND_MODE);
+  mpfr_pow_ui (scale, radix, power, GFC_RND_MODE);
 
   if (neg_flag)
-    mpf_div (result->value.real, x->value.real, scale);
+    mpfr_div (result->value.real, x->value.real, scale, GFC_RND_MODE);
   else
-    mpf_mul (result->value.real, x->value.real, scale);
+    mpfr_mul (result->value.real, x->value.real, scale, GFC_RND_MODE);
 
-  mpf_clear (scale);
-  mpf_clear (radix);
+  mpfr_clear (scale);
+  mpfr_clear (radix);
 
   return range_check (result, "SCALE");
 }
@@ -3195,7 +3146,7 @@ gfc_expr *
 gfc_simplify_set_exponent (gfc_expr * x, gfc_expr * i)
 {
   gfc_expr *result;
-  mpf_t i2, ln2, absv, lnx, pow2, frac;
+  mpfr_t i2, ln2, absv, lnx, pow2, frac, zero;
   unsigned long exp2;
 
   if (x->expr_type != EXPR_CONSTANT || i->expr_type != EXPR_CONSTANT)
@@ -3203,44 +3154,51 @@ gfc_simplify_set_exponent (gfc_expr * x, gfc_expr * i)
 
   result = gfc_constant_result (BT_REAL, x->ts.kind, &x->where);
 
-  if (mpf_cmp (x->value.real, mpf_zero) == 0)
+  gfc_set_model_kind (x->ts.kind);
+  mpfr_init (zero);
+  mpfr_set_ui (zero, 0, GFC_RND_MODE);
+
+  if (mpfr_cmp (x->value.real, zero) == 0)
     {
-      mpf_set (result->value.real, mpf_zero);
+      mpfr_set (result->value.real, zero, GFC_RND_MODE);
+      mpfr_clear (zero);
       return result;
     }
 
-  mpf_init_set_ui (i2, 2);
-  mpf_init (ln2);
-  mpf_init (absv);
-  mpf_init (lnx);
-  mpf_init (pow2);
-  mpf_init (frac);
+  mpfr_init (i2);
+  mpfr_init (ln2);
+  mpfr_init (absv);
+  mpfr_init (lnx);
+  mpfr_init (pow2);
+  mpfr_init (frac);
 
-  natural_logarithm (&i2, &ln2);
+  mpfr_set_ui (i2, 2, GFC_RND_MODE);
+  mpfr_log (ln2, i2, GFC_RND_MODE);
 
-  mpf_abs (absv, x->value.real);
-  natural_logarithm (&absv, &lnx);
+  mpfr_abs (absv, x->value.real, GFC_RND_MODE);
+  mpfr_log (lnx, absv, GFC_RND_MODE);
 
-  mpf_div (lnx, lnx, ln2);
-  mpf_trunc (lnx, lnx);
-  mpf_add_ui (lnx, lnx, 1);
+  mpfr_div (lnx, lnx, ln2, GFC_RND_MODE);
+  mpfr_trunc (lnx, lnx);
+  mpfr_add_ui (lnx, lnx, 1, GFC_RND_MODE);
 
   /* Old exponent value, and fraction.  */
-  exp2 = (unsigned long) mpf_get_d (lnx);
-  mpf_pow_ui (pow2, i2, exp2);
+  exp2 = (unsigned long) mpfr_get_d (lnx, GFC_RND_MODE);
+  mpfr_pow_ui (pow2, i2, exp2, GFC_RND_MODE);
 
-  mpf_div (frac, absv, pow2);
+  mpfr_div (frac, absv, pow2, GFC_RND_MODE);
 
   /* New exponent.  */
   exp2 = (unsigned long) mpz_get_d (i->value.integer);
-  mpf_mul_2exp (result->value.real, frac, exp2);
+  mpfr_mul_2exp (result->value.real, frac, exp2, GFC_RND_MODE);
 
-  mpf_clear (i2);
-  mpf_clear (ln2);
-  mpf_clear (absv);
-  mpf_clear (lnx);
-  mpf_clear (pow2);
-  mpf_clear (frac);
+  mpfr_clear (i2);
+  mpfr_clear (ln2);
+  mpfr_clear (absv);
+  mpfr_clear (lnx);
+  mpfr_clear (pow2);
+  mpfr_clear (frac);
+  mpfr_clear (zero);
 
   return range_check (result, "SET_EXPONENT");
 }
@@ -3352,9 +3310,9 @@ gfc_simplify_sign (gfc_expr * x, gfc_expr * y)
     case BT_REAL:
       /* TODO: Handle -0.0 and +0.0 correctly on machines that support
          it.  */
-      mpf_abs (result->value.real, x->value.real);
-      if (mpf_sgn (y->value.integer) < 0)
-	mpf_neg (result->value.real, result->value.real);
+      mpfr_abs (result->value.real, x->value.real, GFC_RND_MODE);
+      if (mpfr_sgn (y->value.real) < 0)
+	mpfr_neg (result->value.real, result->value.real, GFC_RND_MODE);
 
       break;
 
@@ -3370,7 +3328,7 @@ gfc_expr *
 gfc_simplify_sin (gfc_expr * x)
 {
   gfc_expr *result;
-  mpf_t xp, xq;
+  mpfr_t xp, xq;
 
   if (x->expr_type != EXPR_CONSTANT)
     return NULL;
@@ -3380,23 +3338,24 @@ gfc_simplify_sin (gfc_expr * x)
   switch (x->ts.type)
     {
     case BT_REAL:
-      sine (&x->value.real, &result->value.real);
+      mpfr_sin (result->value.real, x->value.real, GFC_RND_MODE);
       break;
 
     case BT_COMPLEX:
-      mpf_init (xp);
-      mpf_init (xq);
+      gfc_set_model (x->value.real);
+      mpfr_init (xp);
+      mpfr_init (xq);
 
-      sine (&x->value.complex.r, &xp);
-      hypercos (&x->value.complex.i, &xq);
-      mpf_mul (result->value.complex.r, xp, xq);
+      mpfr_sin  (xp, x->value.complex.r, GFC_RND_MODE);
+      mpfr_cosh (xq, x->value.complex.i, GFC_RND_MODE);
+      mpfr_mul (result->value.complex.r, xp, xq, GFC_RND_MODE);
 
-      cosine (&x->value.complex.r, &xp);
-      hypersine (&x->value.complex.i, &xq);
-      mpf_mul (result->value.complex.i, xp, xq);
+      mpfr_cos  (xp, x->value.complex.r, GFC_RND_MODE);
+      mpfr_sinh (xq, x->value.complex.i, GFC_RND_MODE);
+      mpfr_mul (result->value.complex.i, xp, xq, GFC_RND_MODE);
 
-      mpf_clear (xp);
-      mpf_clear (xq);
+      mpfr_clear (xp);
+      mpfr_clear (xq);
       break;
 
     default:
@@ -3417,7 +3376,7 @@ gfc_simplify_sinh (gfc_expr * x)
 
   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
 
-  hypersine (&x->value.real, &result->value.real);
+  mpfr_sinh(result->value.real, x->value.real, GFC_RND_MODE);
 
   return range_check (result, "SINH");
 }
@@ -3443,7 +3402,7 @@ gfc_expr *
 gfc_simplify_spacing (gfc_expr * x)
 {
   gfc_expr *result;
-  mpf_t i1, i2, ln2, absv, lnx;
+  mpfr_t i1, i2, ln2, absv, lnx, zero;
   long diff;
   unsigned long exp2;
   int i, p;
@@ -3459,48 +3418,56 @@ gfc_simplify_spacing (gfc_expr * x)
 
   result = gfc_constant_result (BT_REAL, x->ts.kind, &x->where);
 
-  if (mpf_cmp (x->value.real, mpf_zero) == 0)
+  gfc_set_model_kind (x->ts.kind);
+  mpfr_init (zero);
+  mpfr_set_ui (zero, 0, GFC_RND_MODE);
+
+  if (mpfr_cmp (x->value.real, zero) == 0)
     {
-      mpf_set (result->value.real, gfc_real_kinds[i].tiny);
+      mpfr_set (result->value.real, gfc_real_kinds[i].tiny, GFC_RND_MODE);
+      mpfr_clear (zero);
       return result;
     }
 
-  mpf_init_set_ui (i1, 1);
-  mpf_init_set_ui (i2, 2);
-  mpf_init (ln2);
-  mpf_init (absv);
-  mpf_init (lnx);
+  mpfr_init (i1);
+  mpfr_init (i2);
+  mpfr_init (ln2);
+  mpfr_init (absv);
+  mpfr_init (lnx);
 
-  natural_logarithm (&i2, &ln2);
+  mpfr_set_ui (i1, 1, GFC_RND_MODE);
+  mpfr_set_ui (i2, 2, GFC_RND_MODE);
 
-  mpf_abs (absv, x->value.real);
-  natural_logarithm (&absv, &lnx);
+  mpfr_log (ln2, i2, GFC_RND_MODE);
+  mpfr_abs (absv, x->value.real, GFC_RND_MODE);
+  mpfr_log (lnx, absv, GFC_RND_MODE);
 
-  mpf_div (lnx, lnx, ln2);
-  mpf_trunc (lnx, lnx);
-  mpf_add_ui (lnx, lnx, 1);
+  mpfr_div (lnx, lnx, ln2, GFC_RND_MODE);
+  mpfr_trunc (lnx, lnx);
+  mpfr_add_ui (lnx, lnx, 1, GFC_RND_MODE);
 
-  diff = (long) mpf_get_d (lnx) - (long) p;
+  diff = (long) mpfr_get_d (lnx, GFC_RND_MODE) - (long) p;
   if (diff >= 0)
     {
       exp2 = (unsigned) diff;
-      mpf_mul_2exp (result->value.real, i1, exp2);
+      mpfr_mul_2exp (result->value.real, i1, exp2, GFC_RND_MODE);
     }
   else
     {
       diff = -diff;
       exp2 = (unsigned) diff;
-      mpf_div_2exp (result->value.real, i1, exp2);
+      mpfr_div_2exp (result->value.real, i1, exp2, GFC_RND_MODE);
     }
 
-  mpf_clear (i1);
-  mpf_clear (i2);
-  mpf_clear (ln2);
-  mpf_clear (absv);
-  mpf_clear (lnx);
+  mpfr_clear (i1);
+  mpfr_clear (i2);
+  mpfr_clear (ln2);
+  mpfr_clear (absv);
+  mpfr_clear (lnx);
+  mpfr_clear (zero);
 
-  if (mpf_cmp (result->value.real, gfc_real_kinds[i].tiny) < 0)
-    mpf_set (result->value.real, gfc_real_kinds[i].tiny);
+  if (mpfr_cmp (result->value.real, gfc_real_kinds[i].tiny) < 0)
+    mpfr_set (result->value.real, gfc_real_kinds[i].tiny, GFC_RND_MODE);
 
   return range_check (result, "SPACING");
 }
@@ -3510,7 +3477,7 @@ gfc_expr *
 gfc_simplify_sqrt (gfc_expr * e)
 {
   gfc_expr *result;
-  mpf_t ac, ad, s, t, w;
+  mpfr_t ac, ad, s, t, w;
 
   if (e->expr_type != EXPR_CONSTANT)
     return NULL;
@@ -3520,9 +3487,9 @@ gfc_simplify_sqrt (gfc_expr * e)
   switch (e->ts.type)
     {
     case BT_REAL:
-      if (mpf_cmp_si (e->value.real, 0) < 0)
+      if (mpfr_cmp_si (e->value.real, 0) < 0)
 	goto negative_arg;
-      mpf_sqrt (result->value.real, e->value.real);
+      mpfr_sqrt (result->value.real, e->value.real, GFC_RND_MODE);
 
       break;
 
@@ -3530,82 +3497,83 @@ gfc_simplify_sqrt (gfc_expr * e)
       /* Formula taken from Numerical Recipes to avoid over- and
          underflow.  */
 
-      mpf_init (ac);
-      mpf_init (ad);
-      mpf_init (s);
-      mpf_init (t);
-      mpf_init (w);
+      gfc_set_model (e->value.real);
+      mpfr_init (ac);
+      mpfr_init (ad);
+      mpfr_init (s);
+      mpfr_init (t);
+      mpfr_init (w);
 
-      if (mpf_cmp_ui (e->value.complex.r, 0) == 0
-	  && mpf_cmp_ui (e->value.complex.i, 0) == 0)
+      if (mpfr_cmp_ui (e->value.complex.r, 0) == 0
+	  && mpfr_cmp_ui (e->value.complex.i, 0) == 0)
 	{
 
-	  mpf_set_ui (result->value.complex.r, 0);
-	  mpf_set_ui (result->value.complex.i, 0);
+	  mpfr_set_ui (result->value.complex.r, 0, GFC_RND_MODE);
+	  mpfr_set_ui (result->value.complex.i, 0, GFC_RND_MODE);
 	  break;
 	}
 
-      mpf_abs (ac, e->value.complex.r);
-      mpf_abs (ad, e->value.complex.i);
+      mpfr_abs (ac, e->value.complex.r, GFC_RND_MODE);
+      mpfr_abs (ad, e->value.complex.i, GFC_RND_MODE);
 
-      if (mpf_cmp (ac, ad) >= 0)
+      if (mpfr_cmp (ac, ad) >= 0)
 	{
-	  mpf_div (t, e->value.complex.i, e->value.complex.r);
-	  mpf_mul (t, t, t);
-	  mpf_add_ui (t, t, 1);
-	  mpf_sqrt (t, t);
-	  mpf_add_ui (t, t, 1);
-	  mpf_div_ui (t, t, 2);
-	  mpf_sqrt (t, t);
-	  mpf_sqrt (s, ac);
-	  mpf_mul (w, s, t);
+	  mpfr_div (t, e->value.complex.i, e->value.complex.r, GFC_RND_MODE);
+	  mpfr_mul (t, t, t, GFC_RND_MODE);
+	  mpfr_add_ui (t, t, 1, GFC_RND_MODE);
+	  mpfr_sqrt (t, t, GFC_RND_MODE);
+	  mpfr_add_ui (t, t, 1, GFC_RND_MODE);
+	  mpfr_div_ui (t, t, 2, GFC_RND_MODE);
+	  mpfr_sqrt (t, t, GFC_RND_MODE);
+	  mpfr_sqrt (s, ac, GFC_RND_MODE);
+	  mpfr_mul (w, s, t, GFC_RND_MODE);
 	}
       else
 	{
-	  mpf_div (s, e->value.complex.r, e->value.complex.i);
-	  mpf_mul (t, s, s);
-	  mpf_add_ui (t, t, 1);
-	  mpf_sqrt (t, t);
-	  mpf_abs (s, s);
-	  mpf_add (t, t, s);
-	  mpf_div_ui (t, t, 2);
-	  mpf_sqrt (t, t);
-	  mpf_sqrt (s, ad);
-	  mpf_mul (w, s, t);
+	  mpfr_div (s, e->value.complex.r, e->value.complex.i, GFC_RND_MODE);
+	  mpfr_mul (t, s, s, GFC_RND_MODE);
+	  mpfr_add_ui (t, t, 1, GFC_RND_MODE);
+	  mpfr_sqrt (t, t, GFC_RND_MODE);
+	  mpfr_abs (s, s, GFC_RND_MODE);
+	  mpfr_add (t, t, s, GFC_RND_MODE);
+	  mpfr_div_ui (t, t, 2, GFC_RND_MODE);
+	  mpfr_sqrt (t, t, GFC_RND_MODE);
+	  mpfr_sqrt (s, ad, GFC_RND_MODE);
+	  mpfr_mul (w, s, t, GFC_RND_MODE);
 	}
 
-      if (mpf_cmp_ui (w, 0) != 0 && mpf_cmp_ui (e->value.complex.r, 0) >= 0)
+      if (mpfr_cmp_ui (w, 0) != 0 && mpfr_cmp_ui (e->value.complex.r, 0) >= 0)
 	{
-	  mpf_mul_ui (t, w, 2);
-	  mpf_div (result->value.complex.i, e->value.complex.i, t);
-	  mpf_set (result->value.complex.r, w);
+	  mpfr_mul_ui (t, w, 2, GFC_RND_MODE);
+	  mpfr_div (result->value.complex.i, e->value.complex.i, t, GFC_RND_MODE);
+	  mpfr_set (result->value.complex.r, w, GFC_RND_MODE);
 	}
-      else if (mpf_cmp_ui (w, 0) != 0
-	       && mpf_cmp_ui (e->value.complex.r, 0) < 0
-	       && mpf_cmp_ui (e->value.complex.i, 0) >= 0)
+      else if (mpfr_cmp_ui (w, 0) != 0
+	       && mpfr_cmp_ui (e->value.complex.r, 0) < 0
+	       && mpfr_cmp_ui (e->value.complex.i, 0) >= 0)
 	{
-	  mpf_mul_ui (t, w, 2);
-	  mpf_div (result->value.complex.r, e->value.complex.i, t);
-	  mpf_set (result->value.complex.i, w);
+	  mpfr_mul_ui (t, w, 2, GFC_RND_MODE);
+	  mpfr_div (result->value.complex.r, e->value.complex.i, t, GFC_RND_MODE);
+	  mpfr_set (result->value.complex.i, w, GFC_RND_MODE);
 	}
-      else if (mpf_cmp_ui (w, 0) != 0
-	       && mpf_cmp_ui (e->value.complex.r, 0) < 0
-	       && mpf_cmp_ui (e->value.complex.i, 0) < 0)
+      else if (mpfr_cmp_ui (w, 0) != 0
+	       && mpfr_cmp_ui (e->value.complex.r, 0) < 0
+	       && mpfr_cmp_ui (e->value.complex.i, 0) < 0)
 	{
-	  mpf_mul_ui (t, w, 2);
-	  mpf_div (result->value.complex.r, ad, t);
-	  mpf_neg (w, w);
-	  mpf_set (result->value.complex.i, w);
+	  mpfr_mul_ui (t, w, 2, GFC_RND_MODE);
+	  mpfr_div (result->value.complex.r, ad, t, GFC_RND_MODE);
+	  mpfr_neg (w, w, GFC_RND_MODE);
+	  mpfr_set (result->value.complex.i, w, GFC_RND_MODE);
 	}
       else
 	gfc_internal_error ("invalid complex argument of SQRT at %L",
 			    &e->where);
 
-      mpf_clear (s);
-      mpf_clear (t);
-      mpf_clear (ac);
-      mpf_clear (ad);
-      mpf_clear (w);
+      mpfr_clear (s);
+      mpfr_clear (t);
+      mpfr_clear (ac);
+      mpfr_clear (ad);
+      mpfr_clear (w);
 
       break;
 
@@ -3625,9 +3593,8 @@ negative_arg:
 gfc_expr *
 gfc_simplify_tan (gfc_expr * x)
 {
-  gfc_expr *result;
-  mpf_t mpf_sin, mpf_cos, mag_cos;
   int i;
+  gfc_expr *result;
 
   if (x->expr_type != EXPR_CONSTANT)
     return NULL;
@@ -3638,37 +3605,7 @@ gfc_simplify_tan (gfc_expr * x)
 
   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
 
-  mpf_init (mpf_sin);
-  mpf_init (mpf_cos);
-  mpf_init (mag_cos);
-  sine (&x->value.real, &mpf_sin);
-  cosine (&x->value.real, &mpf_cos);
-  mpf_abs (mag_cos, mpf_cos);
-  if (mpf_cmp_ui (mag_cos, 0) == 0)
-    {
-      gfc_error ("Tangent undefined at %L", &x->where);
-      mpf_clear (mpf_sin);
-      mpf_clear (mpf_cos);
-      mpf_clear (mag_cos);
-      gfc_free_expr (result);
-      return &gfc_bad_expr;
-    }
-  else if (mpf_cmp (mag_cos, gfc_real_kinds[i].tiny) < 0)
-    {
-      gfc_error ("Tangent cannot be accurately evaluated at %L", &x->where);
-      mpf_clear (mpf_sin);
-      mpf_clear (mpf_cos);
-      mpf_clear (mag_cos);
-      gfc_free_expr (result);
-      return &gfc_bad_expr;
-    }
-  else
-    {
-      mpf_div (result->value.real, mpf_sin, mpf_cos);
-      mpf_clear (mpf_sin);
-      mpf_clear (mpf_cos);
-      mpf_clear (mag_cos);
-    }
+  mpfr_tan (result->value.real, x->value.real, GFC_RND_MODE);
 
   return range_check (result, "TAN");
 }
@@ -3678,23 +3615,13 @@ gfc_expr *
 gfc_simplify_tanh (gfc_expr * x)
 {
   gfc_expr *result;
-  mpf_t xp, xq;
 
   if (x->expr_type != EXPR_CONSTANT)
     return NULL;
 
   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
 
-  mpf_init (xp);
-  mpf_init (xq);
-
-  hypersine (&x->value.real, &xq);
-  hypercos (&x->value.real, &xp);
-
-  mpf_div (result->value.real, xq, xp);
-
-  mpf_clear (xp);
-  mpf_clear (xq);
+  mpfr_tanh (result->value.real, x->value.real, GFC_RND_MODE);
 
   return range_check (result, "TANH");
 
@@ -3712,7 +3639,7 @@ gfc_simplify_tiny (gfc_expr * e)
     gfc_internal_error ("gfc_simplify_error(): Bad kind");
 
   result = gfc_constant_result (BT_REAL, e->ts.kind, &e->where);
-  mpf_set (result->value.real, gfc_real_kinds[i].tiny);
+  mpfr_set (result->value.real, gfc_real_kinds[i].tiny, GFC_RND_MODE);
 
   return result;
 }
@@ -3988,21 +3915,5 @@ void
 gfc_simplify_init_1 (void)
 {
 
-  mpf_init_set_str (mpf_zero, "0.0", 10);
-  mpf_init_set_str (mpf_half, "0.5", 10);
-  mpf_init_set_str (mpf_one, "1.0", 10);
-  mpz_init_set_str (mpz_zero, "0", 10);
-
   invert_table (ascii_table, xascii_table);
 }
-
-
-void
-gfc_simplify_done_1 (void)
-{
-
-  mpf_clear (mpf_zero);
-  mpf_clear (mpf_half);
-  mpf_clear (mpf_one);
-  mpz_clear (mpz_zero);
-}
diff --git a/gcc/fortran/trans-const.c b/gcc/fortran/trans-const.c
index 8a716de6c6c..3202a599c39 100644
--- a/gcc/fortran/trans-const.c
+++ b/gcc/fortran/trans-const.c
@@ -234,7 +234,7 @@ gfc_conv_mpz_to_tree (mpz_t i, int kind)
 /* Converts a real constant into backend form.  Uses an intermediate string
    representation.  */
 tree
-gfc_conv_mpf_to_tree (mpf_t f, int kind)
+gfc_conv_mpfr_to_tree (mpfr_t f, int kind)
 {
   tree res;
   tree type;
@@ -251,13 +251,9 @@ gfc_conv_mpf_to_tree (mpf_t f, int kind)
     }
   assert (gfc_real_kinds[n].kind);
 
-  assert (gfc_real_kinds[n].radix == 2);
-
   n = MAX (abs (gfc_real_kinds[n].min_exponent),
 	   abs (gfc_real_kinds[n].max_exponent));
-#if 0
-  edigits = 2 + (int) (log (n) / log (gfc_real_kinds[n].radix));
-#endif
+
   edigits = 1;
   while (n > 0)
     {
@@ -265,8 +261,11 @@ gfc_conv_mpf_to_tree (mpf_t f, int kind)
       edigits += 3;
     }
 
+  if (kind == gfc_default_double_kind())
+    p = mpfr_get_str (NULL, &exp, 10, 17, f, GFC_RND_MODE);
+  else
+    p = mpfr_get_str (NULL, &exp, 10, 8, f, GFC_RND_MODE);
 
-  p = mpf_get_str (NULL, &exp, 10, 0, f);
 
   /* We also have one minus sign, "e", "." and a null terminator.  */
   q = (char *) gfc_getmem (strlen (p) + edigits + 4);
@@ -294,6 +293,7 @@ gfc_conv_mpf_to_tree (mpf_t f, int kind)
 
   type = gfc_get_real_type (kind);
   res = build_real (type, REAL_VALUE_ATOF (q, TYPE_MODE (type)));
+
   gfc_free (q);
   gfc_free (p);
 
@@ -321,16 +321,16 @@ gfc_conv_constant_to_tree (gfc_expr * expr)
       return gfc_conv_mpz_to_tree (expr->value.integer, expr->ts.kind);
 
     case BT_REAL:
-      return gfc_conv_mpf_to_tree (expr->value.real, expr->ts.kind);
+      return gfc_conv_mpfr_to_tree (expr->value.real, expr->ts.kind);
 
     case BT_LOGICAL:
       return build_int_2 (expr->value.logical, 0);
 
     case BT_COMPLEX:
       {
-	tree real = gfc_conv_mpf_to_tree (expr->value.complex.r,
+	tree real = gfc_conv_mpfr_to_tree (expr->value.complex.r,
 					  expr->ts.kind);
-	tree imag = gfc_conv_mpf_to_tree (expr->value.complex.i,
+	tree imag = gfc_conv_mpfr_to_tree (expr->value.complex.i,
 					  expr->ts.kind);
 
 	return build_complex (NULL_TREE, real, imag);
diff --git a/gcc/fortran/trans-const.h b/gcc/fortran/trans-const.h
index 97e831346fe..7a36e9a5bcb 100644
--- a/gcc/fortran/trans-const.h
+++ b/gcc/fortran/trans-const.h
@@ -23,7 +23,7 @@ Software Foundation, 59 Temple Place - Suite 330, Boston, MA
 tree gfc_conv_mpz_to_tree (mpz_t, int);
 
 /* Returns a REAL_CST.  */
-tree gfc_conv_mpf_to_tree (mpf_t, int);
+tree gfc_conv_mpfr_to_tree (mpfr_t, int);
 
 /* Build a tree for a constant.  Must be an EXPR_CONSTANT gfc_expr.
    For CHARACTER literal constants, the caller still has to set the
diff --git a/gcc/fortran/trans-intrinsic.c b/gcc/fortran/trans-intrinsic.c
index 71923fa9c2f..f83e1e177b5 100644
--- a/gcc/fortran/trans-intrinsic.c
+++ b/gcc/fortran/trans-intrinsic.c
@@ -33,9 +33,9 @@ Software Foundation, 59 Temple Place - Suite 330, Boston, MA
 #include "real.h"
 #include "tree-gimple.h"
 #include "flags.h"
-#include <gmp.h>
 #include <assert.h>
 #include "gfortran.h"
+#include "arith.h"
 #include "intrinsic.h"
 #include "trans.h"
 #include "trans-const.h"
@@ -308,7 +308,7 @@ gfc_conv_intrinsic_aint (gfc_se * se, gfc_expr * expr, int op)
   tree arg;
   tree tmp;
   tree cond;
-  mpf_t huge;
+  mpfr_t huge;
   int n;
   int kind;
 
@@ -363,14 +363,15 @@ gfc_conv_intrinsic_aint (gfc_se * se, gfc_expr * expr, int op)
   arg = gfc_evaluate_now (arg, &se->pre);
 
   /* Test if the value is too large to handle sensibly.  */
-  mpf_init (huge);
+  gfc_set_model_kind (kind);
+  mpfr_init (huge);
   n = gfc_validate_kind (BT_INTEGER, kind);
-  mpf_set_z (huge, gfc_integer_kinds[n].huge);
-  tmp = gfc_conv_mpf_to_tree (huge, kind);
+  mpfr_set_z (huge, gfc_integer_kinds[n].huge, GFC_RND_MODE);
+  tmp = gfc_conv_mpfr_to_tree (huge, kind);
   cond = build (LT_EXPR, boolean_type_node, arg, tmp);
 
-  mpf_neg (huge, huge);
-  tmp = gfc_conv_mpf_to_tree (huge, kind);
+  mpfr_neg (huge, huge, GFC_RND_MODE);
+  tmp = gfc_conv_mpfr_to_tree (huge, kind);
   tmp = build (GT_EXPR, boolean_type_node, arg, tmp);
   cond = build (TRUTH_AND_EXPR, boolean_type_node, cond, tmp);
   itype = gfc_get_int_type (kind);
@@ -378,6 +379,7 @@ gfc_conv_intrinsic_aint (gfc_se * se, gfc_expr * expr, int op)
   tmp = build_fix_expr (&se->pre, arg, itype, op);
   tmp = convert (type, tmp);
   se->expr = build (COND_EXPR, type, cond, tmp, arg);
+  mpfr_clear (huge);
 }
 
 
@@ -777,7 +779,7 @@ gfc_conv_intrinsic_mod (gfc_se * se, gfc_expr * expr, int modulo)
   tree zero;
   tree test;
   tree test2;
-  mpf_t huge;
+  mpfr_t huge;
   int n;
 
   arg = gfc_conv_intrinsic_function_args (se, expr);
@@ -799,14 +801,15 @@ gfc_conv_intrinsic_mod (gfc_se * se, gfc_expr * expr, int modulo)
 
       tmp = build (RDIV_EXPR, type, arg, arg2);
       /* Test if the value is too large to handle sensibly.  */
-      mpf_init (huge);
+      gfc_set_model_kind (expr->ts.kind);
+      mpfr_init (huge);
       n = gfc_validate_kind (BT_INTEGER, expr->ts.kind);
-      mpf_set_z (huge, gfc_integer_kinds[n].huge);
-      test = gfc_conv_mpf_to_tree (huge, expr->ts.kind);
+      mpfr_set_z (huge, gfc_integer_kinds[n].huge, GFC_RND_MODE);
+      test = gfc_conv_mpfr_to_tree (huge, expr->ts.kind);
       test2 = build (LT_EXPR, boolean_type_node, tmp, test);
 
-      mpf_neg (huge, huge);
-      test = gfc_conv_mpf_to_tree (huge, expr->ts.kind);
+      mpfr_neg (huge, huge, GFC_RND_MODE);
+      test = gfc_conv_mpfr_to_tree (huge, expr->ts.kind);
       test = build (GT_EXPR, boolean_type_node, tmp, test);
       test2 = build (TRUTH_AND_EXPR, boolean_type_node, test, test2);
 
@@ -816,6 +819,7 @@ gfc_conv_intrinsic_mod (gfc_se * se, gfc_expr * expr, int modulo)
       tmp = build (COND_EXPR, type, test2, tmp, arg);
       tmp = build (MULT_EXPR, type, tmp, arg2);
       se->expr = build (MINUS_EXPR, type, arg, tmp);
+      mpfr_clear (huge);
       break;
 
     default:
@@ -1423,7 +1427,7 @@ gfc_conv_intrinsic_minmaxloc (gfc_se * se, gfc_expr * expr, int op)
   switch (arrayexpr->ts.type)
     {
     case BT_REAL:
-      tmp = gfc_conv_mpf_to_tree (gfc_real_kinds[n].huge, arrayexpr->ts.kind);
+      tmp = gfc_conv_mpfr_to_tree (gfc_real_kinds[n].huge, arrayexpr->ts.kind);
       break;
 
     case BT_INTEGER:
@@ -1564,7 +1568,7 @@ gfc_conv_intrinsic_minmaxval (gfc_se * se, gfc_expr * expr, int op)
   switch (expr->ts.type)
     {
     case BT_REAL:
-      tmp = gfc_conv_mpf_to_tree (gfc_real_kinds[n].huge, expr->ts.kind);
+      tmp = gfc_conv_mpfr_to_tree (gfc_real_kinds[n].huge, expr->ts.kind);
       break;
 
     case BT_INTEGER: