Subversion Repositories scarts

/* 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,

              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;

  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: