OSDN Git Service

2008-11-24 Jerry DeLisle <jvdelisle@gcc.gnu.org>
[pf3gnuchains/gcc-fork.git] / gcc / fortran / arith.c
index 684ae7b..7440be3 100644 (file)
@@ -1,13 +1,13 @@
 /* Compiler arithmetic
-   Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005 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,34 +16,38 @@ 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 "system.h"
 #include "flags.h"
 #include "gfortran.h"
 #include "arith.h"
+#include "target-memory.h"
 
 /* MPFR does not have a direct replacement for mpz_set_f() from GMP.
    It's easily implemented with a few calls though.  */
 
 void
-gfc_mpfr_to_mpz (mpz_t z, mpfr_t x)
+gfc_mpfr_to_mpz (mpz_t z, mpfr_t x, locus *where)
 {
   mp_exp_t e;
 
+  if (mpfr_inf_p (x) || mpfr_nan_p (x))
+    {
+      gfc_error ("Conversion of an Infinity or Not-a-Number at %L "
+                "to INTEGER", where);
+      mpz_set_ui (z, 0);
+      return;
+    }
+
   e = mpfr_get_z_exp (z, x);
-  /* MPFR 2.0.1 (included with GMP 4.1) has a bug whereby mpfr_get_z_exp
-     may set the sign of z incorrectly.  Work around that here.  */
-  if (mpfr_sgn (x) != mpz_sgn (z))
-    mpz_neg (z, z);
 
   if (e > 0)
     mpz_mul_2exp (z, z, e);
@@ -75,57 +79,6 @@ gfc_set_model (mpfr_t x)
   mpfr_set_default_prec (mpfr_get_prec (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 (mpfr_t y, mpfr_t x, mpfr_t result)
-{
-  int i;
-  mpfr_t t;
-
-  gfc_set_model (y);
-  mpfr_init (t);
-
-  i = mpfr_sgn (x);
-
-  if (i > 0)
-    {
-      mpfr_div (t, y, x, GFC_RND_MODE);
-      mpfr_atan (result, t, GFC_RND_MODE);
-    }
-  else if (i < 0)
-    {
-      mpfr_const_pi (result, GFC_RND_MODE);
-      mpfr_div (t, y, x, GFC_RND_MODE);
-      mpfr_abs (t, t, GFC_RND_MODE);
-      mpfr_atan (t, t, GFC_RND_MODE);
-      mpfr_sub (result, result, t, GFC_RND_MODE);
-      if (mpfr_sgn (y) < 0)
-       mpfr_neg (result, result, GFC_RND_MODE);
-    }
-  else
-    {
-      if (mpfr_sgn (y) == 0)
-       mpfr_set_ui (result, 0, GFC_RND_MODE);
-      else
-       {
-          mpfr_const_pi (result, GFC_RND_MODE);
-          mpfr_div_ui (result, result, 2, GFC_RND_MODE);
-         if (mpfr_sgn (y) < 0)
-           mpfr_neg (result, result, GFC_RND_MODE);
-       }
-    }
-
-  mpfr_clear (t);
-
-}
-
 
 /* Given an arithmetic error code, return a pointer to a string that
    explains the error.  */
@@ -138,25 +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_NAN:
-      p = "Arithmetic NaN";
+      p = _("Arithmetic NaN at %L");
       break;
     case ARITH_DIV0:
-      p = "Division by zero";
+      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";
+      p =
+       _("Integer outside symmetric range implied by Standard Fortran at %L");
       break;
     default:
       gfc_internal_error ("gfc_arith_error(): Bad error code");
@@ -173,36 +127,33 @@ gfc_arith_init_1 (void)
 {
   gfc_integer_info *int_info;
   gfc_real_info *real_info;
-  mpfr_t a, b, c;
-  mpz_t r;
+  mpfr_t a, b;
   int i;
 
   mpfr_set_default_prec (128);
   mpfr_init (a);
-  mpz_init (r);
 
-  /* Convert the minimum/maximum values for each kind into their
+  /* Convert the minimum and maximum values for each kind into their
      GNU MP representation.  */
   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");
 
       /* 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.  */
+        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.  */
 
       mpz_init (int_info->pedantic_min_int);
       mpz_neg (int_info->pedantic_min_int, int_info->huge);
@@ -210,16 +161,11 @@ gfc_arith_init_1 (void)
       mpz_init (int_info->min_int);
       mpz_sub_ui (int_info->min_int, int_info->pedantic_min_int, 1);
 
-      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 */
+      /* Range  */
       mpfr_set_z (a, int_info->huge, GFC_RND_MODE);
       mpfr_log10 (a, a, GFC_RND_MODE);
       mpfr_trunc (a, a);
-      gfc_mpfr_to_mpz (r, a);
-      int_info->range = mpz_get_si (r);
+      int_info->range = (int) mpfr_get_si (a, GFC_RND_MODE);
     }
 
   mpfr_clear (a);
@@ -230,83 +176,68 @@ gfc_arith_init_1 (void)
 
       mpfr_init (a);
       mpfr_init (b);
-      mpfr_init (c);
 
       /* huge(x) = (1 - b**(-p)) * b**(emax-1) * b  */
-      /* a = 1 - b**(-p) */
-      mpfr_set_ui (a, 1, GFC_RND_MODE);
-      mpfr_set_ui (b, real_info->radix, GFC_RND_MODE);
-      mpfr_pow_si (b, b, -real_info->digits, GFC_RND_MODE);
-      mpfr_sub (a, a, b, GFC_RND_MODE);
-
-      /* c = b**(emax-1) */
-      mpfr_set_ui (b, real_info->radix, GFC_RND_MODE);
-      mpfr_pow_ui (c, b, real_info->max_exponent - 1, GFC_RND_MODE);
-
-      /* a = a * c = (1 - b**(-p)) * b**(emax-1) */
-      mpfr_mul (a, a, c, GFC_RND_MODE);
-
-      /* a = (1 - b**(-p)) * b**(emax-1) * b */
-      mpfr_mul_ui (a, a, real_info->radix, GFC_RND_MODE);
-
+      /* 1 - b**(-p)  */
       mpfr_init (real_info->huge);
-      mpfr_set (real_info->huge, a, GFC_RND_MODE);
+      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);
 
-      /* tiny(x) = b**(emin-1) */
-      mpfr_set_ui (b, real_info->radix, GFC_RND_MODE);
-      mpfr_pow_si (b, b, real_info->min_exponent - 1, GFC_RND_MODE);
+      /* 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);
 
-      mpfr_init (real_info->tiny);
-      mpfr_set (real_info->tiny, b, GFC_RND_MODE);
+      /* (1 - b**(-p)) * b**(emax-1)  */
+      mpfr_mul (real_info->huge, real_info->huge, a, GFC_RND_MODE);
 
-      /* subnormal (x) = b**(emin - digit) */
-      mpfr_set_ui (b, real_info->radix, GFC_RND_MODE);
-      mpfr_pow_si (b, b, real_info->min_exponent - real_info->digits,
+      /* (1 - b**(-p)) * b**(emax-1) * b  */
+      mpfr_mul_ui (real_info->huge, real_info->huge, real_info->radix,
                   GFC_RND_MODE);
 
-      mpfr_init (real_info->subnormal);
-      mpfr_set (real_info->subnormal, b, 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);
 
-      /* epsilon(x) = b**(1-p) */
-      mpfr_set_ui (b, real_info->radix, GFC_RND_MODE);
-      mpfr_pow_si (b, b, 1 - real_info->digits, GFC_RND_MODE);
+      /* 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 (real_info->epsilon, b, GFC_RND_MODE);
+      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)) */
+      /* 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);
 
-      if (mpfr_cmp (a, b) > 0)
-       mpfr_set (a, b, GFC_RND_MODE);          /* a = min(a, b) */
-
+      /* a = min(a, b)  */
+      mpfr_min (a, a, b, GFC_RND_MODE);
       mpfr_trunc (a, a);
-      gfc_mpfr_to_mpz (r, a);
-      real_info->range = mpz_get_si (r);
+      real_info->range = (int) mpfr_get_si (a, GFC_RND_MODE);
 
-      /* precision(x) = int((p - 1) * log10(b)) + k */
+      /* 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);
-      gfc_mpfr_to_mpz (r, a);
-      real_info->precision = mpz_get_si (r);
+      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.  */
+      /* 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++;
 
-      mpfr_clear (a);
-      mpfr_clear (b);
-      mpfr_clear (c);
+      mpfr_clears (a, b, NULL);
     }
-
-  mpz_clear (r);
 }
 
 
@@ -321,16 +252,29 @@ gfc_arith_done_1 (void)
   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++)
-    {
-      mpfr_clear (rp->epsilon);
-      mpfr_clear (rp->huge);
-      mpfr_clear (rp->tiny);
-    }
+    mpfr_clears (rp->epsilon, rp->huge, rp->tiny, rp->subnormal, NULL);
+}
+
+
+/* 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)
+{
+  /* As wide characters are stored as 32-bit values, they're all
+     representable in UCS=4.  */
+  if (kind == 4)
+    return true;
+
+  if (kind == 1)
+    return c <= 255 ? true : false;
+
+  gcc_unreachable ();
 }
 
 
@@ -338,7 +282,7 @@ gfc_arith_done_1 (void)
    the range of the kind.  Returns ARITH_OK, ARITH_ASYMMETRIC or
    ARITH_OVERFLOW.  */
 
-static arith
+arith
 gfc_check_integer_range (mpz_t p, int kind)
 {
   arith result;
@@ -350,11 +294,15 @@ gfc_check_integer_range (mpz_t p, int kind)
   if (pedantic)
     {
       if (mpz_cmp (p, gfc_integer_kinds[i].pedantic_min_int) < 0)
-        result = ARITH_ASYMMETRIC;
+       result = ARITH_ASYMMETRIC;
     }
 
+
+  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;
@@ -378,47 +326,72 @@ gfc_check_real_range (mpfr_t p, int kind)
   mpfr_init (q);
   mpfr_abs (q, p, GFC_RND_MODE);
 
-  if (mpfr_sgn (q) == 0)
-    retval = ARITH_OK;
+  retval = ARITH_OK;
+
+  if (mpfr_inf_p (p))
+    {
+      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)
-    retval = ARITH_OVERFLOW;
+    {
+      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)
-    retval = ARITH_UNDERFLOW;
+    {
+      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)
     {
-      /* MPFR operates on a numbers with a given precision and enormous
-       exponential range.  To represent subnormal numbers the exponent is
-       allowed to become smaller than emin, but always retains the full
-       precision.  This function resets unused bits to 0 to alleviate
-       rounding problems.  Note, a future version of MPFR will have a
-       mpfr_subnormalize() function, which handles this truncation in a
-       more efficient and robust way.  */
-
-      int j, k;
-      char *bin, *s;
-      mp_exp_t e;
-
-      bin = mpfr_get_str (NULL, &e, gfc_real_kinds[i].radix, 0, q, GMP_RNDN);
-      k = gfc_real_kinds[i].digits - (gfc_real_kinds[i].min_exponent - e);
-      for (j = k; j < gfc_real_kinds[i].digits; j++)
-       bin[j] = '0';
-      /* Need space for '0.', bin, 'E', and e */
-      s = (char *) gfc_getmem (strlen(bin)+10);
-      sprintf (s, "0.%sE%d", bin, (int) e);
-      mpfr_set_str (q, s, gfc_real_kinds[i].radix, GMP_RNDN);
+      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);
-
-      gfc_free (s);
-      gfc_free (bin);
-
-      retval = ARITH_OK;
     }
-  else
-    retval = ARITH_OK;
 
   mpfr_clear (q);
 
@@ -426,17 +399,15 @@ gfc_check_real_range (mpfr_t p, int kind)
 }
 
 
-/* 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 ();
 
@@ -477,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;
 
@@ -490,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;
 
@@ -504,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;
 
@@ -518,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;
 
@@ -532,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;
 
@@ -550,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)
     {
@@ -563,20 +535,32 @@ 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);
+       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_UNDERFLOW)
-        mpfr_set_ui (e->value.complex.r, 0, GFC_RND_MODE);
-      if (rc == ARITH_OK || rc == ARITH_UNDERFLOW)
-        {
-          rc = gfc_check_real_range (e->value.complex.i, e->ts.kind);
-          if (rc == ARITH_UNDERFLOW)
-            mpfr_set_ui (e->value.complex.i, 0, GFC_RND_MODE);
-        }
+       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:
@@ -591,20 +575,20 @@ gfc_range_check (gfc_expr * e)
    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)
+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 ("%s at %L", gfc_arith_error (val), &x->where);
+       gfc_warning (gfc_arith_error (val), &x->where);
       val = ARITH_OK;
     }
 
   if (val == ARITH_ASYMMETRIC)
     {
-      gfc_warning ("%s at %L", gfc_arith_error (val), &x->where);
+      gfc_warning (gfc_arith_error (val), &x->where);
       val = ARITH_OK;
     }
 
@@ -619,10 +603,11 @@ check_result (arith rc, gfc_expr * x, gfc_expr * r, gfc_expr ** rp)
 
 /* 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;
@@ -630,7 +615,7 @@ gfc_arith_uplus (gfc_expr * op1, gfc_expr ** resultp)
 
 
 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;
@@ -663,7 +648,7 @@ gfc_arith_uminus (gfc_expr * op1, gfc_expr ** 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;
@@ -678,15 +663,15 @@ gfc_arith_plus (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
 
     case BT_REAL:
       mpfr_add (result->value.real, op1->value.real, op2->value.real,
-               GFC_RND_MODE);
+              GFC_RND_MODE);
       break;
 
     case BT_COMPLEX:
       mpfr_add (result->value.complex.r, op1->value.complex.r,
-              op2->value.complex.r, GFC_RND_MODE);
+               op2->value.complex.r, GFC_RND_MODE);
 
       mpfr_add (result->value.complex.i, op1->value.complex.i,
-              op2->value.complex.i, GFC_RND_MODE);
+               op2->value.complex.i, GFC_RND_MODE);
       break;
 
     default:
@@ -700,7 +685,7 @@ gfc_arith_plus (gfc_expr * op1, gfc_expr * op2, gfc_expr ** 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;
@@ -715,15 +700,15 @@ gfc_arith_minus (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
 
     case BT_REAL:
       mpfr_sub (result->value.real, op1->value.real, op2->value.real,
-                GFC_RND_MODE);
+               GFC_RND_MODE);
       break;
 
     case BT_COMPLEX:
       mpfr_sub (result->value.complex.r, op1->value.complex.r,
-              op2->value.complex.r, GFC_RND_MODE);
+               op2->value.complex.r, GFC_RND_MODE);
 
       mpfr_sub (result->value.complex.i, op1->value.complex.i,
-              op2->value.complex.i, GFC_RND_MODE);
+               op2->value.complex.i, GFC_RND_MODE);
       break;
 
     default:
@@ -737,7 +722,7 @@ gfc_arith_minus (gfc_expr * op1, gfc_expr * op2, gfc_expr ** 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;
   mpfr_t x, y;
@@ -753,13 +738,10 @@ gfc_arith_times (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
 
     case BT_REAL:
       mpfr_mul (result->value.real, op1->value.real, op2->value.real,
-               GFC_RND_MODE);
+              GFC_RND_MODE);
       break;
 
     case BT_COMPLEX:
-
-      /* FIXME:  possible numericals problem.  */
-
       gfc_set_model (op1->value.complex.r);
       mpfr_init (x);
       mpfr_init (y);
@@ -772,9 +754,7 @@ gfc_arith_times (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
       mpfr_mul (y, op1->value.complex.i, op2->value.complex.r, GFC_RND_MODE);
       mpfr_add (result->value.complex.i, x, y, GFC_RND_MODE);
 
-      mpfr_clear (x);
-      mpfr_clear (y);
-
+      mpfr_clears (x, y, NULL);
       break;
 
     default:
@@ -788,7 +768,7 @@ gfc_arith_times (gfc_expr * op1, gfc_expr * op2, gfc_expr ** 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;
   mpfr_t x, y, div;
@@ -812,21 +792,20 @@ gfc_arith_divide (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
       break;
 
     case BT_REAL:
-      /* FIXME: MPFR correctly generates NaN.  This may not be needed.  */
-      if (mpfr_sgn (op2->value.real) == 0)
+      if (mpfr_sgn (op2->value.real) == 0 && gfc_option.flag_range_check == 1)
        {
          rc = ARITH_DIV0;
          break;
        }
 
       mpfr_div (result->value.real, op1->value.real, op2->value.real,
-               GFC_RND_MODE);
+              GFC_RND_MODE);
       break;
 
     case BT_COMPLEX:
-      /* FIXME: MPFR correctly generates NaN.  This may not be needed.  */
       if (mpfr_sgn (op2->value.complex.r) == 0
-         && mpfr_sgn (op2->value.complex.i) == 0)
+         && mpfr_sgn (op2->value.complex.i) == 0
+         && gfc_option.flag_range_check == 1)
        {
          rc = ARITH_DIV0;
          break;
@@ -837,7 +816,6 @@ gfc_arith_divide (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
       mpfr_init (y);
       mpfr_init (div);
 
-      /* FIXME: possible numerical problems.  */
       mpfr_mul (x, op2->value.complex.r, op2->value.complex.r, GFC_RND_MODE);
       mpfr_mul (y, op2->value.complex.i, op2->value.complex.i, GFC_RND_MODE);
       mpfr_add (div, x, y, GFC_RND_MODE);
@@ -846,18 +824,15 @@ gfc_arith_divide (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
       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);
+               GFC_RND_MODE);
 
       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_clear (x);
-      mpfr_clear (y);
-      mpfr_clear (div);
+               GFC_RND_MODE);
 
+      mpfr_clears (x, y, div, NULL);
       break;
 
     default:
@@ -874,94 +849,105 @@ gfc_arith_divide (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
 /* Compute the reciprocal of a complex number (guaranteed nonzero).  */
 
 static void
-complex_reciprocal (gfc_expr * op)
+complex_reciprocal (gfc_expr *op)
 {
-  mpfr_t mod, a, re, im;
+  mpfr_t mod, tmp;
 
   gfc_set_model (op->value.complex.r);
   mpfr_init (mod);
-  mpfr_init (a);
-  mpfr_init (re);
-  mpfr_init (im);
+  mpfr_init (tmp);
 
-  /* FIXME:  another possible numerical problem.  */
   mpfr_mul (mod, op->value.complex.r, op->value.complex.r, GFC_RND_MODE);
-  mpfr_mul (a, op->value.complex.i, op->value.complex.i, GFC_RND_MODE);
-  mpfr_add (mod, mod, a, GFC_RND_MODE);
-
-  mpfr_div (re, op->value.complex.r, mod, 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);
 
-  mpfr_neg (im, op->value.complex.i, GFC_RND_MODE);
-  mpfr_div (im, im, mod, GFC_RND_MODE);
+  mpfr_div (op->value.complex.r, op->value.complex.r, mod, GFC_RND_MODE);
 
-  mpfr_set (op->value.complex.r, re, GFC_RND_MODE);
-  mpfr_set (op->value.complex.i, im, GFC_RND_MODE);
+  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);
 
-  mpfr_clear (re);
-  mpfr_clear (im);
-  mpfr_clear (mod);
-  mpfr_clear (a);
+  mpfr_clears (tmp, mod, NULL);
 }
 
 
-/* Raise a complex number to positive power.  */
+/* Raise a complex number to positive power (power > 0).
+   This function will modify the content of power.
+
+   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.  */
 
 static void
-complex_pow_ui (gfc_expr * base, int power, gfc_expr * result)
+complex_pow (gfc_expr *result, gfc_expr *base, mpz_t power)
 {
-  mpfr_t re, im, a;
+  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);
-  mpfr_init (a);
 
+  /* res = 1 */
   mpfr_set_ui (result->value.complex.r, 1, GFC_RND_MODE);
   mpfr_set_ui (result->value.complex.i, 0, GFC_RND_MODE);
 
-  for (; power > 0; power--)
+  /* 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))
     {
-      mpfr_mul (re, base->value.complex.r, result->value.complex.r,
-                GFC_RND_MODE);
-      mpfr_mul (a, base->value.complex.i, result->value.complex.i,
-                GFC_RND_MODE);
-      mpfr_sub (re, re, a, GFC_RND_MODE);
-
-      mpfr_mul (im, base->value.complex.r, result->value.complex.i,
-                GFC_RND_MODE);
-      mpfr_mul (a, base->value.complex.i, result->value.complex.r,
-                GFC_RND_MODE);
-      mpfr_add (im, im, a, GFC_RND_MODE);
-
-      mpfr_set (result->value.complex.r, re, GFC_RND_MODE);
-      mpfr_set (result->value.complex.i, im, GFC_RND_MODE);
+      /* 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);
+
+      /* power /= 2; */
+      mpz_fdiv_q_ui (power, power, 2);
     }
 
-  mpfr_clear (re);
-  mpfr_clear (im);
-  mpfr_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;
-  mpfr_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)
+  if (power_sign == 0)
     {
       /* Handle something to the zeroth power.  Since we're dealing
         with integral exponents, there is no ambiguity in the
@@ -987,45 +973,86 @@ gfc_arith_power (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
     }
   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:
-         mpfr_pow_ui (result->value.real, op1->value.real, apower,
-                       GFC_RND_MODE);
-
-         if (power < 0)
-           {
-              gfc_set_model (op1->value.real);
-             mpfr_init (unity_f);
-             mpfr_set_ui (unity_f, 1, GFC_RND_MODE);
-             mpfr_div (result->value.real, unity_f, result->value.real,
-                        GFC_RND_MODE);
-             mpfr_clear (unity_f);
-           }
+         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);
-         break;
+         {
+           mpz_t apower;
+
+           /* 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;
@@ -1042,24 +1069,26 @@ gfc_arith_power (gfc_expr * op1, gfc_expr * op2, gfc_expr ** 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';
 
@@ -1068,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;
 
@@ -1084,11 +1144,11 @@ gfc_compare_expr (gfc_expr * op1, gfc_expr * op2)
       break;
 
     case BT_REAL:
-      rc = mpfr_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:
@@ -1105,40 +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 (mpfr_cmp (op1->value.complex.r, op2->value.complex.r) == 0
-         && mpfr_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)
@@ -1148,7 +1235,6 @@ gfc_compare_string (gfc_expr * a, gfc_expr * b, const int *xcoll_table)
     }
 
   /* Strings are equal */
-
   return 0;
 }
 
@@ -1156,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,
                                &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;
@@ -1171,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,
                                &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;
@@ -1186,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,
                                &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;
@@ -1200,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,
                                &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;
@@ -1214,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,
                                &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;
@@ -1228,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,
                                &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;
@@ -1242,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;
@@ -1257,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;
 
@@ -1286,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;
@@ -1298,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;
 
@@ -1327,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;
@@ -1339,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;
 
@@ -1366,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;
@@ -1380,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)
@@ -1394,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;
 
@@ -1427,8 +1525,7 @@ 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);
@@ -1461,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;
@@ -1470,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;
-
       unary = 1;
       break;
 
-      /* Logical binary operators */
+    /* Logical binary operators  */
     case INTRINSIC_OR:
     case INTRINSIC_AND:
     case INTRINSIC_NEQV:
@@ -1492,24 +1589,33 @@ eval_intrinsic (gfc_intrinsic_op operator,
 
       temp.ts.type = BT_LOGICAL;
       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;
@@ -1517,42 +1623,52 @@ eval_intrinsic (gfc_intrinsic_op operator,
          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;
+
+         /* 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.value.op.operator = operator;
+      temp.value.op.op = op;
 
       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;
@@ -1561,13 +1677,14 @@ eval_intrinsic (gfc_intrinsic_op operator,
       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;
 
@@ -1579,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)
@@ -1601,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;
     }
 
@@ -1611,12 +1726,12 @@ 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->value.op.operator = operator;
+  result->value.op.op = op;
 
   result->value.op.op1 = op1;
   result->value.op.op2 = op2;
@@ -1628,20 +1743,27 @@ 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");
 
-  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;
       break;
@@ -1657,7 +1779,7 @@ 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;
@@ -1671,7 +1793,7 @@ 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))
     {
@@ -1690,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;
@@ -1700,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
@@ -1870,7 +2019,7 @@ 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;
 
@@ -1885,7 +2034,7 @@ 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;
 
@@ -1903,19 +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)
 {
-  gfc_error ("%s converting %s to %s at %L", gfc_arith_error (rc),
-            gfc_typename (from), gfc_typename (to), 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");
+    }
 
-  /* TODO: Do something about the error, ie, throw exception, return
+  /* 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;
@@ -1924,19 +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)
     {
       if (rc == ARITH_ASYMMETRIC)
-        {
-          gfc_warning ("%s at %L", gfc_arith_error (rc), &src->where);
-        }
+       {
+         gfc_warning (gfc_arith_error (rc), &src->where);
+       }
       else
-        {
-          arith_error (rc, &src->ts, &result->ts, &src->where);
-          gfc_free_expr (result);
-          return NULL;
-        }
+       {
+         arith_error (rc, &src->ts, &result->ts, &src->where);
+         gfc_free_expr (result);
+         return NULL;
+       }
     }
 
   return result;
@@ -1946,7 +2130,7 @@ 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;
@@ -1969,7 +2153,7 @@ 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;
@@ -1993,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);
 
-  gfc_mpfr_to_mpz (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);
@@ -2017,7 +2200,7 @@ 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;
@@ -2031,7 +2214,7 @@ gfc_real2real (gfc_expr * src, int kind)
   if (rc == ARITH_UNDERFLOW)
     {
       if (gfc_option.warn_underflow)
-        gfc_warning ("%s at %L", gfc_arith_error (rc), &src->where);
+       gfc_warning (gfc_arith_error (rc), &src->where);
       mpfr_set_ui (result->value.real, 0, GFC_RND_MODE);
     }
   else if (rc != ARITH_OK)
@@ -2048,7 +2231,7 @@ 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;
@@ -2063,7 +2246,7 @@ gfc_real2complex (gfc_expr * src, int kind)
   if (rc == ARITH_UNDERFLOW)
     {
       if (gfc_option.warn_underflow)
-        gfc_warning ("%s at %L", gfc_arith_error (rc), &src->where);
+       gfc_warning (gfc_arith_error (rc), &src->where);
       mpfr_set_ui (result->value.complex.r, 0, GFC_RND_MODE);
     }
   else if (rc != ARITH_OK)
@@ -2080,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);
 
-  gfc_mpfr_to_mpz (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);
@@ -2104,7 +2286,7 @@ 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;
@@ -2118,7 +2300,7 @@ gfc_complex2real (gfc_expr * src, int kind)
   if (rc == ARITH_UNDERFLOW)
     {
       if (gfc_option.warn_underflow)
-        gfc_warning ("%s at %L", gfc_arith_error (rc), &src->where);
+       gfc_warning (gfc_arith_error (rc), &src->where);
       mpfr_set_ui (result->value.real, 0, GFC_RND_MODE);
     }
   if (rc != ARITH_OK)
@@ -2135,7 +2317,7 @@ 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;
@@ -2150,7 +2332,7 @@ gfc_complex2complex (gfc_expr * src, int kind)
   if (rc == ARITH_UNDERFLOW)
     {
       if (gfc_option.warn_underflow)
-        gfc_warning ("%s at %L", gfc_arith_error (rc), &src->where);
+       gfc_warning (gfc_arith_error (rc), &src->where);
       mpfr_set_ui (result->value.complex.r, 0, GFC_RND_MODE);
     }
   else if (rc != ARITH_OK)
@@ -2165,7 +2347,7 @@ gfc_complex2complex (gfc_expr * src, int kind)
   if (rc == ARITH_UNDERFLOW)
     {
       if (gfc_option.warn_underflow)
-        gfc_warning ("%s at %L", gfc_arith_error (rc), &src->where);
+       gfc_warning (gfc_arith_error (rc), &src->where);
       mpfr_set_ui (result->value.complex.i, 0, GFC_RND_MODE);
     }
   else if (rc != ARITH_OK)
@@ -2182,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;
 
@@ -2192,25 +2374,218 @@ 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;
+}