Merge tree-ssa-20020619-branch into mainline.

From-SVN: r81764

1 files changed, 2763 insertions, 0 deletions

diff --git a/gcc/fortran/arith.c b/gcc/fortran/arith.c
new file mode 100644
index 0000000..bd03fba
--- /dev/null
+++ b/gcc/fortran/arith.c
@@ -0,0 +1,2763 @@
+/* Compiler arithmetic
+   Copyright (C) 2000, 2001. 2002, 2003 Free Software Foundation, Inc.
+   Contributed by Andy Vaught
+
+This file is part of GNU G95.
+
+GNU G95 is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2, or (at your option)
+any later version.
+
+GNU G95 is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with GNU G95; see the file COPYING.  If not, write to
+the Free Software Foundation, 59 Temple Place - Suite 330,
+Boston, MA 02111-1307, USA.  */
+
+/* Since target arithmetic must be done on the host, there has to
+   be some way of evaluating arithmetic expressions as the host
+   would evaluate them.  We use the GNU MP library to do arithmetic,
+   and this file provides the interface.  */
+
+#include "config.h"
+
+#include <string.h>
+
+#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.
+   The first entry is the default kind, the second entry of the real
+   structure is the default double kind.  */
+
+#define MPZ_NULL {{0,0,0}}
+#define MPF_NULL {{0,0,0,0}}
+
+#define DEF_GFC_INTEGER_KIND(KIND,RADIX,DIGITS,BIT_SIZE)		\
+	{KIND, RADIX, DIGITS, BIT_SIZE, 0, MPZ_NULL, MPZ_NULL, MPZ_NULL}
+
+#define DEF_GFC_LOGICAL_KIND(KIND,BIT_SIZE)				\
+	{KIND, BIT_SIZE}
+
+#define DEF_GFC_REAL_KIND(KIND,RADIX,DIGITS,MIN_EXP, MAX_EXP)		\
+	{KIND, RADIX, DIGITS, MIN_EXP, MAX_EXP,				\
+	 0, 0, MPF_NULL, MPF_NULL, MPF_NULL}
+
+gfc_integer_info gfc_integer_kinds[] = {
+  DEF_GFC_INTEGER_KIND (4, 2, 31, 32),
+  DEF_GFC_INTEGER_KIND (8, 2, 63, 64),
+  DEF_GFC_INTEGER_KIND (2, 2, 15, 16),
+  DEF_GFC_INTEGER_KIND (1, 2,  7,  8),
+  DEF_GFC_INTEGER_KIND (0, 0,  0,  0)
+};
+
+gfc_logical_info gfc_logical_kinds[] = {
+  DEF_GFC_LOGICAL_KIND (4, 32),
+  DEF_GFC_LOGICAL_KIND (8, 64),
+  DEF_GFC_LOGICAL_KIND (2, 16),
+  DEF_GFC_LOGICAL_KIND (1,  8),
+  DEF_GFC_LOGICAL_KIND (0,  0)
+};
+
+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 (0, 0,  0,     0,    0)
+};
+
+
+/* The integer kind to use for array indices.  This will be set to the
+   proper value based on target information from the backend.  */
+
+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 + (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.  */
+
+void
+natural_logarithm (mpf_t * arg, mpf_t * result)
+{
+  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);
+
+      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);
+  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
+   logaritm 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);
+}
+
+/* Calculate exp(arg).
+
+   We use a reduction of the form
+
+     x = Nln2 + r
+
+   Then we obtain exp(r) from the McLaurin series.
+   exp(x) is then recovered from the identity
+
+     exp(x) = 2^N*exp(r).  */
+
+void
+exponential (mpf_t * arg, mpf_t * result)
+{
+  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++)
+	{
+	  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 McLaurin 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);
+    }
+
+  mpf_clear (x);
+}
+
+
+/* Calculate atan(arg).
+
+   Similar to sine but requires special handling for x near 1.  */
+
+void
+arctangent (mpf_t * arg, mpf_t * result)
+{
+  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)
+    {
+      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);
+    }
+  mpf_clear (x);
+}
+
+
+/* Calculate atan2 (y, x)
+
+atan2(y, x) = atan(y/x)				if x > 0,
+	      sign(y)*(pi - atan(|y/x|))	if x < 0,
+	      0					if x = 0 && y == 0,
+	      sign(y)*pi/2			if x = 0 && y != 0.
+*/
+
+void
+arctangent2 (mpf_t * y, mpf_t * x, mpf_t * result)
+{
+  mpf_t t;
+
+  mpf_init (t);
+
+  switch (mpf_sgn (*x))
+    {
+    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);
+      else
+	{
+	  mpf_set (*result, half_pi);
+	  if (mpf_sgn (*y) == -1)
+	    mpf_neg (*result, *result);
+	}
+       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);
+
+  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);
+}
+
+
+/* Given an arithmetic error code, return a pointer to a string that
+   explains the error.  */
+
+static const char *
+gfc_arith_error (arith code)
+{
+  const char *p;
+
+  switch (code)
+    {
+    case ARITH_OK:
+      p = "Arithmetic OK";
+      break;
+    case ARITH_OVERFLOW:
+      p = "Arithmetic overflow";
+      break;
+    case ARITH_UNDERFLOW:
+      p = "Arithmetic underflow";
+      break;
+    case ARITH_DIV0:
+      p = "Division by zero";
+      break;
+    case ARITH_0TO0:
+      p = "Indeterminate form 0 ** 0";
+      break;
+    case ARITH_INCOMMENSURATE:
+      p = "Array operands are incommensurate";
+      break;
+    default:
+      gfc_internal_error ("gfc_arith_error(): Bad error code");
+    }
+
+  return p;
+}
+
+
+/* Get things ready to do math.  */
+
+void
+gfc_arith_init_1 (void)
+{
+  gfc_integer_info *int_info;
+  gfc_real_info *real_info;
+  mpf_t a, b;
+  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);
+
+      mpf_add (pi, pi, b);
+    }
+
+  mpf_mul_ui (two_pi, pi, 2);
+  mpf_div_ui (half_pi, pi, 2);
+
+  /* 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 */
+      mpz_set_ui (r, int_info->radix);
+      mpz_pow_ui (r, r, int_info->digits);
+
+      mpz_init (int_info->huge);
+      mpz_sub_ui (int_info->huge, r, 1);
+
+      /* These are the numbers that are actually representable by the
+         target.  For bases other than two, this needs to be changed.  */
+      if (int_info->radix != 2)
+	gfc_internal_error ("Fix min_int, max_int calculation");
+
+      mpz_init (int_info->min_int);
+      mpz_neg (int_info->min_int, int_info->huge);
+      /* No -1 here, because the representation is symmetric.  */
+
+      mpz_init (int_info->max_int);
+      mpz_add (int_info->max_int, int_info->huge, int_info->huge);
+      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);
+      int_info->range = mpz_get_si (r);
+    }
+
+  /*  mpf_set_default_prec(GFC_REAL_BITS); */
+  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);
+
+      mpf_pow_ui (a, a, real_info->max_exponent);
+      mpf_pow_ui (b, b, real_info->max_exponent - real_info->digits);
+
+      mpf_init (real_info->huge);
+      mpf_sub (real_info->huge, a, b);
+
+      /* Tiny */
+      mpf_set_ui (b, real_info->radix);
+      mpf_pow_ui (b, b, 1 - real_info->min_exponent);
+
+      mpf_init (real_info->tiny);
+      mpf_ui_div (real_info->tiny, 1, b);
+
+      /* Epsilon */
+      mpf_set_ui (b, real_info->radix);
+      mpf_pow_ui (b, b, real_info->digits - 1);
+
+      mpf_init (real_info->epsilon);
+      mpf_ui_div (real_info->epsilon, 1, b);
+
+      /* Range */
+      common_logarithm (&real_info->huge, &a);
+      common_logarithm (&real_info->tiny, &b);
+      mpf_neg (b, b);
+
+      if (mpf_cmp (a, b) > 0)
+	mpf_set (a, b);		/* a = min(a, b) */
+
+      mpf_trunc (a, a);
+      mpz_set_f (r, a);
+      real_info->range = mpz_get_si (r);
+
+      /* Precision */
+      mpf_set_ui (a, real_info->radix);
+      common_logarithm (&a, &a);
+
+      mpf_mul_ui (a, a, real_info->digits - 1);
+      mpf_trunc (a, a);
+      mpz_set_f (r, a);
+      real_info->precision = mpz_get_si (r);
+
+      /* If the radix is an integral power of 10, add one to the
+         precision.  */
+      for (i = 10; i <= real_info->radix; i *= 10)
+	if (i == real_info->radix)
+	  real_info->precision++;
+    }
+
+  mpz_clear (r);
+  mpf_clear (a);
+  mpf_clear (b);
+}
+
+
+/* Clean up, get rid of numeric constants.  */
+
+void
+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);
+      mpz_clear (ip->max_int);
+      mpz_clear (ip->huge);
+    }
+
+  for (rp = gfc_real_kinds; rp->kind; rp++)
+    {
+      mpf_clear (rp->epsilon);
+      mpf_clear (rp->huge);
+      mpf_clear (rp->tiny);
+    }
+}
+
+
+/* Return default kinds.  */
+
+int
+gfc_default_integer_kind (void)
+{
+  return gfc_integer_kinds[gfc_option.i8 ? 1 : 0].kind;
+}
+
+int
+gfc_default_real_kind (void)
+{
+  return gfc_real_kinds[gfc_option.r8 ? 1 : 0].kind;
+}
+
+int
+gfc_default_double_kind (void)
+{
+  return gfc_real_kinds[1].kind;
+}
+
+int
+gfc_default_character_kind (void)
+{
+  return 1;
+}
+
+int
+gfc_default_logical_kind (void)
+{
+  return gfc_logical_kinds[gfc_option.i8 ? 1 : 0].kind;
+}
+
+int
+gfc_default_complex_kind (void)
+{
+  return gfc_default_real_kind ();
+}
+
+
+/* Make sure that a valid kind is present.  Returns an index into the
+   gfc_integer_kinds array, -1 if the kind is not present.  */
+
+static int
+validate_integer (int kind)
+{
+  int i;
+
+  for (i = 0;; i++)
+    {
+      if (gfc_integer_kinds[i].kind == 0)
+	{
+	  i = -1;
+	  break;
+	}
+      if (gfc_integer_kinds[i].kind == kind)
+	break;
+    }
+
+  return i;
+}
+
+
+static int
+validate_real (int kind)
+{
+  int i;
+
+  for (i = 0;; i++)
+    {
+      if (gfc_real_kinds[i].kind == 0)
+	{
+	  i = -1;
+	  break;
+	}
+      if (gfc_real_kinds[i].kind == kind)
+	break;
+    }
+
+  return i;
+}
+
+
+static int
+validate_logical (int kind)
+{
+  int i;
+
+  for (i = 0;; i++)
+    {
+      if (gfc_logical_kinds[i].kind == 0)
+	{
+	  i = -1;
+	  break;
+	}
+      if (gfc_logical_kinds[i].kind == kind)
+	break;
+    }
+
+  return i;
+}
+
+
+static int
+validate_character (int kind)
+{
+
+  if (kind == gfc_default_character_kind ())
+    return 0;
+  return -1;
+}
+
+
+/* Validate a kind given a basic type.  The return value is the same
+   for the child functions, with -1 indicating nonexistence of the
+   type.  */
+
+int
+gfc_validate_kind (bt type, int kind)
+{
+  int rc;
+
+  switch (type)
+    {
+    case BT_REAL:		/* Fall through */
+    case BT_COMPLEX:
+      rc = validate_real (kind);
+      break;
+    case BT_INTEGER:
+      rc = validate_integer (kind);
+      break;
+    case BT_LOGICAL:
+      rc = validate_logical (kind);
+      break;
+    case BT_CHARACTER:
+      rc = validate_character (kind);
+      break;
+
+    default:
+      gfc_internal_error ("gfc_validate_kind(): Got bad type");
+    }
+
+  return rc;
+}
+
+
+/* Given an integer and a kind, make sure that the integer lies within
+   the range of the kind.  Returns ARITH_OK or ARITH_OVERFLOW.  */
+
+static arith
+gfc_check_integer_range (mpz_t p, int kind)
+{
+  arith result;
+  int i;
+
+  i = validate_integer (kind);
+  if (i == -1)
+    gfc_internal_error ("gfc_check_integer_range(): Bad kind");
+
+  result = ARITH_OK;
+
+  if (mpz_cmp (p, gfc_integer_kinds[i].min_int) < 0
+      || mpz_cmp (p, gfc_integer_kinds[i].max_int) > 0)
+    result = ARITH_OVERFLOW;
+
+  return result;
+}
+
+
+/* Given a real and a kind, make sure that the real lies within the
+   range of the kind.  Returns ARITH_OK, ARITH_OVERFLOW or
+   ARITH_UNDERFLOW.  */
+
+static arith
+gfc_check_real_range (mpf_t p, int kind)
+{
+  arith retval;
+  mpf_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");
+
+  retval = ARITH_OK;
+  if (mpf_sgn (q) == 0)
+    goto done;
+
+  if (mpf_cmp (q, gfc_real_kinds[i].huge) == 1)
+    {
+      retval = ARITH_OVERFLOW;
+      goto done;
+    }
+
+  if (mpf_cmp (q, gfc_real_kinds[i].tiny) == -1)
+    retval = ARITH_UNDERFLOW;
+
+done:
+  mpf_clear (q);
+
+  return retval;
+}
+
+
+/* Function to return a constant expression node of a given type and
+   kind.  */
+
+gfc_expr *
+gfc_constant_result (bt type, int kind, locus * where)
+{
+  gfc_expr *result;
+
+  if (!where)
+    gfc_internal_error
+      ("gfc_constant_result(): locus 'where' cannot be NULL");
+
+  result = gfc_get_expr ();
+
+  result->expr_type = EXPR_CONSTANT;
+  result->ts.type = type;
+  result->ts.kind = kind;
+  result->where = *where;
+
+  switch (type)
+    {
+    case BT_INTEGER:
+      mpz_init (result->value.integer);
+      break;
+
+    case BT_REAL:
+      mpf_init (result->value.real);
+      break;
+
+    case BT_COMPLEX:
+      mpf_init (result->value.complex.r);
+      mpf_init (result->value.complex.i);
+      break;
+
+    default:
+      break;
+    }
+
+  return result;
+}
+
+
+/* Low-level arithmetic functions.  All of these subroutines assume
+   that all operands are of the same type and return an operand of the
+   same type.  The other thing about these subroutines is that they
+   can fail in various ways -- overflow, underflow, division by zero,
+   zero raised to the zero, etc.  */
+
+static arith
+gfc_arith_not (gfc_expr * op1, gfc_expr ** resultp)
+{
+  gfc_expr *result;
+
+  result = gfc_constant_result (BT_LOGICAL, op1->ts.kind, &op1->where);
+  result->value.logical = !op1->value.logical;
+  *resultp = result;
+
+  return ARITH_OK;
+}
+
+
+static arith
+gfc_arith_and (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
+{
+  gfc_expr *result;
+
+  result = gfc_constant_result (BT_LOGICAL, gfc_kind_max (op1, op2),
+				&op1->where);
+  result->value.logical = op1->value.logical && op2->value.logical;
+  *resultp = result;
+
+  return ARITH_OK;
+}
+
+
+static arith
+gfc_arith_or (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
+{
+  gfc_expr *result;
+
+  result = gfc_constant_result (BT_LOGICAL, gfc_kind_max (op1, op2),
+				&op1->where);
+  result->value.logical = op1->value.logical || op2->value.logical;
+  *resultp = result;
+
+  return ARITH_OK;
+}
+
+
+static arith
+gfc_arith_eqv (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
+{
+  gfc_expr *result;
+
+  result = gfc_constant_result (BT_LOGICAL, gfc_kind_max (op1, op2),
+				&op1->where);
+  result->value.logical = op1->value.logical == op2->value.logical;
+  *resultp = result;
+
+  return ARITH_OK;
+}
+
+
+static arith
+gfc_arith_neqv (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
+{
+  gfc_expr *result;
+
+  result = gfc_constant_result (BT_LOGICAL, gfc_kind_max (op1, op2),
+				&op1->where);
+  result->value.logical = op1->value.logical != op2->value.logical;
+  *resultp = result;
+
+  return ARITH_OK;
+}
+
+
+/* Make sure a constant numeric expression is within the range for
+   it's type and kind.  Note that there's also a gfc_check_range(),
+   but that one deals with the intrinsic RANGE function.  */
+
+arith
+gfc_range_check (gfc_expr * e)
+{
+  arith rc;
+
+  switch (e->ts.type)
+    {
+    case BT_INTEGER:
+      rc = gfc_check_integer_range (e->value.integer, e->ts.kind);
+      break;
+
+    case BT_REAL:
+      rc = gfc_check_real_range (e->value.real, e->ts.kind);
+      break;
+
+    case BT_COMPLEX:
+      rc = gfc_check_real_range (e->value.complex.r, e->ts.kind);
+      if (rc != ARITH_OK)
+	rc = gfc_check_real_range (e->value.complex.i, e->ts.kind);
+
+      break;
+
+    default:
+      gfc_internal_error ("gfc_range_check(): Bad type");
+    }
+
+  return rc;
+}
+
+
+/* It may seem silly to have a subroutine that actually computes the
+   unary plus of a constant, but it prevents us from making exceptions
+   in the code elsewhere.  */
+
+static arith
+gfc_arith_uplus (gfc_expr * op1, gfc_expr ** resultp)
+{
+
+  *resultp = gfc_copy_expr (op1);
+  return ARITH_OK;
+}
+
+
+static arith
+gfc_arith_uminus (gfc_expr * op1, gfc_expr ** resultp)
+{
+  gfc_expr *result;
+  arith rc;
+
+  result = gfc_constant_result (op1->ts.type, op1->ts.kind, &op1->where);
+
+  switch (op1->ts.type)
+    {
+    case BT_INTEGER:
+      mpz_neg (result->value.integer, op1->value.integer);
+      break;
+
+    case BT_REAL:
+      mpf_neg (result->value.real, op1->value.real);
+      break;
+
+    case BT_COMPLEX:
+      mpf_neg (result->value.complex.r, op1->value.complex.r);
+      mpf_neg (result->value.complex.i, op1->value.complex.i);
+      break;
+
+    default:
+      gfc_internal_error ("gfc_arith_uminus(): Bad basic type");
+    }
+
+  rc = gfc_range_check (result);
+
+  if (rc != ARITH_OK)
+    gfc_free_expr (result);
+  else
+    *resultp = result;
+
+  return rc;
+}
+
+
+static arith
+gfc_arith_plus (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
+{
+  gfc_expr *result;
+  arith rc;
+
+  result = gfc_constant_result (op1->ts.type, op1->ts.kind, &op1->where);
+
+  switch (op1->ts.type)
+    {
+    case BT_INTEGER:
+      mpz_add (result->value.integer, op1->value.integer, op2->value.integer);
+      break;
+
+    case BT_REAL:
+      mpf_add (result->value.real, op1->value.real, op2->value.real);
+      break;
+
+    case BT_COMPLEX:
+      mpf_add (result->value.complex.r, op1->value.complex.r,
+	       op2->value.complex.r);
+
+      mpf_add (result->value.complex.i, op1->value.complex.i,
+	       op2->value.complex.i);
+      break;
+
+    default:
+      gfc_internal_error ("gfc_arith_plus(): Bad basic type");
+    }
+
+  rc = gfc_range_check (result);
+
+  if (rc != ARITH_OK)
+    gfc_free_expr (result);
+  else
+    *resultp = result;
+
+  return rc;
+}
+
+
+static arith
+gfc_arith_minus (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
+{
+  gfc_expr *result;
+  arith rc;
+
+  result = gfc_constant_result (op1->ts.type, op1->ts.kind, &op1->where);
+
+  switch (op1->ts.type)
+    {
+    case BT_INTEGER:
+      mpz_sub (result->value.integer, op1->value.integer, op2->value.integer);
+      break;
+
+    case BT_REAL:
+      mpf_sub (result->value.real, op1->value.real, op2->value.real);
+      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);
+
+      break;
+
+    default:
+      gfc_internal_error ("gfc_arith_minus(): Bad basic type");
+    }
+
+  rc = gfc_range_check (result);
+
+  if (rc != ARITH_OK)
+    gfc_free_expr (result);
+  else
+    *resultp = result;
+
+  return rc;
+}
+
+
+static arith
+gfc_arith_times (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
+{
+  gfc_expr *result;
+  mpf_t x, y;
+  arith rc;
+
+  result = gfc_constant_result (op1->ts.type, op1->ts.kind, &op1->where);
+
+  switch (op1->ts.type)
+    {
+    case BT_INTEGER:
+      mpz_mul (result->value.integer, op1->value.integer, op2->value.integer);
+      break;
+
+    case BT_REAL:
+      mpf_mul (result->value.real, op1->value.real, op2->value.real);
+      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);
+
+      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);
+
+      mpf_clear (x);
+      mpf_clear (y);
+
+      break;
+
+    default:
+      gfc_internal_error ("gfc_arith_times(): Bad basic type");
+    }
+
+  rc = gfc_range_check (result);
+
+  if (rc != ARITH_OK)
+    gfc_free_expr (result);
+  else
+    *resultp = result;
+
+  return rc;
+}
+
+
+static arith
+gfc_arith_divide (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
+{
+  gfc_expr *result;
+  mpf_t x, y, div;
+  arith rc;
+
+  rc = ARITH_OK;
+
+  result = gfc_constant_result (op1->ts.type, op1->ts.kind, &op1->where);
+
+  switch (op1->ts.type)
+    {
+    case BT_INTEGER:
+      if (mpz_sgn (op2->value.integer) == 0)
+	{
+	  rc = ARITH_DIV0;
+	  break;
+	}
+
+      mpz_tdiv_q (result->value.integer, op1->value.integer,
+		  op2->value.integer);
+      break;
+
+    case BT_REAL:
+      if (mpf_sgn (op2->value.real) == 0)
+	{
+	  rc = ARITH_DIV0;
+	  break;
+	}
+
+      mpf_div (result->value.real, op1->value.real, op2->value.real);
+      break;
+
+    case BT_COMPLEX:
+      if (mpf_sgn (op2->value.complex.r) == 0
+	  && mpf_sgn (op2->value.complex.i) == 0)
+	{
+	  rc = ARITH_DIV0;
+	  break;
+	}
+
+      mpf_init (x);
+      mpf_init (y);
+      mpf_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);
+
+      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);
+
+      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);
+
+      mpf_clear (x);
+      mpf_clear (y);
+      mpf_clear (div);
+
+      break;
+
+    default:
+      gfc_internal_error ("gfc_arith_divide(): Bad basic type");
+    }
+
+  if (rc == ARITH_OK)
+    rc = gfc_range_check (result);
+
+  if (rc != ARITH_OK)
+    gfc_free_expr (result);
+  else
+    *resultp = result;
+
+  return rc;
+}
+
+
+/* Compute the reciprocal of a complex number (guaranteed nonzero).  */
+
+static void
+complex_reciprocal (gfc_expr * op)
+{
+  mpf_t mod, a, result_r, result_i;
+
+  mpf_init (mod);
+  mpf_init (a);
+
+  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);
+
+  mpf_init (result_r);
+  mpf_div (result_r, op->value.complex.r, mod);
+
+  mpf_init (result_i);
+  mpf_neg (result_i, op->value.complex.i);
+  mpf_div (result_i, result_i, mod);
+
+  mpf_set (op->value.complex.r, result_r);
+  mpf_set (op->value.complex.i, result_i);
+
+  mpf_clear (result_r);
+  mpf_clear (result_i);
+
+  mpf_clear (mod);
+  mpf_clear (a);
+}
+
+
+/* Raise a complex number to positive power.  */
+
+static void
+complex_pow_ui (gfc_expr * base, int power, gfc_expr * result)
+{
+  mpf_t temp_r, temp_i, a;
+
+  mpf_set_ui (result->value.complex.r, 1);
+  mpf_set_ui (result->value.complex.i, 0);
+
+  mpf_init (temp_r);
+  mpf_init (temp_i);
+  mpf_init (a);
+
+  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);
+
+      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);
+
+      mpf_set (result->value.complex.r, temp_r);
+      mpf_set (result->value.complex.i, temp_i);
+    }
+
+  mpf_clear (temp_r);
+  mpf_clear (temp_i);
+  mpf_clear (a);
+}
+
+
+/* Raise a number to an integer power.  */
+
+static arith
+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;
+  arith rc;
+
+  rc = ARITH_OK;
+
+  if (gfc_extract_int (op2, &power) != NULL)
+    gfc_internal_error ("gfc_arith_power(): Bad exponent");
+
+  result = gfc_constant_result (op1->ts.type, op1->ts.kind, &op1->where);
+
+  if (power == 0)
+    {				/* Handle something to the zeroth power */
+      switch (op1->ts.type)
+	{
+	case BT_INTEGER:
+	  if (mpz_sgn (op1->value.integer) == 0)
+	    rc = ARITH_0TO0;
+	  else
+	    mpz_set_ui (result->value.integer, 1);
+
+	  break;
+
+	case BT_REAL:
+	  if (mpf_sgn (op1->value.real) == 0)
+	    rc = ARITH_0TO0;
+	  else
+	    mpf_set_ui (result->value.real, 1);
+
+	  break;
+
+	case BT_COMPLEX:
+	  if (mpf_sgn (op1->value.complex.r) == 0
+	      && mpf_sgn (op1->value.complex.i) == 0)
+	    rc = ARITH_0TO0;
+	  else
+	    {
+	      mpf_set_ui (result->value.complex.r, 1);
+	      mpf_set_ui (result->value.complex.r, 0);
+	    }
+
+	  break;
+
+	default:
+	  gfc_internal_error ("gfc_arith_power(): Bad base");
+	}
+    }
+
+  if (power != 0)
+    {
+      apower = power;
+      if (power < 0)
+	apower = -power;
+
+      switch (op1->ts.type)
+	{
+	case BT_INTEGER:
+	  mpz_pow_ui (result->value.integer, op1->value.integer, apower);
+
+	  if (power < 0)
+	    {
+	      mpz_init_set_ui (unity_z, 1);
+	      mpz_tdiv_q (result->value.integer, unity_z,
+			  result->value.integer);
+	      mpz_clear (unity_z);
+	    }
+
+	  break;
+
+	case BT_REAL:
+	  mpf_pow_ui (result->value.real, op1->value.real, apower);
+
+	  if (power < 0)
+	    {
+	      mpf_init_set_ui (unity_f, 1);
+	      mpf_div (result->value.real, unity_f, result->value.real);
+	      mpf_clear (unity_f);
+	    }
+
+	  break;
+
+	case BT_COMPLEX:
+	  complex_pow_ui (op1, apower, result);
+	  if (power < 0)
+	    complex_reciprocal (result);
+
+	  break;
+
+	default:
+	  break;
+	}
+    }
+
+  if (rc == ARITH_OK)
+    rc = gfc_range_check (result);
+
+  if (rc != ARITH_OK)
+    gfc_free_expr (result);
+  else
+    *resultp = result;
+
+  return rc;
+}
+
+
+/* Concatenate two string constants.  */
+
+static arith
+gfc_arith_concat (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
+{
+  gfc_expr *result;
+  int len;
+
+  result = gfc_constant_result (BT_CHARACTER, gfc_default_character_kind (),
+				&op1->where);
+
+  len = op1->value.character.length + op2->value.character.length;
+
+  result->value.character.string = gfc_getmem (len + 1);
+  result->value.character.length = len;
+
+  memcpy (result->value.character.string, op1->value.character.string,
+	  op1->value.character.length);
+
+  memcpy (result->value.character.string + op1->value.character.length,
+	  op2->value.character.string, op2->value.character.length);
+
+  result->value.character.string[len] = '\0';
+
+  *resultp = result;
+
+  return ARITH_OK;
+}
+
+
+/* Comparison operators.  Assumes that the two expression nodes
+   contain two constants of the same type.  */
+
+int
+gfc_compare_expr (gfc_expr * op1, gfc_expr * op2)
+{
+  int rc;
+
+  switch (op1->ts.type)
+    {
+    case BT_INTEGER:
+      rc = mpz_cmp (op1->value.integer, op2->value.integer);
+      break;
+
+    case BT_REAL:
+      rc = mpf_cmp (op1->value.real, op2->value.real);
+      break;
+
+    case BT_CHARACTER:
+      rc = gfc_compare_string (op1, op2, NULL);
+      break;
+
+    case BT_LOGICAL:
+      rc = ((!op1->value.logical && op2->value.logical)
+	    || (op1->value.logical && !op2->value.logical));
+      break;
+
+    default:
+      gfc_internal_error ("gfc_compare_expr(): Bad basic type");
+    }
+
+  return rc;
+}
+
+
+/* Compare a pair of complex numbers.  Naturally, this is only for
+   equality/nonequality.  */
+
+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);
+}
+
+
+/* Given two constant strings and the inverse collating sequence,
+   compare the strings.  We return -1 for a<b, 0 for a==b and 1 for
+   a>b.  If the xcoll_table is NULL, we use the processor's default
+   collating sequence.  */
+
+int
+gfc_compare_string (gfc_expr * a, gfc_expr * b, const int *xcoll_table)
+{
+  int len, alen, blen, i, ac, bc;
+
+  alen = a->value.character.length;
+  blen = b->value.character.length;
+
+  len = (alen > blen) ? alen : blen;
+
+  for (i = 0; i < len; i++)
+    {
+      ac = (i < alen) ? a->value.character.string[i] : ' ';
+      bc = (i < blen) ? b->value.character.string[i] : ' ';
+
+      if (xcoll_table != NULL)
+	{
+	  ac = xcoll_table[ac];
+	  bc = xcoll_table[bc];
+	}
+
+      if (ac < bc)
+	return -1;
+      if (ac > bc)
+	return 1;
+    }
+
+  /* Strings are equal */
+
+  return 0;
+}
+
+
+/* Specific comparison subroutines.  */
+
+static arith
+gfc_arith_eq (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
+{
+  gfc_expr *result;
+
+  result = gfc_constant_result (BT_LOGICAL, gfc_default_logical_kind (),
+				&op1->where);
+  result->value.logical = (op1->ts.type == BT_COMPLEX) ?
+    compare_complex (op1, op2) : (gfc_compare_expr (op1, op2) == 0);
+
+  *resultp = result;
+  return ARITH_OK;
+}
+
+
+static arith
+gfc_arith_ne (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
+{
+  gfc_expr *result;
+
+  result = gfc_constant_result (BT_LOGICAL, gfc_default_logical_kind (),
+				&op1->where);
+  result->value.logical = (op1->ts.type == BT_COMPLEX) ?
+    !compare_complex (op1, op2) : (gfc_compare_expr (op1, op2) != 0);
+
+  *resultp = result;
+  return ARITH_OK;
+}
+
+
+static arith
+gfc_arith_gt (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
+{
+  gfc_expr *result;
+
+  result = gfc_constant_result (BT_LOGICAL, gfc_default_logical_kind (),
+				&op1->where);
+  result->value.logical = (gfc_compare_expr (op1, op2) > 0);
+  *resultp = result;
+
+  return ARITH_OK;
+}
+
+
+static arith
+gfc_arith_ge (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
+{
+  gfc_expr *result;
+
+  result = gfc_constant_result (BT_LOGICAL, gfc_default_logical_kind (),
+				&op1->where);
+  result->value.logical = (gfc_compare_expr (op1, op2) >= 0);
+  *resultp = result;
+
+  return ARITH_OK;
+}
+
+
+static arith
+gfc_arith_lt (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
+{
+  gfc_expr *result;
+
+  result = gfc_constant_result (BT_LOGICAL, gfc_default_logical_kind (),
+				&op1->where);
+  result->value.logical = (gfc_compare_expr (op1, op2) < 0);
+  *resultp = result;
+
+  return ARITH_OK;
+}
+
+
+static arith
+gfc_arith_le (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
+{
+  gfc_expr *result;
+
+  result = gfc_constant_result (BT_LOGICAL, gfc_default_logical_kind (),
+				&op1->where);
+  result->value.logical = (gfc_compare_expr (op1, op2) <= 0);
+  *resultp = result;
+
+  return ARITH_OK;
+}
+
+
+static arith
+reduce_unary (arith (*eval) (gfc_expr *, gfc_expr **), gfc_expr * op,
+	      gfc_expr ** result)
+{
+  gfc_constructor *c, *head;
+  gfc_expr *r;
+  arith rc;
+
+  if (op->expr_type == EXPR_CONSTANT)
+    return eval (op, result);
+
+  rc = ARITH_OK;
+  head = gfc_copy_constructor (op->value.constructor);
+
+  for (c = head; c; c = c->next)
+    {
+      rc = eval (c->expr, &r);
+      if (rc != ARITH_OK)
+	break;
+
+      gfc_replace_expr (c->expr, r);
+    }
+
+  if (rc != ARITH_OK)
+    gfc_free_constructor (head);
+  else
+    {
+      r = gfc_get_expr ();
+      r->expr_type = EXPR_ARRAY;
+      r->value.constructor = head;
+      r->shape = gfc_copy_shape (op->shape, op->rank);
+
+      r->ts = head->expr->ts;
+      r->where = op->where;
+      r->rank = op->rank;
+
+      *result = r;
+    }
+
+  return rc;
+}
+
+
+static arith
+reduce_binary_ac (arith (*eval) (gfc_expr *, gfc_expr *, gfc_expr **),
+		  gfc_expr * op1, gfc_expr * op2,
+		  gfc_expr ** result)
+{
+  gfc_constructor *c, *head;
+  gfc_expr *r;
+  arith rc;
+
+  head = gfc_copy_constructor (op1->value.constructor);
+  rc = ARITH_OK;
+
+  for (c = head; c; c = c->next)
+    {
+      rc = eval (c->expr, op2, &r);
+      if (rc != ARITH_OK)
+	break;
+
+      gfc_replace_expr (c->expr, r);
+    }
+
+  if (rc != ARITH_OK)
+    gfc_free_constructor (head);
+  else
+    {
+      r = gfc_get_expr ();
+      r->expr_type = EXPR_ARRAY;
+      r->value.constructor = head;
+      r->shape = gfc_copy_shape (op1->shape, op1->rank);
+
+      r->ts = head->expr->ts;
+      r->where = op1->where;
+      r->rank = op1->rank;
+
+      *result = r;
+    }
+
+  return rc;
+}
+
+
+static arith
+reduce_binary_ca (arith (*eval) (gfc_expr *, gfc_expr *, gfc_expr **),
+		  gfc_expr * op1, gfc_expr * op2,
+		  gfc_expr ** result)
+{
+  gfc_constructor *c, *head;
+  gfc_expr *r;
+  arith rc;
+
+  head = gfc_copy_constructor (op2->value.constructor);
+  rc = ARITH_OK;
+
+  for (c = head; c; c = c->next)
+    {
+      rc = eval (op1, c->expr, &r);
+      if (rc != ARITH_OK)
+	break;
+
+      gfc_replace_expr (c->expr, r);
+    }
+
+  if (rc != ARITH_OK)
+    gfc_free_constructor (head);
+  else
+    {
+      r = gfc_get_expr ();
+      r->expr_type = EXPR_ARRAY;
+      r->value.constructor = head;
+      r->shape = gfc_copy_shape (op2->shape, op2->rank);
+
+      r->ts = head->expr->ts;
+      r->where = op2->where;
+      r->rank = op2->rank;
+
+      *result = r;
+    }
+
+  return rc;
+}
+
+
+static arith
+reduce_binary_aa (arith (*eval) (gfc_expr *, gfc_expr *, gfc_expr **),
+		  gfc_expr * op1, gfc_expr * op2,
+		  gfc_expr ** result)
+{
+  gfc_constructor *c, *d, *head;
+  gfc_expr *r;
+  arith rc;
+
+  head = gfc_copy_constructor (op1->value.constructor);
+
+  rc = ARITH_OK;
+  d = op2->value.constructor;
+
+  if (gfc_check_conformance ("Elemental binary operation", op1, op2)
+      != SUCCESS)
+    rc = ARITH_INCOMMENSURATE;
+  else
+    {
+
+      for (c = head; c; c = c->next, d = d->next)
+	{
+	  if (d == NULL)
+	    {
+	      rc = ARITH_INCOMMENSURATE;
+	      break;
+	    }
+
+	  rc = eval (c->expr, d->expr, &r);
+	  if (rc != ARITH_OK)
+	    break;
+
+	  gfc_replace_expr (c->expr, r);
+	}
+
+      if (d != NULL)
+	rc = ARITH_INCOMMENSURATE;
+    }
+
+  if (rc != ARITH_OK)
+    gfc_free_constructor (head);
+  else
+    {
+      r = gfc_get_expr ();
+      r->expr_type = EXPR_ARRAY;
+      r->value.constructor = head;
+      r->shape = gfc_copy_shape (op1->shape, op1->rank);
+
+      r->ts = head->expr->ts;
+      r->where = op1->where;
+      r->rank = op1->rank;
+
+      *result = r;
+    }
+
+  return rc;
+}
+
+
+static arith
+reduce_binary (arith (*eval) (gfc_expr *, gfc_expr *, gfc_expr **),
+	       gfc_expr * op1, gfc_expr * op2,
+	       gfc_expr ** result)
+{
+
+  if (op1->expr_type == EXPR_CONSTANT && op2->expr_type == EXPR_CONSTANT)
+    return eval (op1, op2, result);
+
+  if (op1->expr_type == EXPR_CONSTANT && op2->expr_type == EXPR_ARRAY)
+    return reduce_binary_ca (eval, op1, op2, result);
+
+  if (op1->expr_type == EXPR_ARRAY && op2->expr_type == EXPR_CONSTANT)
+    return reduce_binary_ac (eval, op1, op2, result);
+
+  return reduce_binary_aa (eval, op1, op2, result);
+}
+
+
+typedef union
+{
+  arith (*f2)(gfc_expr *, gfc_expr **);
+  arith (*f3)(gfc_expr *, gfc_expr *, gfc_expr **);
+}
+eval_f;
+
+/* High level arithmetic subroutines.  These subroutines go into
+   eval_intrinsic(), which can do one of several things to its
+   operands.  If the operands are incompatible with the intrinsic
+   operation, we return a node pointing to the operands and hope that
+   an operator interface is found during resolution.
+
+   If the operands are compatible and are constants, then we try doing
+   the arithmetic.  We also handle the cases where either or both
+   operands are array constructors.  */
+
+static gfc_expr *
+eval_intrinsic (gfc_intrinsic_op operator,
+		eval_f eval, gfc_expr * op1, gfc_expr * op2)
+{
+  gfc_expr temp, *result;
+  int unary;
+  arith rc;
+
+  gfc_clear_ts (&temp.ts);
+
+  switch (operator)
+    {
+    case INTRINSIC_NOT:	/* Logical unary */
+      if (op1->ts.type != BT_LOGICAL)
+	goto runtime;
+
+      temp.ts.type = BT_LOGICAL;
+      temp.ts.kind = gfc_default_logical_kind ();
+
+      unary = 1;
+      break;
+
+      /* Logical binary operators */
+    case INTRINSIC_OR:
+    case INTRINSIC_AND:
+    case INTRINSIC_NEQV:
+    case INTRINSIC_EQV:
+      if (op1->ts.type != BT_LOGICAL || op2->ts.type != BT_LOGICAL)
+	goto runtime;
+
+      temp.ts.type = BT_LOGICAL;
+      temp.ts.kind = gfc_default_logical_kind ();
+
+      unary = 0;
+      break;
+
+    case INTRINSIC_UPLUS:
+    case INTRINSIC_UMINUS:	/* Numeric unary */
+      if (!gfc_numeric_ts (&op1->ts))
+	goto runtime;
+
+      temp.ts = op1->ts;
+
+      unary = 1;
+      break;
+
+    case INTRINSIC_GE:
+    case INTRINSIC_LT:		/* Additional restrictions  */
+    case INTRINSIC_LE:          /* for ordering relations.  */
+    case INTRINSIC_GT:
+      if (op1->ts.type == BT_COMPLEX || op2->ts.type == BT_COMPLEX)
+	{
+	  temp.ts.type = BT_LOGICAL;
+	  temp.ts.kind = gfc_default_logical_kind();
+	  goto runtime;
+	}
+
+      /* else fall through */
+
+    case INTRINSIC_EQ:
+    case INTRINSIC_NE:
+      if (op1->ts.type == BT_CHARACTER && op2->ts.type == BT_CHARACTER)
+	{
+	  unary = 0;
+	  temp.ts.type = BT_LOGICAL;
+	  temp.ts.kind = gfc_default_logical_kind();
+	  break;
+	}
+
+      /* else fall through */
+
+    case INTRINSIC_PLUS:
+    case INTRINSIC_MINUS:
+    case INTRINSIC_TIMES:
+    case INTRINSIC_DIVIDE:
+    case INTRINSIC_POWER:	/* Numeric binary */
+      if (!gfc_numeric_ts (&op1->ts) || !gfc_numeric_ts (&op2->ts))
+	goto runtime;
+
+      /* Insert any necessary type conversions to make the operands compatible.  */
+
+      temp.expr_type = EXPR_OP;
+      gfc_clear_ts (&temp.ts);
+      temp.operator = operator;
+
+      temp.op1 = op1;
+      temp.op2 = op2;
+
+      gfc_type_convert_binary (&temp);
+
+      if (operator == INTRINSIC_EQ || operator == INTRINSIC_NE
+	  || operator == INTRINSIC_GE || operator == INTRINSIC_GT
+	  || operator == INTRINSIC_LE || operator == INTRINSIC_LT)
+	{
+	  temp.ts.type = BT_LOGICAL;
+	  temp.ts.kind = gfc_default_logical_kind ();
+	}
+
+      unary = 0;
+      break;
+
+    case INTRINSIC_CONCAT:	/* Character binary */
+      if (op1->ts.type != BT_CHARACTER || op2->ts.type != BT_CHARACTER)
+	goto runtime;
+
+      temp.ts.type = BT_CHARACTER;
+      temp.ts.kind = gfc_default_character_kind ();
+
+      unary = 0;
+      break;
+
+    case INTRINSIC_USER:
+      goto runtime;
+
+    default:
+      gfc_internal_error ("eval_intrinsic(): Bad operator");
+    }
+
+  /* Try to combine the operators.  */
+  if (operator == INTRINSIC_POWER && op2->ts.type != BT_INTEGER)
+    goto runtime;
+
+  if (op1->expr_type != EXPR_CONSTANT
+      && (op1->expr_type != EXPR_ARRAY
+	  || !gfc_is_constant_expr (op1)
+	  || !gfc_expanded_ac (op1)))
+    goto runtime;
+
+  if (op2 != NULL
+      && op2->expr_type != EXPR_CONSTANT
+      && (op2->expr_type != EXPR_ARRAY
+	  || !gfc_is_constant_expr (op2)
+	  || !gfc_expanded_ac (op2)))
+    goto runtime;
+
+  if (unary)
+    rc = reduce_unary (eval.f2, op1, &result);
+  else
+    rc = reduce_binary (eval.f3, op1, op2, &result);
+
+  if (rc != ARITH_OK)
+    {				/* Something went wrong */
+      gfc_error ("%s at %L", gfc_arith_error (rc), &op1->where);
+      return NULL;
+    }
+
+  gfc_free_expr (op1);
+  gfc_free_expr (op2);
+  return result;
+
+runtime:
+  /* Create a run-time expression */
+  result = gfc_get_expr ();
+  result->ts = temp.ts;
+
+  result->expr_type = EXPR_OP;
+  result->operator = operator;
+
+  result->op1 = op1;
+  result->op2 = op2;
+
+  result->where = op1->where;
+
+  return result;
+}
+
+
+/* Modify type of expression for zero size array.  */
+static gfc_expr *
+eval_type_intrinsic0 (gfc_intrinsic_op operator, gfc_expr *op)
+{
+  if (op == NULL)
+    gfc_internal_error("eval_type_intrinsic0(): op NULL");
+
+  switch(operator)
+    {
+    case INTRINSIC_GE:
+    case INTRINSIC_LT:
+    case INTRINSIC_LE:
+    case INTRINSIC_GT:
+    case INTRINSIC_EQ:
+    case INTRINSIC_NE:
+      op->ts.type = BT_LOGICAL;
+      op->ts.kind = gfc_default_logical_kind();
+      break;
+
+    default:
+      break;
+    }
+
+  return op;
+}
+
+
+/* Return nonzero if the expression is a zero size array.  */
+
+static int
+gfc_zero_size_array (gfc_expr * e)
+{
+
+  if (e->expr_type != EXPR_ARRAY)
+    return 0;
+
+  return e->value.constructor == NULL;
+}
+
+
+/* Reduce a binary expression where at least one of the operands
+   involves a zero-length array.  Returns NULL if neither of the
+   operands is a zero-length array.  */
+
+static gfc_expr *
+reduce_binary0 (gfc_expr * op1, gfc_expr * op2)
+{
+
+  if (gfc_zero_size_array (op1))
+    {
+      gfc_free_expr (op2);
+      return op1;
+    }
+
+  if (gfc_zero_size_array (op2))
+    {
+      gfc_free_expr (op1);
+      return op2;
+    }
+
+  return NULL;
+}
+
+
+static gfc_expr *
+eval_intrinsic_f2 (gfc_intrinsic_op operator,
+		   arith (*eval) (gfc_expr *, gfc_expr **),
+		   gfc_expr * op1, gfc_expr * op2)
+{
+  gfc_expr *result;
+  eval_f f;
+
+  if (op2 == NULL)
+    {
+      if (gfc_zero_size_array (op1))
+	return eval_type_intrinsic0(operator, op1);
+    }
+  else
+    {
+      result = reduce_binary0 (op1, op2);
+      if (result != NULL)
+	return eval_type_intrinsic0(operator, result);
+    }
+
+  f.f2 = eval;
+  return eval_intrinsic (operator, f, op1, op2);
+}
+
+
+static gfc_expr *
+eval_intrinsic_f3 (gfc_intrinsic_op operator,
+		   arith (*eval) (gfc_expr *, gfc_expr *, gfc_expr **),
+		   gfc_expr * op1, gfc_expr * op2)
+{
+  gfc_expr *result;
+  eval_f f;
+
+  result = reduce_binary0 (op1, op2);
+  if (result != NULL)
+    return eval_type_intrinsic0(operator, result);
+
+  f.f3 = eval;
+  return eval_intrinsic (operator, f, op1, op2);
+}
+
+
+
+gfc_expr *
+gfc_uplus (gfc_expr * op)
+{
+  return eval_intrinsic_f2 (INTRINSIC_UPLUS, gfc_arith_uplus, op, NULL);
+}
+
+gfc_expr *
+gfc_uminus (gfc_expr * op)
+{
+  return eval_intrinsic_f2 (INTRINSIC_UMINUS, gfc_arith_uminus, op, NULL);
+}
+
+gfc_expr *
+gfc_add (gfc_expr * op1, gfc_expr * op2)
+{
+  return eval_intrinsic_f3 (INTRINSIC_PLUS, gfc_arith_plus, op1, op2);
+}
+
+gfc_expr *
+gfc_subtract (gfc_expr * op1, gfc_expr * op2)
+{
+  return eval_intrinsic_f3 (INTRINSIC_MINUS, gfc_arith_minus, op1, op2);
+}
+
+gfc_expr *
+gfc_multiply (gfc_expr * op1, gfc_expr * op2)
+{
+  return eval_intrinsic_f3 (INTRINSIC_TIMES, gfc_arith_times, op1, op2);
+}
+
+gfc_expr *
+gfc_divide (gfc_expr * op1, gfc_expr * op2)
+{
+  return eval_intrinsic_f3 (INTRINSIC_DIVIDE, gfc_arith_divide, op1, op2);
+}
+
+gfc_expr *
+gfc_power (gfc_expr * op1, gfc_expr * op2)
+{
+  return eval_intrinsic_f3 (INTRINSIC_POWER, gfc_arith_power, op1, op2);
+}
+
+gfc_expr *
+gfc_concat (gfc_expr * op1, gfc_expr * op2)
+{
+  return eval_intrinsic_f3 (INTRINSIC_CONCAT, gfc_arith_concat, op1, op2);
+}
+
+gfc_expr *
+gfc_and (gfc_expr * op1, gfc_expr * op2)
+{
+  return eval_intrinsic_f3 (INTRINSIC_AND, gfc_arith_and, op1, op2);
+}
+
+gfc_expr *
+gfc_or (gfc_expr * op1, gfc_expr * op2)
+{
+  return eval_intrinsic_f3 (INTRINSIC_OR, gfc_arith_or, op1, op2);
+}
+
+gfc_expr *
+gfc_not (gfc_expr * op1)
+{
+  return eval_intrinsic_f2 (INTRINSIC_NOT, gfc_arith_not, op1, NULL);
+}
+
+gfc_expr *
+gfc_eqv (gfc_expr * op1, gfc_expr * op2)
+{
+  return eval_intrinsic_f3 (INTRINSIC_EQV, gfc_arith_eqv, op1, op2);
+}
+
+gfc_expr *
+gfc_neqv (gfc_expr * op1, gfc_expr * op2)
+{
+  return eval_intrinsic_f3 (INTRINSIC_NEQV, gfc_arith_neqv, op1, op2);
+}
+
+gfc_expr *
+gfc_eq (gfc_expr * op1, gfc_expr * op2)
+{
+  return eval_intrinsic_f3 (INTRINSIC_EQ, gfc_arith_eq, op1, op2);
+}
+
+gfc_expr *
+gfc_ne (gfc_expr * op1, gfc_expr * op2)
+{
+  return eval_intrinsic_f3 (INTRINSIC_NE, gfc_arith_ne, op1, op2);
+}
+
+gfc_expr *
+gfc_gt (gfc_expr * op1, gfc_expr * op2)
+{
+  return eval_intrinsic_f3 (INTRINSIC_GT, gfc_arith_gt, op1, op2);
+}
+
+gfc_expr *
+gfc_ge (gfc_expr * op1, gfc_expr * op2)
+{
+  return eval_intrinsic_f3 (INTRINSIC_GE, gfc_arith_ge, op1, op2);
+}
+
+gfc_expr *
+gfc_lt (gfc_expr * op1, gfc_expr * op2)
+{
+  return eval_intrinsic_f3 (INTRINSIC_LT, gfc_arith_lt, op1, op2);
+}
+
+gfc_expr *
+gfc_le (gfc_expr * op1, gfc_expr * op2)
+{
+  return eval_intrinsic_f3 (INTRINSIC_LE, gfc_arith_le, op1, op2);
+}
+
+
+/* Convert an integer string to an expression node.  */
+
+gfc_expr *
+gfc_convert_integer (const char *buffer, int kind, int radix, locus * where)
+{
+  gfc_expr *e;
+  const char *t;
+
+  e = gfc_constant_result (BT_INTEGER, kind, where);
+  /* a leading plus is allowed, but not by mpz_set_str */
+  if (buffer[0] == '+')
+    t = buffer + 1;
+  else
+    t = buffer;
+  mpz_set_str (e->value.integer, t, radix);
+
+  return e;
+}
+
+
+/* Convert a real string to an expression node.  */
+
+gfc_expr *
+gfc_convert_real (const char *buffer, int kind, locus * where)
+{
+  gfc_expr *e;
+  const char *t;
+
+  e = gfc_constant_result (BT_REAL, kind, where);
+  /* a leading plus is allowed, but not by mpf_set_str */
+  if (buffer[0] == '+')
+    t = buffer + 1;
+  else
+    t = buffer;
+  mpf_set_str (e->value.real, t, 10);
+
+  return e;
+}
+
+
+/* Convert a pair of real, constant expression nodes to a single
+   complex expression node.  */
+
+gfc_expr *
+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);
+
+  return e;
+}
+
+
+/******* Simplification of intrinsic functions with constant arguments *****/
+
+
+/* Deal with an arithmetic error.  */
+
+static void
+arith_error (arith rc, gfc_typespec * from, gfc_typespec * to, locus * where)
+{
+
+  gfc_error ("%s converting %s to %s at %L", gfc_arith_error (rc),
+	     gfc_typename (from), gfc_typename (to), where);
+
+  /* TODO: Do something about the error, ie underflow rounds to 0,
+     throw exception, return NaN, etc.  */
+}
+
+/* Convert integers to integers.  */
+
+gfc_expr *
+gfc_int2int (gfc_expr * src, int kind)
+{
+  gfc_expr *result;
+  arith rc;
+
+  result = gfc_constant_result (BT_INTEGER, kind, &src->where);
+
+  mpz_set (result->value.integer, src->value.integer);
+
+  if ((rc = gfc_check_integer_range (result->value.integer, kind))
+      != ARITH_OK)
+    {
+      arith_error (rc, &src->ts, &result->ts, &src->where);
+      gfc_free_expr (result);
+      return NULL;
+    }
+
+  return result;
+}
+
+
+/* Convert integers to reals.  */
+
+gfc_expr *
+gfc_int2real (gfc_expr * src, int kind)
+{
+  gfc_expr *result;
+  arith rc;
+
+  result = gfc_constant_result (BT_REAL, kind, &src->where);
+
+  mpf_set_z (result->value.real, src->value.integer);
+
+  if ((rc = gfc_check_real_range (result->value.real, kind)) != ARITH_OK)
+    {
+      arith_error (rc, &src->ts, &result->ts, &src->where);
+      gfc_free_expr (result);
+      return NULL;
+    }
+
+  return result;
+}
+
+
+/* Convert default integer to default complex.  */
+
+gfc_expr *
+gfc_int2complex (gfc_expr * src, int kind)
+{
+  gfc_expr *result;
+  arith rc;
+
+  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);
+
+  if ((rc = gfc_check_real_range (result->value.complex.i, kind)) != ARITH_OK)
+    {
+      arith_error (rc, &src->ts, &result->ts, &src->where);
+      gfc_free_expr (result);
+      return NULL;
+    }
+
+  return result;
+}
+
+
+/* Convert default real to default integer.  */
+
+gfc_expr *
+gfc_real2int (gfc_expr * src, int kind)
+{
+  gfc_expr *result;
+  arith rc;
+
+  result = gfc_constant_result (BT_INTEGER, kind, &src->where);
+
+  mpz_set_f (result->value.integer, src->value.real);
+
+  if ((rc = gfc_check_integer_range (result->value.integer, kind))
+      != ARITH_OK)
+    {
+      arith_error (rc, &src->ts, &result->ts, &src->where);
+      gfc_free_expr (result);
+      return NULL;
+    }
+
+  return result;
+}
+
+
+/* Convert real to real.  */
+
+gfc_expr *
+gfc_real2real (gfc_expr * src, int kind)
+{
+  gfc_expr *result;
+  arith rc;
+
+  result = gfc_constant_result (BT_REAL, kind, &src->where);
+
+  mpf_set (result->value.real, src->value.real);
+
+  if ((rc = gfc_check_real_range (result->value.real, kind)) != ARITH_OK)
+    {
+      arith_error (rc, &src->ts, &result->ts, &src->where);
+      gfc_free_expr (result);
+      return NULL;
+    }
+
+  return result;
+}
+
+
+/* Convert real to complex.  */
+
+gfc_expr *
+gfc_real2complex (gfc_expr * src, int kind)
+{
+  gfc_expr *result;
+  arith rc;
+
+  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);
+
+  if ((rc = gfc_check_real_range (result->value.complex.i, kind)) != ARITH_OK)
+    {
+      arith_error (rc, &src->ts, &result->ts, &src->where);
+      gfc_free_expr (result);
+      return NULL;
+    }
+
+  return result;
+}
+
+
+/* Convert complex to integer.  */
+
+gfc_expr *
+gfc_complex2int (gfc_expr * src, int kind)
+{
+  gfc_expr *result;
+  arith rc;
+
+  result = gfc_constant_result (BT_INTEGER, kind, &src->where);
+
+  mpz_set_f (result->value.integer, src->value.complex.r);
+
+  if ((rc = gfc_check_integer_range (result->value.integer, kind))
+      != ARITH_OK)
+    {
+      arith_error (rc, &src->ts, &result->ts, &src->where);
+      gfc_free_expr (result);
+      return NULL;
+    }
+
+  return result;
+}
+
+
+/* Convert complex to real.  */
+
+gfc_expr *
+gfc_complex2real (gfc_expr * src, int kind)
+{
+  gfc_expr *result;
+  arith rc;
+
+  result = gfc_constant_result (BT_REAL, kind, &src->where);
+
+  mpf_set (result->value.real, src->value.complex.r);
+
+  if ((rc = gfc_check_real_range (result->value.real, kind)) != ARITH_OK)
+    {
+      arith_error (rc, &src->ts, &result->ts, &src->where);
+      gfc_free_expr (result);
+      return NULL;
+    }
+
+  return result;
+}
+
+
+/* Convert complex to complex.  */
+
+gfc_expr *
+gfc_complex2complex (gfc_expr * src, int kind)
+{
+  gfc_expr *result;
+  arith rc;
+
+  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);
+
+  if ((rc = gfc_check_real_range (result->value.complex.r, kind)) != ARITH_OK
+      || (rc =
+	  gfc_check_real_range (result->value.complex.i, kind)) != ARITH_OK)
+    {
+      arith_error (rc, &src->ts, &result->ts, &src->where);
+      gfc_free_expr (result);
+      return NULL;
+    }
+
+  return result;
+}
+
+
+/* Logical kind conversion.  */
+
+gfc_expr *
+gfc_log2log (gfc_expr * src, int kind)
+{
+  gfc_expr *result;
+
+  result = gfc_constant_result (BT_LOGICAL, kind, &src->where);
+  result->value.logical = src->value.logical;
+
+  return result;
+}