arith.c

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

/* Compiler arithmetic
   Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005 Free Software Foundation,
   Inc.
   Contributed by Andy Vaught

This file is part of GCC.

GCC is free software; you can redistribute it and/or modify it under
the terms of the GNU General Public License as published by the Free
Software Foundation; either version 2, or (at your option) any later
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 COPYING.  If not, write to the Free
Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA
02110-1301, USA.  */

/* Since target arithmetic must be done on the host, there has to
   be some way of evaluating arithmetic expressions as the host
   would evaluate them.  We use the GNU MP library to do arithmetic,
   and this file provides the interface.  */

#include "config.h"
#include "system.h"
#include "flags.h"
#include "gfortran.h"
#include "arith.h"

/* MPFR does not have a direct replacement for mpz_set_f() from GMP.
   It's easily implemented with a few calls though.  */

void
gfc_mpfr_to_mpz (mpz_t z, mpfr_t x)
{
  mp_exp_t e;

  e = mpfr_get_z_exp (z, x);
  /* MPFR 2.0.1 (included with GMP 4.1) has a bug whereby mpfr_get_z_exp
     may set the sign of z incorrectly.  Work around that here.  */
  if (mpfr_sgn (x) != mpz_sgn (z))
    mpz_neg (z, z);

  if (e > 0)
    mpz_mul_2exp (z, z, e);
  else
    mpz_tdiv_q_2exp (z, z, -e);
}


/* Set the model number precision by the requested KIND.  */

void
gfc_set_model_kind (int kind)
{
  int index = gfc_validate_kind (BT_REAL, kind, false);
  int base2prec;

  base2prec = gfc_real_kinds[index].digits;
  if (gfc_real_kinds[index].radix != 2)
    base2prec *= gfc_real_kinds[index].radix / 2;
  mpfr_set_default_prec (base2prec);
}


/* Set the model number precision from mpfr_t x.  */

void
gfc_set_model (mpfr_t x)
{
  mpfr_set_default_prec (mpfr_get_prec (x));
}

/* Calculate atan2 (y, x)

atan2(y, x) = atan(y/x)				if x > 0,
	      sign(y)*(pi - atan(|y/x|))	if x < 0,
	      0					if x = 0 && y == 0,
	      sign(y)*pi/2			if x = 0 && y != 0.
*/

void
arctangent2 (mpfr_t y, mpfr_t x, mpfr_t result)
{
  int i;
  mpfr_t t;

  gfc_set_model (y);
  mpfr_init (t);

  i = mpfr_sgn (x);

  if (i > 0)
    {
      mpfr_div (t, y, x, GFC_RND_MODE);
      mpfr_atan (result, t, GFC_RND_MODE);
    }
  else if (i < 0)
    {
      mpfr_const_pi (result, GFC_RND_MODE);
      mpfr_div (t, y, x, GFC_RND_MODE);
      mpfr_abs (t, t, GFC_RND_MODE);
      mpfr_atan (t, t, GFC_RND_MODE);
      mpfr_sub (result, result, t, GFC_RND_MODE);
      if (mpfr_sgn (y) < 0)
	mpfr_neg (result, result, GFC_RND_MODE);
    }
  else
    {
      if (mpfr_sgn (y) == 0)
	mpfr_set_ui (result, 0, GFC_RND_MODE);
      else
	{
          mpfr_const_pi (result, GFC_RND_MODE);
          mpfr_div_ui (result, result, 2, GFC_RND_MODE);
	  if (mpfr_sgn (y) < 0)
	    mpfr_neg (result, result, GFC_RND_MODE);
	}
    }

  mpfr_clear (t);

}


/* Given an arithmetic error code, return a pointer to a string that
   explains the error.  */

static const char *
gfc_arith_error (arith code)
{
  const char *p;

  switch (code)
    {
    case ARITH_OK:
      p = _("Arithmetic OK at %L");
      break;
    case ARITH_OVERFLOW:
      p = _("Arithmetic overflow at %L");
      break;
    case ARITH_UNDERFLOW:
      p = _("Arithmetic underflow at %L");
      break;
    case ARITH_NAN:
      p = _("Arithmetic NaN at %L");
      break;
    case ARITH_DIV0:
      p = _("Division by zero at %L");
      break;
    case ARITH_INCOMMENSURATE:
      p = _("Array operands are incommensurate at %L");
      break;
    case ARITH_ASYMMETRIC:
      p =
	_("Integer outside symmetric range implied by Standard Fortran at %L");
      break;
    default:
      gfc_internal_error ("gfc_arith_error(): Bad error code");
    }

  return p;
}


/* Get things ready to do math.  */

void
gfc_arith_init_1 (void)
{
  gfc_integer_info *int_info;
  gfc_real_info *real_info;
  mpfr_t a, b, c;
  mpz_t r;
  int i;

  mpfr_set_default_prec (128);
  mpfr_init (a);
  mpz_init (r);

  /* Convert the minimum/maximum values for each kind into their
     GNU MP representation.  */
  for (int_info = gfc_integer_kinds; int_info->kind != 0; int_info++)
    {
      /* Huge */
      mpz_set_ui (r, int_info->radix);
      mpz_pow_ui (r, r, int_info->digits);

      mpz_init (int_info->huge);
      mpz_sub_ui (int_info->huge, r, 1);

      /* These are the numbers that are actually representable by the
         target.  For bases other than two, this needs to be changed.  */
      if (int_info->radix != 2)
        gfc_internal_error ("Fix min_int, max_int calculation");

      /* See PRs 13490 and 17912, related to integer ranges.
         The pedantic_min_int exists for range checking when a program
         is compiled with -pedantic, and reflects the belief that
         Standard Fortran requires integers to be symmetrical, i.e.
         every negative integer must have a representable positive
         absolute value, and vice versa.  */

      mpz_init (int_info->pedantic_min_int);
      mpz_neg (int_info->pedantic_min_int, int_info->huge);

      mpz_init (int_info->min_int);
      mpz_sub_ui (int_info->min_int, int_info->pedantic_min_int, 1);

      mpz_init (int_info->max_int);
      mpz_add (int_info->max_int, int_info->huge, int_info->huge);
      mpz_add_ui (int_info->max_int, int_info->max_int, 1);

      /* Range */
      mpfr_set_z (a, int_info->huge, GFC_RND_MODE);
      mpfr_log10 (a, a, GFC_RND_MODE);
      mpfr_trunc (a, a);
      gfc_mpfr_to_mpz (r, a);
      int_info->range = mpz_get_si (r);
    }

  mpfr_clear (a);

  for (real_info = gfc_real_kinds; real_info->kind != 0; real_info++)
    {
      gfc_set_model_kind (real_info->kind);

      mpfr_init (a);
      mpfr_init (b);
      mpfr_init (c);

      /* huge(x) = (1 - b**(-p)) * b**(emax-1) * b  */
      /* a = 1 - b**(-p) */
      mpfr_set_ui (a, 1, GFC_RND_MODE);
      mpfr_set_ui (b, real_info->radix, GFC_RND_MODE);
      mpfr_pow_si (b, b, -real_info->digits, GFC_RND_MODE);
      mpfr_sub (a, a, b, GFC_RND_MODE);

      /* c = b**(emax-1) */
      mpfr_set_ui (b, real_info->radix, GFC_RND_MODE);
      mpfr_pow_ui (c, b, real_info->max_exponent - 1, GFC_RND_MODE);

      /* a = a * c = (1 - b**(-p)) * b**(emax-1) */
      mpfr_mul (a, a, c, GFC_RND_MODE);

      /* a = (1 - b**(-p)) * b**(emax-1) * b */
      mpfr_mul_ui (a, a, real_info->radix, GFC_RND_MODE);

      mpfr_init (real_info->huge);
      mpfr_set (real_info->huge, a, GFC_RND_MODE);

      /* tiny(x) = b**(emin-1) */
      mpfr_set_ui (b, real_info->radix, GFC_RND_MODE);
      mpfr_pow_si (b, b, real_info->min_exponent - 1, GFC_RND_MODE);

      mpfr_init (real_info->tiny);
      mpfr_set (real_info->tiny, b, GFC_RND_MODE);

      /* subnormal (x) = b**(emin - digit) */
      mpfr_set_ui (b, real_info->radix, GFC_RND_MODE);
      mpfr_pow_si (b, b, real_info->min_exponent - real_info->digits,
		   GFC_RND_MODE);

      mpfr_init (real_info->subnormal);
      mpfr_set (real_info->subnormal, b, GFC_RND_MODE);

      /* epsilon(x) = b**(1-p) */
      mpfr_set_ui (b, real_info->radix, GFC_RND_MODE);
      mpfr_pow_si (b, b, 1 - real_info->digits, GFC_RND_MODE);

      mpfr_init (real_info->epsilon);
      mpfr_set (real_info->epsilon, b, GFC_RND_MODE);

      /* range(x) = int(min(log10(huge(x)), -log10(tiny)) */
      mpfr_log10 (a, real_info->huge, GFC_RND_MODE);
      mpfr_log10 (b, real_info->tiny, GFC_RND_MODE);
      mpfr_neg (b, b, GFC_RND_MODE);

      if (mpfr_cmp (a, b) > 0)
	mpfr_set (a, b, GFC_RND_MODE);		/* a = min(a, b) */

      mpfr_trunc (a, a);
      gfc_mpfr_to_mpz (r, a);
      real_info->range = mpz_get_si (r);

      /* precision(x) = int((p - 1) * log10(b)) + k */
      mpfr_set_ui (a, real_info->radix, GFC_RND_MODE);
      mpfr_log10 (a, a, GFC_RND_MODE);

      mpfr_mul_ui (a, a, real_info->digits - 1, GFC_RND_MODE);
      mpfr_trunc (a, a);
      gfc_mpfr_to_mpz (r, a);
      real_info->precision = mpz_get_si (r);

      /* If the radix is an integral power of 10, add one to the
         precision.  */
      for (i = 10; i <= real_info->radix; i *= 10)
	if (i == real_info->radix)
	  real_info->precision++;

      mpfr_clear (a);
      mpfr_clear (b);
      mpfr_clear (c);
    }

  mpz_clear (r);
}


/* Clean up, get rid of numeric constants.  */

void
gfc_arith_done_1 (void)
{
  gfc_integer_info *ip;
  gfc_real_info *rp;

  for (ip = gfc_integer_kinds; ip->kind; ip++)
    {
      mpz_clear (ip->min_int);
      mpz_clear (ip->max_int);
      mpz_clear (ip->huge);
    }

  for (rp = gfc_real_kinds; rp->kind; rp++)
    {
      mpfr_clear (rp->epsilon);
      mpfr_clear (rp->huge);
      mpfr_clear (rp->tiny);
    }
}


/* Given an integer and a kind, make sure that the integer lies within
   the range of the kind.  Returns ARITH_OK, ARITH_ASYMMETRIC or
   ARITH_OVERFLOW.  */

arith
gfc_check_integer_range (mpz_t p, int kind)
{
  arith result;
  int i;

  i = gfc_validate_kind (BT_INTEGER, kind, false);
  result = ARITH_OK;

  if (pedantic)
    {
      if (mpz_cmp (p, gfc_integer_kinds[i].pedantic_min_int) < 0)
        result = ARITH_ASYMMETRIC;
    }

  if (mpz_cmp (p, gfc_integer_kinds[i].min_int) < 0
      || mpz_cmp (p, gfc_integer_kinds[i].max_int) > 0)
    result = ARITH_OVERFLOW;

  return result;
}


/* Given a real and a kind, make sure that the real lies within the
   range of the kind.  Returns ARITH_OK, ARITH_OVERFLOW or
   ARITH_UNDERFLOW.  */

static arith
gfc_check_real_range (mpfr_t p, int kind)
{
  arith retval;
  mpfr_t q;
  int i;

  i = gfc_validate_kind (BT_REAL, kind, false);

  gfc_set_model (p);
  mpfr_init (q);
  mpfr_abs (q, p, GFC_RND_MODE);

  if (mpfr_sgn (q) == 0)
    retval = ARITH_OK;
  else if (mpfr_cmp (q, gfc_real_kinds[i].huge) > 0)
    retval = ARITH_OVERFLOW;
  else if (mpfr_cmp (q, gfc_real_kinds[i].subnormal) < 0)
    retval = ARITH_UNDERFLOW;
  else if (mpfr_cmp (q, gfc_real_kinds[i].tiny) < 0)
    {
      /* MPFR operates on a numbers with a given precision and enormous
	exponential range.  To represent subnormal numbers the exponent is
	allowed to become smaller than emin, but always retains the full
	precision.  This function resets unused bits to 0 to alleviate
	rounding problems.  Note, a future version of MPFR will have a
 	mpfr_subnormalize() function, which handles this truncation in a
	more efficient and robust way.  */

      int j, k;
      char *bin, *s;
      mp_exp_t e;

      bin = mpfr_get_str (NULL, &e, gfc_real_kinds[i].radix, 0, q, GMP_RNDN);
      k = gfc_real_kinds[i].digits - (gfc_real_kinds[i].min_exponent - e);
      for (j = k; j < gfc_real_kinds[i].digits; j++)
	bin[j] = '0';
      /* Need space for '0.', bin, 'E', and e */
      s = (char *) gfc_getmem (strlen(bin)+10);
      sprintf (s, "0.%sE%d", bin, (int) e);
      mpfr_set_str (q, s, gfc_real_kinds[i].radix, GMP_RNDN);

      if (mpfr_sgn (p) < 0)
	mpfr_neg (p, q, GMP_RNDN);
      else
	mpfr_set (p, q, GMP_RNDN);

      gfc_free (s);
      gfc_free (bin);

      retval = ARITH_OK;
    }
  else
    retval = ARITH_OK;

  mpfr_clear (q);

  return retval;
}


/* Function to return a constant expression node of a given type and
   kind.  */

gfc_expr *
gfc_constant_result (bt type, int kind, locus * where)
{
  gfc_expr *result;

  if (!where)
    gfc_internal_error
      ("gfc_constant_result(): locus 'where' cannot be NULL");

  result = gfc_get_expr ();

  result->expr_type = EXPR_CONSTANT;
  result->ts.type = type;
  result->ts.kind = kind;
  result->where = *where;

  switch (type)
    {
    case BT_INTEGER:
      mpz_init (result->value.integer);
      break;

    case BT_REAL:
      gfc_set_model_kind (kind);
      mpfr_init (result->value.real);
      break;

    case BT_COMPLEX:
      gfc_set_model_kind (kind);
      mpfr_init (result->value.complex.r);
      mpfr_init (result->value.complex.i);
      break;

    default:
      break;
    }

  return result;
}


/* Low-level arithmetic functions.  All of these subroutines assume
   that all operands are of the same type and return an operand of the
   same type.  The other thing about these subroutines is that they
   can fail in various ways -- overflow, underflow, division by zero,
   zero raised to the zero, etc.  */

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

  result = gfc_constant_result (BT_LOGICAL, op1->ts.kind, &op1->where);
  result->value.logical = !op1->value.logical;
  *resultp = result;

  return ARITH_OK;
}


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

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

  return ARITH_OK;
}


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

  result = gfc_constant_result (BT_LOGICAL, gfc_kind_max (op1, op2),
				&op1->where);
  result->value.logical = op1->value.logical || op2->value.logical;
  *resultp = result;

  return ARITH_OK;
}


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

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

  return ARITH_OK;
}


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

  result = gfc_constant_result (BT_LOGICAL, gfc_kind_max (op1, op2),
				&op1->where);
  result->value.logical = op1->value.logical != op2->value.logical;
  *resultp = result;

  return ARITH_OK;
}


/* Make sure a constant numeric expression is within the range for
   its type and kind.  Note that there's also a gfc_check_range(),
   but that one deals with the intrinsic RANGE function.  */

arith
gfc_range_check (gfc_expr * e)
{
  arith rc;

  switch (e->ts.type)
    {
    case BT_INTEGER:
      rc = gfc_check_integer_range (e->value.integer, e->ts.kind);
      break;

    case BT_REAL:
      rc = gfc_check_real_range (e->value.real, e->ts.kind);
      if (rc == ARITH_UNDERFLOW)
        mpfr_set_ui (e->value.real, 0, GFC_RND_MODE);
      break;

    case BT_COMPLEX:
      rc = gfc_check_real_range (e->value.complex.r, e->ts.kind);
      if (rc == ARITH_UNDERFLOW)
        mpfr_set_ui (e->value.complex.r, 0, GFC_RND_MODE);
      if (rc == ARITH_OK || rc == ARITH_UNDERFLOW)
        {
          rc = gfc_check_real_range (e->value.complex.i, e->ts.kind);
          if (rc == ARITH_UNDERFLOW)
            mpfr_set_ui (e->value.complex.i, 0, GFC_RND_MODE);
        }

      break;

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

  return rc;
}


/* Several of the following routines use the same set of statements to
   check the validity of the result.  Encapsulate the checking here.  */

static arith
check_result (arith rc, gfc_expr * x, gfc_expr * r, gfc_expr ** rp)
{
  arith val = rc;

  if (val == ARITH_UNDERFLOW)
    {
      if (gfc_option.warn_underflow)
	gfc_warning (gfc_arith_error (val), &x->where);
      val = ARITH_OK;
    }

  if (val == ARITH_ASYMMETRIC)
    {
      gfc_warning (gfc_arith_error (val), &x->where);
      val = ARITH_OK;
    }

  if (val != ARITH_OK)
    gfc_free_expr (r);
  else
    *rp = r;