gcc/fortran/simplify.c

/* Simplify intrinsic functions at compile-time.
   Copyright (C) 2000-2018 Free Software Foundation, Inc.
   Contributed by Andy Vaught & Katherine Holcomb

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 3, or (at your option) any later
version.

GCC is distributed in the hope that it will be useful, but WITHOUT ANY
WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
for more details.

You should have received a copy of the GNU General Public License
along with GCC; see the file COPYING3.  If not see
<http://www.gnu.org/licenses/>.  */

#include "config.h"
#include "system.h"
#include "coretypes.h"
#include "tm.h"		/* For BITS_PER_UNIT.  */
#include "gfortran.h"
#include "arith.h"
#include "intrinsic.h"
#include "match.h"
#include "target-memory.h"
#include "constructor.h"
#include "version.h"	/* For version_string.  */


gfc_expr gfc_bad_expr;

static gfc_expr *simplify_size (gfc_expr *, gfc_expr *, int);


/* Note that 'simplification' is not just transforming expressions.
   For functions that are not simplified at compile time, range
   checking is done if possible.

   The return convention is that each simplification function returns:

     A new expression node corresponding to the simplified arguments.
     The original arguments are destroyed by the caller, and must not
     be a part of the new expression.

     NULL pointer indicating that no simplification was possible and
     the original expression should remain intact.

     An expression pointer to gfc_bad_expr (a static placeholder)
     indicating that some error has prevented simplification.  The
     error is generated within the function and should be propagated
     upwards

   By the time a simplification function gets control, it has been
   decided that the function call is really supposed to be the
   intrinsic.  No type checking is strictly necessary, since only
   valid types will be passed on.  On the other hand, a simplification
   subroutine may have to look at the type of an argument as part of
   its processing.

   Array arguments are only passed to these subroutines that implement
   the simplification of transformational intrinsics.

   The functions in this file don't have much comment with them, but
   everything is reasonably straight-forward.  The Standard, chapter 13
   is the best comment you'll find for this file anyway.  */

/* Range checks an expression node.  If all goes well, returns the
   node, otherwise returns &gfc_bad_expr and frees the node.  */

static gfc_expr *
range_check (gfc_expr *result, const char *name)
{
  if (result == NULL)
    return &gfc_bad_expr;

  if (result->expr_type != EXPR_CONSTANT)
    return result;

  switch (gfc_range_check (result))
    {
      case ARITH_OK:
	return result;

      case ARITH_OVERFLOW:
	gfc_error ("Result of %s overflows its kind at %L", name,
		   &result->where);
	break;

      case ARITH_UNDERFLOW:
	gfc_error ("Result of %s underflows its kind at %L", name,
		   &result->where);
	break;

      case ARITH_NAN:
	gfc_error ("Result of %s is NaN at %L", name, &result->where);
	break;

      default:
	gfc_error ("Result of %s gives range error for its kind at %L", name,
		   &result->where);
	break;
    }

  gfc_free_expr (result);
  return &gfc_bad_expr;
}


/* A helper function that gets an optional and possibly missing
   kind parameter.  Returns the kind, -1 if something went wrong.  */

static int
get_kind (bt type, gfc_expr *k, const char *name, int default_kind)
{
  int kind;

  if (k == NULL)
    return default_kind;

  if (k->expr_type != EXPR_CONSTANT)
    {
      gfc_error ("KIND parameter of %s at %L must be an initialization "
		 "expression", name, &k->where);
      return -1;
    }

  if (gfc_extract_int (k, &kind)
      || gfc_validate_kind (type, kind, true) < 0)
    {
      gfc_error ("Invalid KIND parameter of %s at %L", name, &k->where);
      return -1;
    }

  return kind;
}


/* Converts an mpz_t signed variable into an unsigned one, assuming
   two's complement representations and a binary width of bitsize.
   The conversion is a no-op unless x is negative; otherwise, it can
   be accomplished by masking out the high bits.  */

static void
convert_mpz_to_unsigned (mpz_t x, int bitsize)
{
  mpz_t mask;

  if (mpz_sgn (x) < 0)
    {
      /* Confirm that no bits above the signed range are unset if we
	 are doing range checking.  */
      if (flag_range_check != 0)
	gcc_assert (mpz_scan0 (x, bitsize-1) == ULONG_MAX);

      mpz_init_set_ui (mask, 1);
      mpz_mul_2exp (mask, mask, bitsize);
      mpz_sub_ui (mask, mask, 1);

      mpz_and (x, x, mask);

      mpz_clear (mask);
    }
  else
    {
      /* Confirm that no bits above the signed range are set.  */
      gcc_assert (mpz_scan1 (x, bitsize-1) == ULONG_MAX);
    }
}


/* Converts an mpz_t unsigned variable into a signed one, assuming
   two's complement representations and a binary width of bitsize.
   If the bitsize-1 bit is set, this is taken as a sign bit and
   the number is converted to the corresponding negative number.  */

void
gfc_convert_mpz_to_signed (mpz_t x, int bitsize)
{
  mpz_t mask;

  /* Confirm that no bits above the unsigned range are set if we are
     doing range checking.  */
  if (flag_range_check != 0)
    gcc_assert (mpz_scan1 (x, bitsize) == ULONG_MAX);

  if (mpz_tstbit (x, bitsize - 1) == 1)
    {
      mpz_init_set_ui (mask, 1);
      mpz_mul_2exp (mask, mask, bitsize);
      mpz_sub_ui (mask, mask, 1);

      /* We negate the number by hand, zeroing the high bits, that is
	 make it the corresponding positive number, and then have it
	 negated by GMP, giving the correct representation of the
	 negative number.  */
      mpz_com (x, x);
      mpz_add_ui (x, x, 1);
      mpz_and (x, x, mask);

      mpz_neg (x, x);

      mpz_clear (mask);
    }
}


/* In-place convert BOZ to REAL of the specified kind.  */

static gfc_expr *
convert_boz (gfc_expr *x, int kind)
{
  if (x && x->ts.type == BT_INTEGER && x->is_boz)
    {
      gfc_typespec ts;
      gfc_clear_ts (&ts);
      ts.type = BT_REAL;
      ts.kind = kind;

      if (!gfc_convert_boz (x, &ts))
	return &gfc_bad_expr;
    }

  return x;
}


/* 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 = gfc_constructor_first (e->value.constructor);
       c; c = gfc_constructor_next (c))
    if (c->expr->expr_type != EXPR_CONSTANT
	  && c->expr->expr_type != EXPR_STRUCTURE)
      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 = gfc_constructor_first (e->value.constructor);
      while (ctor)
	{
	  init_result_expr (ctor->expr, init, array);
	  ctor = gfc_constructor_next (ctor);
	}
    }
  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;
   if conj_a is true, the matrix_a is complex conjugated.  */

static gfc_expr *
compute_dot_product (gfc_expr *matrix_a, int stride_a, int offset_a,
		     gfc_expr *matrix_b, int stride_b, int offset_b,
		     bool conj_a)
{
  gfc_expr *result, *a, *b, *c;

  result = gfc_get_constant_expr (matrix_a->ts.type, matrix_a->ts.kind,
				  &matrix_a->where);
  init_result_expr (result, 0, NULL);

  a = gfc_constructor_lookup_expr (matrix_a->value.constructor, offset_a);
  b = gfc_constructor_lookup_expr (matrix_b->value.constructor, offset_b);
  while (a && 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 (a),
				      gfc_copy_expr (b)));
	    break;

	  case BT_INTEGER:
	  case BT_REAL:
	  case BT_COMPLEX:
	    if (conj_a && a->ts.type == BT_COMPLEX)
	      c = gfc_simplify_conjg (a);
	    else
	      c = gfc_copy_expr (a);
	    result = gfc_add (result, gfc_multiply (c, gfc_copy_expr (b)));
	    break;

	  default:
	    gcc_unreachable();
	}

      offset_a += stride_a;
      a = gfc_constructor_lookup_expr (matrix_a->value.constructor, offset_a);

      offset_b += stride_b;
      b = gfc_constructor_lookup_expr (matrix_b->value.constructor, offset_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_get_constant_expr (type, kind, where);

  result = gfc_get_array_expr (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_constructor_append_expr (&result->value.constructor,
				   gfc_get_constant_expr (type, kind, where),
				   NULL);
    }

  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 implementation 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 = gfc_constructor_first (array->value.constructor);
  mask_ctor = NULL;
  if (mask && mask->expr_type == EXPR_ARRAY)
    mask_ctor = gfc_constructor_first (mask->value.constructor);

  while (array_ctor)
    {
      a = array_ctor->expr;
      array_ctor = gfc_constructor_next (array_ctor);

      /* A constant MASK equals .TRUE. here and can be ignored.  */
      if (mask_ctor)
	{
	  m = mask_ctor->expr;
	  mask_ctor = gfc_constructor_next (mask_ctor);
	  if (!m->value.logical)
	    continue;
	}

      result = op (result, gfc_copy_expr (a));
      if (!result)
	return result;
    }

  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().
   The result might be post processed using post_op.  */

static gfc_expr *
simplify_transformation_to_array (gfc_expr *result, gfc_expr *array, gfc_expr *dim,
				  gfc_expr *mask, transformational_op op,
				  transformational_op post_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);
  mpz_clear (size);

  arrayvec = XCNEWVEC (gfc_expr*, arraysize);

  array_ctor = gfc_constructor_first (array->value.constructor);
  mask_ctor = NULL;
  if (mask && mask->expr_type == EXPR_ARRAY)
    mask_ctor = gfc_constructor_first (mask->value.constructor);

  for (i = 0; i < arraysize; ++i)
    {
      arrayvec[i] = array_ctor->expr;
      array_ctor = gfc_constructor_next (array_ctor);

      if (mask_ctor)
	{
	  if (!mask_ctor->expr->value.logical)
	    arrayvec[i] = NULL;

	  mask_ctor = gfc_constructor_next (mask_ctor);
	}
    }

  /* Same for the result expression.  */
  gfc_array_size (result, &size);
  resultsize = mpz_get_ui (size);
  mpz_clear (size);

  resultvec = XCNEWVEC (gfc_expr*, resultsize);
  result_ctor = gfc_constructor_first (result->value.constructor);
  for (i = 0; i < resultsize; ++i)
    {
      resultvec[i] = result_ctor->expr;
      result_ctor = gfc_constructor_next (result_ctor);
    }

  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)
	    {
	      /* If the nested loop is unrolled GFC_MAX_DIMENSIONS
		 times, we'd warn for the last iteration, because the
		 array index will have already been incremented to the
		 array sizes, and we can't tell that this must make
		 the test against result->rank false, because ranks
		 must not exceed GFC_MAX_DIMENSIONS.  */
	      GCC_DIAGNOSTIC_PUSH_IGNORED (-Warray-bounds)
	      count[n]++;
	      base += sstride[n];
	      dest += dstride[n];
	      GCC_DIAGNOSTIC_POP
	    }
	  else
	    done = true;
       }
    }

  /* Place updated expression in result constructor.  */
  result_ctor = gfc_constructor_first (result->value.constructor);
  for (i = 0; i < resultsize; ++i)
    {
      if (post_op)
	result_ctor->expr = post_op (result_ctor->expr, resultvec[i]);
      else
	result_ctor->expr = resultvec[i];
      result_ctor = gfc_constructor_next (result_ctor);
    }

  free (arrayvec);
  free (resultvec);
  return result;
}


static gfc_expr *
simplify_transformation (gfc_expr *array, gfc_expr *dim, gfc_expr *mask,
			 int init_val, transformational_op op)
{
  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, init_val, NULL);

  return !dim || array->rank == 1 ?
    simplify_transformation_to_scalar (result, array, mask, op) :
    simplify_transformation_to_array (result, array, dim, mask, op, NULL);
}


/********************** Simplification functions *****************************/

gfc_expr *
gfc_simplify_abs (gfc_expr *e)
{
  gfc_expr *result;

  if (e->expr_type != EXPR_CONSTANT)
    return NULL;

  switch (e->ts.type)
    {
      case BT_INTEGER:
	result = gfc_get_constant_expr (BT_INTEGER, e->ts.kind, &e->where);
	mpz_abs (result->value.integer, e->value.integer);
	return range_check (result, "IABS");

      case BT_REAL:
	result = gfc_get_constant_expr (BT_REAL, e->ts.kind, &e->where);
	mpfr_abs (result->value.real, e->value.real, GFC_RND_MODE);
	return range_check (result, "ABS");

      case BT_COMPLEX:
	gfc_set_model_kind (e->ts.kind);
	result = gfc_get_constant_expr (BT_REAL, e->ts.kind, &e->where);
	mpc_abs (result->value.real, e->value.complex, GFC_RND_MODE);
	return range_check (result, "CABS");

      default:
	gfc_internal_error ("gfc_simplify_abs(): Bad type");
    }
}


static gfc_expr *
simplify_achar_char (gfc_expr *e, gfc_expr *k, const char *name, bool ascii)
{
  gfc_expr *result;
  int kind;
  bool too_large = false;

  if (e->expr_type != EXPR_CONSTANT)
    return NULL;

  kind = get_kind (BT_CHARACTER, k, name, gfc_default_character_kind);
  if (kind == -1)
    return &gfc_bad_expr;

  if (mpz_cmp_si (e->value.integer, 0) < 0)
    {
      gfc_error ("Argument of %s function at %L is negative", name,
		 &e->where);
      return &gfc_bad_expr;
    }

  if (ascii && warn_surprising && mpz_cmp_si (e->value.integer, 127) > 0)
    gfc_warning (OPT_Wsurprising,
		 "Argument of %s function at %L outside of range [0,127]",
		 name, &e->where);

  if (kind == 1 && mpz_cmp_si (e->value.integer, 255) > 0)
    too_large = true;
  else if (kind == 4)
    {
      mpz_t t;
      mpz_init_set_ui (t, 2);
      mpz_pow_ui (t, t, 32);
      mpz_sub_ui (t, t, 1);
      if (mpz_cmp (e->value.integer, t) > 0)
	too_large = true;
      mpz_clear (t);
    }

  if (too_large)
    {
      gfc_error ("Argument of %s function at %L is too large for the "
		 "collating sequence of kind %d", name, &e->where, kind);
      return &gfc_bad_expr;
    }

  result = gfc_get_character_expr (kind, &e->where, NULL, 1);
  result->value.character.string[0] = mpz_get_ui (e->value.integer);

  return result;
}



/* We use the processor's collating sequence, because all
   systems that gfortran currently works on are ASCII.  */

gfc_expr *
gfc_simplify_achar (gfc_expr *e, gfc_expr *k)
{
  return simplify_achar_char (e, k, "ACHAR", true);
}


gfc_expr *
gfc_simplify_acos (gfc_expr *x)
{
  gfc_expr *result;

  if (x->expr_type != EXPR_CONSTANT)
    return NULL;

  switch (x->ts.type)
    {
      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_get_constant_expr (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_get_constant_expr (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");
    }

  return range_check (result, "ACOS");
}

gfc_expr *
gfc_simplify_acosh (gfc_expr *x)
{
  gfc_expr *result;

  if (x->expr_type != EXPR_CONSTANT)
    return NULL;

  switch (x->ts.type)
    {
      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;
	  }

	result = gfc_get_constant_expr (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_get_constant_expr (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");
}

gfc_expr *
gfc_simplify_adjustl (gfc_expr *e)
{
  gfc_expr *result;
  int count, i, len;
  gfc_char_t ch;

  if (e->expr_type != EXPR_CONSTANT)
    return NULL;

  len = e->value.character.length;

  for (count = 0, i = 0; i < len; ++i)
    {
      ch = e->value.character.string[i];
      if (ch != ' ')
	break;
      ++count;
    }

  result = gfc_get_character_expr (e->ts.kind, &e->where, NULL, len);
  for (i = 0; i < len - count; ++i)
    result->value.character.string[i] = e->value.character.string[count + i];

  return result;
}


gfc_expr *
gfc_simplify_adjustr (gfc_expr *e)
{
  gfc_expr *result;
  int count, i, len;
  gfc_char_t ch;

  if (e->expr_type != EXPR_CONSTANT)
    return NULL;

  len = e->value.character.length;

  for (count = 0, i = len - 1; i >= 0; --i)
    {
      ch = e->value.character.string[i];
      if (ch != ' ')
	break;
      ++count;
    }

  result = gfc_get_character_expr (e->ts.kind, &e->where, NULL, len);
  for (i = 0; i < count; ++i)
    result->value.character.string[i] = ' ';

  for (i = count; i < len; ++i)
    result->value.character.string[i] = e->value.character.string[i - count];

  return result;
}


gfc_expr *
gfc_simplify_aimag (gfc_expr *e)
{
  gfc_expr *result;

  if (e->expr_type != EXPR_CONSTANT)
    return NULL;

  result = gfc_get_constant_expr (BT_REAL, e->ts.kind, &e->where);
  mpfr_set (result->value.real, mpc_imagref (e->value.complex), GFC_RND_MODE);

  return range_check (result, "AIMAG");
}


gfc_expr *
gfc_simplify_aint (gfc_expr *e, gfc_expr *k)
{
  gfc_expr *rtrunc, *result;
  int kind;

  kind = get_kind (BT_REAL, k, "AINT", e->ts.kind);
  if (kind == -1)
    return &gfc_bad_expr;

  if (e->expr_type != EXPR_CONSTANT)
    return NULL;

  rtrunc = gfc_copy_expr (e);
  mpfr_trunc (rtrunc->value.real, e->value.real);

  result = gfc_real2real (rtrunc, kind);

  gfc_free_expr (rtrunc);

  return range_check (result, "AINT");
}


gfc_expr *
gfc_simplify_all (gfc_expr *mask, gfc_expr *dim)
{
  return simplify_transformation (mask, dim, NULL, true, gfc_and);
}


gfc_expr *
gfc_simplify_dint (gfc_expr *e)
{
  gfc_expr *rtrunc, *result;

  if (e->expr_type != EXPR_CONSTANT)
    return NULL;

  rtrunc = gfc_copy_expr (e);
  mpfr_trunc (rtrunc->value.real, e->value.real);

  result = gfc_real2real (rtrunc, gfc_default_double_kind);

  gfc_free_expr (rtrunc);

  return range_check (result, "DINT");
}


gfc_expr *
gfc_simplify_dreal (gfc_expr *e)
{
  gfc_expr *result = NULL;

  if (e->expr_type != EXPR_CONSTANT)
    return NULL;

  result = gfc_get_constant_expr (BT_REAL, e->ts.kind, &e->where);
  mpc_real (result->value.real, e->value.complex, GFC_RND_MODE);

  return range_check (result, "DREAL");
}


gfc_expr *
gfc_simplify_anint (gfc_expr *e, gfc_expr *k)
{
  gfc_expr *result;
  int kind;

  kind = get_kind (BT_REAL, k, "ANINT", e->ts.kind);
  if (kind == -1)
    return &gfc_bad_expr;

  if (e->expr_type != EXPR_CONSTANT)
    return NULL;

  result = gfc_get_constant_expr (e->ts.type, kind, &e->where);
  mpfr_round (result->value.real, e->value.real);

  return range_check (result, "ANINT");
}


gfc_expr *
gfc_simplify_and (gfc_expr *x, gfc_expr *y)
{
  gfc_expr *result;
  int kind;

  if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
    return NULL;

  kind = x->ts.kind > y->ts.kind ? x->ts.kind : y->ts.kind;

  switch (x->ts.type)
    {
      case BT_INTEGER:
	result = gfc_get_constant_expr (BT_INTEGER, kind, &x->where);
	mpz_and (result->value.integer, x->value.integer, y->value.integer);
	return range_check (result, "AND");

      case BT_LOGICAL:
	return gfc_get_logical_expr (kind, &x->where,
				     x->value.logical && y->value.logical);

      default:
	gcc_unreachable ();
    }
}


gfc_expr *
gfc_simplify_any (gfc_expr *mask, gfc_expr *dim)
{
  return simplify_transformation (mask, dim, NULL, false, gfc_or);
}


gfc_expr *
gfc_simplify_dnint (gfc_expr *e)
{
  gfc_expr *result;

  if (e->expr_type != EXPR_CONSTANT)
    return NULL;

  result = gfc_get_constant_expr (BT_REAL, gfc_default_double_kind, &e->where);
  mpfr_round (result->value.real, e->value.real);

  return range_check (result, "DNINT");
}


gfc_expr *
gfc_simplify_asin (gfc_expr *x)
{
  gfc_expr *result;

  if (x->expr_type != EXPR_CONSTANT)
    return NULL;

  switch (x->ts.type)
    {
      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_get_constant_expr (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_get_constant_expr (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");
    }

  return range_check (result, "ASIN");
}


gfc_expr *
gfc_simplify_asinh (gfc_expr *x)
{
  gfc_expr *result;

  if (x->expr_type != EXPR_CONSTANT)
    return NULL;

  result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);

  switch (x->ts.type)
    {
      case BT_REAL:
	mpfr_asinh (result->value.real, x->value.real, GFC_RND_MODE);
	break;

      case BT_COMPLEX:
	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");
}


gfc_expr *
gfc_simplify_atan (gfc_expr *x)
{
  gfc_expr *result;

  if (x->expr_type != EXPR_CONSTANT)
    return NULL;

  result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);

  switch (x->ts.type)
    {
      case BT_REAL:
	mpfr_atan (result->value.real, x->value.real, GFC_RND_MODE);
	break;

      case BT_COMPLEX:
	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");
}


gfc_expr *
gfc_simplify_atanh (gfc_expr *x)
{
  gfc_expr *result;

  if (x->expr_type != EXPR_CONSTANT)
    return NULL;

  switch (x->ts.type)
    {
      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;
	  }
	result = gfc_get_constant_expr (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_get_constant_expr (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");
}


gfc_expr *
gfc_simplify_atan2 (gfc_expr *y, gfc_expr *x)
{
  gfc_expr *result;

  if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
    return NULL;

  if (mpfr_zero_p (y->value.real) && mpfr_zero_p (x->value.real))
    {
      gfc_error ("If first argument of ATAN2 %L is zero, then the "
		 "second argument must not be zero", &x->where);
      return &gfc_bad_expr;
    }

  result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
  mpfr_atan2 (result->value.real, y->value.real, x->value.real, GFC_RND_MODE);

  return range_check (result, "ATAN2");
}


gfc_expr *
gfc_simplify_bessel_j0 (gfc_expr *x)
{
  gfc_expr *result;

  if (x->expr_type != EXPR_CONSTANT)
    return NULL;

  result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
  mpfr_j0 (result->value.real, x->value.real, GFC_RND_MODE);

  return range_check (result, "BESSEL_J0");
}


gfc_expr *
gfc_simplify_bessel_j1 (gfc_expr *x)
{
  gfc_expr *result;

  if (x->expr_type != EXPR_CONSTANT)
    return NULL;

  result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
  mpfr_j1 (result->value.real, x->value.real, GFC_RND_MODE);

  return range_check (result, "BESSEL_J1");
}


gfc_expr *
gfc_simplify_bessel_jn (gfc_expr *order, gfc_expr *x)
{
  gfc_expr *result;
  long n;

  if (x->expr_type != EXPR_CONSTANT || order->expr_type != EXPR_CONSTANT)
    return NULL;

  n = mpz_get_si (order->value.integer);
  result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
  mpfr_jn (result->value.real, n, x->value.real, GFC_RND_MODE);

  return range_check (result, "BESSEL_JN");
}


/* Simplify transformational form of JN and YN.  */

static gfc_expr *
gfc_simplify_bessel_n2 (gfc_expr *order1, gfc_expr *order2, gfc_expr *x,
			bool jn)
{
  gfc_expr *result;
  gfc_expr *e;
  long n1, n2;
  int i;
  mpfr_t x2rev, last1, last2;

  if (x->expr_type != EXPR_CONSTANT || order1->expr_type != EXPR_CONSTANT
      || order2->expr_type != EXPR_CONSTANT)
    return NULL;

  n1 = mpz_get_si (order1->value.integer);
  n2 = mpz_get_si (order2->value.integer);
  result = gfc_get_array_expr (x->ts.type, x->ts.kind, &x->where);
  result->rank = 1;
  result->shape = gfc_get_shape (1);
  mpz_init_set_ui (result->shape[0], MAX (n2-n1+1, 0));

  if (n2 < n1)
    return result;

  /* Special case: x == 0; it is J0(0.0) == 1, JN(N > 0, 0.0) == 0; and
     YN(N, 0.0) = -Inf.  */

  if (mpfr_cmp_ui (x->value.real, 0.0) == 0)
    {
      if (!jn && flag_range_check)
	{
	  gfc_error ("Result of BESSEL_YN is -INF at %L", &result->where);
 	  gfc_free_expr (result);
	  return &gfc_bad_expr;
	}

      if (jn && n1 == 0)
	{
	  e = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
	  mpfr_set_ui (e->value.real, 1, GFC_RND_MODE);
	  gfc_constructor_append_expr (&result->value.constructor, e,
				       &x->where);
	  n1++;
	}

      for (i = n1; i <= n2; i++)
	{
	  e = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
	  if (jn)
	    mpfr_set_ui (e->value.real, 0, GFC_RND_MODE);
	  else
	    mpfr_set_inf (e->value.real, -1);
	  gfc_constructor_append_expr (&result->value.constructor, e,
				       &x->where);
	}

      return result;
    }

  /* Use the faster but more verbose recurrence algorithm. Bessel functions
     are stable for downward recursion and Neumann functions are stable
     for upward recursion. It is
       x2rev = 2.0/x,
       J(N-1, x) = x2rev * N * J(N, x) - J(N+1, x),
       Y(N+1, x) = x2rev * N * Y(N, x) - Y(N-1, x).
     Cf. http://dlmf.nist.gov/10.74#iv and http://dlmf.nist.gov/10.6#E1  */

  gfc_set_model_kind (x->ts.kind);

  /* Get first recursion anchor.  */

  mpfr_init (last1);
  if (jn)
    mpfr_jn (last1, n2, x->value.real, GFC_RND_MODE);
  else
    mpfr_yn (last1, n1, x->value.real, GFC_RND_MODE);

  e = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
  mpfr_set (e->value.real, last1, GFC_RND_MODE);
  if (range_check (e, jn ? "BESSEL_JN" : "BESSEL_YN") == &gfc_bad_expr)
    {
      mpfr_clear (last1);
      gfc_free_expr (e);
      gfc_free_expr (result);
      return &gfc_bad_expr;
    }
  gfc_constructor_append_expr (&result->value.constructor, e, &x->where);

  if (n1 == n2)
    {
      mpfr_clear (last1);
      return result;
    }

  /* Get second recursion anchor.  */

  mpfr_init (last2);
  if (jn)
    mpfr_jn (last2, n2-1, x->value.real, GFC_RND_MODE);
  else
    mpfr_yn (last2, n1+1, x->value.real, GFC_RND_MODE);

  e = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
  mpfr_set (e->value.real, last2, GFC_RND_MODE);
  if (range_check (e, jn ? "BESSEL_JN" : "BESSEL_YN") == &gfc_bad_expr)
    {
      mpfr_clear (last1);
      mpfr_clear (last2);
      gfc_free_expr (e);
      gfc_free_expr (result);
      return &gfc_bad_expr;
    }
  if (jn)
    gfc_constructor_insert_expr (&result->value.constructor, e, &x->where, -2);
  else
    gfc_constructor_append_expr (&result->value.constructor, e, &x->where);

  if (n1 + 1 == n2)
    {
      mpfr_clear (last1);
      mpfr_clear (last2);
      return result;
    }

  /* Start actual recursion.  */

  mpfr_init (x2rev);
  mpfr_ui_div (x2rev, 2, x->value.real, GFC_RND_MODE);

  for (i = 2; i <= n2-n1; i++)
    {
      e = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);

      /* Special case: For YN, if the previous N gave -INF, set
	 also N+1 to -INF.  */
      if (!jn && !flag_range_check && mpfr_inf_p (last2))
	{
	  mpfr_set_inf (e->value.real, -1);
	  gfc_constructor_append_expr (&result->value.constructor, e,
				       &x->where);
	  continue;
	}

      mpfr_mul_si (e->value.real, x2rev, jn ? (n2-i+1) : (n1+i-1),
		   GFC_RND_MODE);
      mpfr_mul (e->value.real, e->value.real, last2, GFC_RND_MODE);
      mpfr_sub (e->value.real, e->value.real, last1, GFC_RND_MODE);

      if (range_check (e, jn ? "BESSEL_JN" : "BESSEL_YN") == &gfc_bad_expr)
	{
	  /* Range_check frees "e" in that case.  */
	  e = NULL;
	  goto error;
	}

      if (jn)
	gfc_constructor_insert_expr (&result->value.constructor, e, &x->where,
				     -i-1);
      else
	gfc_constructor_append_expr (&result->value.constructor, e, &x->where);

      mpfr_set (last1, last2, GFC_RND_MODE);
      mpfr_set (last2, e->value.real, GFC_RND_MODE);
    }

  mpfr_clear (last1);
  mpfr_clear (last2);
  mpfr_clear (x2rev);
  return result;

error:
  mpfr_clear (last1);
  mpfr_clear (last2);
  mpfr_clear (x2rev);
  gfc_free_expr (e);
  gfc_free_expr (result);
  return &gfc_bad_expr;
}


gfc_expr *
gfc_simplify_bessel_jn2 (gfc_expr *order1, gfc_expr *order2, gfc_expr *x)
{
  return gfc_simplify_bessel_n2 (order1, order2, x, true);
}


gfc_expr *
gfc_simplify_bessel_y0 (gfc_expr *x)
{
  gfc_expr *result;

  if (x->expr_type != EXPR_CONSTANT)
    return NULL;

  result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
  mpfr_y0 (result->value.real, x->value.real, GFC_RND_MODE);

  return range_check (result, "BESSEL_Y0");
}


gfc_expr *
gfc_simplify_bessel_y1 (gfc_expr *x)
{
  gfc_expr *result;

  if (x->expr_type != EXPR_CONSTANT)
    return NULL;

  result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
  mpfr_y1 (result->value.real, x->value.real, GFC_RND_MODE);

  return range_check (result, "BESSEL_Y1");
}


gfc_expr *
gfc_simplify_bessel_yn (gfc_expr *order, gfc_expr *x)
{
  gfc_expr *result;
  long n;

  if (x->expr_type != EXPR_CONSTANT || order->expr_type != EXPR_CONSTANT)
    return NULL;

  n = mpz_get_si (order->value.integer);
  result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
  mpfr_yn (result->value.real, n, x->value.real, GFC_RND_MODE);

  return range_check (result, "BESSEL_YN");
}


gfc_expr *
gfc_simplify_bessel_yn2 (gfc_expr *order1, gfc_expr *order2, gfc_expr *x)
{
  return gfc_simplify_bessel_n2 (order1, order2, x, false);
}


gfc_expr *
gfc_simplify_bit_size (gfc_expr *e)
{
  int i = gfc_validate_kind (e->ts.type, e->ts.kind, false);
  return gfc_get_int_expr (e->ts.kind, &e->where,
			   gfc_integer_kinds[i].bit_size);
}


gfc_expr *
gfc_simplify_btest (gfc_expr *e, gfc_expr *bit)
{
  int b;

  if (e->expr_type != EXPR_CONSTANT || bit->expr_type != EXPR_CONSTANT)
    return NULL;

  if (gfc_extract_int (bit, &b) || b < 0)
    return gfc_get_logical_expr (gfc_default_logical_kind, &e->where, false);

  return gfc_get_logical_expr (gfc_default_logical_kind, &e->where,
			       mpz_tstbit (e->value.integer, b));
}


static int
compare_bitwise (gfc_expr *i, gfc_expr *j)
{
  mpz_t x, y;
  int k, res;

  gcc_assert (i->ts.type == BT_INTEGER);
  gcc_assert (j->ts.type == BT_INTEGER);

  mpz_init_set (x, i->value.integer);
  k = gfc_validate_kind (i->ts.type, i->ts.kind, false);
  convert_mpz_to_unsigned (x, gfc_integer_kinds[k].bit_size);

  mpz_init_set (y, j->value.integer);
  k = gfc_validate_kind (j->ts.type, j->ts.kind, false);
  convert_mpz_to_unsigned (y, gfc_integer_kinds[k].bit_size);

  res = mpz_cmp (x, y);
  mpz_clear (x);
  mpz_clear (y);
  return res;
}


gfc_expr *
gfc_simplify_bge (gfc_expr *i, gfc_expr *j)
{
  if (i->expr_type != EXPR_CONSTANT || j->expr_type != EXPR_CONSTANT)
    return NULL;

  return gfc_get_logical_expr (gfc_default_logical_kind, &i->where,
			       compare_bitwise (i, j) >= 0);
}


gfc_expr *
gfc_simplify_bgt (gfc_expr *i, gfc_expr *j)
{
  if (i->expr_type != EXPR_CONSTANT || j->expr_type != EXPR_CONSTANT)
    return NULL;

  return gfc_get_logical_expr (gfc_default_logical_kind, &i->where,
			       compare_bitwise (i, j) > 0);
}


gfc_expr *
gfc_simplify_ble (gfc_expr *i, gfc_expr *j)
{
  if (i->expr_type != EXPR_CONSTANT || j->expr_type != EXPR_CONSTANT)
    return NULL;

  return gfc_get_logical_expr (gfc_default_logical_kind, &i->where,
			       compare_bitwise (i, j) <= 0);
}


gfc_expr *
gfc_simplify_blt (gfc_expr *i, gfc_expr *j)
{
  if (i->expr_type != EXPR_CONSTANT || j->expr_type != EXPR_CONSTANT)
    return NULL;

  return gfc_get_logical_expr (gfc_default_logical_kind, &i->where,
			       compare_bitwise (i, j) < 0);
}


gfc_expr *
gfc_simplify_ceiling (gfc_expr *e, gfc_expr *k)
{
  gfc_expr *ceil, *result;
  int kind;

  kind = get_kind (BT_INTEGER, k, "CEILING", gfc_default_integer_kind);
  if (kind == -1)
    return &gfc_bad_expr;

  if (e->expr_type != EXPR_CONSTANT)
    return NULL;

  ceil = gfc_copy_expr (e);
  mpfr_ceil (ceil->value.real, e->value.real);

  result = gfc_get_constant_expr (BT_INTEGER, kind, &e->where);
  gfc_mpfr_to_mpz (result->value.integer, ceil->value.real, &e->where);

  gfc_free_expr (ceil);

  return range_check (result, "CEILING");
}


gfc_expr *
gfc_simplify_char (gfc_expr *e, gfc_expr *k)
{
  return simplify_achar_char (e, k, "CHAR", false);
}


/* Common subroutine for simplifying CMPLX, COMPLEX and DCMPLX.  */

static gfc_expr *
simplify_cmplx (const char *name, gfc_expr *x, gfc_expr *y, int kind)
{
  gfc_expr *result;

  if (convert_boz (x, kind) == &gfc_bad_expr)
    return &gfc_bad_expr;

  if (convert_boz (y, kind) == &gfc_bad_expr)
    return &gfc_bad_expr;

  if (x->expr_type != EXPR_CONSTANT
      || (y != NULL && y->expr_type != EXPR_CONSTANT))
    return NULL;

  result = gfc_get_constant_expr (BT_COMPLEX, kind, &x->where);

  switch (x->ts.type)
    {
      case BT_INTEGER:
	mpc_set_z (result->value.complex, x->value.integer, GFC_MPC_RND_MODE);
	break;

      case BT_REAL:
	mpc_set_fr (result->value.complex, x->value.real, GFC_RND_MODE);
	break;

      case BT_COMPLEX:
	mpc_set (result->value.complex, x->value.complex, GFC_MPC_RND_MODE);
	break;

      default:
	gfc_internal_error ("gfc_simplify_dcmplx(): Bad type (x)");
    }

  if (!y)
    return range_check (result, name);

  switch (y->ts.type)
    {
      case BT_INTEGER:
	mpfr_set_z (mpc_imagref (result->value.complex),
		    y->value.integer, GFC_RND_MODE);
	break;

      case BT_REAL:
	mpfr_set (mpc_imagref (result->value.complex),
		  y->value.real, GFC_RND_MODE);
	break;

      default:
	gfc_internal_error ("gfc_simplify_dcmplx(): Bad type (y)");
    }

  return range_check (result, name);
}


gfc_expr *
gfc_simplify_cmplx (gfc_expr *x, gfc_expr *y, gfc_expr *k)
{
  int kind;

  kind = get_kind (BT_REAL, k, "CMPLX", gfc_default_complex_kind);
  if (kind == -1)
    return &gfc_bad_expr;

  return simplify_cmplx ("CMPLX", x, y, kind);
}


gfc_expr *
gfc_simplify_complex (gfc_expr *x, gfc_expr *y)
{
  int kind;

  if (x->ts.type == BT_INTEGER && y->ts.type == BT_INTEGER)
    kind = gfc_default_complex_kind;
  else if (x->ts.type == BT_REAL || y->ts.type == BT_INTEGER)
    kind = x->ts.kind;
  else if (x->ts.type == BT_INTEGER || y->ts.type == BT_REAL)
    kind = y->ts.kind;
  else if (x->ts.type == BT_REAL && y->ts.type == BT_REAL)
    kind = (x->ts.kind > y->ts.kind) ? x->ts.kind : y->ts.kind;
  else
    gcc_unreachable ();

  return simplify_cmplx ("COMPLEX", x, y, kind);
}


gfc_expr *
gfc_simplify_conjg (gfc_expr *e)
{
  gfc_expr *result;

  if (e->expr_type != EXPR_CONSTANT)
    return NULL;

  result = gfc_copy_expr (e);
  mpc_conj (result->value.complex, result->value.complex, GFC_MPC_RND_MODE);

  return range_check (result, "CONJG");
}

/* Return the simplification of the constant expression in icall, or NULL
   if the expression is not constant.  */

static gfc_expr *
simplify_trig_call (gfc_expr *icall)
{
  gfc_isym_id func = icall->value.function.isym->id;
  gfc_expr *x = icall->value.function.actual->expr;

  /* The actual simplifiers will return NULL for non-constant x.  */
  switch (func)
    {
    case GFC_ISYM_ACOS:
      return gfc_simplify_acos (x);
    case GFC_ISYM_ASIN:
      return gfc_simplify_asin (x);
    case GFC_ISYM_ATAN:
      return gfc_simplify_atan (x);
    case GFC_ISYM_COS:
      return gfc_simplify_cos (x);
    case GFC_ISYM_COTAN:
      return gfc_simplify_cotan (x);
    case GFC_ISYM_SIN:
      return gfc_simplify_sin (x);
    case GFC_ISYM_TAN:
      return gfc_simplify_tan (x);
    default:
      gfc_internal_error ("in simplify_trig_call(): Bad intrinsic");
    }
}

/* Convert a floating-point number from radians to degrees.  */

static void
degrees_f (mpfr_t x, mp_rnd_t rnd_mode)
{
  mpfr_t tmp;
  mpfr_init (tmp);

  /* Set x = x % 2pi to avoid offsets with large angles.  */
  mpfr_const_pi (tmp, rnd_mode);
  mpfr_mul_ui (tmp, tmp, 2, rnd_mode);
  mpfr_fmod (tmp, x, tmp, rnd_mode);

  /* Set x = x * 180.  */
  mpfr_mul_ui (x, x, 180, rnd_mode);

  /* Set x = x / pi.  */
  mpfr_const_pi (tmp, rnd_mode);
  mpfr_div (x, x, tmp, rnd_mode);

  mpfr_clear (tmp);
}

/* Convert a floating-point number from degrees to radians.  */

static void
radians_f (mpfr_t x, mp_rnd_t rnd_mode)
{
  mpfr_t tmp;
  mpfr_init (tmp);

  /* Set x = x % 360 to avoid offsets with large angles.  */
  mpfr_set_ui (tmp, 360, rnd_mode);
  mpfr_fmod (tmp, x, tmp, rnd_mode);

  /* Set x = x * pi.  */
  mpfr_const_pi (tmp, rnd_mode);
  mpfr_mul (x, x, tmp, rnd_mode);

  /* Set x = x / 180.  */
  mpfr_div_ui (x, x, 180, rnd_mode);

  mpfr_clear (tmp);
}


/* Convert argument to radians before calling a trig function.  */

gfc_expr *
gfc_simplify_trigd (gfc_expr *icall)
{
  gfc_expr *arg;

  arg = icall->value.function.actual->expr;

  if (arg->ts.type != BT_REAL)
    gfc_internal_error ("in gfc_simplify_trigd(): Bad type");

  if (arg->expr_type == EXPR_CONSTANT)
    /* Convert constant to radians before passing off to simplifier.  */
    radians_f (arg->value.real, GFC_RND_MODE);

  /* Let the usual simplifier take over - we just simplified the arg.  */
  return simplify_trig_call (icall);
}

/* Convert result of an inverse trig function to degrees.  */

gfc_expr *
gfc_simplify_atrigd (gfc_expr *icall)
{
  gfc_expr *result;

  if (icall->value.function.actual->expr->ts.type != BT_REAL)
    gfc_internal_error ("in gfc_simplify_atrigd(): Bad type");

  /* See if another simplifier has work to do first.  */
  result = simplify_trig_call (icall);

  if (result && result->expr_type == EXPR_CONSTANT)
    {
      /* Convert constant to degrees after passing off to actual simplifier.  */
      degrees_f (result->value.real, GFC_RND_MODE);
      return result;
    }

  /* Let gfc_resolve_atrigd take care of the non-constant case.  */
  return NULL;
}

/* Convert the result of atan2 to degrees.  */

gfc_expr *
gfc_simplify_atan2d (gfc_expr *y, gfc_expr *x)
{
  gfc_expr *result;

  if (x->ts.type != BT_REAL || y->ts.type != BT_REAL)
    gfc_internal_error ("in gfc_simplify_atan2d(): Bad type");

  if (x->expr_type == EXPR_CONSTANT && y->expr_type == EXPR_CONSTANT)
    {
      result = gfc_simplify_atan2 (y, x);
      if (result != NULL)
	{
	  degrees_f (result->value.real, GFC_RND_MODE);
	  return result;
	}
    }

  /* Let gfc_resolve_atan2d take care of the non-constant case.  */
  return NULL;
}

gfc_expr *
gfc_simplify_cos (gfc_expr *x)
{
  gfc_expr *result;

  if (x->expr_type != EXPR_CONSTANT)
    return NULL;

  result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);

  switch (x->ts.type)
    {
      case BT_REAL:
	mpfr_cos (result->value.real, x->value.real, GFC_RND_MODE);
	break;

      case BT_COMPLEX:
	gfc_set_model_kind (x->ts.kind);
	mpc_cos (result->value.complex, x->value.complex, GFC_MPC_RND_MODE);
	break;

      default:
	gfc_internal_error ("in gfc_simplify_cos(): Bad type");
    }

  return range_check (result, "COS");
}


gfc_expr *
gfc_simplify_cosh (gfc_expr *x)
{
  gfc_expr *result;

  if (x->expr_type != EXPR_CONSTANT)
    return NULL;

  result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);

  switch (x->ts.type)
    {
      case BT_REAL:
	mpfr_cosh (result->value.real, x->value.real, GFC_RND_MODE);
	break;

      case BT_COMPLEX:
	mpc_cosh (result->value.complex, x->value.complex, GFC_MPC_RND_MODE);
	break;

      default:
	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, NULL);
}


gfc_expr *
gfc_simplify_cshift (gfc_expr *array, gfc_expr *shift, gfc_expr *dim)
{
  gfc_expr *a, *result;
  int dm;

  /* DIM is only useful for rank > 1, but deal with it here as one can
     set DIM = 1 for rank = 1.  */
  if (dim)
    {
      if (!gfc_is_constant_expr (dim))
	return NULL;
      dm = mpz_get_si (dim->value.integer);
    }
  else
    dm = 1;

  /* Copy array into 'a', simplify it, and then test for a constant array.  */
  a = gfc_copy_expr (array);
  gfc_simplify_expr (a, 0);
  if (!is_constant_array_expr (a))
    {
      gfc_free_expr (a);
      return NULL;
    }

  if (a->rank == 1)
    {
      gfc_constructor *ca, *cr;
      mpz_t size;
      int i, j, shft, sz;

      if (!gfc_is_constant_expr (shift))
	{
	  gfc_free_expr (a);
	  return NULL;
	}

      shft = mpz_get_si (shift->value.integer);

      /*  Case (i):  If ARRAY has rank one, element i of the result is
	  ARRAY (1 + MODULO (i + SHIFT - 1, SIZE (ARRAY))).  */

      mpz_init (size);
      gfc_array_size (a, &size);
      sz = mpz_get_si (size);
      mpz_clear (size);

      /* Special case: Zero-sized array.  */
      if (sz == 0)
	return a;

      /* Adjust shft to deal with right or left shifts. */
      shft = shft % sz;
      if (shft < 0)
	shft += sz;

      /* Special case: Shift to the original order!  */
      if (shft % sz == 0)
	return a;

      result = gfc_copy_expr (a);
      cr = gfc_constructor_first (result->value.constructor);
      for (i = 0; i < sz; i++, cr = gfc_constructor_next (cr))
	{
	  j = (i + shft) % sz;
	  ca = gfc_constructor_first (a->value.constructor);
	  while (j-- > 0)
	    ca = gfc_constructor_next (ca);
	  cr->expr = gfc_copy_expr (ca->expr);
	}

      gfc_free_expr (a);
      return result;
    }
  else
    {
      /* FIXME: Deal with rank > 1 arrays.  For now, don't leak memory.  */

      /* GCC bootstrap is too stupid to realize that the above code for dm
	 is correct.  First, dim can be specified for a rank 1 array.  It is
	 not needed in this nor used here.  Second, the code is simply waiting
	 for someone to implement rank > 1 simplification.   For now, add a
	 pessimization to the code that has a zero valid reason to be here.  */
      if (dm > array->rank)
	gcc_unreachable ();

      gfc_free_expr (a);
    }

  return NULL;
}


gfc_expr *
gfc_simplify_dcmplx (gfc_expr *x, gfc_expr *y)
{
  return simplify_cmplx ("DCMPLX", x, y, gfc_default_double_kind);
}


gfc_expr *
gfc_simplify_dble (gfc_expr *e)
{
  gfc_expr *result = NULL;

  if (e->expr_type != EXPR_CONSTANT)
    return NULL;

  if (convert_boz (e, gfc_default_double_kind) == &gfc_bad_expr)
    return &gfc_bad_expr;

  result = gfc_convert_constant (e, BT_REAL, gfc_default_double_kind);
  if (result == &gfc_bad_expr)
    return &gfc_bad_expr;

  return range_check (result, "DBLE");
}


gfc_expr *
gfc_simplify_digits (gfc_expr *x)
{
  int i, digits;

  i = gfc_validate_kind (x->ts.type, x->ts.kind, false);

  switch (x->ts.type)
    {
      case BT_INTEGER:
	digits = gfc_integer_kinds[i].digits;
	break;

      case BT_REAL:
      case BT_COMPLEX:
	digits = gfc_real_kinds[i].digits;
	break;

      default:
	gcc_unreachable ();
    }

  return gfc_get_int_expr (gfc_default_integer_kind, NULL, digits);
}


gfc_expr *
gfc_simplify_dim (gfc_expr *x, gfc_expr *y)
{
  gfc_expr *result;
  int kind;

  if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
    return NULL;

  kind = x->ts.kind > y->ts.kind ? x->ts.kind : y->ts.kind;
  result = gfc_get_constant_expr (x->ts.type, kind, &x->where);

  switch (x->ts.type)
    {
      case BT_INTEGER:
	if (mpz_cmp (x->value.integer, y->value.integer) > 0)
	  mpz_sub (result->value.integer, x->value.integer, y->value.integer);
	else
	  mpz_set_ui (result->value.integer, 0);

	break;

      case BT_REAL:
	if (mpfr_cmp (x->value.real, y->value.real) > 0)
	  mpfr_sub (result->value.real, x->value.real, y->value.real,
		    GFC_RND_MODE);
	else
	  mpfr_set_ui (result->value.real, 0, GFC_RND_MODE);

	break;

      default:
	gfc_internal_error ("gfc_simplify_dim(): Bad type");
    }

  return range_check (result, "DIM");
}


gfc_expr*
gfc_simplify_dot_product (gfc_expr *vector_a, gfc_expr *vector_b)
{

  gfc_expr temp;

  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);

  temp.expr_type = EXPR_OP;
  gfc_clear_ts (&temp.ts);
  temp.value.op.op = INTRINSIC_NONE;
  temp.value.op.op1 = vector_a;
  temp.value.op.op2 = vector_b;
  gfc_type_convert_binary (&temp, 1);

  return compute_dot_product (vector_a, 1, 0, vector_b, 1, 0, true);
}


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;

  a1 = gfc_real2real (x, gfc_default_double_kind);
  a2 = gfc_real2real (y, gfc_default_double_kind);

  result = gfc_get_constant_expr (BT_REAL, gfc_default_double_kind, &x->where);
  mpfr_mul (result->value.real, a1->value.real, a2->value.real, GFC_RND_MODE);

  gfc_free_expr (a2);
  gfc_free_expr (a1);

  return range_check (result, "DPROD");
}


static gfc_expr *
simplify_dshift (gfc_expr *arg1, gfc_expr *arg2, gfc_expr *shiftarg,
		      bool right)
{
  gfc_expr *result;
  int i, k, size, shift;

  if (arg1->expr_type != EXPR_CONSTANT || arg2->expr_type != EXPR_CONSTANT
      || shiftarg->expr_type != EXPR_CONSTANT)
    return NULL;

  k = gfc_validate_kind (BT_INTEGER, arg1->ts.kind, false);
  size = gfc_integer_kinds[k].bit_size;

  gfc_extract_int (shiftarg, &shift);

  /* DSHIFTR(I,J,SHIFT) = DSHIFTL(I,J,SIZE-SHIFT).  */
  if (right)
    shift = size - shift;

  result = gfc_get_constant_expr (BT_INTEGER, arg1->ts.kind, &arg1->where);
  mpz_set_ui (result->value.integer, 0);

  for (i = 0; i < shift; i++)
    if (mpz_tstbit (arg2->value.integer, size - shift + i))
      mpz_setbit (result->value.integer, i);

  for (i = 0; i < size - shift; i++)
    if (mpz_tstbit (arg1->value.integer, i))
      mpz_setbit (result->value.integer, shift + i);

  /* Convert to a signed value.  */
  gfc_convert_mpz_to_signed (result->value.integer, size);

  return result;
}


gfc_expr *
gfc_simplify_dshiftr (gfc_expr *arg1, gfc_expr *arg2, gfc_expr *shiftarg)
{
  return simplify_dshift (arg1, arg2, shiftarg, true);
}


gfc_expr *
gfc_simplify_dshiftl (gfc_expr *arg1, gfc_expr *arg2, gfc_expr *shiftarg)
{
  return simplify_dshift (arg1, arg2, shiftarg, false);
}


gfc_expr *
gfc_simplify_erf (gfc_expr *x)
{
  gfc_expr *result;

  if (x->expr_type != EXPR_CONSTANT)
    return NULL;

  result = gfc_get_constant_expr (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_get_constant_expr (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:

    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)
{
  mpfr_t sum, x, u, v, w, oldsum, sumtrunc;
  mpz_t num;
  mp_prec_t prec;
  unsigned i;

  prec = mpfr_get_default_prec ();
  mpfr_set_default_prec (2 * prec);

  mpfr_init (sum);
  mpfr_init (x);
  mpfr_init (u);
  mpfr_init (v);
  mpfr_init (w);
  mpz_init (num);

  mpfr_init (oldsum);
  mpfr_init (sumtrunc);
  mpfr_set_prec (oldsum, prec);
  mpfr_set_prec (sumtrunc, prec);

  mpfr_set (x, arg, GFC_RND_MODE);
  mpfr_set_ui (sum, 1, GFC_RND_MODE);
  mpz_set_ui (num, 1);

  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);

  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);

    mpfr_set (w, u, GFC_RND_MODE);
    mpfr_pow_ui (w, w, i, GFC_RND_MODE);

    mpfr_set_z (v, num, GFC_RND_MODE);
    mpfr_mul (v, v, w, GFC_RND_MODE);

    mpfr_add (sum, sum, v, GFC_RND_MODE);

    mpfr_set (sumtrunc, sum, GFC_RND_MODE);
    if (mpfr_cmp (sumtrunc, oldsum) == 0)
      break;
  }

  /* 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_scaled (gfc_expr *x)
{
  gfc_expr *result;

  if (x->expr_type != EXPR_CONSTANT)
    return NULL;

  result = gfc_get_constant_expr (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);

  return range_check (result, "ERFC_SCALED");
}

#undef MAX_ITER
#undef ARG_LIMIT


gfc_expr *
gfc_simplify_epsilon (gfc_expr *e)
{
  gfc_expr *result;
  int i;

  i = gfc_validate_kind (e->ts.type, e->ts.kind, false);

  result = gfc_get_constant_expr (BT_REAL, e->ts.kind, &e->where);
  mpfr_set (result->value.real, gfc_real_kinds[i].epsilon, GFC_RND_MODE);

  return range_check (result, "EPSILON");
}


gfc_expr *
gfc_simplify_exp (gfc_expr *x)
{
  gfc_expr *result;

  if (x->expr_type != EXPR_CONSTANT)
    return NULL;

  result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);

  switch (x->ts.type)
    {
      case BT_REAL:
	mpfr_exp (result->value.real, x->value.real, GFC_RND_MODE);
	break;

      case BT_COMPLEX:
	gfc_set_model_kind (x->ts.kind);
	mpc_exp (result->value.complex, x->value.complex, GFC_MPC_RND_MODE);
	break;

      default:
	gfc_internal_error ("in gfc_simplify_exp(): Bad type");
    }

  return range_check (result, "EXP");
}


gfc_expr *
gfc_simplify_exponent (gfc_expr *x)
{
  long int val;
  gfc_expr *result;

  if (x->expr_type != EXPR_CONSTANT)
    return NULL;

  result = gfc_get_constant_expr (BT_INTEGER, gfc_default_integer_kind,
				  &x->where);

  /* EXPONENT(inf) = EXPONENT(nan) = HUGE(0) */
  if (mpfr_inf_p (x->value.real) || mpfr_nan_p (x->value.real))
    {
      int i = gfc_validate_kind (BT_INTEGER, gfc_default_integer_kind, false);
      mpz_set (result->value.integer, gfc_integer_kinds[i].huge);
      return result;
    }

  /* EXPONENT(+/- 0.0) = 0  */
  if (mpfr_zero_p (x->value.real))
    {
      mpz_set_ui (result->value.integer, 0);
      return result;
    }

  gfc_set_model (x->value.real);

  val = (long int) mpfr_get_exp (x->value.real);
  mpz_set_si (result->value.integer, val);

  return range_check (result, "EXPONENT");
}


gfc_expr *
gfc_simplify_failed_or_stopped_images (gfc_expr *team ATTRIBUTE_UNUSED,
				       gfc_expr *kind)
{
  if (flag_coarray == GFC_FCOARRAY_NONE)
    {
      gfc_current_locus = *gfc_current_intrinsic_where;
      gfc_fatal_error ("Coarrays disabled at %C, use %<-fcoarray=%> to enable");
      return &gfc_bad_expr;
    }

  if (flag_coarray == GFC_FCOARRAY_SINGLE)
    {
      gfc_expr *result;
      int actual_kind;
      if (kind)
	gfc_extract_int (kind, &actual_kind);
      else
	actual_kind = gfc_default_integer_kind;

      result = gfc_get_array_expr (BT_INTEGER, actual_kind, &gfc_current_locus);
      result->rank = 1;
      return result;
    }

  /* For fcoarray = lib no simplification is possible, because it is not known
     what images failed or are stopped at compile time.  */
  return NULL;
}


gfc_expr *
gfc_simplify_float (gfc_expr *a)
{
  gfc_expr *result;

  if (a->expr_type != EXPR_CONSTANT)
    return NULL;

  if (a->is_boz)
    {
      if (convert_boz (a, gfc_default_real_kind) == &gfc_bad_expr)
	return &gfc_bad_expr;

      result = gfc_copy_expr (a);
    }
  else
    result = gfc_int2real (a, gfc_default_real_kind);

  return range_check (result, "FLOAT");
}


static bool
is_last_ref_vtab (gfc_expr *e)
{
  gfc_ref *ref;
  gfc_component *comp = NULL;

  if (e->expr_type != EXPR_VARIABLE)
    return false;

  for (ref = e->ref; ref; ref = ref->next)
    if (ref->type == REF_COMPONENT)
      comp = ref->u.c.component;

  if (!e->ref || !comp)
    return e->symtree->n.sym->attr.vtab;

  if (comp->name[0] == '_' && strcmp (comp->name, "_vptr") == 0)
    return true;

  return false;
}


gfc_expr *
gfc_simplify_extends_type_of (gfc_expr *a, gfc_expr *mold)
{
  /* Avoid simplification of resolved symbols.  */
  if (is_last_ref_vtab (a) || is_last_ref_vtab (mold))
    return NULL;

  if (a->ts.type == BT_DERIVED && mold->ts.type == BT_DERIVED)
    return gfc_get_logical_expr (gfc_default_logical_kind, &a->where,
				 gfc_type_is_extension_of (mold->ts.u.derived,
							   a->ts.u.derived));

  if (UNLIMITED_POLY (a) || UNLIMITED_POLY (mold))
    return NULL;

  /* Return .false. if the dynamic type can never be an extension.  */
  if ((a->ts.type == BT_CLASS && mold->ts.type == BT_CLASS
       && !gfc_type_is_extension_of
			(mold->ts.u.derived->components->ts.u.derived,
			 a->ts.u.derived->components->ts.u.derived)
       && !gfc_type_is_extension_of
			(a->ts.u.derived->components->ts.u.derived,
			 mold->ts.u.derived->components->ts.u.derived))
      || (a->ts.type == BT_DERIVED && mold->ts.type == BT_CLASS
	  && !gfc_type_is_extension_of
			(mold->ts.u.derived->components->ts.u.derived,
			 a->ts.u.derived))
      || (a->ts.type == BT_CLASS && mold->ts.type == BT_DERIVED
	  && !gfc_type_is_extension_of
			(mold->ts.u.derived,
			 a->ts.u.derived->components->ts.u.derived)
	  && !gfc_type_is_extension_of
			(a->ts.u.derived->components->ts.u.derived,
			 mold->ts.u.derived)))
    return gfc_get_logical_expr (gfc_default_logical_kind, &a->where, false);

  /* Return .true. if the dynamic type is guaranteed to be an extension.  */
  if (a->ts.type == BT_CLASS && mold->ts.type == BT_DERIVED
      && gfc_type_is_extension_of (mold->ts.u.derived,
				   a->ts.u.derived->components->ts.u.derived))
    return gfc_get_logical_expr (gfc_default_logical_kind, &a->where, true);

  return NULL;
}


gfc_expr *
gfc_simplify_same_type_as (gfc_expr *a, gfc_expr *b)
{
  /* Avoid simplification of resolved symbols.  */
  if (is_last_ref_vtab (a) || is_last_ref_vtab (b))
    return NULL;

  /* Return .false. if the dynamic type can never be the
     same.  */
  if (((a->ts.type == BT_CLASS && gfc_expr_attr (a).class_ok)
       || (b->ts.type == BT_CLASS && gfc_expr_attr (b).class_ok))
      && !gfc_type_compatible (&a->ts, &b->ts)
      && !gfc_type_compatible (&b->ts, &a->ts))
    return gfc_get_logical_expr (gfc_default_logical_kind, &a->where, false);

  if (a->ts.type != BT_DERIVED || b->ts.type != BT_DERIVED)
     return NULL;

  return gfc_get_logical_expr (gfc_default_logical_kind, &a->where,
			       gfc_compare_derived_types (a->ts.u.derived,
							  b->ts.u.derived));
}


gfc_expr *
gfc_simplify_floor (gfc_expr *e, gfc_expr *k)
{
  gfc_expr *result;
  mpfr_t floor;
  int kind;

  kind = get_kind (BT_INTEGER, k, "FLOOR", gfc_default_integer_kind);
  if (kind == -1)
    gfc_internal_error ("gfc_simplify_floor(): Bad kind");

  if (e->expr_type != EXPR_CONSTANT)
    return NULL;

  mpfr_init2 (floor, mpfr_get_prec (e->value.real));
  mpfr_floor (floor, e->value.real);

  result = gfc_get_constant_expr (BT_INTEGER, kind, &e->where);
  gfc_mpfr_to_mpz (result->value.integer, floor, &e->where);

  mpfr_clear (floor);

  return range_check (result, "FLOOR");
}


gfc_expr *
gfc_simplify_fraction (gfc_expr *x)
{
  gfc_expr *result;

#if MPFR_VERSION < MPFR_VERSION_NUM(3,1,0)
  mpfr_t absv, exp, pow2;
#else
  mpfr_exp_t e;
#endif

  if (x->expr_type != EXPR_CONSTANT)
    return NULL;

  result = gfc_get_constant_expr (BT_REAL, x->ts.kind, &x->where);

  /* FRACTION(inf) = NaN.  */
  if (mpfr_inf_p (x->value.real))
    {
      mpfr_set_nan (result->value.real);
      return result;
    }

#if MPFR_VERSION < MPFR_VERSION_NUM(3,1,0)

  /* MPFR versions before 3.1.0 do not include mpfr_frexp.  
     TODO: remove the kludge when MPFR 3.1.0 or newer will be required */

  if (mpfr_sgn (x->value.real) == 0)
    {
      mpfr_set (result->value.real, x->value.real, GFC_RND_MODE);
      return result;
    }

  gfc_set_model_kind (x->ts.kind);
  mpfr_init (exp);
  mpfr_init (absv);
  mpfr_init (pow2);

  mpfr_abs (absv, x->value.real, GFC_RND_MODE);
  mpfr_log2 (exp, absv, GFC_RND_MODE);

  mpfr_trunc (exp, exp);
  mpfr_add_ui (exp, exp, 1, GFC_RND_MODE);

  mpfr_ui_pow (pow2, 2, exp, GFC_RND_MODE);

  mpfr_div (result->value.real, x->value.real, pow2, GFC_RND_MODE);

  mpfr_clears (exp, absv, pow2, NULL);

#else

  /* mpfr_frexp() correctly handles zeros and NaNs.  */
  mpfr_frexp (&e, result->value.real, x->value.real, GFC_RND_MODE);

#endif

  return range_check (result, "FRACTION");
}


gfc_expr *
gfc_simplify_gamma (gfc_expr *x)
{
  gfc_expr *result;

  if (x->expr_type != EXPR_CONSTANT)
    return NULL;

  result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
  mpfr_gamma (result->value.real, x->value.real, GFC_RND_MODE);

  return range_check (result, "GAMMA");
}


gfc_expr *
gfc_simplify_huge (gfc_expr *e)
{
  gfc_expr *result;
  int i;

  i = gfc_validate_kind (e->ts.type, e->ts.kind, false);
  result = gfc_get_constant_expr (e->ts.type, e->ts.kind, &e->where);

  switch (e->ts.type)
    {
      case BT_INTEGER:
	mpz_set (result->value.integer, gfc_integer_kinds[i].huge);
	break;

      case BT_REAL:
	mpfr_set (result->value.real, gfc_real_kinds[i].huge, GFC_RND_MODE);
	break;

      default:
	gcc_unreachable ();
    }

  return result;
}


gfc_expr *
gfc_simplify_hypot (gfc_expr *x, gfc_expr *y)
{
  gfc_expr *result;

  if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
    return NULL;

  result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
  mpfr_hypot (result->value.real, x->value.real, y->value.real, GFC_RND_MODE);
  return range_check (result, "HYPOT");
}


/* We use the processor's collating sequence, because all
   systems that gfortran currently works on are ASCII.  */

gfc_expr *
gfc_simplify_iachar (gfc_expr *e, gfc_expr *kind)
{
  gfc_expr *result;
  gfc_char_t index;
  int k;

  if (e->expr_type != EXPR_CONSTANT)
    return NULL;

  if (e->value.character.length != 1)
    {
      gfc_error ("Argument of IACHAR at %L must be of length one", &e->where);
      return &gfc_bad_expr;
    }

  index = e->value.character.string[0];

  if (warn_surprising && index > 127)
    gfc_warning (OPT_Wsurprising,
		 "Argument of IACHAR function at %L outside of range 0..127",
		 &e->where);

  k = get_kind (BT_INTEGER, kind, "IACHAR", gfc_default_integer_kind);
  if (k == -1)
    return &gfc_bad_expr;

  result = gfc_get_int_expr (k, &e->where, index);

  return range_check (result, "IACHAR");
}


static gfc_expr *
do_bit_and (gfc_expr *result, gfc_expr *e)
{
  gcc_assert (e->ts.type == BT_INTEGER && e->expr_type == EXPR_CONSTANT);
  gcc_assert (result->ts.type == BT_INTEGER
	      && result->expr_type == EXPR_CONSTANT);

  mpz_and (result->value.integer, result->value.integer, e->value.integer);
  return result;
}


gfc_expr *
gfc_simplify_iall (gfc_expr *array, gfc_expr *dim, gfc_expr *mask)
{
  return simplify_transformation (array, dim, mask, -1, do_bit_and);
}


static gfc_expr *
do_bit_ior (gfc_expr *result, gfc_expr *e)
{
  gcc_assert (e->ts.type == BT_INTEGER && e->expr_type == EXPR_CONSTANT);
  gcc_assert (result->ts.type == BT_INTEGER
	      && result->expr_type == EXPR_CONSTANT);

  mpz_ior (result->value.integer, result->value.integer, e->value.integer);
  return result;
}


gfc_expr *
gfc_simplify_iany (gfc_expr *array, gfc_expr *dim, gfc_expr *mask)
{
  return simplify_transformation (array, dim, mask, 0, do_bit_ior);
}


gfc_expr *
gfc_simplify_iand (gfc_expr *x, gfc_expr *y)
{
  gfc_expr *result;

  if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
    return NULL;

  result = gfc_get_constant_expr (BT_INTEGER, x->ts.kind, &x->where);
  mpz_and (result->value.integer, x->value.integer, y->value.integer);

  return range_check (result, "IAND");
}


gfc_expr *
gfc_simplify_ibclr (gfc_expr *x, gfc_expr *y)
{
  gfc_expr *result;
  int k, pos;

  if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
    return NULL;

  gfc_extract_int (y, &pos);

  k = gfc_validate_kind (x->ts.type, x->ts.kind, false);

  result = gfc_copy_expr (x);

  convert_mpz_to_unsigned (result->value.integer,
			   gfc_integer_kinds[k].bit_size);

  mpz_clrbit (result->value.integer, pos);

  gfc_convert_mpz_to_signed (result->value.integer,
			 gfc_integer_kinds[k].bit_size);

  return result;
}


gfc_expr *
gfc_simplify_ibits (gfc_expr *x, gfc_expr *y, gfc_expr *z)
{
  gfc_expr *result;
  int pos, len;
  int i, k, bitsize;
  int *bits;

  if (x->expr_type != EXPR_CONSTANT
      || y->expr_type != EXPR_CONSTANT
      || z->expr_type != EXPR_CONSTANT)
    return NULL;

  gfc_extract_int (y, &pos);
  gfc_extract_int (z, &len);

  k = gfc_validate_kind (BT_INTEGER, x->ts.kind, false);

  bitsize = gfc_integer_kinds[k].bit_size;

  if (pos + len > bitsize)
    {
      gfc_error ("Sum of second and third arguments of IBITS exceeds "
		 "bit size at %L", &y->where);
      return &gfc_bad_expr;
    }

  result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
  convert_mpz_to_unsigned (result->value.integer,
			   gfc_integer_kinds[k].bit_size);

  bits = XCNEWVEC (int, bitsize);

  for (i = 0; i < bitsize; i++)
    bits[i] = 0;

  for (i = 0; i < len; i++)
    bits[i] = mpz_tstbit (x->value.integer, i + pos);

  for (i = 0; i < bitsize; i++)
    {
      if (bits[i] == 0)
	mpz_clrbit (result->value.integer, i);
      else if (bits[i] == 1)
	mpz_setbit (result->value.integer, i);
      else
	gfc_internal_error ("IBITS: Bad bit");
    }

  free (bits);

  gfc_convert_mpz_to_signed (result->value.integer,
			 gfc_integer_kinds[k].bit_size);

  return result;
}


gfc_expr *
gfc_simplify_ibset (gfc_expr *x, gfc_expr *y)
{
  gfc_expr *result;
  int k, pos;

  if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
    return NULL;

  gfc_extract_int (y, &pos);

  k = gfc_validate_kind (x->ts.type, x->ts.kind, false);

  result = gfc_copy_expr (x);

  convert_mpz_to_unsigned (result->value.integer,
			   gfc_integer_kinds[k].bit_size);

  mpz_setbit (result->value.integer, pos);

  gfc_convert_mpz_to_signed (result->value.integer,
			 gfc_integer_kinds[k].bit_size);

  return result;
}


gfc_expr *
gfc_simplify_ichar (gfc_expr *e, gfc_expr *kind)
{
  gfc_expr *result;
  gfc_char_t index;
  int k;

  if (e->expr_type != EXPR_CONSTANT)
    return NULL;

  if (e->value.character.length != 1)
    {
      gfc_error ("Argument of ICHAR at %L must be of length one", &e->where);
      return &gfc_bad_expr;
    }

  index = e->value.character.string[0];

  k = get_kind (BT_INTEGER, kind, "ICHAR", gfc_default_integer_kind);
  if (k == -1)
    return &gfc_bad_expr;

  result = gfc_get_int_expr (k, &e->where, index);

  return range_check (result, "ICHAR");
}


gfc_expr *
gfc_simplify_ieor (gfc_expr *x, gfc_expr *y)
{
  gfc_expr *result;

  if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
    return NULL;

  result = gfc_get_constant_expr (BT_INTEGER, x->ts.kind, &x->where);
  mpz_xor (result->value.integer, x->value.integer, y->value.integer);

  return range_check (result, "IEOR");
}


gfc_expr *
gfc_simplify_index (gfc_expr *x, gfc_expr *y, gfc_expr *b, gfc_expr *kind)
{
  gfc_expr *result;
  int back, len, lensub;
  int i, j, k, count, index = 0, start;

  if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT
      || ( b != NULL && b->expr_type !=  EXPR_CONSTANT))
    return NULL;

  if (b != NULL && b->value.logical != 0)
    back = 1;
  else
    back = 0;

  k = get_kind (BT_INTEGER, kind, "INDEX", gfc_default_integer_kind);
  if (k == -1)
    return &gfc_bad_expr;

  result = gfc_get_constant_expr (BT_INTEGER, k, &x->where);

  len = x->value.character.length;
  lensub = y->value.character.length;

  if (len < lensub)
    {
      mpz_set_si (result->value.integer, 0);
      return result;
    }

  if (back == 0)
    {
      if (lensub == 0)
	{
	  mpz_set_si (result->value.integer, 1);
	  return result;
	}
      else if (lensub == 1)
	{
	  for (i = 0; i < len; i++)
	    {
	      for (j = 0; j < lensub; j++)
		{
		  if (y->value.character.string[j]
		      == x->value.character.string[i])
		    {
		      index = i + 1;
		      goto done;
		    }
		}
	    }
	}
      else
	{
	  for (i = 0; i < len; i++)
	    {
	      for (j = 0; j < lensub; j++)
		{
		  if (y->value.character.string[j]
		      == x->value.character.string[i])
		    {
		      start = i;
		      count = 0;

		      for (k = 0; k < lensub; k++)
			{
			  if (y->value.character.string[k]
			      == x->value.character.string[k + start])
			    count++;
			}

		      if (count == lensub)
			{
			  index = start + 1;
			  goto done;