OSDN Git Service

2010-02-20 Paul Thomas <pault@gcc.gnu.org>
[pf3gnuchains/gcc-fork.git] / gcc / fortran / simplify.c
index 429c515..8768cb6 100644 (file)
@@ -1,5 +1,5 @@
 /* Simplify intrinsic functions at compile-time.
-   Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008
+   Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009
    Free Software Foundation, Inc.
    Contributed by Andy Vaught & Katherine Holcomb
 
@@ -27,6 +27,10 @@ along with GCC; see the file COPYING3.  If not see
 #include "intrinsic.h"
 #include "target-memory.h"
 
+/* Savely advance an array constructor by 'n' elements.
+   Mainly used by simplifiers of transformational intrinsics.  */
+#define ADVANCE(ctor, n) do { int i; for (i = 0; i < n && ctor; ++i) ctor = ctor->next; } while (0)
+
 gfc_expr gfc_bad_expr;
 
 
@@ -210,6 +214,397 @@ convert_mpz_to_signed (mpz_t x, int bitsize)
     }
 }
 
+/* Test that the expression is an constant array.  */
+
+static bool
+is_constant_array_expr (gfc_expr *e)
+{
+  gfc_constructor *c;
+
+  if (e == NULL)
+    return true;
+
+  if (e->expr_type != EXPR_ARRAY || !gfc_is_constant_expr (e))
+    return false;
+
+  for (c = e->value.constructor; c; c = c->next)
+    if (c->expr->expr_type != EXPR_CONSTANT)
+      return false;
+
+  return true;
+}
+
+
+/* Initialize a transformational result expression with a given value.  */
+
+static void
+init_result_expr (gfc_expr *e, int init, gfc_expr *array)
+{
+  if (e && e->expr_type == EXPR_ARRAY)
+    {
+      gfc_constructor *ctor = e->value.constructor;
+      while (ctor)
+       {
+         init_result_expr (ctor->expr, init, array);
+         ctor = ctor->next;
+       }
+    }
+  else if (e && e->expr_type == EXPR_CONSTANT)
+    {
+      int i = gfc_validate_kind (e->ts.type, e->ts.kind, false);
+      int length;
+      gfc_char_t *string;
+
+      switch (e->ts.type)
+       {
+         case BT_LOGICAL:
+           e->value.logical = (init ? 1 : 0);
+           break;
+
+         case BT_INTEGER:
+           if (init == INT_MIN)
+             mpz_set (e->value.integer, gfc_integer_kinds[i].min_int);
+           else if (init == INT_MAX)
+             mpz_set (e->value.integer, gfc_integer_kinds[i].huge);
+           else
+             mpz_set_si (e->value.integer, init);
+           break;
+
+         case BT_REAL:
+           if (init == INT_MIN)
+             {
+               mpfr_set (e->value.real, gfc_real_kinds[i].huge, GFC_RND_MODE);
+               mpfr_neg (e->value.real, e->value.real, GFC_RND_MODE);
+             }
+           else if (init == INT_MAX)
+             mpfr_set (e->value.real, gfc_real_kinds[i].huge, GFC_RND_MODE);
+           else
+             mpfr_set_si (e->value.real, init, GFC_RND_MODE);
+           break;
+
+         case BT_COMPLEX:
+           mpc_set_si (e->value.complex, init, GFC_MPC_RND_MODE);
+           break;
+
+         case BT_CHARACTER:
+           if (init == INT_MIN)
+             {
+               gfc_expr *len = gfc_simplify_len (array, NULL);
+               gfc_extract_int (len, &length);
+               string = gfc_get_wide_string (length + 1);
+               gfc_wide_memset (string, 0, length);
+             }
+           else if (init == INT_MAX)
+             {
+               gfc_expr *len = gfc_simplify_len (array, NULL);
+               gfc_extract_int (len, &length);
+               string = gfc_get_wide_string (length + 1);
+               gfc_wide_memset (string, 255, length);
+             }
+           else
+             {
+               length = 0;
+               string = gfc_get_wide_string (1);
+             }
+
+           string[length] = '\0';
+           e->value.character.length = length;
+           e->value.character.string = string;
+           break;
+
+         default:
+           gcc_unreachable();
+       }
+    }
+  else
+    gcc_unreachable();
+}
+
+
+/* Helper function for gfc_simplify_dot_product() and gfc_simplify_matmul.  */
+
+static gfc_expr *
+compute_dot_product (gfc_constructor *ctor_a, int stride_a,
+                    gfc_constructor *ctor_b, int stride_b)
+{
+  gfc_expr *result;
+  gfc_expr *a = ctor_a->expr, *b = ctor_b->expr;
+
+  gcc_assert (gfc_compare_types (&a->ts, &b->ts));
+
+  result = gfc_constant_result (a->ts.type, a->ts.kind, &a->where);
+  init_result_expr (result, 0, NULL);
+
+  while (ctor_a && ctor_b)
+    {
+      /* Copying of expressions is required as operands are free'd
+        by the gfc_arith routines.  */
+      switch (result->ts.type)
+       {
+         case BT_LOGICAL:
+           result = gfc_or (result,
+                            gfc_and (gfc_copy_expr (ctor_a->expr),
+                                     gfc_copy_expr (ctor_b->expr)));
+           break;
+
+         case BT_INTEGER:
+         case BT_REAL:
+         case BT_COMPLEX:
+           result = gfc_add (result,
+                             gfc_multiply (gfc_copy_expr (ctor_a->expr),
+                                           gfc_copy_expr (ctor_b->expr)));
+           break;
+
+         default:
+           gcc_unreachable();
+       }
+
+      ADVANCE (ctor_a, stride_a);
+      ADVANCE (ctor_b, stride_b);
+    }
+
+  return result;
+}
+
+
+/* Build a result expression for transformational intrinsics, 
+   depending on DIM. */
+
+static gfc_expr *
+transformational_result (gfc_expr *array, gfc_expr *dim, bt type,
+                        int kind, locus* where)
+{
+  gfc_expr *result;
+  int i, nelem;
+
+  if (!dim || array->rank == 1)
+    return gfc_constant_result (type, kind, where);
+
+  result = gfc_start_constructor (type, kind, where);
+  result->shape = gfc_copy_shape_excluding (array->shape, array->rank, dim);
+  result->rank = array->rank - 1;
+
+  /* gfc_array_size() would count the number of elements in the constructor,
+     we have not built those yet.  */
+  nelem = 1;
+  for  (i = 0; i < result->rank; ++i)
+    nelem *= mpz_get_ui (result->shape[i]);
+
+  for (i = 0; i < nelem; ++i)
+    {
+      gfc_expr *e = gfc_constant_result (type, kind, where);
+      gfc_append_constructor (result, e);
+    }
+
+  return result;
+}
+
+
+typedef gfc_expr* (*transformational_op)(gfc_expr*, gfc_expr*);
+
+/* Wrapper function, implements 'op1 += 1'. Only called if MASK
+   of COUNT intrinsic is .TRUE..
+
+   Interface and implimentation mimics arith functions as
+   gfc_add, gfc_multiply, etc.  */
+
+static gfc_expr* gfc_count (gfc_expr *op1, gfc_expr *op2)
+{
+  gfc_expr *result;
+
+  gcc_assert (op1->ts.type == BT_INTEGER);
+  gcc_assert (op2->ts.type == BT_LOGICAL);
+  gcc_assert (op2->value.logical);
+
+  result = gfc_copy_expr (op1);
+  mpz_add_ui (result->value.integer, result->value.integer, 1);
+
+  gfc_free_expr (op1);
+  gfc_free_expr (op2);
+  return result;
+}
+
+
+/* Transforms an ARRAY with operation OP, according to MASK, to a
+   scalar RESULT. E.g. called if
+
+     REAL, PARAMETER :: array(n, m) = ...
+     REAL, PARAMETER :: s = SUM(array)
+
+  where OP == gfc_add().  */
+
+static gfc_expr *
+simplify_transformation_to_scalar (gfc_expr *result, gfc_expr *array, gfc_expr *mask,
+                                  transformational_op op)
+{
+  gfc_expr *a, *m;
+  gfc_constructor *array_ctor, *mask_ctor;
+
+  /* Shortcut for constant .FALSE. MASK.  */
+  if (mask
+      && mask->expr_type == EXPR_CONSTANT
+      && !mask->value.logical)
+    return result;
+
+  array_ctor = array->value.constructor;
+  mask_ctor = NULL;
+  if (mask && mask->expr_type == EXPR_ARRAY)
+    mask_ctor = mask->value.constructor;
+
+  while (array_ctor)
+    {
+      a = array_ctor->expr;
+      array_ctor = array_ctor->next;
+
+      /* A constant MASK equals .TRUE. here and can be ignored.  */
+      if (mask_ctor)
+       {
+         m = mask_ctor->expr;
+         mask_ctor = mask_ctor->next;
+         if (!m->value.logical)
+           continue;
+       }
+
+      result = op (result, gfc_copy_expr (a));
+    }
+
+  return result;
+}
+
+/* Transforms an ARRAY with operation OP, according to MASK, to an
+   array RESULT. E.g. called if
+
+     REAL, PARAMETER :: array(n, m) = ...
+     REAL, PARAMETER :: s(n) = PROD(array, DIM=1)
+
+  where OP == gfc_multiply().  */
+
+static gfc_expr *
+simplify_transformation_to_array (gfc_expr *result, gfc_expr *array, gfc_expr *dim,
+                                 gfc_expr *mask, transformational_op op)
+{
+  mpz_t size;
+  int done, i, n, arraysize, resultsize, dim_index, dim_extent, dim_stride;
+  gfc_expr **arrayvec, **resultvec, **base, **src, **dest;
+  gfc_constructor *array_ctor, *mask_ctor, *result_ctor;
+
+  int count[GFC_MAX_DIMENSIONS], extent[GFC_MAX_DIMENSIONS],
+      sstride[GFC_MAX_DIMENSIONS], dstride[GFC_MAX_DIMENSIONS],
+      tmpstride[GFC_MAX_DIMENSIONS];
+
+  /* Shortcut for constant .FALSE. MASK.  */
+  if (mask
+      && mask->expr_type == EXPR_CONSTANT
+      && !mask->value.logical)
+    return result;
+
+  /* Build an indexed table for array element expressions to minimize
+     linked-list traversal. Masked elements are set to NULL.  */
+  gfc_array_size (array, &size);
+  arraysize = mpz_get_ui (size);
+
+  arrayvec = (gfc_expr**) gfc_getmem (sizeof (gfc_expr*) * arraysize);
+
+  array_ctor = array->value.constructor;
+  mask_ctor = NULL;
+  if (mask && mask->expr_type == EXPR_ARRAY)
+    mask_ctor = mask->value.constructor;
+
+  for (i = 0; i < arraysize; ++i)
+    {
+      arrayvec[i] = array_ctor->expr;
+      array_ctor = array_ctor->next;
+
+      if (mask_ctor)
+       {
+         if (!mask_ctor->expr->value.logical)
+           arrayvec[i] = NULL;
+
+         mask_ctor = mask_ctor->next;
+       }
+    }
+
+  /* Same for the result expression.  */
+  gfc_array_size (result, &size);
+  resultsize = mpz_get_ui (size);
+  mpz_clear (size);
+
+  resultvec = (gfc_expr**) gfc_getmem (sizeof (gfc_expr*) * resultsize);
+  result_ctor = result->value.constructor;
+  for (i = 0; i < resultsize; ++i)
+    {
+      resultvec[i] = result_ctor->expr;
+      result_ctor = result_ctor->next;
+    }
+
+  gfc_extract_int (dim, &dim_index);
+  dim_index -= 1;               /* zero-base index */
+  dim_extent = 0;
+  dim_stride = 0;
+
+  for (i = 0, n = 0; i < array->rank; ++i)
+    {
+      count[i] = 0;
+      tmpstride[i] = (i == 0) ? 1 : tmpstride[i-1] * mpz_get_si (array->shape[i-1]);
+      if (i == dim_index)
+       {
+         dim_extent = mpz_get_si (array->shape[i]);
+         dim_stride = tmpstride[i];
+         continue;
+       }
+
+      extent[n] = mpz_get_si (array->shape[i]);
+      sstride[n] = tmpstride[i];
+      dstride[n] = (n == 0) ? 1 : dstride[n-1] * extent[n-1];
+      n += 1;
+    }
+
+  done = false;
+  base = arrayvec;
+  dest = resultvec;
+  while (!done)
+    {
+      for (src = base, n = 0; n < dim_extent; src += dim_stride, ++n)
+       if (*src)
+         *dest = op (*dest, gfc_copy_expr (*src));
+
+      count[0]++;
+      base += sstride[0];
+      dest += dstride[0];
+
+      n = 0;
+      while (!done && count[n] == extent[n])
+       {
+         count[n] = 0;
+         base -= sstride[n] * extent[n];
+         dest -= dstride[n] * extent[n];
+
+         n++;
+         if (n < result->rank)
+           {
+             count [n]++;
+             base += sstride[n];
+             dest += dstride[n];
+           }
+         else
+           done = true;
+       }
+    }
+
+  /* Place updated expression in result constructor.  */
+  result_ctor = result->value.constructor;
+  for (i = 0; i < resultsize; ++i)
+    {
+      result_ctor->expr = resultvec[i];
+      result_ctor = result_ctor->next;
+    }
+
+  gfc_free (arrayvec);
+  gfc_free (resultvec);
+  return result;
+}
+
+
 
 /********************** Simplification functions *****************************/
 
@@ -244,8 +639,7 @@ gfc_simplify_abs (gfc_expr *e)
 
       gfc_set_model_kind (e->ts.kind);
 
-      mpfr_hypot (result->value.real, e->value.complex.r, 
-                 e->value.complex.i, GFC_RND_MODE);
+      mpc_abs (result->value.real, e->value.complex, GFC_RND_MODE);
       result = range_check (result, "CABS");
       break;
 
@@ -331,17 +725,27 @@ gfc_simplify_acos (gfc_expr *x)
   if (x->expr_type != EXPR_CONSTANT)
     return NULL;
 
-  if (mpfr_cmp_si (x->value.real, 1) > 0
-      || mpfr_cmp_si (x->value.real, -1) < 0)
+  switch (x->ts.type)
     {
-      gfc_error ("Argument of ACOS at %L must be between -1 and 1",
-                &x->where);
-      return &gfc_bad_expr;
+      case BT_REAL:
+       if (mpfr_cmp_si (x->value.real, 1) > 0
+           || mpfr_cmp_si (x->value.real, -1) < 0)
+         {
+           gfc_error ("Argument of ACOS at %L must be between -1 and 1",
+                      &x->where);
+           return &gfc_bad_expr;
+         }
+       result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
+       mpfr_acos (result->value.real, x->value.real, GFC_RND_MODE);
+       break;
+      case BT_COMPLEX:
+       result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
+       mpc_acos (result->value.complex, x->value.complex, GFC_MPC_RND_MODE);
+       break;
+      default:
+       gfc_internal_error ("in gfc_simplify_acos(): Bad type");
     }
 
-  result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
-
-  mpfr_acos (result->value.real, x->value.real, GFC_RND_MODE);
 
   return range_check (result, "ACOS");
 }
@@ -354,16 +758,26 @@ gfc_simplify_acosh (gfc_expr *x)
   if (x->expr_type != EXPR_CONSTANT)
     return NULL;
 
-  if (mpfr_cmp_si (x->value.real, 1) < 0)
+  switch (x->ts.type)
     {
-      gfc_error ("Argument of ACOSH at %L must not be less than 1",
-                &x->where);
-      return &gfc_bad_expr;
-    }
-
-  result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
+      case BT_REAL:
+       if (mpfr_cmp_si (x->value.real, 1) < 0)
+         {
+           gfc_error ("Argument of ACOSH at %L must not be less than 1",
+                      &x->where);
+           return &gfc_bad_expr;
+         }
 
-  mpfr_acosh (result->value.real, x->value.real, GFC_RND_MODE);
+       result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
+       mpfr_acosh (result->value.real, x->value.real, GFC_RND_MODE);
+       break;
+      case BT_COMPLEX:
+       result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
+       mpc_acosh (result->value.complex, x->value.complex, GFC_MPC_RND_MODE);
+       break;
+      default:
+       gfc_internal_error ("in gfc_simplify_acosh(): Bad type");
+    }
 
   return range_check (result, "ACOSH");
 }
@@ -451,7 +865,7 @@ gfc_simplify_aimag (gfc_expr *e)
     return NULL;
 
   result = gfc_constant_result (BT_REAL, e->ts.kind, &e->where);
-  mpfr_set (result->value.real, e->value.complex.i, GFC_RND_MODE);
+  mpfr_set (result->value.real, mpc_imagref (e->value.complex), GFC_RND_MODE);
 
   return range_check (result, "AIMAG");
 }
@@ -482,6 +896,25 @@ gfc_simplify_aint (gfc_expr *e, gfc_expr *k)
 
 
 gfc_expr *
+gfc_simplify_all (gfc_expr *mask, gfc_expr *dim)
+{
+  gfc_expr *result;
+
+  if (!is_constant_array_expr (mask)
+      || !gfc_is_constant_expr (dim))
+    return NULL;
+
+  result = transformational_result (mask, dim, mask->ts.type,
+                                   mask->ts.kind, &mask->where);
+  init_result_expr (result, true, NULL);
+
+  return !dim || mask->rank == 1 ?
+    simplify_transformation_to_scalar (result, mask, NULL, gfc_and) :
+    simplify_transformation_to_array (result, mask, dim, NULL, gfc_and);
+}
+
+
+gfc_expr *
 gfc_simplify_dint (gfc_expr *e)
 {
   gfc_expr *rtrunc, *result;
@@ -547,6 +980,25 @@ gfc_simplify_and (gfc_expr *x, gfc_expr *y)
 
 
 gfc_expr *
+gfc_simplify_any (gfc_expr *mask, gfc_expr *dim)
+{
+  gfc_expr *result;
+
+  if (!is_constant_array_expr (mask)
+      || !gfc_is_constant_expr (dim))
+    return NULL;
+
+  result = transformational_result (mask, dim, mask->ts.type,
+                                   mask->ts.kind, &mask->where);
+  init_result_expr (result, false, NULL);
+
+  return !dim || mask->rank == 1 ?
+    simplify_transformation_to_scalar (result, mask, NULL, gfc_or) :
+    simplify_transformation_to_array (result, mask, dim, NULL, gfc_or);
+}
+
+
+gfc_expr *
 gfc_simplify_dnint (gfc_expr *e)
 {
   gfc_expr *result;
@@ -570,18 +1022,27 @@ gfc_simplify_asin (gfc_expr *x)
   if (x->expr_type != EXPR_CONSTANT)
     return NULL;
 
-  if (mpfr_cmp_si (x->value.real, 1) > 0
-      || mpfr_cmp_si (x->value.real, -1) < 0)
+  switch (x->ts.type)
     {
-      gfc_error ("Argument of ASIN at %L must be between -1 and 1",
-                &x->where);
-      return &gfc_bad_expr;
+      case BT_REAL:
+       if (mpfr_cmp_si (x->value.real, 1) > 0
+           || mpfr_cmp_si (x->value.real, -1) < 0)
+         {
+           gfc_error ("Argument of ASIN at %L must be between -1 and 1",
+                      &x->where);
+           return &gfc_bad_expr;
+         }
+       result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
+       mpfr_asin (result->value.real, x->value.real, GFC_RND_MODE);
+       break;
+      case BT_COMPLEX:
+       result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
+       mpc_asin (result->value.complex, x->value.complex, GFC_MPC_RND_MODE);
+       break;
+      default:
+       gfc_internal_error ("in gfc_simplify_asin(): Bad type");
     }
 
-  result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
-
-  mpfr_asin (result->value.real, x->value.real, GFC_RND_MODE);
-
   return range_check (result, "ASIN");
 }
 
@@ -594,9 +1055,19 @@ gfc_simplify_asinh (gfc_expr *x)
   if (x->expr_type != EXPR_CONSTANT)
     return NULL;
 
-  result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
-
-  mpfr_asinh (result->value.real, x->value.real, GFC_RND_MODE);
+  switch (x->ts.type)
+    {
+      case BT_REAL:
+       result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
+       mpfr_asinh (result->value.real, x->value.real, GFC_RND_MODE);
+       break;
+      case BT_COMPLEX:
+       result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
+       mpc_asinh (result->value.complex, x->value.complex, GFC_MPC_RND_MODE);
+       break;
+      default:
+       gfc_internal_error ("in gfc_simplify_asinh(): Bad type");
+    }
 
   return range_check (result, "ASINH");
 }
@@ -610,9 +1081,19 @@ gfc_simplify_atan (gfc_expr *x)
   if (x->expr_type != EXPR_CONSTANT)
     return NULL;
     
-  result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
-
-  mpfr_atan (result->value.real, x->value.real, GFC_RND_MODE);
+  switch (x->ts.type)
+    {
+      case BT_REAL:
+       result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
+       mpfr_atan (result->value.real, x->value.real, GFC_RND_MODE);
+       break;
+      case BT_COMPLEX:
+       result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
+       mpc_atan (result->value.complex, x->value.complex, GFC_MPC_RND_MODE);
+       break;
+      default:
+       gfc_internal_error ("in gfc_simplify_atan(): Bad type");
+    }
 
   return range_check (result, "ATAN");
 }
@@ -626,17 +1107,27 @@ gfc_simplify_atanh (gfc_expr *x)
   if (x->expr_type != EXPR_CONSTANT)
     return NULL;
 
-  if (mpfr_cmp_si (x->value.real, 1) >= 0
-      || mpfr_cmp_si (x->value.real, -1) <= 0)
+  switch (x->ts.type)
     {
-      gfc_error ("Argument of ATANH at %L must be inside the range -1 to 1",
-                &x->where);
-      return &gfc_bad_expr;
-    }
-
-  result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
+      case BT_REAL:
+       if (mpfr_cmp_si (x->value.real, 1) >= 0
+           || mpfr_cmp_si (x->value.real, -1) <= 0)
+         {
+           gfc_error ("Argument of ATANH at %L must be inside the range -1 "
+                      "to 1", &x->where);
+           return &gfc_bad_expr;
+         }
 
-  mpfr_atanh (result->value.real, x->value.real, GFC_RND_MODE);
+       result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
+       mpfr_atanh (result->value.real, x->value.real, GFC_RND_MODE);
+       break;
+      case BT_COMPLEX:
+       result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
+       mpc_atanh (result->value.complex, x->value.complex, GFC_MPC_RND_MODE);
+       break;
+      default:
+       gfc_internal_error ("in gfc_simplify_atanh(): Bad type");
+    }
 
   return range_check (result, "ATANH");
 }
@@ -668,7 +1159,6 @@ gfc_simplify_atan2 (gfc_expr *y, gfc_expr *x)
 gfc_expr *
 gfc_simplify_bessel_j0 (gfc_expr *x ATTRIBUTE_UNUSED)
 {
-#if MPFR_VERSION >= MPFR_VERSION_NUM(2,3,0)
   gfc_expr *result;
 
   if (x->expr_type != EXPR_CONSTANT)
@@ -678,16 +1168,12 @@ gfc_simplify_bessel_j0 (gfc_expr *x ATTRIBUTE_UNUSED)
   mpfr_j0 (result->value.real, x->value.real, GFC_RND_MODE);
 
   return range_check (result, "BESSEL_J0");
-#else
-  return NULL;
-#endif
 }
 
 
 gfc_expr *
 gfc_simplify_bessel_j1 (gfc_expr *x ATTRIBUTE_UNUSED)
 {
-#if MPFR_VERSION >= MPFR_VERSION_NUM(2,3,0)
   gfc_expr *result;
 
   if (x->expr_type != EXPR_CONSTANT)
@@ -697,9 +1183,6 @@ gfc_simplify_bessel_j1 (gfc_expr *x ATTRIBUTE_UNUSED)
   mpfr_j1 (result->value.real, x->value.real, GFC_RND_MODE);
 
   return range_check (result, "BESSEL_J1");
-#else
-  return NULL;
-#endif
 }
 
 
@@ -707,7 +1190,6 @@ gfc_expr *
 gfc_simplify_bessel_jn (gfc_expr *order ATTRIBUTE_UNUSED,
                        gfc_expr *x ATTRIBUTE_UNUSED)
 {
-#if MPFR_VERSION >= MPFR_VERSION_NUM(2,3,0)
   gfc_expr *result;
   long n;
 
@@ -719,16 +1201,12 @@ gfc_simplify_bessel_jn (gfc_expr *order ATTRIBUTE_UNUSED,
   mpfr_jn (result->value.real, n, x->value.real, GFC_RND_MODE);
 
   return range_check (result, "BESSEL_JN");
-#else
-  return NULL;
-#endif
 }
 
 
 gfc_expr *
 gfc_simplify_bessel_y0 (gfc_expr *x ATTRIBUTE_UNUSED)
 {
-#if MPFR_VERSION >= MPFR_VERSION_NUM(2,3,0)
   gfc_expr *result;
 
   if (x->expr_type != EXPR_CONSTANT)
@@ -738,16 +1216,12 @@ gfc_simplify_bessel_y0 (gfc_expr *x ATTRIBUTE_UNUSED)
   mpfr_y0 (result->value.real, x->value.real, GFC_RND_MODE);
 
   return range_check (result, "BESSEL_Y0");
-#else
-  return NULL;
-#endif
 }
 
 
 gfc_expr *
 gfc_simplify_bessel_y1 (gfc_expr *x ATTRIBUTE_UNUSED)
 {
-#if MPFR_VERSION >= MPFR_VERSION_NUM(2,3,0)
   gfc_expr *result;
 
   if (x->expr_type != EXPR_CONSTANT)
@@ -757,9 +1231,6 @@ gfc_simplify_bessel_y1 (gfc_expr *x ATTRIBUTE_UNUSED)
   mpfr_y1 (result->value.real, x->value.real, GFC_RND_MODE);
 
   return range_check (result, "BESSEL_Y1");
-#else
-  return NULL;
-#endif
 }
 
 
@@ -767,7 +1238,6 @@ gfc_expr *
 gfc_simplify_bessel_yn (gfc_expr *order ATTRIBUTE_UNUSED,
                        gfc_expr *x ATTRIBUTE_UNUSED)
 {
-#if MPFR_VERSION >= MPFR_VERSION_NUM(2,3,0)
   gfc_expr *result;
   long n;
 
@@ -779,9 +1249,6 @@ gfc_simplify_bessel_yn (gfc_expr *order ATTRIBUTE_UNUSED,
   mpfr_yn (result->value.real, n, x->value.real, GFC_RND_MODE);
 
   return range_check (result, "BESSEL_YN");
-#else
-  return NULL;
-#endif
 }
 
 
@@ -832,7 +1299,7 @@ gfc_simplify_ceiling (gfc_expr *e, gfc_expr *k)
   ceil = gfc_copy_expr (e);
 
   mpfr_ceil (ceil->value.real, e->value.real);
-  gfc_mpfr_to_mpz (result->value.integer, ceil->value.real);
+  gfc_mpfr_to_mpz (result->value.integer, ceil->value.real, &e->where);
 
   gfc_free_expr (ceil);
 
@@ -856,22 +1323,19 @@ simplify_cmplx (const char *name, gfc_expr *x, gfc_expr *y, int kind)
 
   result = gfc_constant_result (BT_COMPLEX, kind, &x->where);
 
-  mpfr_set_ui (result->value.complex.i, 0, GFC_RND_MODE);
-
   switch (x->ts.type)
     {
     case BT_INTEGER:
       if (!x->is_boz)
-       mpfr_set_z (result->value.complex.r, x->value.integer, GFC_RND_MODE);
+       mpc_set_z (result->value.complex, x->value.integer, GFC_MPC_RND_MODE);
       break;
 
     case BT_REAL:
-      mpfr_set (result->value.complex.r, x->value.real, GFC_RND_MODE);
+      mpc_set_fr (result->value.complex, x->value.real, GFC_RND_MODE);
       break;
 
     case BT_COMPLEX:
-      mpfr_set (result->value.complex.r, x->value.complex.r, GFC_RND_MODE);
-      mpfr_set (result->value.complex.i, x->value.complex.i, GFC_RND_MODE);
+      mpc_set (result->value.complex, x->value.complex, GFC_MPC_RND_MODE);
       break;
 
     default:
@@ -884,12 +1348,13 @@ simplify_cmplx (const char *name, gfc_expr *x, gfc_expr *y, int kind)
        {
        case BT_INTEGER:
          if (!y->is_boz)
-           mpfr_set_z (result->value.complex.i, y->value.integer,
-                       GFC_RND_MODE);
+           mpfr_set_z (mpc_imagref (result->value.complex),
+                       y->value.integer, GFC_RND_MODE);
          break;
 
        case BT_REAL:
-         mpfr_set (result->value.complex.i, y->value.real, GFC_RND_MODE);
+         mpfr_set (mpc_imagref (result->value.complex),
+                   y->value.real, GFC_RND_MODE);
          break;
 
        default:
@@ -906,7 +1371,8 @@ simplify_cmplx (const char *name, gfc_expr *x, gfc_expr *y, int kind)
       ts.type = BT_REAL;
       if (!gfc_convert_boz (x, &ts))
        return &gfc_bad_expr;
-      mpfr_set (result->value.complex.r, x->value.real, GFC_RND_MODE);
+      mpfr_set (mpc_realref (result->value.complex),
+               x->value.real, GFC_RND_MODE);
     }
 
   if (y && y->is_boz)
@@ -917,7 +1383,8 @@ simplify_cmplx (const char *name, gfc_expr *x, gfc_expr *y, int kind)
       ts.type = BT_REAL;
       if (!gfc_convert_boz (y, &ts))
        return &gfc_bad_expr;
-      mpfr_set (result->value.complex.i, y->value.real, GFC_RND_MODE);
+      mpfr_set (mpc_imagref (result->value.complex),
+               y->value.real, GFC_RND_MODE);
     }
 
   return range_check (result, name);
@@ -999,8 +1466,7 @@ gfc_simplify_conjg (gfc_expr *e)
     return NULL;
 
   result = gfc_copy_expr (e);
-  mpfr_neg (result->value.complex.i, result->value.complex.i, GFC_RND_MODE);
-
+  mpc_conj (result->value.complex, result->value.complex, GFC_MPC_RND_MODE);
   return range_check (result, "CONJG");
 }
 
@@ -1009,7 +1475,6 @@ gfc_expr *
 gfc_simplify_cos (gfc_expr *x)
 {
   gfc_expr *result;
-  mpfr_t xp, xq;
 
   if (x->expr_type != EXPR_CONSTANT)
     return NULL;
@@ -1023,19 +1488,7 @@ gfc_simplify_cos (gfc_expr *x)
       break;
     case BT_COMPLEX:
       gfc_set_model_kind (x->ts.kind);
-      mpfr_init (xp);
-      mpfr_init (xq);
-
-      mpfr_cos  (xp, x->value.complex.r, GFC_RND_MODE);
-      mpfr_cosh (xq, x->value.complex.i, GFC_RND_MODE);
-      mpfr_mul (result->value.complex.r, xp, xq, GFC_RND_MODE);
-
-      mpfr_sin  (xp, x->value.complex.r, GFC_RND_MODE);
-      mpfr_sinh (xq, x->value.complex.i, GFC_RND_MODE);
-      mpfr_mul (xp, xp, xq, GFC_RND_MODE);
-      mpfr_neg (result->value.complex.i, xp, GFC_RND_MODE );
-
-      mpfr_clears (xp, xq, NULL);
+      mpc_cos (result->value.complex, x->value.complex, GFC_MPC_RND_MODE);
       break;
     default:
       gfc_internal_error ("in gfc_simplify_cos(): Bad type");
@@ -1056,13 +1509,44 @@ gfc_simplify_cosh (gfc_expr *x)
 
   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
 
-  mpfr_cosh (result->value.real, x->value.real, GFC_RND_MODE);
+  if (x->ts.type == BT_REAL)
+    mpfr_cosh (result->value.real, x->value.real, GFC_RND_MODE);
+  else if (x->ts.type == BT_COMPLEX)
+    mpc_cosh (result->value.complex, x->value.complex, GFC_MPC_RND_MODE);
+  else
+    gcc_unreachable ();
 
   return range_check (result, "COSH");
 }
 
 
 gfc_expr *
+gfc_simplify_count (gfc_expr *mask, gfc_expr *dim, gfc_expr *kind)
+{
+  gfc_expr *result;
+
+  if (!is_constant_array_expr (mask)
+      || !gfc_is_constant_expr (dim)
+      || !gfc_is_constant_expr (kind))
+    return NULL;
+
+  result = transformational_result (mask, dim,
+                                   BT_INTEGER,
+                                   get_kind (BT_INTEGER, kind, "COUNT",
+                                             gfc_default_integer_kind),
+                                   &mask->where);
+
+  init_result_expr (result, 0, NULL);
+
+  /* Passing MASK twice, once as data array, once as mask.
+     Whenever gfc_count is called, '1' is added to the result.  */
+  return !dim || mask->rank == 1 ?
+    simplify_transformation_to_scalar (result, mask, mask, gfc_count) :
+    simplify_transformation_to_array (result, mask, dim, mask, gfc_count);
+}
+
+
+gfc_expr *
 gfc_simplify_dcmplx (gfc_expr *x, gfc_expr *y)
 {
 
@@ -1179,50 +1663,208 @@ gfc_simplify_dim (gfc_expr *x, gfc_expr *y)
       gfc_internal_error ("gfc_simplify_dim(): Bad type");
     }
 
-  return range_check (result, "DIM");
+  return range_check (result, "DIM");
+}
+
+
+gfc_expr*
+gfc_simplify_dot_product (gfc_expr *vector_a, gfc_expr *vector_b)
+{
+  gfc_expr *result;
+
+  if (!is_constant_array_expr (vector_a)
+      || !is_constant_array_expr (vector_b))
+    return NULL;
+
+  gcc_assert (vector_a->rank == 1);
+  gcc_assert (vector_b->rank == 1);
+  gcc_assert (gfc_compare_types (&vector_a->ts, &vector_b->ts));
+
+  if (vector_a->value.constructor && vector_b->value.constructor)
+    return compute_dot_product (vector_a->value.constructor, 1,
+                               vector_b->value.constructor, 1);
+
+  /* Zero sized array ...  */
+  result = gfc_constant_result (vector_a->ts.type,
+                               vector_a->ts.kind,
+                               &vector_a->where);
+  init_result_expr (result, 0, NULL);
+  return result;
+}
+
+
+gfc_expr *
+gfc_simplify_dprod (gfc_expr *x, gfc_expr *y)
+{
+  gfc_expr *a1, *a2, *result;
+
+  if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
+    return NULL;
+
+  result = gfc_constant_result (BT_REAL, gfc_default_double_kind, &x->where);
+
+  a1 = gfc_real2real (x, gfc_default_double_kind);
+  a2 = gfc_real2real (y, gfc_default_double_kind);
+
+  mpfr_mul (result->value.real, a1->value.real, a2->value.real, GFC_RND_MODE);
+
+  gfc_free_expr (a1);
+  gfc_free_expr (a2);
+
+  return range_check (result, "DPROD");
+}
+
+
+gfc_expr *
+gfc_simplify_erf (gfc_expr *x)
+{
+  gfc_expr *result;
+
+  if (x->expr_type != EXPR_CONSTANT)
+    return NULL;
+
+  result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
+
+  mpfr_erf (result->value.real, x->value.real, GFC_RND_MODE);
+
+  return range_check (result, "ERF");
+}
+
+
+gfc_expr *
+gfc_simplify_erfc (gfc_expr *x)
+{
+  gfc_expr *result;
+
+  if (x->expr_type != EXPR_CONSTANT)
+    return NULL;
+
+  result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
+
+  mpfr_erfc (result->value.real, x->value.real, GFC_RND_MODE);
+
+  return range_check (result, "ERFC");
+}
+
+
+/* Helper functions to simplify ERFC_SCALED(x) = ERFC(x) * EXP(X**2).  */
+
+#define MAX_ITER 200
+#define ARG_LIMIT 12
+
+/* Calculate ERFC_SCALED directly by its definition:
+
+     ERFC_SCALED(x) = ERFC(x) * EXP(X**2)
+
+   using a large precision for intermediate results.  This is used for all
+   but large values of the argument.  */
+static void
+fullprec_erfc_scaled (mpfr_t res, mpfr_t arg)
+{
+  mp_prec_t prec;
+  mpfr_t a, b;
+
+  prec = mpfr_get_default_prec ();
+  mpfr_set_default_prec (10 * prec);
+
+  mpfr_init (a);
+  mpfr_init (b);
+
+  mpfr_set (a, arg, GFC_RND_MODE);
+  mpfr_sqr (b, a, GFC_RND_MODE);
+  mpfr_exp (b, b, GFC_RND_MODE);
+  mpfr_erfc (a, a, GFC_RND_MODE);
+  mpfr_mul (a, a, b, GFC_RND_MODE);
+
+  mpfr_set (res, a, GFC_RND_MODE);
+  mpfr_set_default_prec (prec);
+
+  mpfr_clear (a);
+  mpfr_clear (b);
 }
 
+/* Calculate ERFC_SCALED using a power series expansion in 1/arg:
 
-gfc_expr *
-gfc_simplify_dprod (gfc_expr *x, gfc_expr *y)
+    ERFC_SCALED(x) = 1 / (x * sqrt(pi))
+                     * (1 + Sum_n (-1)**n * (1 * 3 * 5 * ... * (2n-1))
+                                          / (2 * x**2)**n)
+
+  This is used for large values of the argument.  Intermediate calculations
+  are performed with twice the precision.  We don't do a fixed number of
+  iterations of the sum, but stop when it has converged to the required
+  precision.  */
+static void
+asympt_erfc_scaled (mpfr_t res, mpfr_t arg)
 {
-  gfc_expr *a1, *a2, *result;
+  mpfr_t sum, x, u, v, w, oldsum, sumtrunc;
+  mpz_t num;
+  mp_prec_t prec;
+  unsigned i;
 
-  if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
-    return NULL;
+  prec = mpfr_get_default_prec ();
+  mpfr_set_default_prec (2 * prec);
 
-  result = gfc_constant_result (BT_REAL, gfc_default_double_kind, &x->where);
+  mpfr_init (sum);
+  mpfr_init (x);
+  mpfr_init (u);
+  mpfr_init (v);
+  mpfr_init (w);
+  mpz_init (num);
 
-  a1 = gfc_real2real (x, gfc_default_double_kind);
-  a2 = gfc_real2real (y, gfc_default_double_kind);
+  mpfr_init (oldsum);
+  mpfr_init (sumtrunc);
+  mpfr_set_prec (oldsum, prec);
+  mpfr_set_prec (sumtrunc, prec);
 
-  mpfr_mul (result->value.real, a1->value.real, a2->value.real, GFC_RND_MODE);
+  mpfr_set (x, arg, GFC_RND_MODE);
+  mpfr_set_ui (sum, 1, GFC_RND_MODE);
+  mpz_set_ui (num, 1);
 
-  gfc_free_expr (a1);
-  gfc_free_expr (a2);
+  mpfr_set (u, x, GFC_RND_MODE);
+  mpfr_sqr (u, u, GFC_RND_MODE);
+  mpfr_mul_ui (u, u, 2, GFC_RND_MODE);
+  mpfr_pow_si (u, u, -1, GFC_RND_MODE);
 
-  return range_check (result, "DPROD");
-}
+  for (i = 1; i < MAX_ITER; i++)
+  {
+    mpfr_set (oldsum, sum, GFC_RND_MODE);
 
+    mpz_mul_ui (num, num, 2 * i - 1);
+    mpz_neg (num, num);
 
-gfc_expr *
-gfc_simplify_erf (gfc_expr *x)
-{
-  gfc_expr *result;
+    mpfr_set (w, u, GFC_RND_MODE);
+    mpfr_pow_ui (w, w, i, GFC_RND_MODE);
 
-  if (x->expr_type != EXPR_CONSTANT)
-    return NULL;
+    mpfr_set_z (v, num, GFC_RND_MODE);
+    mpfr_mul (v, v, w, GFC_RND_MODE);
 
-  result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
+    mpfr_add (sum, sum, v, GFC_RND_MODE);
 
-  mpfr_erf (result->value.real, x->value.real, GFC_RND_MODE);
+    mpfr_set (sumtrunc, sum, GFC_RND_MODE);
+    if (mpfr_cmp (sumtrunc, oldsum) == 0)
+      break;
+  }
 
-  return range_check (result, "ERF");
+  /* We should have converged by now; otherwise, ARG_LIMIT is probably
+     set too low.  */
+  gcc_assert (i < MAX_ITER);
+
+  /* Divide by x * sqrt(Pi).  */
+  mpfr_const_pi (u, GFC_RND_MODE);
+  mpfr_sqrt (u, u, GFC_RND_MODE);
+  mpfr_mul (u, u, x, GFC_RND_MODE);
+  mpfr_div (sum, sum, u, GFC_RND_MODE);
+
+  mpfr_set (res, sum, GFC_RND_MODE);
+  mpfr_set_default_prec (prec);
+
+  mpfr_clears (sum, x, u, v, w, oldsum, sumtrunc, NULL);
+  mpz_clear (num);
 }
 
 
 gfc_expr *
-gfc_simplify_erfc (gfc_expr *x)
+gfc_simplify_erfc_scaled (gfc_expr *x)
 {
   gfc_expr *result;
 
@@ -1230,12 +1872,17 @@ gfc_simplify_erfc (gfc_expr *x)
     return NULL;
 
   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
+  if (mpfr_cmp_d (x->value.real, ARG_LIMIT) >= 0)
+    asympt_erfc_scaled (result->value.real, x->value.real);
+  else
+    fullprec_erfc_scaled (result->value.real, x->value.real);
 
-  mpfr_erfc (result->value.real, x->value.real, GFC_RND_MODE);
-
-  return range_check (result, "ERFC");
+  return range_check (result, "ERFC_SCALED");
 }
 
+#undef MAX_ITER
+#undef ARG_LIMIT
+
 
 gfc_expr *
 gfc_simplify_epsilon (gfc_expr *e)
@@ -1257,7 +1904,6 @@ gfc_expr *
 gfc_simplify_exp (gfc_expr *x)
 {
   gfc_expr *result;
-  mpfr_t xp, xq;
 
   if (x->expr_type != EXPR_CONSTANT)
     return NULL;
@@ -1272,14 +1918,7 @@ gfc_simplify_exp (gfc_expr *x)
 
     case BT_COMPLEX:
       gfc_set_model_kind (x->ts.kind);
-      mpfr_init (xp);
-      mpfr_init (xq);
-      mpfr_exp (xq, x->value.complex.r, GFC_RND_MODE);
-      mpfr_cos (xp, x->value.complex.i, GFC_RND_MODE);
-      mpfr_mul (result->value.complex.r, xq, xp, GFC_RND_MODE);
-      mpfr_sin (xp, x->value.complex.i, GFC_RND_MODE);
-      mpfr_mul (result->value.complex.i, xq, xp, GFC_RND_MODE);
-      mpfr_clears (xp, xq, NULL);
+      mpc_exp (result->value.complex, x->value.complex, GFC_MPC_RND_MODE);
       break;
 
     default:
@@ -1365,7 +2004,7 @@ gfc_simplify_floor (gfc_expr *e, gfc_expr *k)
   mpfr_init (floor);
   mpfr_floor (floor, e->value.real);
 
-  gfc_mpfr_to_mpz (result->value.integer, floor);
+  gfc_mpfr_to_mpz (result->value.integer, floor, &e->where);
 
   mpfr_clear (floor);
 
@@ -1949,7 +2588,7 @@ gfc_simplify_ifix (gfc_expr *e)
   rtrunc = gfc_copy_expr (e);
 
   mpfr_trunc (rtrunc->value.real, e->value.real);
-  gfc_mpfr_to_mpz (result->value.integer, rtrunc->value.real);
+  gfc_mpfr_to_mpz (result->value.integer, rtrunc->value.real, &e->where);
 
   gfc_free_expr (rtrunc);
   return range_check (result, "IFIX");
@@ -1970,7 +2609,7 @@ gfc_simplify_idint (gfc_expr *e)
   rtrunc = gfc_copy_expr (e);
 
   mpfr_trunc (rtrunc->value.real, e->value.real);
-  gfc_mpfr_to_mpz (result->value.integer, rtrunc->value.real);
+  gfc_mpfr_to_mpz (result->value.integer, rtrunc->value.real, &e->where);
 
   gfc_free_expr (rtrunc);
   return range_check (result, "IDINT");
@@ -1993,6 +2632,54 @@ gfc_simplify_ior (gfc_expr *x, gfc_expr *y)
 
 
 gfc_expr *
+gfc_simplify_is_iostat_end (gfc_expr *x)
+{
+  gfc_expr *result;
+
+  if (x->expr_type != EXPR_CONSTANT)
+    return NULL;
+
+  result = gfc_constant_result (BT_LOGICAL, gfc_default_logical_kind,
+                               &x->where);
+  result->value.logical = (mpz_cmp_si (x->value.integer, LIBERROR_END) == 0);
+
+  return result;
+}
+
+
+gfc_expr *
+gfc_simplify_is_iostat_eor (gfc_expr *x)
+{
+  gfc_expr *result;
+
+  if (x->expr_type != EXPR_CONSTANT)
+    return NULL;
+
+  result = gfc_constant_result (BT_LOGICAL, gfc_default_logical_kind,
+                               &x->where);
+  result->value.logical = (mpz_cmp_si (x->value.integer, LIBERROR_EOR) == 0);
+
+  return result;
+}
+
+
+gfc_expr *
+gfc_simplify_isnan (gfc_expr *x)
+{
+  gfc_expr *result;
+
+  if (x->expr_type != EXPR_CONSTANT)
+    return NULL;
+
+  result = gfc_constant_result (BT_LOGICAL, gfc_default_logical_kind,
+                               &x->where);
+  result->value.logical = mpfr_nan_p (x->value.real);
+
+  return result;
+}
+
+
+gfc_expr *
 gfc_simplify_ishft (gfc_expr *e, gfc_expr *s)
 {
   gfc_expr *result;
@@ -2202,7 +2889,7 @@ gfc_simplify_kind (gfc_expr *e)
 
 static gfc_expr *
 simplify_bound_dim (gfc_expr *array, gfc_expr *kind, int d, int upper,
-                   gfc_array_spec *as)
+                   gfc_array_spec *as, gfc_ref *ref)
 {
   gfc_expr *l, *u, *result;
   int k;
@@ -2216,13 +2903,6 @@ simplify_bound_dim (gfc_expr *array, gfc_expr *kind, int d, int upper,
        return NULL;
     }
 
-  /* Then, we need to know the extent of the given dimension.  */
-  l = as->lower[d-1];
-  u = as->upper[d-1];
-
-  if (l->expr_type != EXPR_CONSTANT || u->expr_type != EXPR_CONSTANT)
-    return NULL;
-
   k = get_kind (BT_INTEGER, kind, upper ? "UBOUND" : "LBOUND",
                gfc_default_integer_kind); 
   if (k == -1)
@@ -2230,21 +2910,43 @@ simplify_bound_dim (gfc_expr *array, gfc_expr *kind, int d, int upper,
 
   result = gfc_constant_result (BT_INTEGER, k, &array->where);
 
-  if (mpz_cmp (l->value.integer, u->value.integer) > 0)
+
+  /* Then, we need to know the extent of the given dimension.  */
+  if (ref->u.ar.type == AR_FULL)
     {
-      /* Zero extent.  */
-      if (upper)
-       mpz_set_si (result->value.integer, 0);
+      l = as->lower[d-1];
+      u = as->upper[d-1];
+
+      if (l->expr_type != EXPR_CONSTANT || u->expr_type != EXPR_CONSTANT)
+       return NULL;
+
+      if (mpz_cmp (l->value.integer, u->value.integer) > 0)
+       {
+         /* Zero extent.  */
+         if (upper)
+           mpz_set_si (result->value.integer, 0);
+         else
+           mpz_set_si (result->value.integer, 1);
+       }
       else
-       mpz_set_si (result->value.integer, 1);
+       {
+         /* Nonzero extent.  */
+         if (upper)
+           mpz_set (result->value.integer, u->value.integer);
+         else
+           mpz_set (result->value.integer, l->value.integer);
+       }
     }
   else
     {
-      /* Nonzero extent.  */
       if (upper)
-       mpz_set (result->value.integer, u->value.integer);
+       {
+         if (gfc_ref_dimen_size (&ref->u.ar, d-1, &result->value.integer)
+             != SUCCESS)
+           return NULL;
+       }
       else
-       mpz_set (result->value.integer, l->value.integer);
+       mpz_set_si (result->value.integer, (long int) 1);
     }
 
   return range_check (result, upper ? "UBOUND" : "LBOUND");
@@ -2277,11 +2979,17 @@ simplify_bound (gfc_expr *array, gfc_expr *dim, gfc_expr *kind, int upper)
            case AR_FULL:
              /* We're done because 'as' has already been set in the
                 previous iteration.  */
-             goto done;
+             if (!ref->next)
+               goto done;
+
+           /* Fall through.  */
 
-           case AR_SECTION:
            case AR_UNKNOWN:
              return NULL;
+
+           case AR_SECTION:
+             as = ref->u.ar.as;
+             goto done;
            }
 
          gcc_unreachable ();
@@ -2321,7 +3029,7 @@ simplify_bound (gfc_expr *array, gfc_expr *dim, gfc_expr *kind, int upper)
       /* Simplify the bounds for each dimension.  */
       for (d = 0; d < array->rank; d++)
        {
-         bounds[d] = simplify_bound_dim (array, kind, d + 1, upper, as);
+         bounds[d] = simplify_bound_dim (array, kind, d + 1, upper, as, ref);
          if (bounds[d] == NULL || bounds[d] == &gfc_bad_expr)
            {
              int j;
@@ -2387,7 +3095,7 @@ simplify_bound (gfc_expr *array, gfc_expr *dim, gfc_expr *kind, int upper)
          return &gfc_bad_expr;
        }
 
-      return simplify_bound_dim (array, kind, d, upper, as);
+      return simplify_bound_dim (array, kind, d, upper, as, ref);
     }
 }
 
@@ -2413,10 +3121,13 @@ gfc_simplify_leadz (gfc_expr *e)
   bs = gfc_integer_kinds[i].bit_size;
   if (mpz_cmp_si (e->value.integer, 0) == 0)
     lz = bs;
+  else if (mpz_cmp_si (e->value.integer, 0) < 0)
+    lz = 0;
   else
     lz = bs - mpz_sizeinbase (e->value.integer, 2);
 
-  result = gfc_constant_result (BT_INTEGER, gfc_default_integer_kind, &e->where);
+  result = gfc_constant_result (BT_INTEGER, gfc_default_integer_kind,
+                               &e->where);
   mpz_set_ui (result->value.integer, lz);
 
   return result;
@@ -2436,16 +3147,28 @@ gfc_simplify_len (gfc_expr *e, gfc_expr *kind)
     {
       result = gfc_constant_result (BT_INTEGER, k, &e->where);
       mpz_set_si (result->value.integer, e->value.character.length);
-      return range_check (result, "LEN");
+      if (gfc_range_check (result) == ARITH_OK)
+       return result;
+      else
+       {
+         gfc_free_expr (result);
+         return NULL;
+       }
     }
 
-  if (e->ts.cl != NULL && e->ts.cl->length != NULL
-      && e->ts.cl->length->expr_type == EXPR_CONSTANT
-      && e->ts.cl->length->ts.type == BT_INTEGER)
+  if (e->ts.u.cl != NULL && e->ts.u.cl->length != NULL
+      && e->ts.u.cl->length->expr_type == EXPR_CONSTANT
+      && e->ts.u.cl->length->ts.type == BT_INTEGER)
     {
       result = gfc_constant_result (BT_INTEGER, k, &e->where);
-      mpz_set (result->value.integer, e->ts.cl->length->value.integer);
-      return range_check (result, "LEN");
+      mpz_set (result->value.integer, e->ts.u.cl->length->value.integer);
+      if (gfc_range_check (result) == ARITH_OK)
+       return result;
+      else
+       {
+         gfc_free_expr (result);
+         return NULL;
+       }
     }
 
   return NULL;
@@ -2483,7 +3206,6 @@ gfc_simplify_len_trim (gfc_expr *e, gfc_expr *kind)
 gfc_expr *
 gfc_simplify_lgamma (gfc_expr *x ATTRIBUTE_UNUSED)
 {
-#if MPFR_VERSION >= MPFR_VERSION_NUM(2,3,0)
   gfc_expr *result;
   int sg;
 
@@ -2495,9 +3217,6 @@ gfc_simplify_lgamma (gfc_expr *x ATTRIBUTE_UNUSED)
   mpfr_lgamma (result->value.real, &sg, x->value.real, GFC_RND_MODE);
 
   return range_check (result, "LGAMMA");
-#else
-  return NULL;
-#endif
 }
 
 
@@ -2546,7 +3265,6 @@ gfc_expr *
 gfc_simplify_log (gfc_expr *x)
 {
   gfc_expr *result;
-  mpfr_t xr, xi;
 
   if (x->expr_type != EXPR_CONSTANT)
     return NULL;
@@ -2569,8 +3287,8 @@ gfc_simplify_log (gfc_expr *x)
       break;
 
     case BT_COMPLEX:
-      if ((mpfr_sgn (x->value.complex.r) == 0)
-         && (mpfr_sgn (x->value.complex.i) == 0))
+      if ((mpfr_sgn (mpc_realref (x->value.complex)) == 0)
+         && (mpfr_sgn (mpc_imagref (x->value.complex)) == 0))
        {
          gfc_error ("Complex argument of LOG at %L cannot be zero",
                     &x->where);
@@ -2579,20 +3297,7 @@ gfc_simplify_log (gfc_expr *x)
        }
 
       gfc_set_model_kind (x->ts.kind);
-      mpfr_init (xr);
-      mpfr_init (xi);
-
-      mpfr_atan2 (result->value.complex.i, x->value.complex.i,
-                 x->value.complex.r, GFC_RND_MODE);
-
-      mpfr_mul (xr, x->value.complex.r, x->value.complex.r, GFC_RND_MODE);
-      mpfr_mul (xi, x->value.complex.i, x->value.complex.i, GFC_RND_MODE);
-      mpfr_add (xr, xr, xi, GFC_RND_MODE);
-      mpfr_sqrt (xr, xr, GFC_RND_MODE);
-      mpfr_log (result->value.complex.r, xr, GFC_RND_MODE);
-
-      mpfr_clears (xr, xi, NULL);
-
+      mpc_log (result->value.complex, x->value.complex, GFC_MPC_RND_MODE);
       break;
 
     default:
@@ -2647,6 +3352,156 @@ gfc_simplify_logical (gfc_expr *e, gfc_expr *k)
 }
 
 
+gfc_expr*
+gfc_simplify_matmul (gfc_expr *matrix_a, gfc_expr *matrix_b)
+{
+  gfc_expr *result;
+  gfc_constructor *ma_ctor, *mb_ctor;
+  int row, result_rows, col, result_columns, stride_a, stride_b;
+
+  if (!is_constant_array_expr (matrix_a)
+      || !is_constant_array_expr (matrix_b))
+    return NULL;
+
+  gcc_assert (gfc_compare_types (&matrix_a->ts, &matrix_b->ts));
+  result = gfc_start_constructor (matrix_a->ts.type,
+                                 matrix_a->ts.kind,
+                                 &matrix_a->where);
+
+  if (matrix_a->rank == 1 && matrix_b->rank == 2)
+    {
+      result_rows = 1;
+      result_columns = mpz_get_si (matrix_b->shape[0]);
+      stride_a = 1;
+      stride_b = mpz_get_si (matrix_b->shape[0]);
+
+      result->rank = 1;
+      result->shape = gfc_get_shape (result->rank);
+      mpz_init_set_si (result->shape[0], result_columns);
+    }
+  else if (matrix_a->rank == 2 && matrix_b->rank == 1)
+    {
+      result_rows = mpz_get_si (matrix_b->shape[0]);
+      result_columns = 1;
+      stride_a = mpz_get_si (matrix_a->shape[0]);
+      stride_b = 1;
+
+      result->rank = 1;
+      result->shape = gfc_get_shape (result->rank);
+      mpz_init_set_si (result->shape[0], result_rows);
+    }
+  else if (matrix_a->rank == 2 && matrix_b->rank == 2)
+    {
+      result_rows = mpz_get_si (matrix_a->shape[0]);
+      result_columns = mpz_get_si (matrix_b->shape[1]);
+      stride_a = mpz_get_si (matrix_a->shape[1]);
+      stride_b = mpz_get_si (matrix_b->shape[0]);
+
+      result->rank = 2;
+      result->shape = gfc_get_shape (result->rank);
+      mpz_init_set_si (result->shape[0], result_rows);
+      mpz_init_set_si (result->shape[1], result_columns);
+    }
+  else
+    gcc_unreachable();
+
+  ma_ctor = matrix_a->value.constructor;
+  mb_ctor = matrix_b->value.constructor;
+
+  for (col = 0; col < result_columns; ++col)
+    {
+      ma_ctor = matrix_a->value.constructor;
+
+      for (row = 0; row < result_rows; ++row)
+       {
+         gfc_expr *e;
+         e = compute_dot_product (ma_ctor, stride_a,
+                                  mb_ctor, 1);
+
+         gfc_append_constructor (result, e);
+
+         ADVANCE (ma_ctor, 1);
+       }
+
+      ADVANCE (mb_ctor, stride_b);
+    }
+
+  return result;
+}
+
+
+gfc_expr *
+gfc_simplify_merge (gfc_expr *tsource, gfc_expr *fsource, gfc_expr *mask)
+{
+  if (tsource->expr_type != EXPR_CONSTANT
+      || fsource->expr_type != EXPR_CONSTANT
+      || mask->expr_type != EXPR_CONSTANT)
+    return NULL;
+
+  return gfc_copy_expr (mask->value.logical ? tsource : fsource);
+}
+
+
+/* Selects bewteen current value and extremum for simplify_min_max
+   and simplify_minval_maxval.  */
+static void
+min_max_choose (gfc_expr *arg, gfc_expr *extremum, int sign)
+{
+  switch (arg->ts.type)
+    {
+      case BT_INTEGER:
+       if (mpz_cmp (arg->value.integer,
+                       extremum->value.integer) * sign > 0)
+       mpz_set (extremum->value.integer, arg->value.integer);
+       break;
+
+      case BT_REAL:
+       /* We need to use mpfr_min and mpfr_max to treat NaN properly.  */
+       if (sign > 0)
+         mpfr_max (extremum->value.real, extremum->value.real,
+                     arg->value.real, GFC_RND_MODE);
+       else
+         mpfr_min (extremum->value.real, extremum->value.real,
+                     arg->value.real, GFC_RND_MODE);
+       break;
+
+      case BT_CHARACTER:
+#define LENGTH(x) ((x)->value.character.length)
+#define STRING(x) ((x)->value.character.string)
+       if (LENGTH(extremum) < LENGTH(arg))
+         {
+           gfc_char_t *tmp = STRING(extremum);
+
+           STRING(extremum) = gfc_get_wide_string (LENGTH(arg) + 1);
+           memcpy (STRING(extremum), tmp,
+                     LENGTH(extremum) * sizeof (gfc_char_t));
+           gfc_wide_memset (&STRING(extremum)[LENGTH(extremum)], ' ',
+                              LENGTH(arg) - LENGTH(extremum));
+           STRING(extremum)[LENGTH(arg)] = '\0';  /* For debugger  */
+           LENGTH(extremum) = LENGTH(arg);
+           gfc_free (tmp);
+         }
+
+       if (gfc_compare_string (arg, extremum) * sign > 0)
+         {
+           gfc_free (STRING(extremum));
+           STRING(extremum) = gfc_get_wide_string (LENGTH(extremum) + 1);
+           memcpy (STRING(extremum), STRING(arg),
+                     LENGTH(arg) * sizeof (gfc_char_t));
+           gfc_wide_memset (&STRING(extremum)[LENGTH(arg)], ' ',
+                              LENGTH(extremum) - LENGTH(arg));
+           STRING(extremum)[LENGTH(extremum)] = '\0';  /* For debugger  */
+         }
+#undef LENGTH
+#undef STRING
+       break;
+             
+      default:
+       gfc_internal_error ("simplify_min_max(): Bad type in arglist");
+    }
+}
+
+
 /* This function is special since MAX() can take any number of
    arguments.  The simplified expression is a rewritten version of the
    argument list containing at most one constant element.  Other
@@ -2677,59 +3532,7 @@ simplify_min_max (gfc_expr *expr, int sign)
          continue;
        }
 
-      switch (arg->expr->ts.type)
-       {
-       case BT_INTEGER:
-         if (mpz_cmp (arg->expr->value.integer,
-                      extremum->expr->value.integer) * sign > 0)
-           mpz_set (extremum->expr->value.integer, arg->expr->value.integer);
-         break;
-
-       case BT_REAL:
-         /* We need to use mpfr_min and mpfr_max to treat NaN properly.  */
-         if (sign > 0)
-           mpfr_max (extremum->expr->value.real, extremum->expr->value.real,
-                     arg->expr->value.real, GFC_RND_MODE);
-         else
-           mpfr_min (extremum->expr->value.real, extremum->expr->value.real,
-                     arg->expr->value.real, GFC_RND_MODE);
-         break;
-
-       case BT_CHARACTER:
-#define LENGTH(x) ((x)->expr->value.character.length)
-#define STRING(x) ((x)->expr->value.character.string)
-         if (LENGTH(extremum) < LENGTH(arg))
-           {
-             gfc_char_t *tmp = STRING(extremum);
-
-             STRING(extremum) = gfc_get_wide_string (LENGTH(arg) + 1);
-             memcpy (STRING(extremum), tmp,
-                     LENGTH(extremum) * sizeof (gfc_char_t));
-             gfc_wide_memset (&STRING(extremum)[LENGTH(extremum)], ' ',
-                              LENGTH(arg) - LENGTH(extremum));
-             STRING(extremum)[LENGTH(arg)] = '\0';  /* For debugger  */
-             LENGTH(extremum) = LENGTH(arg);
-             gfc_free (tmp);
-           }
-
-         if (gfc_compare_string (arg->expr, extremum->expr) * sign > 0)
-           {
-             gfc_free (STRING(extremum));
-             STRING(extremum) = gfc_get_wide_string (LENGTH(extremum) + 1);
-             memcpy (STRING(extremum), STRING(arg),
-                     LENGTH(arg) * sizeof (gfc_char_t));
-             gfc_wide_memset (&STRING(extremum)[LENGTH(arg)], ' ',
-                              LENGTH(extremum) - LENGTH(arg));
-             STRING(extremum)[LENGTH(extremum)] = '\0';  /* For debugger  */
-           }
-#undef LENGTH
-#undef STRING
-         break;
-             
-
-       default:
-         gfc_internal_error ("simplify_min_max(): Bad type in arglist");
-       }
+      min_max_choose (arg->expr, extremum->expr, sign);
 
       /* Delete the extra constant argument.  */
       if (last == NULL)
@@ -2742,35 +3545,98 @@ simplify_min_max (gfc_expr *expr, int sign)
       arg = last;
     }
 
-  /* If there is one value left, replace the function call with the
-     expression.  */
-  if (expr->value.function.actual->next != NULL)
+  /* If there is one value left, replace the function call with the
+     expression.  */
+  if (expr->value.function.actual->next != NULL)
+    return NULL;
+
+  /* Convert to the correct type and kind.  */
+  if (expr->ts.type != BT_UNKNOWN) 
+    return gfc_convert_constant (expr->value.function.actual->expr,
+       expr->ts.type, expr->ts.kind);
+
+  if (specific->ts.type != BT_UNKNOWN) 
+    return gfc_convert_constant (expr->value.function.actual->expr,
+       specific->ts.type, specific->ts.kind); 
+  return gfc_copy_expr (expr->value.function.actual->expr);
+}
+
+
+gfc_expr *
+gfc_simplify_min (gfc_expr *e)
+{
+  return simplify_min_max (e, -1);
+}
+
+
+gfc_expr *
+gfc_simplify_max (gfc_expr *e)
+{
+  return simplify_min_max (e, 1);
+}
+
+
+/* This is a simplified version of simplify_min_max to provide
+   simplification of minval and maxval for a vector.  */
+
+static gfc_expr *
+simplify_minval_maxval (gfc_expr *expr, int sign)
+{
+  gfc_constructor *ctr, *extremum;
+  gfc_intrinsic_sym * specific;
+
+  extremum = NULL;
+  specific = expr->value.function.isym;
+
+  ctr = expr->value.constructor;
+
+  for (; ctr; ctr = ctr->next)
+    {
+      if (ctr->expr->expr_type != EXPR_CONSTANT)
+       return NULL;
+
+      if (extremum == NULL)
+       {
+         extremum = ctr;
+         continue;
+       }
+
+      min_max_choose (ctr->expr, extremum->expr, sign);
+     }
+
+  if (extremum == NULL)
     return NULL;
 
   /* Convert to the correct type and kind.  */
   if (expr->ts.type != BT_UNKNOWN) 
-    return gfc_convert_constant (expr->value.function.actual->expr,
+    return gfc_convert_constant (extremum->expr,
        expr->ts.type, expr->ts.kind);
 
   if (specific->ts.type != BT_UNKNOWN) 
-    return gfc_convert_constant (expr->value.function.actual->expr,
+    return gfc_convert_constant (extremum->expr,
        specific->ts.type, specific->ts.kind); 
  
-  return gfc_copy_expr (expr->value.function.actual->expr);
+  return gfc_copy_expr (extremum->expr);
 }
 
 
 gfc_expr *
-gfc_simplify_min (gfc_expr *e)
+gfc_simplify_minval (gfc_expr *array, gfc_expr* dim, gfc_expr *mask)
 {
-  return simplify_min_max (e, -1);
+  if (array->expr_type != EXPR_ARRAY || array->rank != 1 || dim || mask)
+    return NULL;
+  
+  return simplify_minval_maxval (array, -1);
 }
 
 
 gfc_expr *
-gfc_simplify_max (gfc_expr *e)
+gfc_simplify_maxval (gfc_expr *array, gfc_expr* dim, gfc_expr *mask)
 {
-  return simplify_min_max (e, 1);
+  if (array->expr_type != EXPR_ARRAY || array->rank != 1 || dim || mask)
+    return NULL;
+  return simplify_minval_maxval (array, 1);
 }
 
 
@@ -2950,6 +3816,7 @@ gfc_simplify_nearest (gfc_expr *x, gfc_expr *s)
   mpfr_set_emin ((mp_exp_t) gfc_real_kinds[kind].min_exponent -
                mpfr_get_prec(result->value.real) + 1);
   mpfr_set_emax ((mp_exp_t) gfc_real_kinds[kind].max_exponent - 1);
+  mpfr_check_range (result->value.real, 0, GMP_RNDU);
 
   if (mpfr_sgn (s->value.real) > 0)
     {
@@ -2997,7 +3864,7 @@ simplify_nint (const char *name, gfc_expr *e, gfc_expr *k)
 
   mpfr_round (itrunc->value.real, e->value.real);
 
-  gfc_mpfr_to_mpz (result->value.integer, itrunc->value.real);
+  gfc_mpfr_to_mpz (result->value.integer, itrunc->value.real, &e->where);
 
   gfc_free_expr (itrunc);
 
@@ -3093,6 +3960,75 @@ gfc_simplify_or (gfc_expr *x, gfc_expr *y)
 
 
 gfc_expr *
+gfc_simplify_pack (gfc_expr *array, gfc_expr *mask, gfc_expr *vector)
+{
+  gfc_expr *result;
+  gfc_constructor *array_ctor, *mask_ctor, *vector_ctor;
+
+  if (!is_constant_array_expr(array)
+      || !is_constant_array_expr(vector)
+      || (!gfc_is_constant_expr (mask)
+          && !is_constant_array_expr(mask)))
+    return NULL;
+
+  result = gfc_start_constructor (array->ts.type, 
+                                 array->ts.kind,
+                                 &array->where);
+
+  array_ctor = array->value.constructor;
+  vector_ctor = vector ? vector->value.constructor : NULL;
+
+  if (mask->expr_type == EXPR_CONSTANT
+      && mask->value.logical)
+    {
+      /* Copy all elements of ARRAY to RESULT.  */
+      while (array_ctor)
+       {
+         gfc_append_constructor (result, 
+                                 gfc_copy_expr (array_ctor->expr));
+
+         ADVANCE (array_ctor, 1);
+         ADVANCE (vector_ctor, 1);
+       }
+    }
+  else if (mask->expr_type == EXPR_ARRAY)
+    {
+      /* Copy only those elements of ARRAY to RESULT whose 
+        MASK equals .TRUE..  */
+      mask_ctor = mask->value.constructor;
+      while (mask_ctor)
+       {
+         if (mask_ctor->expr->value.logical)
+           {
+             gfc_append_constructor (result, 
+                                     gfc_copy_expr (array_ctor->expr)); 
+             ADVANCE (vector_ctor, 1);
+           }
+
+         ADVANCE (array_ctor, 1);
+         ADVANCE (mask_ctor, 1);
+       }
+    }
+
+  /* Append any left-over elements from VECTOR to RESULT.  */
+  while (vector_ctor)
+    {
+      gfc_append_constructor (result, 
+                             gfc_copy_expr (vector_ctor->expr));
+      ADVANCE (vector_ctor, 1);
+    }
+
+  result->shape = gfc_get_shape (1);
+  gfc_array_size (result, &result->shape[0]);
+
+  if (array->ts.type == BT_CHARACTER)
+    result->ts.u.cl = array->ts.u.cl;
+
+  return result;
+}
+
+
+gfc_expr *
 gfc_simplify_precision (gfc_expr *e)
 {
   gfc_expr *result;
@@ -3108,6 +4044,30 @@ gfc_simplify_precision (gfc_expr *e)
 
 
 gfc_expr *
+gfc_simplify_product (gfc_expr *array, gfc_expr *dim, gfc_expr *mask)
+{
+  gfc_expr *result;
+
+  if (!is_constant_array_expr (array)
+      || !gfc_is_constant_expr (dim))
+    return NULL;
+
+  if (mask
+      && !is_constant_array_expr (mask)
+      && mask->expr_type != EXPR_CONSTANT)
+    return NULL;
+
+  result = transformational_result (array, dim, array->ts.type,
+                                   array->ts.kind, &array->where);
+  init_result_expr (result, 1, NULL);
+
+  return !dim || array->rank == 1 ?
+    simplify_transformation_to_scalar (result, array, mask, gfc_multiply) :
+    simplify_transformation_to_array (result, array, dim, mask, gfc_multiply);
+}
+
+
+gfc_expr *
 gfc_simplify_radix (gfc_expr *e)
 {
   gfc_expr *result;
@@ -3230,8 +4190,7 @@ gfc_simplify_realpart (gfc_expr *e)
     return NULL;
 
   result = gfc_constant_result (BT_REAL, e->ts.kind, &e->where);
-  mpfr_set (result->value.real, e->value.complex.r, GFC_RND_MODE);
-
+  mpc_real (result->value.real, e->value.complex, GFC_RND_MODE);
   return range_check (result, "REALPART");
 }
 
@@ -3256,14 +4215,14 @@ gfc_simplify_repeat (gfc_expr *e, gfc_expr *n)
     }
 
   /* If we don't know the character length, we can do no more.  */
-  if (e->ts.cl && e->ts.cl->length
-       && e->ts.cl->length->expr_type == EXPR_CONSTANT)
+  if (e->ts.u.cl && e->ts.u.cl->length
+       && e->ts.u.cl->length->expr_type == EXPR_CONSTANT)
     {
-      len = mpz_get_si (e->ts.cl->length->value.integer);
+      len = mpz_get_si (e->ts.u.cl->length->value.integer);
       have_length = true;
     }
   else if (e->expr_type == EXPR_CONSTANT
-            && (e->ts.cl == NULL || e->ts.cl->length == NULL))
+            && (e->ts.u.cl == NULL || e->ts.u.cl->length == NULL))
     {
       len = e->value.character.length;
     }
@@ -3291,7 +4250,7 @@ gfc_simplify_repeat (gfc_expr *e, gfc_expr *n)
       if (have_length)
        {
          mpz_tdiv_q (max, gfc_integer_kinds[i].huge,
-                     e->ts.cl->length->value.integer);
+                     e->ts.u.cl->length->value.integer);
        }
       else
        {
@@ -3320,8 +4279,8 @@ gfc_simplify_repeat (gfc_expr *e, gfc_expr *n)
     return NULL;
 
   if (len || 
-      (e->ts.cl->length && 
-       mpz_sgn (e->ts.cl->length->value.integer)) != 0)
+      (e->ts.u.cl->length && 
+       mpz_sgn (e->ts.u.cl->length->value.integer)) != 0)
     {
       const char *res = gfc_extract_int (n, &ncop);
       gcc_assert (res == NULL);
@@ -3354,30 +4313,6 @@ gfc_simplify_repeat (gfc_expr *e, gfc_expr *n)
 }
 
 
-/* Test that the expression is an constant array.  */
-
-static bool
-is_constant_array_expr (gfc_expr *e)
-{
-  gfc_constructor *c;
-
-  if (e == NULL)
-    return true;
-
-  if (e->expr_type != EXPR_ARRAY || !gfc_is_constant_expr (e))
-    return false;
-  
-  if (e->value.constructor == NULL)
-    return false;
-  
-  for (c = e->value.constructor; c; c = c->next)
-    if (c->expr->expr_type != EXPR_CONSTANT)
-      return false;
-
-  return true;
-}
-
-
 /* This one is a bear, but mainly has to do with shuffling elements.  */
 
 gfc_expr *
@@ -3393,16 +4328,10 @@ gfc_simplify_reshape (gfc_expr *source, gfc_expr *shape_exp,
   gfc_expr *e;
 
   /* Check that argument expression types are OK.  */
-  if (!is_constant_array_expr (source))
-    return NULL;
-
-  if (!is_constant_array_expr (shape_exp))
-    return NULL;
-
-  if (!is_constant_array_expr (pad))
-    return NULL;
-
-  if (!is_constant_array_expr (order_exp))
+  if (!is_constant_array_expr (source)
+      || !is_constant_array_expr (shape_exp)
+      || !is_constant_array_expr (pad)
+      || !is_constant_array_expr (order_exp))
     return NULL;
 
   /* Proceed with simplification, unpacking the array.  */
@@ -3417,40 +4346,16 @@ gfc_simplify_reshape (gfc_expr *source, gfc_expr *shape_exp,
       if (e == NULL)
        break;
 
-      if (gfc_extract_int (e, &shape[rank]) != NULL)
-       {
-         gfc_error ("Integer too large in shape specification at %L",
-                    &e->where);
-         gfc_free_expr (e);
-         goto bad_reshape;
-       }
-
-      if (rank >= GFC_MAX_DIMENSIONS)
-       {
-         gfc_error ("Too many dimensions in shape specification for RESHAPE "
-                    "at %L", &e->where);
-         gfc_free_expr (e);
-         goto bad_reshape;
-       }
+      gfc_extract_int (e, &shape[rank]);
 
-      if (shape[rank] < 0)
-       {
-         gfc_error ("Shape specification at %L cannot be negative",
-                    &e->where);
-         gfc_free_expr (e);
-         goto bad_reshape;
-       }
+      gcc_assert (rank >= 0 && rank < GFC_MAX_DIMENSIONS);
+      gcc_assert (shape[rank] >= 0);
 
       gfc_free_expr (e);
       rank++;
     }
 
-  if (rank == 0)
-    {
-      gfc_error ("Shape specification at %L cannot be the null array",
-                &shape_exp->where);
-      goto bad_reshape;
-    }
+  gcc_assert (rank > 0);
 
   /* Now unpack the order array if present.  */
   if (order_exp == NULL)
@@ -3466,41 +4371,14 @@ gfc_simplify_reshape (gfc_expr *source, gfc_expr *shape_exp,
       for (i = 0; i < rank; i++)
        {
          e = gfc_get_array_element (order_exp, i);
-         if (e == NULL)
-           {
-             gfc_error ("ORDER parameter of RESHAPE at %L is not the same "
-                        "size as SHAPE parameter", &order_exp->where);
-             goto bad_reshape;
-           }
-
-         if (gfc_extract_int (e, &order[i]) != NULL)
-           {
-             gfc_error ("Error in ORDER parameter of RESHAPE at %L",
-                        &e->where);
-             gfc_free_expr (e);
-             goto bad_reshape;
-           }
-
-         if (order[i] < 1 || order[i] > rank)
-           {
-             gfc_error ("ORDER parameter of RESHAPE at %L is out of range",
-                        &e->where);
-             gfc_free_expr (e);
-             goto bad_reshape;
-           }
-
-         order[i]--;
-
-         if (x[order[i]])
-           {
-             gfc_error ("Invalid permutation in ORDER parameter at %L",
-                        &e->where);
-             gfc_free_expr (e);
-             goto bad_reshape;
-           }
+         gcc_assert (e);
 
+         gfc_extract_int (e, &order[i]);
          gfc_free_expr (e);
 
+         gcc_assert (order[i] >= 1 && order[i] <= rank);
+         order[i]--;
+         gcc_assert (x[order[i]] == 0);
          x[order[i]] = 1;
        }
     }
@@ -3527,7 +4405,7 @@ gfc_simplify_reshape (gfc_expr *source, gfc_expr *shape_exp,
   for (i = 0; i < rank; i++)
     x[i] = 0;
 
-  for (;;)
+  while (nsource > 0 || npad > 0)
     {
       /* Figure out which element to extract.  */
       mpz_set_ui (index, 0);
@@ -3548,18 +4426,13 @@ gfc_simplify_reshape (gfc_expr *source, gfc_expr *shape_exp,
        e = gfc_get_array_element (source, j);
       else
        {
-         j = j - nsource;
-
-         if (npad == 0)
-           {
-             gfc_error ("PAD parameter required for short SOURCE parameter "
-                        "at %L", &source->where);
-             goto bad_reshape;
-           }
+         gcc_assert (npad > 0);
 
+         j = j - nsource;
          j = j % npad;
          e = gfc_get_array_element (pad, j);
        }
+      gcc_assert (e);
 
       if (head == NULL)
        head = tail = gfc_get_constructor ();
@@ -3569,9 +4442,6 @@ gfc_simplify_reshape (gfc_expr *source, gfc_expr *shape_exp,
          tail = tail->next;
        }
 
-      if (e == NULL)
-       goto bad_reshape;
-
       tail->where = e->where;
       tail->expr = e;
 
@@ -3603,11 +4473,6 @@ inc:
   e->rank = rank;
 
   return e;
-
-bad_reshape:
-  gfc_free_constructor (head);
-  mpz_clear (index);
-  return &gfc_bad_expr;
 }
 
 
@@ -4067,16 +4932,15 @@ gfc_simplify_sign (gfc_expr *x, gfc_expr *y)
       mpz_abs (result->value.integer, x->value.integer);
       if (mpz_sgn (y->value.integer) < 0)
        mpz_neg (result->value.integer, result->value.integer);
-
       break;
 
     case BT_REAL:
-      /* TODO: Handle -0.0 and +0.0 correctly on machines that support
-        it.  */
-      mpfr_abs (result->value.real, x->value.real, GFC_RND_MODE);
-      if (mpfr_sgn (y->value.real) < 0)
-       mpfr_neg (result->value.real, result->value.real, GFC_RND_MODE);
-
+      if (gfc_option.flag_sign_zero)
+       mpfr_copysign (result->value.real, x->value.real, y->value.real,
+                      GFC_RND_MODE);
+      else
+       mpfr_setsign (result->value.real, x->value.real,
+                     mpfr_sgn (y->value.real) < 0 ? 1 : 0, GFC_RND_MODE);
       break;
 
     default:
@@ -4091,7 +4955,6 @@ gfc_expr *
 gfc_simplify_sin (gfc_expr *x)
 {
   gfc_expr *result;
-  mpfr_t xp, xq;
 
   if (x->expr_type != EXPR_CONSTANT)
     return NULL;
@@ -4106,18 +4969,7 @@ gfc_simplify_sin (gfc_expr *x)
 
     case BT_COMPLEX:
       gfc_set_model (x->value.real);
-      mpfr_init (xp);
-      mpfr_init (xq);
-
-      mpfr_sin  (xp, x->value.complex.r, GFC_RND_MODE);
-      mpfr_cosh (xq, x->value.complex.i, GFC_RND_MODE);
-      mpfr_mul (result->value.complex.r, xp, xq, GFC_RND_MODE);
-
-      mpfr_cos  (xp, x->value.complex.r, GFC_RND_MODE);
-      mpfr_sinh (xq, x->value.complex.i, GFC_RND_MODE);
-      mpfr_mul (result->value.complex.i, xp, xq, GFC_RND_MODE);
-
-      mpfr_clears (xp, xq, NULL);
+      mpc_sin (result->value.complex, x->value.complex, GFC_MPC_RND_MODE);
       break;
 
     default:
@@ -4138,7 +4990,13 @@ gfc_simplify_sinh (gfc_expr *x)
 
   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
 
-  mpfr_sinh (result->value.real, x->value.real, GFC_RND_MODE);
+  if (x->ts.type == BT_REAL)
+    mpfr_sinh (result->value.real, x->value.real, GFC_RND_MODE);
+  else if (x->ts.type == BT_COMPLEX)
+    mpc_sinh (result->value.complex, x->value.complex, GFC_MPC_RND_MODE);
+  else
+    gcc_unreachable ();
+
 
   return range_check (result, "SINH");
 }
@@ -4199,10 +5057,116 @@ gfc_simplify_spacing (gfc_expr *x)
 
 
 gfc_expr *
+gfc_simplify_spread (gfc_expr *source, gfc_expr *dim_expr, gfc_expr *ncopies_expr)
+{
+  gfc_expr *result = 0L;
+  int i, j, dim, ncopies;
+  mpz_t size;
+
+  if ((!gfc_is_constant_expr (source)
+       && !is_constant_array_expr (source))
+      || !gfc_is_constant_expr (dim_expr)
+      || !gfc_is_constant_expr (ncopies_expr))
+    return NULL;
+
+  gcc_assert (dim_expr->ts.type == BT_INTEGER);
+  gfc_extract_int (dim_expr, &dim);
+  dim -= 1;   /* zero-base DIM */
+
+  gcc_assert (ncopies_expr->ts.type == BT_INTEGER);
+  gfc_extract_int (ncopies_expr, &ncopies);
+  ncopies = MAX (ncopies, 0);
+
+  /* Do not allow the array size to exceed the limit for an array
+     constructor.  */
+  if (source->expr_type == EXPR_ARRAY)
+    {
+      if (gfc_array_size (source, &size) == FAILURE)
+       gfc_internal_error ("Failure getting length of a constant array.");
+    }
+  else
+    mpz_init_set_ui (size, 1);
+
+  if (mpz_get_si (size)*ncopies > gfc_option.flag_max_array_constructor)
+    return NULL;
+
+  if (source->expr_type == EXPR_CONSTANT)
+    {
+      gcc_assert (dim == 0);
+
+      result = gfc_start_constructor (source->ts.type,
+                                     source->ts.kind,
+                                     &source->where);
+      result->rank = 1;
+      result->shape = gfc_get_shape (result->rank);
+      mpz_init_set_si (result->shape[0], ncopies);
+
+      for (i = 0; i < ncopies; ++i)
+        gfc_append_constructor (result, gfc_copy_expr (source));
+    }
+  else if (source->expr_type == EXPR_ARRAY)
+    {
+      int result_size, rstride[GFC_MAX_DIMENSIONS], extent[GFC_MAX_DIMENSIONS];
+      gfc_constructor *ctor, *source_ctor, *result_ctor;
+
+      gcc_assert (source->rank < GFC_MAX_DIMENSIONS);
+      gcc_assert (dim >= 0 && dim <= source->rank);
+
+      result = gfc_start_constructor (source->ts.type,
+                                     source->ts.kind,
+                                     &source->where);
+      result->rank = source->rank + 1;
+      result->shape = gfc_get_shape (result->rank);
+
+      result_size = 1;
+      for (i = 0, j = 0; i < result->rank; ++i)
+       {
+         if (i != dim)
+           mpz_init_set (result->shape[i], source->shape[j++]);
+         else
+           mpz_init_set_si (result->shape[i], ncopies);
+
+         extent[i] = mpz_get_si (result->shape[i]);
+         rstride[i] = (i == 0) ? 1 : rstride[i-1] * extent[i-1];
+         result_size *= extent[i];
+       }
+
+      for (i = 0; i < result_size; ++i)
+       gfc_append_constructor (result, NULL);
+
+      source_ctor = source->value.constructor;
+      result_ctor = result->value.constructor;
+      while (source_ctor)
+       {
+         ctor = result_ctor;
+
+         for (i = 0; i < ncopies; ++i)
+         {
+           ctor->expr = gfc_copy_expr (source_ctor->expr);
+           ADVANCE (ctor, rstride[dim]);
+         }
+
+         ADVANCE (result_ctor, (dim == 0 ? ncopies : 1));
+         ADVANCE (source_ctor, 1);
+       }
+    }
+  else
+    /* FIXME: Returning here avoids a regression in array_simplify_1.f90.
+       Replace NULL with gcc_unreachable() after implementing
+       gfc_simplify_cshift(). */
+    return NULL;
+
+  if (source->ts.type == BT_CHARACTER)
+    result->ts.u.cl = source->ts.u.cl;
+
+  return result;
+}
+
+
+gfc_expr *
 gfc_simplify_sqrt (gfc_expr *e)
 {
   gfc_expr *result;
-  mpfr_t ac, ad, s, t, w;
 
   if (e->expr_type != EXPR_CONSTANT)
     return NULL;
@@ -4219,82 +5183,8 @@ gfc_simplify_sqrt (gfc_expr *e)
       break;
 
     case BT_COMPLEX:
-      /* Formula taken from Numerical Recipes to avoid over- and
-        underflow.  */
-
       gfc_set_model (e->value.real);
-      mpfr_init (ac);
-      mpfr_init (ad);
-      mpfr_init (s);
-      mpfr_init (t);
-      mpfr_init (w);
-
-      if (mpfr_cmp_ui (e->value.complex.r, 0) == 0
-         && mpfr_cmp_ui (e->value.complex.i, 0) == 0)
-       {
-         mpfr_set_ui (result->value.complex.r, 0, GFC_RND_MODE);
-         mpfr_set_ui (result->value.complex.i, 0, GFC_RND_MODE);
-         break;
-       }
-
-      mpfr_abs (ac, e->value.complex.r, GFC_RND_MODE);
-      mpfr_abs (ad, e->value.complex.i, GFC_RND_MODE);
-
-      if (mpfr_cmp (ac, ad) >= 0)
-       {
-         mpfr_div (t, e->value.complex.i, e->value.complex.r, GFC_RND_MODE);
-         mpfr_mul (t, t, t, GFC_RND_MODE);
-         mpfr_add_ui (t, t, 1, GFC_RND_MODE);
-         mpfr_sqrt (t, t, GFC_RND_MODE);
-         mpfr_add_ui (t, t, 1, GFC_RND_MODE);
-         mpfr_div_ui (t, t, 2, GFC_RND_MODE);
-         mpfr_sqrt (t, t, GFC_RND_MODE);
-         mpfr_sqrt (s, ac, GFC_RND_MODE);
-         mpfr_mul (w, s, t, GFC_RND_MODE);
-       }
-      else
-       {
-         mpfr_div (s, e->value.complex.r, e->value.complex.i, GFC_RND_MODE);
-         mpfr_mul (t, s, s, GFC_RND_MODE);
-         mpfr_add_ui (t, t, 1, GFC_RND_MODE);
-         mpfr_sqrt (t, t, GFC_RND_MODE);
-         mpfr_abs (s, s, GFC_RND_MODE);
-         mpfr_add (t, t, s, GFC_RND_MODE);
-         mpfr_div_ui (t, t, 2, GFC_RND_MODE);
-         mpfr_sqrt (t, t, GFC_RND_MODE);
-         mpfr_sqrt (s, ad, GFC_RND_MODE);
-         mpfr_mul (w, s, t, GFC_RND_MODE);
-       }
-
-      if (mpfr_cmp_ui (w, 0) != 0 && mpfr_cmp_ui (e->value.complex.r, 0) >= 0)
-       {
-         mpfr_mul_ui (t, w, 2, GFC_RND_MODE);
-         mpfr_div (result->value.complex.i, e->value.complex.i, t, GFC_RND_MODE);
-         mpfr_set (result->value.complex.r, w, GFC_RND_MODE);
-       }
-      else if (mpfr_cmp_ui (w, 0) != 0
-              && mpfr_cmp_ui (e->value.complex.r, 0) < 0
-              && mpfr_cmp_ui (e->value.complex.i, 0) >= 0)
-       {
-         mpfr_mul_ui (t, w, 2, GFC_RND_MODE);
-         mpfr_div (result->value.complex.r, e->value.complex.i, t, GFC_RND_MODE);
-         mpfr_set (result->value.complex.i, w, GFC_RND_MODE);
-       }
-      else if (mpfr_cmp_ui (w, 0) != 0
-              && mpfr_cmp_ui (e->value.complex.r, 0) < 0
-              && mpfr_cmp_ui (e->value.complex.i, 0) < 0)
-       {
-         mpfr_mul_ui (t, w, 2, GFC_RND_MODE);
-         mpfr_div (result->value.complex.r, ad, t, GFC_RND_MODE);
-         mpfr_neg (w, w, GFC_RND_MODE);
-         mpfr_set (result->value.complex.i, w, GFC_RND_MODE);
-       }
-      else
-       gfc_internal_error ("invalid complex argument of SQRT at %L",
-                           &e->where);
-
-      mpfr_clears (s, t, ac, ad, w, NULL);
-
+      mpc_sqrt (result->value.complex, e->value.complex, GFC_MPC_RND_MODE);
       break;
 
     default:
@@ -4311,19 +5201,45 @@ negative_arg:
 
 
 gfc_expr *
+gfc_simplify_sum (gfc_expr *array, gfc_expr *dim, gfc_expr *mask)
+{
+  gfc_expr *result;
+
+  if (!is_constant_array_expr (array)
+      || !gfc_is_constant_expr (dim))
+    return NULL;
+
+  if (mask
+      && !is_constant_array_expr (mask)
+      && mask->expr_type != EXPR_CONSTANT)
+    return NULL;
+
+  result = transformational_result (array, dim, array->ts.type,
+                                   array->ts.kind, &array->where);
+  init_result_expr (result, 0, NULL);
+
+  return !dim || array->rank == 1 ?
+    simplify_transformation_to_scalar (result, array, mask, gfc_add) :
+    simplify_transformation_to_array (result, array, dim, mask, gfc_add);
+}
+
+
+gfc_expr *
 gfc_simplify_tan (gfc_expr *x)
 {
-  int i;
   gfc_expr *result;
 
   if (x->expr_type != EXPR_CONSTANT)
     return NULL;
 
-  i = gfc_validate_kind (BT_REAL, x->ts.kind, false);
-
   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
 
-  mpfr_tan (result->value.real, x->value.real, GFC_RND_MODE);
+  if (x->ts.type == BT_REAL)
+    mpfr_tan (result->value.real, x->value.real, GFC_RND_MODE);
+  else if (x->ts.type == BT_COMPLEX)
+    mpc_tan (result->value.complex, x->value.complex, GFC_MPC_RND_MODE);
+  else
+    gcc_unreachable ();
 
   return range_check (result, "TAN");
 }
@@ -4339,7 +5255,12 @@ gfc_simplify_tanh (gfc_expr *x)
 
   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
 
-  mpfr_tanh (result->value.real, x->value.real, GFC_RND_MODE);
+  if (x->ts.type == BT_REAL)
+    mpfr_tanh (result->value.real, x->value.real, GFC_RND_MODE);
+  else if (x->ts.type == BT_COMPLEX)
+    mpc_tanh (result->value.complex, x->value.complex, GFC_MPC_RND_MODE);
+  else
+    gcc_unreachable ();
 
   return range_check (result, "TANH");
 
@@ -4467,6 +5388,7 @@ gfc_simplify_transfer (gfc_expr *source, gfc_expr *mold, gfc_expr *size)
   /* Allocate the buffer to store the binary version of the source.  */
   buffer_size = MAX (source_size, result_size);
   buffer = (unsigned char*)alloca (buffer_size);
+  memset (buffer, 0, buffer_size);
 
   /* Now write source to the buffer.  */
   gfc_target_encode_expr (source, buffer, buffer_size);
@@ -4479,6 +5401,47 @@ gfc_simplify_transfer (gfc_expr *source, gfc_expr *mold, gfc_expr *size)
 
 
 gfc_expr *
+gfc_simplify_transpose (gfc_expr *matrix)
+{
+  int i, matrix_rows;
+  gfc_expr *result;
+  gfc_constructor *matrix_ctor;
+
+  if (!is_constant_array_expr (matrix))
+    return NULL;
+
+  gcc_assert (matrix->rank == 2);
+
+  result = gfc_start_constructor (matrix->ts.type, matrix->ts.kind, &matrix->where);
+  result->rank = 2;
+  result->shape = gfc_get_shape (result->rank);
+  mpz_set (result->shape[0], matrix->shape[1]);
+  mpz_set (result->shape[1], matrix->shape[0]);
+
+  if (matrix->ts.type == BT_CHARACTER)
+    result->ts.u.cl = matrix->ts.u.cl;
+
+  matrix_rows = mpz_get_si (matrix->shape[0]);
+  matrix_ctor = matrix->value.constructor;
+  for (i = 0; i < matrix_rows; ++i)
+    {
+      gfc_constructor *column_ctor = matrix_ctor;
+      while (column_ctor)
+       {
+         gfc_append_constructor (result, 
+                                 gfc_copy_expr (column_ctor->expr));
+
+         ADVANCE (column_ctor, matrix_rows);
+       }
+
+      ADVANCE (matrix_ctor, 1);
+    }
+
+  return result;
+}
+
+
+gfc_expr *
 gfc_simplify_trim (gfc_expr *e)
 {
   gfc_expr *result;
@@ -4521,6 +5484,54 @@ gfc_simplify_ubound (gfc_expr *array, gfc_expr *dim, gfc_expr *kind)
 
 
 gfc_expr *
+gfc_simplify_unpack (gfc_expr *vector, gfc_expr *mask, gfc_expr *field)
+{
+  gfc_expr *result, *e;
+  gfc_constructor *vector_ctor, *mask_ctor, *field_ctor;
+
+  if (!is_constant_array_expr (vector)
+      || !is_constant_array_expr (mask)
+      || (!gfc_is_constant_expr (field)
+         && !is_constant_array_expr(field)))
+    return NULL;
+
+  result = gfc_start_constructor (vector->ts.type,
+                                 vector->ts.kind,
+                                 &vector->where);
+  result->rank = mask->rank;
+  result->shape = gfc_copy_shape (mask->shape, mask->rank);
+
+  if (vector->ts.type == BT_CHARACTER)
+    result->ts.u.cl = vector->ts.u.cl;
+
+  vector_ctor = vector->value.constructor;
+  mask_ctor = mask->value.constructor;
+  field_ctor = field->expr_type == EXPR_ARRAY ? field->value.constructor : NULL;
+
+  while (mask_ctor)
+    {
+      if (mask_ctor->expr->value.logical)
+       {
+         gcc_assert (vector_ctor);
+         e = gfc_copy_expr (vector_ctor->expr);
+         ADVANCE (vector_ctor, 1);
+       }
+      else if (field->expr_type == EXPR_ARRAY)
+       e = gfc_copy_expr (field_ctor->expr);
+      else
+       e = gfc_copy_expr (field);
+
+      gfc_append_constructor (result, e);
+
+      ADVANCE (mask_ctor, 1);
+      ADVANCE (field_ctor, 1);
+    }
+
+  return result;
+}
+
+
+gfc_expr *
 gfc_simplify_verify (gfc_expr *s, gfc_expr *set, gfc_expr *b, gfc_expr *kind)
 {
   gfc_expr *result;
@@ -4872,7 +5883,7 @@ gfc_convert_char_constant (gfc_expr *e, bt type ATTRIBUTE_UNUSED, int kind)
       result->shape = gfc_copy_shape (e->shape, e->rank);
       result->where = e->where;
       result->rank = e->rank;
-      result->ts.cl = e->ts.cl;
+      result->ts.u.cl = e->ts.u.cl;
 
       return result;
     }