OSDN Git Service

2008-11-24 Jerry DeLisle <jvdelisle@gcc.gnu.org>
[pf3gnuchains/gcc-fork.git] / gcc / fortran / arith.c
index 4c036ae..7440be3 100644 (file)
@@ -1,13 +1,13 @@
 /* Compiler arithmetic
-   Copyright (C) 2000, 2001, 2002, 2003, 2004 Free Software Foundation,
-   Inc.
+   Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008
+   Free Software Foundation, Inc.
    Contributed by Andy Vaught
 
 This file is part of GCC.
 
 GCC 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
+Software Foundation; either version 3, or (at your option) any later
 version.
 
 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
@@ -16,606 +16,67 @@ 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 GCC; see the file COPYING.  If not, write to the Free
-Software Foundation, 59 Temple Place - Suite 330, Boston, MA
-02111-1307, USA.  */
+along with GCC; see the file COPYING3.  If not see
+<http://www.gnu.org/licenses/>.  */
 
 /* 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.  */
+   would evaluate them.  We use the GNU MP library and the MPFR
+   library to do arithmetic, and this file provides the interface.  */
 
 #include "config.h"
-
-#include <string.h>
-
+#include "system.h"
+#include "flags.h"
 #include "gfortran.h"
 #include "arith.h"
+#include "target-memory.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/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, locus *where)
 {
-  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);
+  mp_exp_t e;
 
-  for (i = 1; i < GFC_REAL_BITS; i++)
+  if (mpfr_inf_p (x) || mpfr_nan_p (x))
     {
-      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);
+      gfc_error ("Conversion of an Infinity or Not-a-Number at %L "
+                "to INTEGER", where);
+      mpz_set_ui (z, 0);
+      return;
     }
 
-  /* 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 Maclaurin 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);
+  e = mpfr_get_z_exp (z, x);
 
-  if (mpf_cmp_ui (x, 0) == 0)
-    {
-      mpf_set_ui (*result, 1);
-    }
-  else if (mpf_cmp_ui (x, 1) == 0)
-    {
-      mpf_set (*result, e);
-    }
+  if (e > 0)
+    mpz_mul_2exp (z, z, 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);
+    mpz_tdiv_q_2exp (z, z, -e);
 }
 
 
-/* Calculate sin(arg).
-
-   We use a reduction of the form
-
-     x= N*2pi + r
-
-   Then we obtain sin(r) from the Maclaurin series.  */
+/* Set the model number precision by the requested KIND.  */
 
 void
-sine (mpf_t * arg, mpf_t * result)
+gfc_set_model_kind (int kind)
 {
-  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);
-    }
+  int index = gfc_validate_kind (BT_REAL, kind, false);
+  int base2prec;
 
-  mpf_clear (x);
+  base2prec = gfc_real_kinds[index].digits;
+  if (gfc_real_kinds[index].radix != 2)
+    base2prec *= gfc_real_kinds[index].radix / 2;
+  mpfr_set_default_prec (base2prec);
 }
 
 
-/* Calculate cos(arg).
-
-   Similar to sine.  */
+/* Set the model number precision from mpfr_t x.  */
 
 void
-cosine (mpf_t * arg, mpf_t * result)
+gfc_set_model (mpfr_t x)
 {
-  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);
+  mpfr_set_default_prec (mpfr_get_prec (x));
 }
 
 
@@ -630,22 +91,26 @@ gfc_arith_error (arith code)
   switch (code)
     {
     case ARITH_OK:
-      p = "Arithmetic OK";
+      p = _("Arithmetic OK at %L");
       break;
     case ARITH_OVERFLOW:
-      p = "Arithmetic overflow";
+      p = _("Arithmetic overflow at %L");
       break;
     case ARITH_UNDERFLOW:
-      p = "Arithmetic underflow";
+      p = _("Arithmetic underflow at %L");
       break;
-    case ARITH_DIV0:
-      p = "Division by zero";
+    case ARITH_NAN:
+      p = _("Arithmetic NaN at %L");
       break;
-    case ARITH_0TO0:
-      p = "Indeterminate form 0 ** 0";
+    case ARITH_DIV0:
+      p = _("Division by zero at %L");
       break;
     case ARITH_INCOMMENSURATE:
-      p = "Array operands are incommensurate";
+      p = _("Array operands are incommensurate at %L");
+      break;
+    case ARITH_ASYMMETRIC:
+      p =
+       _("Integer outside symmetric range implied by Standard Fortran at %L");
       break;
     default:
       gfc_internal_error ("gfc_arith_error(): Bad error code");
@@ -662,160 +127,117 @@ 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);
-    }
+  mpfr_t a, b;
+  int i;
 
-  mpf_mul_ui (two_pi, pi, 2);
-  mpf_div_ui (half_pi, pi, 2);
+  mpfr_set_default_prec (128);
+  mpfr_init (a);
 
-  /* Convert the minimum/maximum values for each kind into their
+  /* Convert the minimum and 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);
-
+      /* Huge  */
       mpz_init (int_info->huge);
-      mpz_sub_ui (int_info->huge, r, 1);
+      mpz_set_ui (int_info->huge, int_info->radix);
+      mpz_pow_ui (int_info->huge, int_info->huge, int_info->digits);
+      mpz_sub_ui (int_info->huge, int_info->huge, 1);
 
       /* These are the numbers that are actually representable by the
-         target.  For bases other than two, this needs to be changed.  */
+        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");
+       gfc_internal_error ("Fix min_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);
+      /* See PRs 13490 and 17912, related to integer ranges.
+        The pedantic_min_int exists for range checking when a program
+        is compiled with -pedantic, and reflects the belief that
+        Standard Fortran requires integers to be symmetrical, i.e.
+        every negative integer must have a representable positive
+        absolute value, and vice versa.  */
 
-      /* Tiny */
-      mpf_set_ui (b, real_info->radix);
-      mpf_pow_ui (b, b, 1 - real_info->min_exponent);
+      mpz_init (int_info->pedantic_min_int);
+      mpz_neg (int_info->pedantic_min_int, int_info->huge);
 
-      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);
+      mpz_init (int_info->min_int);
+      mpz_sub_ui (int_info->min_int, int_info->pedantic_min_int, 1);
 
-      /* Precision */
-      mpf_set_ui (a, real_info->radix);
-      common_logarithm (&a, &a);
+      /* Range  */
+      mpfr_set_z (a, int_info->huge, GFC_RND_MODE);
+      mpfr_log10 (a, a, GFC_RND_MODE);
+      mpfr_trunc (a, a);
+      int_info->range = (int) mpfr_get_si (a, GFC_RND_MODE);
+    }
 
-      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);
+  mpfr_clear (a);
 
-      /* If the radix is an integral power of 10, add one to the
-         precision.  */
+  for (real_info = gfc_real_kinds; real_info->kind != 0; real_info++)
+    {
+      gfc_set_model_kind (real_info->kind);
+
+      mpfr_init (a);
+      mpfr_init (b);
+
+      /* huge(x) = (1 - b**(-p)) * b**(emax-1) * b  */
+      /* 1 - b**(-p)  */
+      mpfr_init (real_info->huge);
+      mpfr_set_ui (real_info->huge, 1, GFC_RND_MODE);
+      mpfr_set_ui (a, real_info->radix, GFC_RND_MODE);
+      mpfr_pow_si (a, a, -real_info->digits, GFC_RND_MODE);
+      mpfr_sub (real_info->huge, real_info->huge, a, GFC_RND_MODE);
+
+      /* b**(emax-1)  */
+      mpfr_set_ui (a, real_info->radix, GFC_RND_MODE);
+      mpfr_pow_ui (a, a, real_info->max_exponent - 1, GFC_RND_MODE);
+
+      /* (1 - b**(-p)) * b**(emax-1)  */
+      mpfr_mul (real_info->huge, real_info->huge, a, GFC_RND_MODE);
+
+      /* (1 - b**(-p)) * b**(emax-1) * b  */
+      mpfr_mul_ui (real_info->huge, real_info->huge, real_info->radix,
+                  GFC_RND_MODE);
+
+      /* tiny(x) = b**(emin-1)  */
+      mpfr_init (real_info->tiny);
+      mpfr_set_ui (real_info->tiny, real_info->radix, GFC_RND_MODE);
+      mpfr_pow_si (real_info->tiny, real_info->tiny,
+                  real_info->min_exponent - 1, GFC_RND_MODE);
+
+      /* subnormal (x) = b**(emin - digit)  */
+      mpfr_init (real_info->subnormal);
+      mpfr_set_ui (real_info->subnormal, real_info->radix, GFC_RND_MODE);
+      mpfr_pow_si (real_info->subnormal, real_info->subnormal,
+                  real_info->min_exponent - real_info->digits, GFC_RND_MODE);
+
+      /* epsilon(x) = b**(1-p)  */
+      mpfr_init (real_info->epsilon);
+      mpfr_set_ui (real_info->epsilon, real_info->radix, GFC_RND_MODE);
+      mpfr_pow_si (real_info->epsilon, real_info->epsilon,
+                  1 - real_info->digits, GFC_RND_MODE);
+
+      /* 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);
+
+      /* a = min(a, b)  */
+      mpfr_min (a, a, b, GFC_RND_MODE);
+      mpfr_trunc (a, a);
+      real_info->range = (int) mpfr_get_si (a, GFC_RND_MODE);
+
+      /* 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);
+      mpfr_mul_ui (a, a, real_info->digits - 1, GFC_RND_MODE);
+      mpfr_trunc (a, a);
+      real_info->precision = (int) mpfr_get_si (a, GFC_RND_MODE);
+
+      /* 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);
+      mpfr_clears (a, b, NULL);
+    }
 }
 
 
@@ -827,190 +249,60 @@ 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->pedantic_min_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 ();
+    mpfr_clears (rp->epsilon, rp->huge, rp->tiny, rp->subnormal, NULL);
 }
 
 
-/* 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)
+/* Given a wide character value and a character kind, determine whether
+   the character is representable for that kind.  */
+bool
+gfc_check_character_range (gfc_char_t c, int kind)
 {
-  int i;
+  /* As wide characters are stored as 32-bit values, they're all
+     representable in UCS=4.  */
+  if (kind == 4)
+    return true;
 
-  for (i = 0;; i++)
-    {
-      if (gfc_integer_kinds[i].kind == 0)
-       {
-         i = -1;
-         break;
-       }
-      if (gfc_integer_kinds[i].kind == kind)
-       break;
-    }
+  if (kind == 1)
+    return c <= 255 ? true : false;
 
-  return i;
+  gcc_unreachable ();
 }
 
 
-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;
-}
-
+/* Given an integer and a kind, make sure that the integer lies within
+   the range of the kind.  Returns ARITH_OK, ARITH_ASYMMETRIC or
+   ARITH_OVERFLOW.  */
 
-static int
-validate_logical (int kind)
+arith
+gfc_check_integer_range (mpz_t p, int kind)
 {
+  arith result;
   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;
+  i = gfc_validate_kind (BT_INTEGER, kind, false);
+  result = ARITH_OK;
 
-  switch (type)
+  if (pedantic)
     {
-    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");
+      if (mpz_cmp (p, gfc_integer_kinds[i].pedantic_min_int) < 0)
+       result = ARITH_ASYMMETRIC;
     }
 
-  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 (gfc_option.flag_range_check == 0)
+    return result;
 
   if (mpz_cmp (p, gfc_integer_kinds[i].min_int) < 0
-      || mpz_cmp (p, gfc_integer_kinds[i].max_int) > 0)
+      || mpz_cmp (p, gfc_integer_kinds[i].huge) > 0)
     result = ARITH_OVERFLOW;
 
   return result;
@@ -1022,50 +314,100 @@ 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 = gfc_validate_kind (BT_REAL, kind, false);
 
-  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)
-    goto done;
 
-  if (mpf_cmp (q, gfc_real_kinds[i].huge) == 1)
+  if (mpfr_inf_p (p))
     {
-      retval = ARITH_OVERFLOW;
-      goto done;
+      if (gfc_option.flag_range_check != 0)
+       retval = ARITH_OVERFLOW;
+    }
+  else if (mpfr_nan_p (p))
+    {
+      if (gfc_option.flag_range_check != 0)
+       retval = ARITH_NAN;
+    }
+  else if (mpfr_sgn (q) == 0)
+    {
+      mpfr_clear (q);
+      return retval;
+    }
+  else if (mpfr_cmp (q, gfc_real_kinds[i].huge) > 0)
+    {
+      if (gfc_option.flag_range_check == 0)
+       mpfr_set_inf (p, mpfr_sgn (p));
+      else
+       retval = ARITH_OVERFLOW;
+    }
+  else if (mpfr_cmp (q, gfc_real_kinds[i].subnormal) < 0)
+    {
+      if (gfc_option.flag_range_check == 0)
+       {
+         if (mpfr_sgn (p) < 0)
+           {
+             mpfr_set_ui (p, 0, GFC_RND_MODE);
+             mpfr_set_si (q, -1, GFC_RND_MODE);
+             mpfr_copysign (p, p, q, GFC_RND_MODE);
+           }
+         else
+           mpfr_set_ui (p, 0, GFC_RND_MODE);
+       }
+      else
+       retval = ARITH_UNDERFLOW;
+    }
+  else if (mpfr_cmp (q, gfc_real_kinds[i].tiny) < 0)
+    {
+      mp_exp_t emin, emax;
+      int en;
+
+      /* Save current values of emin and emax.  */
+      emin = mpfr_get_emin ();
+      emax = mpfr_get_emax ();
+
+      /* Set emin and emax for the current model number.  */
+      en = gfc_real_kinds[i].min_exponent - gfc_real_kinds[i].digits + 1;
+      mpfr_set_emin ((mp_exp_t) en);
+      mpfr_set_emax ((mp_exp_t) gfc_real_kinds[i].max_exponent);
+      mpfr_check_range (q, 0, GFC_RND_MODE);
+      mpfr_subnormalize (q, 0, GFC_RND_MODE);
+
+      /* Reset emin and emax.  */
+      mpfr_set_emin (emin);
+      mpfr_set_emax (emax);
+
+      /* Copy sign if needed.  */
+      if (mpfr_sgn (p) < 0)
+       mpfr_neg (p, q, GMP_RNDN);
+      else
+       mpfr_set (p, q, GMP_RNDN);
     }
 
-  if (mpf_cmp (q, gfc_real_kinds[i].tiny) == -1)
-    retval = ARITH_UNDERFLOW;
-
-done:
-  mpf_clear (q);
+  mpfr_clear (q);
 
   return retval;
 }
 
 
-/* Function to return a constant expression node of a given type and
-   kind.  */
+/* 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_constant_result (bt type, int kind, locus *where)
 {
   gfc_expr *result;
 
   if (!where)
-    gfc_internal_error
-      ("gfc_constant_result(): locus 'where' cannot be NULL");
+    gfc_internal_error ("gfc_constant_result(): locus 'where' cannot be NULL");
 
   result = gfc_get_expr ();
 
@@ -1081,12 +423,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:
@@ -1104,7 +448,7 @@ gfc_constant_result (bt type, int kind, locus * where)
    zero raised to the zero, etc.  */
 
 static arith
-gfc_arith_not (gfc_expr * op1, gfc_expr ** resultp)
+gfc_arith_not (gfc_expr *op1, gfc_expr **resultp)
 {
   gfc_expr *result;
 
@@ -1117,7 +461,7 @@ gfc_arith_not (gfc_expr * op1, gfc_expr ** resultp)
 
 
 static arith
-gfc_arith_and (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
+gfc_arith_and (gfc_expr *op1, gfc_expr *op2, gfc_expr **resultp)
 {
   gfc_expr *result;
 
@@ -1131,7 +475,7 @@ gfc_arith_and (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
 
 
 static arith
-gfc_arith_or (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
+gfc_arith_or (gfc_expr *op1, gfc_expr *op2, gfc_expr **resultp)
 {
   gfc_expr *result;
 
@@ -1145,7 +489,7 @@ gfc_arith_or (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
 
 
 static arith
-gfc_arith_eqv (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
+gfc_arith_eqv (gfc_expr *op1, gfc_expr *op2, gfc_expr **resultp)
 {
   gfc_expr *result;
 
@@ -1159,7 +503,7 @@ gfc_arith_eqv (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
 
 
 static arith
-gfc_arith_neqv (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
+gfc_arith_neqv (gfc_expr *op1, gfc_expr *op2, gfc_expr **resultp)
 {
   gfc_expr *result;
 
@@ -1177,9 +521,10 @@ gfc_arith_neqv (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
    but that one deals with the intrinsic RANGE function.  */
 
 arith
-gfc_range_check (gfc_expr * e)
+gfc_range_check (gfc_expr *e)
 {
   arith rc;
+  arith rc2;
 
   switch (e->ts.type)
     {
@@ -1189,13 +534,33 @@ gfc_range_check (gfc_expr * e)
 
     case BT_REAL:
       rc = gfc_check_real_range (e->value.real, e->ts.kind);
+      if (rc == ARITH_UNDERFLOW)
+       mpfr_set_ui (e->value.real, 0, GFC_RND_MODE);
+      if (rc == ARITH_OVERFLOW)
+       mpfr_set_inf (e->value.real, mpfr_sgn (e->value.real));
+      if (rc == ARITH_NAN)
+       mpfr_set_nan (e->value.real);
       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);
+      if (rc == ARITH_UNDERFLOW)
+       mpfr_set_ui (e->value.complex.r, 0, GFC_RND_MODE);
+      if (rc == ARITH_OVERFLOW)
+       mpfr_set_inf (e->value.complex.r, mpfr_sgn (e->value.complex.r));
+      if (rc == ARITH_NAN)
+       mpfr_set_nan (e->value.complex.r);
+
+      rc2 = gfc_check_real_range (e->value.complex.i, e->ts.kind);
+      if (rc == ARITH_UNDERFLOW)
+       mpfr_set_ui (e->value.complex.i, 0, GFC_RND_MODE);
+      if (rc == ARITH_OVERFLOW)
+       mpfr_set_inf (e->value.complex.i, mpfr_sgn (e->value.complex.i));
+      if (rc == ARITH_NAN)
+       mpfr_set_nan (e->value.complex.i);
 
+      if (rc == ARITH_OK)
+       rc = rc2;
       break;
 
     default:
@@ -1206,21 +571,51 @@ gfc_range_check (gfc_expr * e)
 }
 
 
+/* Several of the following routines use the same set of statements to
+   check the validity of the result.  Encapsulate the checking here.  */
+
+static arith
+check_result (arith rc, gfc_expr *x, gfc_expr *r, gfc_expr **rp)
+{
+  arith val = rc;
+
+  if (val == ARITH_UNDERFLOW)
+    {
+      if (gfc_option.warn_underflow)
+       gfc_warning (gfc_arith_error (val), &x->where);
+      val = ARITH_OK;
+    }
+
+  if (val == ARITH_ASYMMETRIC)
+    {
+      gfc_warning (gfc_arith_error (val), &x->where);
+      val = ARITH_OK;
+    }
+
+  if (val != ARITH_OK)
+    gfc_free_expr (r);
+  else
+    *rp = r;
+
+  return val;
+}
+
+
 /* 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.  */
+   in the code elsewhere.  Used for unary plus and parenthesized
+   expressions.  */
 
 static arith
-gfc_arith_uplus (gfc_expr * op1, gfc_expr ** resultp)
+gfc_arith_identity (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_arith_uminus (gfc_expr *op1, gfc_expr **resultp)
 {
   gfc_expr *result;
   arith rc;
@@ -1234,12 +629,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:
@@ -1248,17 +643,12 @@ gfc_arith_uminus (gfc_expr * op1, gfc_expr ** resultp)
 
   rc = gfc_range_check (result);
 
-  if (rc != ARITH_OK)
-    gfc_free_expr (result);
-  else
-    *resultp = result;
-
-  return rc;
+  return check_result (rc, op1, result, resultp);
 }
 
 
 static arith
-gfc_arith_plus (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
+gfc_arith_plus (gfc_expr *op1, gfc_expr *op2, gfc_expr **resultp)
 {
   gfc_expr *result;
   arith rc;
@@ -1272,15 +662,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:
@@ -1289,17 +680,12 @@ gfc_arith_plus (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
 
   rc = gfc_range_check (result);
 
-  if (rc != ARITH_OK)
-    gfc_free_expr (result);
-  else
-    *resultp = result;
-
-  return rc;
+  return check_result (rc, op1, result, resultp);
 }
 
 
 static arith
-gfc_arith_minus (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
+gfc_arith_minus (gfc_expr *op1, gfc_expr *op2, gfc_expr **resultp)
 {
   gfc_expr *result;
   arith rc;
@@ -1313,16 +699,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:
@@ -1331,20 +717,15 @@ gfc_arith_minus (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
 
   rc = gfc_range_check (result);
 
-  if (rc != ARITH_OK)
-    gfc_free_expr (result);
-  else
-    *resultp = result;
-
-  return rc;
+  return check_result (rc, op1, result, resultp);
 }
 
 
 static arith
-gfc_arith_times (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
+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);
@@ -1356,24 +737,24 @@ 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);
+      gfc_set_model (op1->value.complex.r);
+      mpfr_init (x);
+      mpfr_init (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);
+      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);
 
-      mpf_clear (x);
-      mpf_clear (y);
+      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_clears (x, y, NULL);
       break;
 
     default:
@@ -1382,20 +763,15 @@ gfc_arith_times (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
 
   rc = gfc_range_check (result);
 
-  if (rc != ARITH_OK)
-    gfc_free_expr (result);
-  else
-    *resultp = result;
-
-  return rc;
+  return check_result (rc, op1, result, resultp);
 }
 
 
 static arith
-gfc_arith_divide (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
+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;
@@ -1416,45 +792,47 @@ gfc_arith_divide (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
       break;
 
     case BT_REAL:
-      if (mpf_sgn (op2->value.real) == 0)
+      if (mpfr_sgn (op2->value.real) == 0 && gfc_option.flag_range_check == 1)
        {
          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)
+      if (mpfr_sgn (op2->value.complex.r) == 0
+         && mpfr_sgn (op2->value.complex.i) == 0
+         && gfc_option.flag_range_check == 1)
        {
          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);
+      gfc_set_model (op1->value.complex.r);
+      mpfr_init (x);
+      mpfr_init (y);
+      mpfr_init (div);
 
-      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, 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.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.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_clear (x);
-      mpf_clear (y);
-      mpf_clear (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);
 
+      mpfr_clears (x, y, div, NULL);
       break;
 
     default:
@@ -1464,175 +842,217 @@ gfc_arith_divide (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
   if (rc == ARITH_OK)
     rc = gfc_range_check (result);
 
-  if (rc != ARITH_OK)
-    gfc_free_expr (result);
-  else
-    *resultp = result;
-
-  return rc;
+  return check_result (rc, op1, result, resultp);
 }
 
 
 /* Compute the reciprocal of a complex number (guaranteed nonzero).  */
 
 static void
-complex_reciprocal (gfc_expr * op)
+complex_reciprocal (gfc_expr *op)
 {
-  mpf_t mod, a, result_r, result_i;
-
-  mpf_init (mod);
-  mpf_init (a);
+  mpfr_t mod, tmp;
 
-  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 (tmp);
 
-  mpf_init (result_r);
-  mpf_div (result_r, op->value.complex.r, mod);
+  mpfr_mul (mod, op->value.complex.r, op->value.complex.r, GFC_RND_MODE);
+  mpfr_mul (tmp, op->value.complex.i, op->value.complex.i, GFC_RND_MODE);
+  mpfr_add (mod, mod, tmp, GFC_RND_MODE);
 
-  mpf_init (result_i);
-  mpf_neg (result_i, op->value.complex.i);
-  mpf_div (result_i, result_i, mod);
+  mpfr_div (op->value.complex.r, 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 (op->value.complex.i, op->value.complex.i, GFC_RND_MODE);
+  mpfr_div (op->value.complex.i, op->value.complex.i, mod, GFC_RND_MODE);
 
-  mpf_clear (result_r);
-  mpf_clear (result_i);
-
-  mpf_clear (mod);
-  mpf_clear (a);
+  mpfr_clears (tmp, mod, NULL);
 }
 
 
-/* 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);
+/* Raise a complex number to positive power (power > 0).
+   This function will modify the content of power.
 
-  mpf_init (temp_r);
-  mpf_init (temp_i);
-  mpf_init (a);
+   Use Binary Method, which is not an optimal but a simple and reasonable
+   arithmetic. See section 4.6.3, "Evaluation of Powers" of Donald E. Knuth,
+   "Seminumerical Algorithms", Vol. 2, "The Art of Computer Programming",
+   3rd Edition, 1998.  */
 
-  for (; power > 0; power--)
+static void
+complex_pow (gfc_expr *result, gfc_expr *base, mpz_t power)
+{
+  mpfr_t x_r, x_i, tmp, re, im;
+
+  gfc_set_model (base->value.complex.r);
+  mpfr_init (x_r);
+  mpfr_init (x_i);
+  mpfr_init (tmp);
+  mpfr_init (re);
+  mpfr_init (im);
+
+  /* res = 1 */
+  mpfr_set_ui (result->value.complex.r, 1, GFC_RND_MODE);
+  mpfr_set_ui (result->value.complex.i, 0, GFC_RND_MODE);
+
+  /* x = base */
+  mpfr_set (x_r, base->value.complex.r, GFC_RND_MODE);
+  mpfr_set (x_i, base->value.complex.i, GFC_RND_MODE);
+
+  /* Macro for complex multiplication. We have to take care that
+     res_r/res_i and a_r/a_i can (and will) be the same variable.  */
+#define CMULT(res_r,res_i,a_r,a_i,b_r,b_i) \
+    mpfr_mul (re, a_r, b_r, GFC_RND_MODE), \
+    mpfr_mul (tmp, a_i, b_i, GFC_RND_MODE), \
+    mpfr_sub (re, re, tmp, GFC_RND_MODE), \
+    \
+    mpfr_mul (im, a_r, b_i, GFC_RND_MODE), \
+    mpfr_mul (tmp, a_i, b_r, GFC_RND_MODE), \
+    mpfr_add (res_i, im, tmp, GFC_RND_MODE), \
+    mpfr_set (res_r, re, GFC_RND_MODE)
+  
+#define res_r result->value.complex.r
+#define res_i result->value.complex.i
+
+  /* for (; power > 0; x *= x) */
+  for (; mpz_cmp_si (power, 0) > 0; CMULT(x_r,x_i,x_r,x_i,x_r,x_i))
     {
-      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);
+      /* if (power & 1) res = res * x; */
+      if (mpz_congruent_ui_p (power, 1, 2))
+       CMULT(res_r,res_i,res_r,res_i,x_r,x_i);
 
-      mpf_set (result->value.complex.r, temp_r);
-      mpf_set (result->value.complex.i, temp_i);
+      /* power /= 2; */
+      mpz_fdiv_q_ui (power, power, 2);
     }
 
-  mpf_clear (temp_r);
-  mpf_clear (temp_i);
-  mpf_clear (a);
+#undef res_r
+#undef res_i
+#undef CMULT
+
+  mpfr_clears (x_r, x_i, tmp, re, im, NULL);
 }
 
 
 /* Raise a number to an integer power.  */
 
 static arith
-gfc_arith_power (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
+gfc_arith_power (gfc_expr *op1, gfc_expr *op2, gfc_expr **resultp)
 {
-  int power, apower;
+  int power_sign;
   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");
+  gcc_assert (op2->expr_type == EXPR_CONSTANT && op2->ts.type == BT_INTEGER);
 
+  rc = ARITH_OK;
   result = gfc_constant_result (op1->ts.type, op1->ts.kind, &op1->where);
+  power_sign = mpz_sgn (op2->value.integer);
 
-  if (power == 0)
-    {                          /* Handle something to the zeroth power */
+  if (power_sign == 0)
+    {
+      /* Handle something to the zeroth power.  Since we're dealing
+        with integral exponents, there is no ambiguity in the
+        limiting procedure used to determine the value of 0**0.  */
       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);
-
+         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);
-
+         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)
-           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;
 
        default:
          gfc_internal_error ("gfc_arith_power(): Bad base");
        }
     }
-
-  if (power != 0)
+  else
     {
-      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);
-           }
-
+         {
+           int power;
+
+           /* First, we simplify the cases of op1 == 1, 0 or -1.  */
+           if (mpz_cmp_si (op1->value.integer, 1) == 0)
+             {
+               /* 1**op2 == 1 */
+               mpz_set_si (result->value.integer, 1);
+             }
+           else if (mpz_cmp_si (op1->value.integer, 0) == 0)
+             {
+               /* 0**op2 == 0, if op2 > 0
+                  0**op2 overflow, if op2 < 0 ; in that case, we
+                  set the result to 0 and return ARITH_DIV0.  */
+               mpz_set_si (result->value.integer, 0);
+               if (mpz_cmp_si (op2->value.integer, 0) < 0)
+                 rc = ARITH_DIV0;
+             }
+           else if (mpz_cmp_si (op1->value.integer, -1) == 0)
+             {
+               /* (-1)**op2 == (-1)**(mod(op2,2)) */
+               unsigned int odd = mpz_fdiv_ui (op2->value.integer, 2);
+               if (odd)
+                 mpz_set_si (result->value.integer, -1);
+               else
+                 mpz_set_si (result->value.integer, 1);
+             }
+           /* Then, we take care of op2 < 0.  */
+           else if (mpz_cmp_si (op2->value.integer, 0) < 0)
+             {
+               /* if op2 < 0, op1**op2 == 0  because abs(op1) > 1.  */
+               mpz_set_si (result->value.integer, 0);
+             }
+           else if (gfc_extract_int (op2, &power) != NULL)
+             {
+               /* If op2 doesn't fit in an int, the exponentiation will
+                  overflow, because op2 > 0 and abs(op1) > 1.  */
+               mpz_t max;
+               int i = gfc_validate_kind (BT_INTEGER, result->ts.kind, false);
+
+               if (gfc_option.flag_range_check)
+                 rc = ARITH_OVERFLOW;
+
+               /* Still, we want to give the same value as the processor.  */
+               mpz_init (max);
+               mpz_add_ui (max, gfc_integer_kinds[i].huge, 1);
+               mpz_mul_ui (max, max, 2);
+               mpz_powm (result->value.integer, op1->value.integer,
+                         op2->value.integer, max);
+               mpz_clear (max);
+             }
+           else
+             mpz_pow_ui (result->value.integer, op1->value.integer, power);
+         }
          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);
-           }
-
+         mpfr_pow_z (result->value.real, op1->value.real, op2->value.integer,
+                     GFC_RND_MODE);
          break;
 
        case BT_COMPLEX:
-         complex_pow_ui (op1, apower, result);
-         if (power < 0)
-           complex_reciprocal (result);
+         {
+           mpz_t apower;
 
-         break;
+           /* Compute op1**abs(op2)  */
+           mpz_init (apower);
+           mpz_abs (apower, op2->value.integer);
+           complex_pow (result, op1, apower);
+           mpz_clear (apower);
+
+           /* If (op2 < 0), compute the inverse.  */
+           if (power_sign < 0)
+             complex_reciprocal (result);
+
+           break;
+         }
 
        default:
          break;
@@ -1642,36 +1062,33 @@ gfc_arith_power (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
   if (rc == ARITH_OK)
     rc = gfc_range_check (result);
 
-  if (rc != ARITH_OK)
-    gfc_free_expr (result);
-  else
-    *resultp = result;
-
-  return rc;
+  return check_result (rc, op1, result, resultp);
 }
 
 
 /* Concatenate two string constants.  */
 
 static arith
-gfc_arith_concat (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
+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 (),
+  gcc_assert (op1->ts.kind == op2->ts.kind);
+  result = gfc_constant_result (BT_CHARACTER, op1->ts.kind,
                                &op1->where);
 
   len = op1->value.character.length + op2->value.character.length;
 
-  result->value.character.string = gfc_getmem (len + 1);
+  result->value.character.string = gfc_get_wide_string (len + 1);
   result->value.character.length = len;
 
   memcpy (result->value.character.string, op1->value.character.string,
-         op1->value.character.length);
+         op1->value.character.length * sizeof (gfc_char_t));
 
-  memcpy (result->value.character.string + op1->value.character.length,
-         op2->value.character.string, op2->value.character.length);
+  memcpy (&result->value.character.string[op1->value.character.length],
+         op2->value.character.string,
+         op2->value.character.length * sizeof (gfc_char_t));
 
   result->value.character.string[len] = '\0';
 
@@ -1680,12 +1097,43 @@ gfc_arith_concat (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
   return ARITH_OK;
 }
 
+/* Comparison between real values; returns 0 if (op1 .op. op2) is true.
+   This function mimics mpfr_cmp but takes NaN into account.  */
+
+static int
+compare_real (gfc_expr *op1, gfc_expr *op2, gfc_intrinsic_op op)
+{
+  int rc;
+  switch (op)
+    {
+      case INTRINSIC_EQ:
+       rc = mpfr_equal_p (op1->value.real, op2->value.real) ? 0 : 1;
+       break;
+      case INTRINSIC_GT:
+       rc = mpfr_greater_p (op1->value.real, op2->value.real) ? 1 : -1;
+       break;
+      case INTRINSIC_GE:
+       rc = mpfr_greaterequal_p (op1->value.real, op2->value.real) ? 1 : -1;
+       break;
+      case INTRINSIC_LT:
+       rc = mpfr_less_p (op1->value.real, op2->value.real) ? -1 : 1;
+       break;
+      case INTRINSIC_LE:
+       rc = mpfr_lessequal_p (op1->value.real, op2->value.real) ? -1 : 1;
+       break;
+      default:
+       gfc_internal_error ("compare_real(): Bad operator");
+    }
+
+  return rc;
+}
 
 /* Comparison operators.  Assumes that the two expression nodes
-   contain two constants of the same type.  */
+   contain two constants of the same type. The op argument is
+   needed to handle NaN correctly.  */
 
 int
-gfc_compare_expr (gfc_expr * op1, gfc_expr * op2)
+gfc_compare_expr (gfc_expr *op1, gfc_expr *op2, gfc_intrinsic_op op)
 {
   int rc;
 
@@ -1696,11 +1144,11 @@ gfc_compare_expr (gfc_expr * op1, gfc_expr * op2)
       break;
 
     case BT_REAL:
-      rc = mpf_cmp (op1->value.real, op2->value.real);
+      rc = compare_real (op1, op2, op);
       break;
 
     case BT_CHARACTER:
-      rc = gfc_compare_string (op1, op2, NULL);
+      rc = gfc_compare_string (op1, op2);
       break;
 
     case BT_LOGICAL:
@@ -1717,41 +1165,67 @@ gfc_compare_expr (gfc_expr * op1, gfc_expr * op2)
 
 
 /* Compare a pair of complex numbers.  Naturally, this is only for
-   equality/nonequality.  */
+   equality and inequality.  */
 
 static int
-compare_complex (gfc_expr * op1, gfc_expr * op2)
+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_equal_p (op1->value.complex.r, op2->value.complex.r)
+         && mpfr_equal_p (op1->value.complex.i, op2->value.complex.i));
 }
 
 
-/* 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.  */
+/* 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. 
+   We use the processor's default collating sequence.  */
 
 int
-gfc_compare_string (gfc_expr * a, gfc_expr * b, const int *xcoll_table)
+gfc_compare_string (gfc_expr *a, gfc_expr *b)
 {
-  int len, alen, blen, i, ac, bc;
+  int len, alen, blen, i;
+  gfc_char_t ac, bc;
 
   alen = a->value.character.length;
   blen = b->value.character.length;
 
-  len = (alen > blen) ? alen : blen;
+  len = MAX(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 (ac < bc)
+       return -1;
+      if (ac > bc)
+       return 1;
+    }
+
+  /* Strings are equal */
+  return 0;
+}
+
+
+int
+gfc_compare_with_Cstring (gfc_expr *a, const char *b, bool case_sensitive)
+{
+  int len, alen, blen, i;
+  gfc_char_t ac, bc;
+
+  alen = a->value.character.length;
+  blen = strlen (b);
+
+  len = MAX(alen, blen);
 
   for (i = 0; i < len; i++)
     {
-      ac = (i < alen) ? a->value.character.string[i] : ' ';
-      bc = (i < blen) ? b->value.character.string[i] : ' ';
+      ac = ((i < alen) ? a->value.character.string[i] : ' ');
+      bc = ((i < blen) ? b[i] : ' ');
 
-      if (xcoll_table != NULL)
+      if (!case_sensitive)
        {
-         ac = xcoll_table[ac];
-         bc = xcoll_table[bc];
+         ac = TOLOWER (ac);
+         bc = TOLOWER (bc);
        }
 
       if (ac < bc)
@@ -1761,7 +1235,6 @@ gfc_compare_string (gfc_expr * a, gfc_expr * b, const int *xcoll_table)
     }
 
   /* Strings are equal */
-
   return 0;
 }
 
@@ -1769,14 +1242,15 @@ gfc_compare_string (gfc_expr * a, gfc_expr * b, const int *xcoll_table)
 /* Specific comparison subroutines.  */
 
 static arith
-gfc_arith_eq (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
+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 (),
+  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);
+  result->value.logical = (op1->ts.type == BT_COMPLEX)
+                       ? compare_complex (op1, op2)
+                       : (gfc_compare_expr (op1, op2, INTRINSIC_EQ) == 0);
 
   *resultp = result;
   return ARITH_OK;
@@ -1784,14 +1258,15 @@ gfc_arith_eq (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
 
 
 static arith
-gfc_arith_ne (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
+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 (),
+  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);
+  result->value.logical = (op1->ts.type == BT_COMPLEX)
+                       ? !compare_complex (op1, op2)
+                       : (gfc_compare_expr (op1, op2, INTRINSIC_EQ) != 0);
 
   *resultp = result;
   return ARITH_OK;
@@ -1799,13 +1274,13 @@ gfc_arith_ne (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
 
 
 static arith
-gfc_arith_gt (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
+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 (),
+  result = gfc_constant_result (BT_LOGICAL, gfc_default_logical_kind,
                                &op1->where);
-  result->value.logical = (gfc_compare_expr (op1, op2) > 0);
+  result->value.logical = (gfc_compare_expr (op1, op2, INTRINSIC_GT) > 0);
   *resultp = result;
 
   return ARITH_OK;
@@ -1813,13 +1288,13 @@ gfc_arith_gt (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
 
 
 static arith
-gfc_arith_ge (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
+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 (),
+  result = gfc_constant_result (BT_LOGICAL, gfc_default_logical_kind,
                                &op1->where);
-  result->value.logical = (gfc_compare_expr (op1, op2) >= 0);
+  result->value.logical = (gfc_compare_expr (op1, op2, INTRINSIC_GE) >= 0);
   *resultp = result;
 
   return ARITH_OK;
@@ -1827,13 +1302,13 @@ gfc_arith_ge (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
 
 
 static arith
-gfc_arith_lt (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
+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 (),
+  result = gfc_constant_result (BT_LOGICAL, gfc_default_logical_kind,
                                &op1->where);
-  result->value.logical = (gfc_compare_expr (op1, op2) < 0);
+  result->value.logical = (gfc_compare_expr (op1, op2, INTRINSIC_LT) < 0);
   *resultp = result;
 
   return ARITH_OK;
@@ -1841,13 +1316,13 @@ gfc_arith_lt (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
 
 
 static arith
-gfc_arith_le (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
+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 (),
+  result = gfc_constant_result (BT_LOGICAL, gfc_default_logical_kind,
                                &op1->where);
-  result->value.logical = (gfc_compare_expr (op1, op2) <= 0);
+  result->value.logical = (gfc_compare_expr (op1, op2, INTRINSIC_LE) <= 0);
   *resultp = result;
 
   return ARITH_OK;
@@ -1855,8 +1330,8 @@ gfc_arith_le (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
 
 
 static arith
-reduce_unary (arith (*eval) (gfc_expr *, gfc_expr **), gfc_expr * op,
-             gfc_expr ** result)
+reduce_unary (arith (*eval) (gfc_expr *, gfc_expr **), gfc_expr *op,
+             gfc_expr **result)
 {
   gfc_constructor *c, *head;
   gfc_expr *r;
@@ -1870,7 +1345,8 @@ reduce_unary (arith (*eval) (gfc_expr *, gfc_expr **), gfc_expr * op,
 
   for (c = head; c; c = c->next)
     {
-      rc = eval (c->expr, &r);
+      rc = reduce_unary (eval, c->expr, &r);
+
       if (rc != ARITH_OK)
        break;
 
@@ -1899,8 +1375,7 @@ reduce_unary (arith (*eval) (gfc_expr *, gfc_expr **), gfc_expr * op,
 
 static arith
 reduce_binary_ac (arith (*eval) (gfc_expr *, gfc_expr *, gfc_expr **),
-                 gfc_expr * op1, gfc_expr * op2,
-                 gfc_expr ** result)
+                 gfc_expr *op1, gfc_expr *op2, gfc_expr **result)
 {
   gfc_constructor *c, *head;
   gfc_expr *r;
@@ -1911,7 +1386,11 @@ reduce_binary_ac (arith (*eval) (gfc_expr *, gfc_expr *, gfc_expr **),
 
   for (c = head; c; c = c->next)
     {
-      rc = eval (c->expr, op2, &r);
+      if (c->expr->expr_type == EXPR_CONSTANT)
+        rc = eval (c->expr, op2, &r);
+      else
+       rc = reduce_binary_ac (eval, c->expr, op2, &r);
+
       if (rc != ARITH_OK)
        break;
 
@@ -1940,8 +1419,7 @@ reduce_binary_ac (arith (*eval) (gfc_expr *, gfc_expr *, gfc_expr **),
 
 static arith
 reduce_binary_ca (arith (*eval) (gfc_expr *, gfc_expr *, gfc_expr **),
-                 gfc_expr * op1, gfc_expr * op2,
-                 gfc_expr ** result)
+                 gfc_expr *op1, gfc_expr *op2, gfc_expr **result)
 {
   gfc_constructor *c, *head;
   gfc_expr *r;
@@ -1952,7 +1430,11 @@ reduce_binary_ca (arith (*eval) (gfc_expr *, gfc_expr *, gfc_expr **),
 
   for (c = head; c; c = c->next)
     {
-      rc = eval (op1, c->expr, &r);
+      if (c->expr->expr_type == EXPR_CONSTANT)
+       rc = eval (op1, c->expr, &r);
+      else
+       rc = reduce_binary_ca (eval, op1, c->expr, &r);
+
       if (rc != ARITH_OK)
        break;
 
@@ -1979,10 +1461,14 @@ reduce_binary_ca (arith (*eval) (gfc_expr *, gfc_expr *, gfc_expr **),
 }
 
 
+/* We need a forward declaration of reduce_binary.  */
+static arith reduce_binary (arith (*eval) (gfc_expr *, gfc_expr *, gfc_expr **),
+                           gfc_expr *op1, gfc_expr *op2, gfc_expr **result);
+
+
 static arith
 reduce_binary_aa (arith (*eval) (gfc_expr *, gfc_expr *, gfc_expr **),
-                 gfc_expr * op1, gfc_expr * op2,
-                 gfc_expr ** result)
+                 gfc_expr *op1, gfc_expr *op2, gfc_expr **result)
 {
   gfc_constructor *c, *d, *head;
   gfc_expr *r;
@@ -1993,12 +1479,11 @@ reduce_binary_aa (arith (*eval) (gfc_expr *, gfc_expr *, gfc_expr **),
   rc = ARITH_OK;
   d = op2->value.constructor;
 
-  if (gfc_check_conformance ("Elemental binary operation", op1, op2)
+  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)
@@ -2007,7 +1492,7 @@ reduce_binary_aa (arith (*eval) (gfc_expr *, gfc_expr *, gfc_expr **),
              break;
            }
 
-         rc = eval (c->expr, d->expr, &r);
+         rc = reduce_binary (eval, c->expr, d->expr, &r);
          if (rc != ARITH_OK)
            break;
 
@@ -2040,10 +1525,8 @@ reduce_binary_aa (arith (*eval) (gfc_expr *, gfc_expr *, gfc_expr **),
 
 static arith
 reduce_binary (arith (*eval) (gfc_expr *, gfc_expr *, gfc_expr **),
-              gfc_expr * op1, gfc_expr * op2,
-              gfc_expr ** result)
+              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);
 
@@ -2075,8 +1558,8 @@ eval_f;
    operands are array constructors.  */
 
 static gfc_expr *
-eval_intrinsic (gfc_intrinsic_op operator,
-               eval_f eval, gfc_expr * op1, gfc_expr * op2)
+eval_intrinsic (gfc_intrinsic_op op,
+               eval_f eval, gfc_expr *op1, gfc_expr *op2)
 {
   gfc_expr temp, *result;
   int unary;
@@ -2084,19 +1567,19 @@ eval_intrinsic (gfc_intrinsic_op operator,
 
   gfc_clear_ts (&temp.ts);
 
-  switch (operator)
+  switch (op)
     {
-    case INTRINSIC_NOT:        /* Logical unary */
+    /* Logical unary  */
+    case INTRINSIC_NOT:
       if (op1->ts.type != BT_LOGICAL)
        goto runtime;
 
       temp.ts.type = BT_LOGICAL;
-      temp.ts.kind = gfc_default_logical_kind ();
-
+      temp.ts.kind = gfc_default_logical_kind;
       unary = 1;
       break;
 
-      /* Logical binary operators */
+    /* Logical binary operators  */
     case INTRINSIC_OR:
     case INTRINSIC_AND:
     case INTRINSIC_NEQV:
@@ -2105,83 +1588,103 @@ eval_intrinsic (gfc_intrinsic_op operator,
        goto runtime;
 
       temp.ts.type = BT_LOGICAL;
-      temp.ts.kind = gfc_default_logical_kind ();
-
+      temp.ts.kind = gfc_default_logical_kind;
       unary = 0;
       break;
 
+    /* Numeric unary  */
     case INTRINSIC_UPLUS:
-    case INTRINSIC_UMINUS:     /* Numeric unary */
+    case INTRINSIC_UMINUS:
       if (!gfc_numeric_ts (&op1->ts))
        goto runtime;
 
       temp.ts = op1->ts;
+      unary = 1;
+      break;
 
+    case INTRINSIC_PARENTHESES:
+      temp.ts = op1->ts;
       unary = 1;
       break;
 
+    /* Additional restrictions for ordering relations.  */
     case INTRINSIC_GE:
-    case INTRINSIC_LT:         /* Additional restrictions  */
-    case INTRINSIC_LE:          /* for ordering relations.  */
+    case INTRINSIC_GE_OS:
+    case INTRINSIC_LT:
+    case INTRINSIC_LT_OS:
+    case INTRINSIC_LE:
+    case INTRINSIC_LE_OS:
     case INTRINSIC_GT:
+    case INTRINSIC_GT_OS:
       if (op1->ts.type == BT_COMPLEX || op2->ts.type == BT_COMPLEX)
        {
          temp.ts.type = BT_LOGICAL;
-         temp.ts.kind = gfc_default_logical_kind();
+         temp.ts.kind = gfc_default_logical_kind;
          goto runtime;
        }
 
-      /* else fall through */
-
+    /* Fall through  */
     case INTRINSIC_EQ:
+    case INTRINSIC_EQ_OS:
     case INTRINSIC_NE:
+    case INTRINSIC_NE_OS:
       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();
+         temp.ts.kind = gfc_default_logical_kind;
+
+         /* If kind mismatch, exit and we'll error out later.  */
+         if (op1->ts.kind != op2->ts.kind)
+           goto runtime;
+
          break;
        }
 
-      /* else fall through */
-
+    /* Fall through  */
+    /* Numeric binary  */
     case INTRINSIC_PLUS:
     case INTRINSIC_MINUS:
     case INTRINSIC_TIMES:
     case INTRINSIC_DIVIDE:
-    case INTRINSIC_POWER:      /* Numeric binary */
+    case INTRINSIC_POWER:
       if (!gfc_numeric_ts (&op1->ts) || !gfc_numeric_ts (&op2->ts))
        goto runtime;
 
-      /* Insert any necessary type conversions to make the operands compatible.  */
+      /* Insert any necessary type conversions to make the operands
+        compatible.  */
 
       temp.expr_type = EXPR_OP;
       gfc_clear_ts (&temp.ts);
-      temp.operator = operator;
+      temp.value.op.op = op;
 
-      temp.op1 = op1;
-      temp.op2 = op2;
+      temp.value.op.op1 = op1;
+      temp.value.op.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)
+      if (op == INTRINSIC_EQ || op == INTRINSIC_NE
+         || op == INTRINSIC_GE || op == INTRINSIC_GT
+         || op == INTRINSIC_LE || op == INTRINSIC_LT
+         || op == INTRINSIC_EQ_OS || op == INTRINSIC_NE_OS
+         || op == INTRINSIC_GE_OS || op == INTRINSIC_GT_OS
+         || op == INTRINSIC_LE_OS || op == INTRINSIC_LT_OS)
        {
          temp.ts.type = BT_LOGICAL;
-         temp.ts.kind = gfc_default_logical_kind ();
+         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)
+    /* Character binary  */
+    case INTRINSIC_CONCAT:
+      if (op1->ts.type != BT_CHARACTER || op2->ts.type != BT_CHARACTER
+         || op1->ts.kind != op2->ts.kind)
        goto runtime;
 
       temp.ts.type = BT_CHARACTER;
-      temp.ts.kind = gfc_default_character_kind ();
-
+      temp.ts.kind = op1->ts.kind;
       unary = 0;
       break;
 
@@ -2193,20 +1696,18 @@ eval_intrinsic (gfc_intrinsic_op operator,
     }
 
   /* Try to combine the operators.  */
-  if (operator == INTRINSIC_POWER && op2->ts.type != BT_INTEGER)
+  if (op == 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)))
+         || !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)))
+        && (op2->expr_type != EXPR_ARRAY
+            || !gfc_is_constant_expr (op2) || !gfc_expanded_ac (op2)))
     goto runtime;
 
   if (unary)
@@ -2215,8 +1716,8 @@ eval_intrinsic (gfc_intrinsic_op operator,
     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);
+    { /* Something went wrong.  */
+      gfc_error (gfc_arith_error (rc), &op1->where);
       return NULL;
     }
 
@@ -2225,15 +1726,15 @@ eval_intrinsic (gfc_intrinsic_op operator,
   return result;
 
 runtime:
-  /* Create a run-time expression */
+  /* Create a run-time expression */
   result = gfc_get_expr ();
   result->ts = temp.ts;
 
   result->expr_type = EXPR_OP;
-  result->operator = operator;
+  result->value.op.op = op;
 
-  result->op1 = op1;
-  result->op2 = op2;
+  result->value.op.op1 = op1;
+  result->value.op.op2 = op2;
 
   result->where = op1->where;
 
@@ -2242,22 +1743,29 @@ runtime:
 
 
 /* Modify type of expression for zero size array.  */
+
 static gfc_expr *
-eval_type_intrinsic0 (gfc_intrinsic_op operator, gfc_expr *op)
+eval_type_intrinsic0 (gfc_intrinsic_op iop, gfc_expr *op)
 {
   if (op == NULL)
-    gfc_internal_error("eval_type_intrinsic0(): op NULL");
+    gfc_internal_error ("eval_type_intrinsic0(): op NULL");
 
-  switch(operator)
+  switch (iop)
     {
     case INTRINSIC_GE:
+    case INTRINSIC_GE_OS:
     case INTRINSIC_LT:
+    case INTRINSIC_LT_OS:
     case INTRINSIC_LE:
+    case INTRINSIC_LE_OS:
     case INTRINSIC_GT:
+    case INTRINSIC_GT_OS:
     case INTRINSIC_EQ:
+    case INTRINSIC_EQ_OS:
     case INTRINSIC_NE:
+    case INTRINSIC_NE_OS:
       op->ts.type = BT_LOGICAL;
-      op->ts.kind = gfc_default_logical_kind();
+      op->ts.kind = gfc_default_logical_kind;
       break;
 
     default:
@@ -2271,9 +1779,8 @@ eval_type_intrinsic0 (gfc_intrinsic_op operator, gfc_expr *op)
 /* Return nonzero if the expression is a zero size array.  */
 
 static int
-gfc_zero_size_array (gfc_expr * e)
+gfc_zero_size_array (gfc_expr *e)
 {
-
   if (e->expr_type != EXPR_ARRAY)
     return 0;
 
@@ -2286,9 +1793,8 @@ gfc_zero_size_array (gfc_expr * e)
    operands is a zero-length array.  */
 
 static gfc_expr *
-reduce_binary0 (gfc_expr * op1, gfc_expr * op2)
+reduce_binary0 (gfc_expr *op1, gfc_expr *op2)
 {
-
   if (gfc_zero_size_array (op1))
     {
       gfc_free_expr (op2);
@@ -2306,9 +1812,9 @@ reduce_binary0 (gfc_expr * op1, gfc_expr * op2)
 
 
 static gfc_expr *
-eval_intrinsic_f2 (gfc_intrinsic_op operator,
+eval_intrinsic_f2 (gfc_intrinsic_op op,
                   arith (*eval) (gfc_expr *, gfc_expr **),
-                  gfc_expr * op1, gfc_expr * op2)
+                  gfc_expr *op1, gfc_expr *op2)
 {
   gfc_expr *result;
   eval_f f;
@@ -2316,163 +1822,190 @@ eval_intrinsic_f2 (gfc_intrinsic_op operator,
   if (op2 == NULL)
     {
       if (gfc_zero_size_array (op1))
-       return eval_type_intrinsic0(operator, op1);
+       return eval_type_intrinsic0 (op, op1);
     }
   else
     {
       result = reduce_binary0 (op1, op2);
       if (result != NULL)
-       return eval_type_intrinsic0(operator, result);
+       return eval_type_intrinsic0 (op, result);
     }
 
   f.f2 = eval;
-  return eval_intrinsic (operator, f, op1, op2);
+  return eval_intrinsic (op, f, op1, op2);
 }
 
 
 static gfc_expr *
-eval_intrinsic_f3 (gfc_intrinsic_op operator,
+eval_intrinsic_f3 (gfc_intrinsic_op op,
                   arith (*eval) (gfc_expr *, gfc_expr *, gfc_expr **),
-                  gfc_expr * op1, gfc_expr * op2)
+                  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);
+    return eval_type_intrinsic0(op, result);
 
   f.f3 = eval;
-  return eval_intrinsic (operator, f, op1, op2);
+  return eval_intrinsic (op, f, op1, op2);
 }
 
 
+gfc_expr *
+gfc_parentheses (gfc_expr *op)
+{
+  if (gfc_is_constant_expr (op))
+    return op;
+
+  return eval_intrinsic_f2 (INTRINSIC_PARENTHESES, gfc_arith_identity,
+                           op, NULL);
+}
 
 gfc_expr *
-gfc_uplus (gfc_expr * op)
+gfc_uplus (gfc_expr *op)
 {
-  return eval_intrinsic_f2 (INTRINSIC_UPLUS, gfc_arith_uplus, op, NULL);
+  return eval_intrinsic_f2 (INTRINSIC_UPLUS, gfc_arith_identity, op, NULL);
 }
 
+
 gfc_expr *
-gfc_uminus (gfc_expr * op)
+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)
+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)
+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)
+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)
+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)
+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)
+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)
+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)
+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)
+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)
+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)
+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)
+gfc_eq (gfc_expr *op1, gfc_expr *op2, gfc_intrinsic_op op)
 {
-  return eval_intrinsic_f3 (INTRINSIC_EQ, gfc_arith_eq, op1, op2);
+  return eval_intrinsic_f3 (op, gfc_arith_eq, op1, op2);
 }
 
+
 gfc_expr *
-gfc_ne (gfc_expr * op1, gfc_expr * op2)
+gfc_ne (gfc_expr *op1, gfc_expr *op2, gfc_intrinsic_op op)
 {
-  return eval_intrinsic_f3 (INTRINSIC_NE, gfc_arith_ne, op1, op2);
+  return eval_intrinsic_f3 (op, gfc_arith_ne, op1, op2);
 }
 
+
 gfc_expr *
-gfc_gt (gfc_expr * op1, gfc_expr * op2)
+gfc_gt (gfc_expr *op1, gfc_expr *op2, gfc_intrinsic_op op)
 {
-  return eval_intrinsic_f3 (INTRINSIC_GT, gfc_arith_gt, op1, op2);
+  return eval_intrinsic_f3 (op, gfc_arith_gt, op1, op2);
 }
 
+
 gfc_expr *
-gfc_ge (gfc_expr * op1, gfc_expr * op2)
+gfc_ge (gfc_expr *op1, gfc_expr *op2, gfc_intrinsic_op op)
 {
-  return eval_intrinsic_f3 (INTRINSIC_GE, gfc_arith_ge, op1, op2);
+  return eval_intrinsic_f3 (op, gfc_arith_ge, op1, op2);
 }
 
+
 gfc_expr *
-gfc_lt (gfc_expr * op1, gfc_expr * op2)
+gfc_lt (gfc_expr *op1, gfc_expr *op2, gfc_intrinsic_op op)
 {
-  return eval_intrinsic_f3 (INTRINSIC_LT, gfc_arith_lt, op1, op2);
+  return eval_intrinsic_f3 (op, gfc_arith_lt, op1, op2);
 }
 
+
 gfc_expr *
-gfc_le (gfc_expr * op1, gfc_expr * op2)
+gfc_le (gfc_expr *op1, gfc_expr *op2, gfc_intrinsic_op op)
 {
-  return eval_intrinsic_f3 (INTRINSIC_LE, gfc_arith_le, op1, op2);
+  return eval_intrinsic_f3 (op, 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_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 */
+  /* A leading plus is allowed, but not by mpz_set_str.  */
   if (buffer[0] == '+')
     t = buffer + 1;
   else
@@ -2486,18 +2019,12 @@ gfc_convert_integer (const char *buffer, int kind, int radix, locus * where)
 /* Convert a real string to an expression node.  */
 
 gfc_expr *
-gfc_convert_real (const char *buffer, int kind, locus * where)
+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);
+  mpfr_set_str (e->value.real, buffer, 10, GFC_RND_MODE);
 
   return e;
 }
@@ -2507,13 +2034,13 @@ gfc_convert_real (const char *buffer, int kind, locus * where)
    complex expression node.  */
 
 gfc_expr *
-gfc_convert_complex (gfc_expr * real, gfc_expr * imag, int kind)
+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;
 }
@@ -2525,20 +2052,55 @@ gfc_convert_complex (gfc_expr * real, gfc_expr * imag, int kind)
 /* Deal with an arithmetic error.  */
 
 static void
-arith_error (arith rc, gfc_typespec * from, gfc_typespec * to, locus * where)
+arith_error (arith rc, gfc_typespec *from, gfc_typespec *to, locus *where)
 {
+  switch (rc)
+    {
+    case ARITH_OK:
+      gfc_error ("Arithmetic OK converting %s to %s at %L",
+                gfc_typename (from), gfc_typename (to), where);
+      break;
+    case ARITH_OVERFLOW:
+      gfc_error ("Arithmetic overflow converting %s to %s at %L. This check "
+                "can be disabled with the option -fno-range-check",
+                gfc_typename (from), gfc_typename (to), where);
+      break;
+    case ARITH_UNDERFLOW:
+      gfc_error ("Arithmetic underflow converting %s to %s at %L. This check "
+                "can be disabled with the option -fno-range-check",
+                gfc_typename (from), gfc_typename (to), where);
+      break;
+    case ARITH_NAN:
+      gfc_error ("Arithmetic NaN converting %s to %s at %L. This check "
+                "can be disabled with the option -fno-range-check",
+                gfc_typename (from), gfc_typename (to), where);
+      break;
+    case ARITH_DIV0:
+      gfc_error ("Division by zero converting %s to %s at %L",
+                gfc_typename (from), gfc_typename (to), where);
+      break;
+    case ARITH_INCOMMENSURATE:
+      gfc_error ("Array operands are incommensurate converting %s to %s at %L",
+                gfc_typename (from), gfc_typename (to), where);
+      break;
+    case ARITH_ASYMMETRIC:
+      gfc_error ("Integer outside symmetric range implied by Standard Fortran"
+                " converting %s to %s at %L",
+                gfc_typename (from), gfc_typename (to), where);
+      break;
+    default:
+      gfc_internal_error ("gfc_arith_error(): Bad error code");
+    }
 
-  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.  */
+  /* TODO: Do something about the error, i.e., throw exception, return
+     NaN, etc.  */
 }
 
+
 /* Convert integers to integers.  */
 
 gfc_expr *
-gfc_int2int (gfc_expr * src, int kind)
+gfc_int2int (gfc_expr *src, int kind)
 {
   gfc_expr *result;
   arith rc;
@@ -2547,12 +2109,18 @@ gfc_int2int (gfc_expr * src, int kind)
 
   mpz_set (result->value.integer, src->value.integer);
 
-  if ((rc = gfc_check_integer_range (result->value.integer, kind))
-      != ARITH_OK)
+  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;
+      if (rc == ARITH_ASYMMETRIC)
+       {
+         gfc_warning (gfc_arith_error (rc), &src->where);
+       }
+      else
+       {
+         arith_error (rc, &src->ts, &result->ts, &src->where);
+         gfc_free_expr (result);
+         return NULL;
+       }
     }
 
   return result;
@@ -2562,14 +2130,14 @@ gfc_int2int (gfc_expr * src, int kind)
 /* Convert integers to reals.  */
 
 gfc_expr *
-gfc_int2real (gfc_expr * src, int kind)
+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);
+  mpfr_set_z (result->value.real, src->value.integer, GFC_RND_MODE);
 
   if ((rc = gfc_check_real_range (result->value.real, kind)) != ARITH_OK)
     {
@@ -2585,15 +2153,15 @@ gfc_int2real (gfc_expr * src, int kind)
 /* Convert default integer to default complex.  */
 
 gfc_expr *
-gfc_int2complex (gfc_expr * src, int kind)
+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);
+  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)
     {
@@ -2609,17 +2177,16 @@ gfc_int2complex (gfc_expr * src, int kind)
 /* Convert default real to default integer.  */
 
 gfc_expr *
-gfc_real2int (gfc_expr * src, int kind)
+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);
+  gfc_mpfr_to_mpz (result->value.integer, src->value.real, &src->where);
 
-  if ((rc = gfc_check_integer_range (result->value.integer, kind))
-      != ARITH_OK)
+  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);
@@ -2633,16 +2200,24 @@ gfc_real2int (gfc_expr * src, int kind)
 /* Convert real to real.  */
 
 gfc_expr *
-gfc_real2real (gfc_expr * src, int kind)
+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);
+  mpfr_set (result->value.real, src->value.real, GFC_RND_MODE);
 
-  if ((rc = gfc_check_real_range (result->value.real, kind)) != ARITH_OK)
+  rc = gfc_check_real_range (result->value.real, kind);
+
+  if (rc == ARITH_UNDERFLOW)
+    {
+      if (gfc_option.warn_underflow)
+       gfc_warning (gfc_arith_error (rc), &src->where);
+      mpfr_set_ui (result->value.real, 0, GFC_RND_MODE);
+    }
+  else if (rc != ARITH_OK)
     {
       arith_error (rc, &src->ts, &result->ts, &src->where);
       gfc_free_expr (result);
@@ -2656,17 +2231,25 @@ gfc_real2real (gfc_expr * src, int kind)
 /* Convert real to complex.  */
 
 gfc_expr *
-gfc_real2complex (gfc_expr * src, int kind)
+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);
+  mpfr_set (result->value.complex.r, src->value.real, 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)
+  rc = gfc_check_real_range (result->value.complex.r, kind);
+
+  if (rc == ARITH_UNDERFLOW)
+    {
+      if (gfc_option.warn_underflow)
+       gfc_warning (gfc_arith_error (rc), &src->where);
+      mpfr_set_ui (result->value.complex.r, 0, GFC_RND_MODE);
+    }
+  else if (rc != ARITH_OK)
     {
       arith_error (rc, &src->ts, &result->ts, &src->where);
       gfc_free_expr (result);
@@ -2680,17 +2263,16 @@ gfc_real2complex (gfc_expr * src, int kind)
 /* Convert complex to integer.  */
 
 gfc_expr *
-gfc_complex2int (gfc_expr * src, int kind)
+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);
+  gfc_mpfr_to_mpz (result->value.integer, src->value.complex.r, &src->where);
 
-  if ((rc = gfc_check_integer_range (result->value.integer, kind))
-      != ARITH_OK)
+  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);
@@ -2704,16 +2286,24 @@ gfc_complex2int (gfc_expr * src, int kind)
 /* Convert complex to real.  */
 
 gfc_expr *
-gfc_complex2real (gfc_expr * src, int kind)
+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);
+  mpfr_set (result->value.real, src->value.complex.r, GFC_RND_MODE);
 
-  if ((rc = gfc_check_real_range (result->value.real, kind)) != ARITH_OK)
+  rc = gfc_check_real_range (result->value.real, kind);
+
+  if (rc == ARITH_UNDERFLOW)
+    {
+      if (gfc_option.warn_underflow)
+       gfc_warning (gfc_arith_error (rc), &src->where);
+      mpfr_set_ui (result->value.real, 0, GFC_RND_MODE);
+    }
+  if (rc != ARITH_OK)
     {
       arith_error (rc, &src->ts, &result->ts, &src->where);
       gfc_free_expr (result);
@@ -2727,19 +2317,40 @@ gfc_complex2real (gfc_expr * src, int kind)
 /* Convert complex to complex.  */
 
 gfc_expr *
-gfc_complex2complex (gfc_expr * src, int kind)
+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);
+  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);
 
-  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)
+  rc = gfc_check_real_range (result->value.complex.r, kind);
+
+  if (rc == ARITH_UNDERFLOW)
+    {
+      if (gfc_option.warn_underflow)
+       gfc_warning (gfc_arith_error (rc), &src->where);
+      mpfr_set_ui (result->value.complex.r, 0, GFC_RND_MODE);
+    }
+  else if (rc != ARITH_OK)
+    {
+      arith_error (rc, &src->ts, &result->ts, &src->where);
+      gfc_free_expr (result);
+      return NULL;
+    }
+
+  rc = gfc_check_real_range (result->value.complex.i, kind);
+
+  if (rc == ARITH_UNDERFLOW)
+    {
+      if (gfc_option.warn_underflow)
+       gfc_warning (gfc_arith_error (rc), &src->where);
+      mpfr_set_ui (result->value.complex.i, 0, GFC_RND_MODE);
+    }
+  else if (rc != ARITH_OK)
     {
       arith_error (rc, &src->ts, &result->ts, &src->where);
       gfc_free_expr (result);
@@ -2753,7 +2364,7 @@ gfc_complex2complex (gfc_expr * src, int kind)
 /* Logical kind conversion.  */
 
 gfc_expr *
-gfc_log2log (gfc_expr * src, int kind)
+gfc_log2log (gfc_expr *src, int kind)
 {
   gfc_expr *result;
 
@@ -2762,3 +2373,219 @@ gfc_log2log (gfc_expr * src, int kind)
 
   return result;
 }
+
+
+/* Convert logical to integer.  */
+
+gfc_expr *
+gfc_log2int (gfc_expr *src, int kind)
+{
+  gfc_expr *result;
+
+  result = gfc_constant_result (BT_INTEGER, kind, &src->where);
+  mpz_set_si (result->value.integer, src->value.logical);
+
+  return result;
+}
+
+
+/* Convert integer to logical.  */
+
+gfc_expr *
+gfc_int2log (gfc_expr *src, int kind)
+{
+  gfc_expr *result;
+
+  result = gfc_constant_result (BT_LOGICAL, kind, &src->where);
+  result->value.logical = (mpz_cmp_si (src->value.integer, 0) != 0);
+
+  return result;
+}
+
+
+/* Helper function to set the representation in a Hollerith conversion.  
+   This assumes that the ts.type and ts.kind of the result have already
+   been set.  */
+
+static void
+hollerith2representation (gfc_expr *result, gfc_expr *src)
+{
+  int src_len, result_len;
+
+  src_len = src->representation.length;
+  result_len = gfc_target_expr_size (result);
+
+  if (src_len > result_len)
+    {
+      gfc_warning ("The Hollerith constant at %L is too long to convert to %s",
+                  &src->where, gfc_typename(&result->ts));
+    }
+
+  result->representation.string = XCNEWVEC (char, result_len + 1);
+  memcpy (result->representation.string, src->representation.string,
+         MIN (result_len, src_len));
+
+  if (src_len < result_len)
+    memset (&result->representation.string[src_len], ' ', result_len - src_len);
+
+  result->representation.string[result_len] = '\0'; /* For debugger  */
+  result->representation.length = result_len;
+}
+
+
+/* Convert Hollerith to integer. The constant will be padded or truncated.  */
+
+gfc_expr *
+gfc_hollerith2int (gfc_expr *src, int kind)
+{
+  gfc_expr *result;
+
+  result = gfc_get_expr ();
+  result->expr_type = EXPR_CONSTANT;
+  result->ts.type = BT_INTEGER;
+  result->ts.kind = kind;
+  result->where = src->where;
+
+  hollerith2representation (result, src);
+  gfc_interpret_integer (kind, (unsigned char *) result->representation.string,
+                        result->representation.length, result->value.integer);
+
+  return result;
+}
+
+
+/* Convert Hollerith to real. The constant will be padded or truncated.  */
+
+gfc_expr *
+gfc_hollerith2real (gfc_expr *src, int kind)
+{
+  gfc_expr *result;
+  int len;
+
+  len = src->value.character.length;
+
+  result = gfc_get_expr ();
+  result->expr_type = EXPR_CONSTANT;
+  result->ts.type = BT_REAL;
+  result->ts.kind = kind;
+  result->where = src->where;
+
+  hollerith2representation (result, src);
+  gfc_interpret_float (kind, (unsigned char *) result->representation.string,
+                      result->representation.length, result->value.real);
+
+  return result;
+}
+
+
+/* Convert Hollerith to complex. The constant will be padded or truncated.  */
+
+gfc_expr *
+gfc_hollerith2complex (gfc_expr *src, int kind)
+{
+  gfc_expr *result;
+  int len;
+
+  len = src->value.character.length;
+
+  result = gfc_get_expr ();
+  result->expr_type = EXPR_CONSTANT;
+  result->ts.type = BT_COMPLEX;
+  result->ts.kind = kind;
+  result->where = src->where;
+
+  hollerith2representation (result, src);
+  gfc_interpret_complex (kind, (unsigned char *) result->representation.string,
+                        result->representation.length, result->value.complex.r,
+                        result->value.complex.i);
+
+  return result;
+}
+
+
+/* Convert Hollerith to character. */
+
+gfc_expr *
+gfc_hollerith2character (gfc_expr *src, int kind)
+{
+  gfc_expr *result;
+
+  result = gfc_copy_expr (src);
+  result->ts.type = BT_CHARACTER;
+  result->ts.kind = kind;
+
+  result->value.character.length = result->representation.length;
+  result->value.character.string
+    = gfc_char_to_widechar (result->representation.string);
+
+  return result;
+}
+
+
+/* Convert Hollerith to logical. The constant will be padded or truncated.  */
+
+gfc_expr *
+gfc_hollerith2logical (gfc_expr *src, int kind)
+{
+  gfc_expr *result;
+  int len;
+
+  len = src->value.character.length;
+
+  result = gfc_get_expr ();
+  result->expr_type = EXPR_CONSTANT;
+  result->ts.type = BT_LOGICAL;
+  result->ts.kind = kind;
+  result->where = src->where;
+
+  hollerith2representation (result, src);
+  gfc_interpret_logical (kind, (unsigned char *) result->representation.string,
+                        result->representation.length, &result->value.logical);
+
+  return result;
+}
+
+
+/* Returns an initializer whose value is one higher than the value of the
+   LAST_INITIALIZER argument.  If the argument is NULL, the
+   initializers value will be set to zero.  The initializer's kind
+   will be set to gfc_c_int_kind.
+
+   If -fshort-enums is given, the appropriate kind will be selected
+   later after all enumerators have been parsed.  A warning is issued
+   here if an initializer exceeds gfc_c_int_kind.  */
+
+gfc_expr *
+gfc_enum_initializer (gfc_expr *last_initializer, locus where)
+{
+  gfc_expr *result;
+
+  result = gfc_get_expr ();
+  result->expr_type = EXPR_CONSTANT;
+  result->ts.type = BT_INTEGER;
+  result->ts.kind = gfc_c_int_kind;
+  result->where = where;
+
+  mpz_init (result->value.integer);
+
+  if (last_initializer != NULL)
+    {
+      mpz_add_ui (result->value.integer, last_initializer->value.integer, 1);
+      result->where = last_initializer->where;
+
+      if (gfc_check_integer_range (result->value.integer,
+            gfc_c_int_kind) != ARITH_OK)
+       {
+         gfc_error ("Enumerator exceeds the C integer type at %C");
+         return NULL;
+       }
+    }
+  else
+    {
+      /* Control comes here, if it's the very first enumerator and no
+        initializer has been given.  It will be initialized to zero.  */
+      mpz_set_si (result->value.integer, 0);
+    }
+
+  return result;
+}