arith.c

gcc-fortran,linux使用fortran的编译软件。很好用的。

  return val;
}


/* It may seem silly to have a subroutine that actually computes the
   unary plus of a constant, but it prevents us from making exceptions
   in the code elsewhere.  */

static arith
gfc_arith_uplus (gfc_expr * op1, gfc_expr ** resultp)
{
  *resultp = gfc_copy_expr (op1);
  return ARITH_OK;
}


static arith
gfc_arith_uminus (gfc_expr * op1, gfc_expr ** resultp)
{
  gfc_expr *result;
  arith rc;

  result = gfc_constant_result (op1->ts.type, op1->ts.kind, &op1->where);

  switch (op1->ts.type)
    {
    case BT_INTEGER:
      mpz_neg (result->value.integer, op1->value.integer);
      break;

    case BT_REAL:
      mpfr_neg (result->value.real, op1->value.real, GFC_RND_MODE);
      break;

    case BT_COMPLEX:
      mpfr_neg (result->value.complex.r, op1->value.complex.r, GFC_RND_MODE);
      mpfr_neg (result->value.complex.i, op1->value.complex.i, GFC_RND_MODE);
      break;

    default:
      gfc_internal_error ("gfc_arith_uminus(): Bad basic type");
    }

  rc = gfc_range_check (result);

  return check_result (rc, op1, result, resultp);
}


static arith
gfc_arith_plus (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
{
  gfc_expr *result;
  arith rc;

  result = gfc_constant_result (op1->ts.type, op1->ts.kind, &op1->where);

  switch (op1->ts.type)
    {
    case BT_INTEGER:
      mpz_add (result->value.integer, op1->value.integer, op2->value.integer);
      break;

    case BT_REAL:
      mpfr_add (result->value.real, op1->value.real, op2->value.real,
               GFC_RND_MODE);
      break;

    case BT_COMPLEX:
      mpfr_add (result->value.complex.r, op1->value.complex.r,
	       op2->value.complex.r, GFC_RND_MODE);

      mpfr_add (result->value.complex.i, op1->value.complex.i,
	       op2->value.complex.i, GFC_RND_MODE);
      break;

    default:
      gfc_internal_error ("gfc_arith_plus(): Bad basic type");
    }

  rc = gfc_range_check (result);

  return check_result (rc, op1, result, resultp);
}


static arith
gfc_arith_minus (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
{
  gfc_expr *result;
  arith rc;

  result = gfc_constant_result (op1->ts.type, op1->ts.kind, &op1->where);

  switch (op1->ts.type)
    {
    case BT_INTEGER:
      mpz_sub (result->value.integer, op1->value.integer, op2->value.integer);
      break;

    case BT_REAL:
      mpfr_sub (result->value.real, op1->value.real, op2->value.real,
                GFC_RND_MODE);
      break;

    case BT_COMPLEX:
      mpfr_sub (result->value.complex.r, op1->value.complex.r,
	       op2->value.complex.r, GFC_RND_MODE);

      mpfr_sub (result->value.complex.i, op1->value.complex.i,
	       op2->value.complex.i, GFC_RND_MODE);
      break;

    default:
      gfc_internal_error ("gfc_arith_minus(): Bad basic type");
    }

  rc = gfc_range_check (result);

  return check_result (rc, op1, result, resultp);
}


static arith
gfc_arith_times (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
{
  gfc_expr *result;
  mpfr_t x, y;
  arith rc;

  result = gfc_constant_result (op1->ts.type, op1->ts.kind, &op1->where);

  switch (op1->ts.type)
    {
    case BT_INTEGER:
      mpz_mul (result->value.integer, op1->value.integer, op2->value.integer);
      break;

    case BT_REAL:
      mpfr_mul (result->value.real, op1->value.real, op2->value.real,
               GFC_RND_MODE);
      break;

    case BT_COMPLEX:

      /* FIXME:  possible numericals problem.  */

      gfc_set_model (op1->value.complex.r);
      mpfr_init (x);
      mpfr_init (y);

      mpfr_mul (x, op1->value.complex.r, op2->value.complex.r, GFC_RND_MODE);
      mpfr_mul (y, op1->value.complex.i, op2->value.complex.i, GFC_RND_MODE);
      mpfr_sub (result->value.complex.r, x, y, GFC_RND_MODE);

      mpfr_mul (x, op1->value.complex.r, op2->value.complex.i, GFC_RND_MODE);
      mpfr_mul (y, op1->value.complex.i, op2->value.complex.r, GFC_RND_MODE);
      mpfr_add (result->value.complex.i, x, y, GFC_RND_MODE);

      mpfr_clear (x);
      mpfr_clear (y);

      break;

    default:
      gfc_internal_error ("gfc_arith_times(): Bad basic type");
    }

  rc = gfc_range_check (result);

  return check_result (rc, op1, result, resultp);
}


static arith
gfc_arith_divide (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
{
  gfc_expr *result;
  mpfr_t x, y, div;
  arith rc;

  rc = ARITH_OK;

  result = gfc_constant_result (op1->ts.type, op1->ts.kind, &op1->where);

  switch (op1->ts.type)
    {
    case BT_INTEGER:
      if (mpz_sgn (op2->value.integer) == 0)
	{
	  rc = ARITH_DIV0;
	  break;
	}

      mpz_tdiv_q (result->value.integer, op1->value.integer,
		  op2->value.integer);
      break;

    case BT_REAL:
      /* FIXME: MPFR correctly generates NaN.  This may not be needed.  */
      if (mpfr_sgn (op2->value.real) == 0)
	{
	  rc = ARITH_DIV0;
	  break;
	}

      mpfr_div (result->value.real, op1->value.real, op2->value.real,
               GFC_RND_MODE);
      break;

    case BT_COMPLEX:
      /* FIXME: MPFR correctly generates NaN.  This may not be needed.  */
      if (mpfr_sgn (op2->value.complex.r) == 0
	  && mpfr_sgn (op2->value.complex.i) == 0)
	{
	  rc = ARITH_DIV0;
	  break;
	}

      gfc_set_model (op1->value.complex.r);
      mpfr_init (x);
      mpfr_init (y);
      mpfr_init (div);

      /* FIXME: possible numerical problems.  */
      mpfr_mul (x, op2->value.complex.r, op2->value.complex.r, GFC_RND_MODE);
      mpfr_mul (y, op2->value.complex.i, op2->value.complex.i, GFC_RND_MODE);
      mpfr_add (div, x, y, GFC_RND_MODE);

      mpfr_mul (x, op1->value.complex.r, op2->value.complex.r, GFC_RND_MODE);
      mpfr_mul (y, op1->value.complex.i, op2->value.complex.i, GFC_RND_MODE);
      mpfr_add (result->value.complex.r, x, y, GFC_RND_MODE);
      mpfr_div (result->value.complex.r, result->value.complex.r, div,
                GFC_RND_MODE);

      mpfr_mul (x, op1->value.complex.i, op2->value.complex.r, GFC_RND_MODE);
      mpfr_mul (y, op1->value.complex.r, op2->value.complex.i, GFC_RND_MODE);
      mpfr_sub (result->value.complex.i, x, y, GFC_RND_MODE);
      mpfr_div (result->value.complex.i, result->value.complex.i, div,
                GFC_RND_MODE);

      mpfr_clear (x);
      mpfr_clear (y);
      mpfr_clear (div);

      break;

    default:
      gfc_internal_error ("gfc_arith_divide(): Bad basic type");
    }

  if (rc == ARITH_OK)
    rc = gfc_range_check (result);

  return check_result (rc, op1, result, resultp);
}


/* Compute the reciprocal of a complex number (guaranteed nonzero).  */

static void
complex_reciprocal (gfc_expr * op)
{
  mpfr_t mod, a, re, im;

  gfc_set_model (op->value.complex.r);
  mpfr_init (mod);
  mpfr_init (a);
  mpfr_init (re);
  mpfr_init (im);

  /* FIXME:  another possible numerical problem.  */
  mpfr_mul (mod, op->value.complex.r, op->value.complex.r, GFC_RND_MODE);
  mpfr_mul (a, op->value.complex.i, op->value.complex.i, GFC_RND_MODE);
  mpfr_add (mod, mod, a, GFC_RND_MODE);

  mpfr_div (re, op->value.complex.r, mod, GFC_RND_MODE);

  mpfr_neg (im, op->value.complex.i, GFC_RND_MODE);
  mpfr_div (im, im, mod, GFC_RND_MODE);

  mpfr_set (op->value.complex.r, re, GFC_RND_MODE);
  mpfr_set (op->value.complex.i, im, GFC_RND_MODE);

  mpfr_clear (re);
  mpfr_clear (im);
  mpfr_clear (mod);
  mpfr_clear (a);
}


/* Raise a complex number to positive power.  */

static void
complex_pow_ui (gfc_expr * base, int power, gfc_expr * result)
{
  mpfr_t re, im, a;

  gfc_set_model (base->value.complex.r);
  mpfr_init (re);
  mpfr_init (im);
  mpfr_init (a);

  mpfr_set_ui (result->value.complex.r, 1, GFC_RND_MODE);
  mpfr_set_ui (result->value.complex.i, 0, GFC_RND_MODE);

  for (; power > 0; power--)
    {
      mpfr_mul (re, base->value.complex.r, result->value.complex.r,
                GFC_RND_MODE);
      mpfr_mul (a, base->value.complex.i, result->value.complex.i,
                GFC_RND_MODE);
      mpfr_sub (re, re, a, GFC_RND_MODE);

      mpfr_mul (im, base->value.complex.r, result->value.complex.i,
                GFC_RND_MODE);
      mpfr_mul (a, base->value.complex.i, result->value.complex.r,
                GFC_RND_MODE);
      mpfr_add (im, im, a, GFC_RND_MODE);

      mpfr_set (result->value.complex.r, re, GFC_RND_MODE);
      mpfr_set (result->value.complex.i, im, GFC_RND_MODE);
    }

  mpfr_clear (re);
  mpfr_clear (im);
  mpfr_clear (a);
}


/* Raise a number to an integer power.  */

static arith
gfc_arith_power (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
{
  int power, apower;
  gfc_expr *result;
  mpz_t unity_z;
  mpfr_t unity_f;
  arith rc;

  rc = ARITH_OK;

  if (gfc_extract_int (op2, &power) != NULL)
    gfc_internal_error ("gfc_arith_power(): Bad exponent");

  result = gfc_constant_result (op1->ts.type, op1->ts.kind, &op1->where);

  if (power == 0)
    {
      /* Handle something to the zeroth power.  Since we're dealing
	 with integral exponents, there is no ambiguity in the
	 limiting procedure used to determine the value of 0**0.  */
      switch (op1->ts.type)
	{
	case BT_INTEGER:
	  mpz_set_ui (result->value.integer, 1);
	  break;

	case BT_REAL:
	  mpfr_set_ui (result->value.real, 1, GFC_RND_MODE);
	  break;

	case BT_COMPLEX:
	  mpfr_set_ui (result->value.complex.r, 1, GFC_RND_MODE);
	  mpfr_set_ui (result->value.complex.i, 0, GFC_RND_MODE);
	  break;

	default:
	  gfc_internal_error ("gfc_arith_power(): Bad base");
	}
    }
  else
    {
      apower = power;
      if (power < 0)
	apower = -power;

      switch (op1->ts.type)
	{
	case BT_INTEGER:
	  mpz_pow_ui (result->value.integer, op1->value.integer, apower);

	  if (power < 0)
	    {
	      mpz_init_set_ui (unity_z, 1);
	      mpz_tdiv_q (result->value.integer, unity_z,
			  result->value.integer);
	      mpz_clear (unity_z);
	    }

	  break;

	case BT_REAL:
	  mpfr_pow_ui (result->value.real, op1->value.real, apower,
                       GFC_RND_MODE);

	  if (power < 0)
	    {
              gfc_set_model (op1->value.real);
	      mpfr_init (unity_f);
	      mpfr_set_ui (unity_f, 1, GFC_RND_MODE);
	      mpfr_div (result->value.real, unity_f, result->value.real,
                        GFC_RND_MODE);
	      mpfr_clear (unity_f);
	    }
	  break;

	case BT_COMPLEX:
	  complex_pow_ui (op1, apower, result);
	  if (power < 0)
	    complex_reciprocal (result);
	  break;

	default:
	  break;
	}
    }

  if (rc == ARITH_OK)
    rc = gfc_range_check (result);

  return check_result (rc, op1, result, resultp);
}


/* Concatenate two string constants.  */

static arith
gfc_arith_concat (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
{
  gfc_expr *result;
  int len;

  result = gfc_constant_result (BT_CHARACTER, gfc_default_character_kind,
				&op1->where);

  len = op1->value.character.length + op2->value.character.length;

  result->value.character.string = gfc_getmem (len + 1);
  result->value.character.length = len;

  memcpy (result->value.character.string, op1->value.character.string,
	  op1->value.character.length);

  memcpy (result->value.character.string + op1->value.character.length,
	  op2->value.character.string, op2->value.character.length);

  result->value.character.string[len] = '\0';

  *resultp = result;

  return ARITH_OK;
}


/* Comparison operators.  Assumes that the two expression nodes
   contain two constants of the same type.  */

int
gfc_compare_expr (gfc_expr * op1, gfc_expr * op2)
{
  int rc;

  switch (op1->ts.type)
    {
    case BT_INTEGER:
      rc = mpz_cmp (op1->value.integer, op2->value.integer);
      break;

    case BT_REAL:
      rc = mpfr_cmp (op1->value.real, op2->value.real);
      break;

    case BT_CHARACTER:
      rc = gfc_compare_string (op1, op2, NULL);
      break;

    case BT_LOGICAL:
      rc = ((!op1->value.logical && op2->value.logical)
	    || (op1->value.logical && !op2->value.logical));
      break;

    default:
      gfc_internal_error ("gfc_compare_expr(): Bad basic type");
    }

  return rc;
}


/* Compare a pair of complex numbers.  Naturally, this is only for
   equality/nonequality.  */

static int
compare_complex (gfc_expr * op1, gfc_expr * op2)
{
  return (mpfr_cmp (op1->value.complex.r, op2->value.complex.r) == 0
	  && mpfr_cmp (op1->value.complex.i, op2->value.complex.i) == 0);
}


/* Given two constant strings and the inverse collating sequence,
   compare the strings.  We return -1 for a<b, 0 for a==b and 1 for
   a>b.  If the xcoll_table is NULL, we use the processor's default
   collating sequence.  */

int
gfc_compare_string (gfc_expr * a, gfc_expr * b, const int *xcoll_table)
{
  int len, alen, blen, i, ac, bc;

  alen = a->value.character.length;
  blen = b->value.character.length;

  len = (alen > blen) ? alen : blen;

  for (i = 0; i < len; i++)
    {
      ac = (i < alen) ? a->value.character.string[i] : ' ';
      bc = (i < blen) ? b->value.character.string[i] : ' ';

      if (xcoll_table != NULL)
	{
	  ac = xcoll_table[ac];
	  bc = xcoll_table[bc];
	}

      if (ac < bc)
	return -1;
      if (ac > bc)
	return 1;
    }

  /* Strings are equal */

  return 0;
}


/* Specific comparison subroutines.  */

static arith
gfc_arith_eq (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
{
  gfc_expr *result;

  result = gfc_constant_result (BT_LOGICAL, gfc_default_logical_kind,
				&op1->where);
  result->value.logical = (op1->ts.type == BT_COMPLEX) ?
    compare_complex (op1, op2) : (gfc_compare_expr (op1, op2) == 0);

  *resultp = result;
  return ARITH_OK;
}


static arith
gfc_arith_ne (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
{
  gfc_expr *result;

  result = gfc_constant_result (BT_LOGICAL, gfc_default_logical_kind,
				&op1->where);
  result->value.logical = (op1->ts.type == BT_COMPLEX) ?
    !compare_complex (op1, op2) : (gfc_compare_expr (op1, op2) != 0);

  *resultp = result;
  return ARITH_OK;
}


static arith
gfc_arith_gt (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
{
  gfc_expr *result;

  result = gfc_constant_result (BT_LOGICAL, gfc_default_logical_kind,
				&op1->where);
  result->value.logical = (gfc_compare_expr (op1, op2) > 0);
  *resultp = result;

  return ARITH_OK;
}


static arith
gfc_arith_ge (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
{
  gfc_expr *result;

  result = gfc_constant_result (BT_LOGICAL, gfc_default_logical_kind,
				&op1->where);
  result->value.logical = (gfc_compare_expr (op1, op2) >= 0);
  *resultp = result;

  return ARITH_OK;
}


static arith
gfc_arith_lt (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
{
  gfc_expr *result;

  result = gfc_constant_result (BT_LOGICAL, gfc_default_logical_kind,
				&op1->where);
  result->value.logical = (gfc_compare_expr (op1, op2) < 0);
  *resultp = result;

  return ARITH_OK;
}