gcc/fortran/arith.c

   1 /* Compiler arithmetic
   2    Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006
   3    Free Software Foundation, Inc.
   4    Contributed by Andy Vaught
   5
   6 This file is part of GCC.
   7
   8 GCC is free software; you can redistribute it and/or modify it under
   9 the terms of the GNU General Public License as published by the Free
  10 Software Foundation; either version 2, or (at your option) any later
  11 version.
  12
  13 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
  14 WARRANTY; without even the implied warranty of MERCHANTABILITY or
  15 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  16 for more details.
  17
  18 You should have received a copy of the GNU General Public License
  19 along with GCC; see the file COPYING.  If not, write to the Free
  20 Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA
  21 02110-1301, USA.  */
  22
  23 /* Since target arithmetic must be done on the host, there has to
  24    be some way of evaluating arithmetic expressions as the host
  25    would evaluate them.  We use the GNU MP library and the MPFR
  26    library to do arithmetic, and this file provides the interface.  */
  27
  28 #include "config.h"
  29 #include "system.h"
  30 #include "flags.h"
  31 #include "gfortran.h"
  32 #include "arith.h"
  33
  34 /* MPFR does not have a direct replacement for mpz_set_f() from GMP.
  35    It's easily implemented with a few calls though.  */
  36
  37 void
  38 gfc_mpfr_to_mpz (mpz_t z, mpfr_t x)
  39 {
  40   mp_exp_t e;
  41
  42   e = mpfr_get_z_exp (z, x);
  43   /* MPFR 2.0.1 (included with GMP 4.1) has a bug whereby mpfr_get_z_exp
  44      may set the sign of z incorrectly.  Work around that here.  */
  45   if (mpfr_sgn (x) != mpz_sgn (z))
  46     mpz_neg (z, z);
  47
  48   if (e > 0)
  49     mpz_mul_2exp (z, z, e);
  50   else
  51     mpz_tdiv_q_2exp (z, z, -e);
  52 }
  53
  54
  55 /* Set the model number precision by the requested KIND.  */
  56
  57 void
  58 gfc_set_model_kind (int kind)
  59 {
  60   int index = gfc_validate_kind (BT_REAL, kind, false);
  61   int base2prec;
  62
  63   base2prec = gfc_real_kinds[index].digits;
  64   if (gfc_real_kinds[index].radix != 2)
  65     base2prec *= gfc_real_kinds[index].radix / 2;
  66   mpfr_set_default_prec (base2prec);
  67 }
  68
  69
  70 /* Set the model number precision from mpfr_t x.  */
  71
  72 void
  73 gfc_set_model (mpfr_t x)
  74 {
  75   mpfr_set_default_prec (mpfr_get_prec (x));
  76 }
  77
  78 #if MPFR_VERSION_MAJOR < 2 || (MPFR_VERSION_MAJOR == 2 && MPFR_VERSION_MINOR < 2)
  79 /* Calculate atan2 (y, x)
  80
  81 atan2(y, x) = atan(y/x)                         if x > 0,
  82               sign(y)*(pi - atan(|y/x|))        if x < 0,
  83               0                                 if x = 0 && y == 0,
  84               sign(y)*pi/2                      if x = 0 && y != 0.
  85 */
  86
  87 void
  88 arctangent2 (mpfr_t y, mpfr_t x, mpfr_t result)
  89 {
  90   int i;
  91   mpfr_t t;
  92
  93   gfc_set_model (y);
  94   mpfr_init (t);
  95
  96   i = mpfr_sgn (x);
  97
  98   if (i > 0)
  99     {
 100       mpfr_div (t, y, x, GFC_RND_MODE);
 101       mpfr_atan (result, t, GFC_RND_MODE);
 102     }
 103   else if (i < 0)
 104     {
 105       mpfr_const_pi (result, GFC_RND_MODE);
 106       mpfr_div (t, y, x, GFC_RND_MODE);
 107       mpfr_abs (t, t, GFC_RND_MODE);
 108       mpfr_atan (t, t, GFC_RND_MODE);
 109       mpfr_sub (result, result, t, GFC_RND_MODE);
 110       if (mpfr_sgn (y) < 0)
 111         mpfr_neg (result, result, GFC_RND_MODE);
 112     }
 113   else
 114     {
 115       if (mpfr_sgn (y) == 0)
 116         mpfr_set_ui (result, 0, GFC_RND_MODE);
 117       else
 118         {
 119           mpfr_const_pi (result, GFC_RND_MODE);
 120           mpfr_div_ui (result, result, 2, GFC_RND_MODE);
 121           if (mpfr_sgn (y) < 0)
 122             mpfr_neg (result, result, GFC_RND_MODE);
 123         }
 124     }
 125
 126   mpfr_clear (t);
 127 }
 128 #endif
 129
 130 /* Given an arithmetic error code, return a pointer to a string that
 131    explains the error.  */
 132
 133 static const char *
 134 gfc_arith_error (arith code)
 135 {
 136   const char *p;
 137
 138   switch (code)
 139     {
 140     case ARITH_OK:
 141       p = _("Arithmetic OK at %L");
 142       break;
 143     case ARITH_OVERFLOW:
 144       p = _("Arithmetic overflow at %L");
 145       break;
 146     case ARITH_UNDERFLOW:
 147       p = _("Arithmetic underflow at %L");
 148       break;
 149     case ARITH_NAN:
 150       p = _("Arithmetic NaN at %L");
 151       break;
 152     case ARITH_DIV0:
 153       p = _("Division by zero at %L");
 154       break;
 155     case ARITH_INCOMMENSURATE:
 156       p = _("Array operands are incommensurate at %L");
 157       break;
 158     case ARITH_ASYMMETRIC:
 159       p =
 160         _("Integer outside symmetric range implied by Standard Fortran at %L");
 161       break;
 162     default:
 163       gfc_internal_error ("gfc_arith_error(): Bad error code");
 164     }
 165
 166   return p;
 167 }
 168
 169
 170 /* Get things ready to do math.  */
 171
 172 void
 173 gfc_arith_init_1 (void)
 174 {
 175   gfc_integer_info *int_info;
 176   gfc_real_info *real_info;
 177   mpfr_t a, b, c;
 178   mpz_t r;
 179   int i;
 180
 181   mpfr_set_default_prec (128);
 182   mpfr_init (a);
 183   mpz_init (r);
 184
 185   /* Convert the minimum and maximum values for each kind into their
 186      GNU MP representation.  */
 187   for (int_info = gfc_integer_kinds; int_info->kind != 0; int_info++)
 188     {
 189       /* Huge  */
 190       mpz_set_ui (r, int_info->radix);
 191       mpz_pow_ui (r, r, int_info->digits);
 192
 193       mpz_init (int_info->huge);
 194       mpz_sub_ui (int_info->huge, r, 1);
 195
 196       /* These are the numbers that are actually representable by the
 197          target.  For bases other than two, this needs to be changed.  */
 198       if (int_info->radix != 2)
 199         gfc_internal_error ("Fix min_int calculation");
 200
 201       /* See PRs 13490 and 17912, related to integer ranges.
 202          The pedantic_min_int exists for range checking when a program
 203          is compiled with -pedantic, and reflects the belief that
 204          Standard Fortran requires integers to be symmetrical, i.e.
 205          every negative integer must have a representable positive
 206          absolute value, and vice versa.  */
 207
 208       mpz_init (int_info->pedantic_min_int);
 209       mpz_neg (int_info->pedantic_min_int, int_info->huge);
 210
 211       mpz_init (int_info->min_int);
 212       mpz_sub_ui (int_info->min_int, int_info->pedantic_min_int, 1);
 213
 214       /* Range  */
 215       mpfr_set_z (a, int_info->huge, GFC_RND_MODE);
 216       mpfr_log10 (a, a, GFC_RND_MODE);
 217       mpfr_trunc (a, a);
 218       gfc_mpfr_to_mpz (r, a);
 219       int_info->range = mpz_get_si (r);
 220     }
 221
 222   mpfr_clear (a);
 223
 224   for (real_info = gfc_real_kinds; real_info->kind != 0; real_info++)
 225     {
 226       gfc_set_model_kind (real_info->kind);
 227
 228       mpfr_init (a);
 229       mpfr_init (b);
 230       mpfr_init (c);
 231
 232       /* huge(x) = (1 - b**(-p)) * b**(emax-1) * b  */
 233       /* a = 1 - b**(-p)  */
 234       mpfr_set_ui (a, 1, GFC_RND_MODE);
 235       mpfr_set_ui (b, real_info->radix, GFC_RND_MODE);
 236       mpfr_pow_si (b, b, -real_info->digits, GFC_RND_MODE);
 237       mpfr_sub (a, a, b, GFC_RND_MODE);
 238
 239       /* c = b**(emax-1)  */
 240       mpfr_set_ui (b, real_info->radix, GFC_RND_MODE);
 241       mpfr_pow_ui (c, b, real_info->max_exponent - 1, GFC_RND_MODE);
 242
 243       /* a = a * c = (1 - b**(-p)) * b**(emax-1)  */
 244       mpfr_mul (a, a, c, GFC_RND_MODE);
 245
 246       /* a = (1 - b**(-p)) * b**(emax-1) * b  */
 247       mpfr_mul_ui (a, a, real_info->radix, GFC_RND_MODE);
 248
 249       mpfr_init (real_info->huge);
 250       mpfr_set (real_info->huge, a, GFC_RND_MODE);
 251
 252       /* tiny(x) = b**(emin-1)  */
 253       mpfr_set_ui (b, real_info->radix, GFC_RND_MODE);
 254       mpfr_pow_si (b, b, real_info->min_exponent - 1, GFC_RND_MODE);
 255
 256       mpfr_init (real_info->tiny);
 257       mpfr_set (real_info->tiny, b, GFC_RND_MODE);
 258
 259       /* subnormal (x) = b**(emin - digit)  */
 260       mpfr_set_ui (b, real_info->radix, GFC_RND_MODE);
 261       mpfr_pow_si (b, b, real_info->min_exponent - real_info->digits,
 262                    GFC_RND_MODE);
 263
 264       mpfr_init (real_info->subnormal);
 265       mpfr_set (real_info->subnormal, b, GFC_RND_MODE);
 266
 267       /* epsilon(x) = b**(1-p)  */
 268       mpfr_set_ui (b, real_info->radix, GFC_RND_MODE);
 269       mpfr_pow_si (b, b, 1 - real_info->digits, GFC_RND_MODE);
 270
 271       mpfr_init (real_info->epsilon);
 272       mpfr_set (real_info->epsilon, b, GFC_RND_MODE);
 273
 274       /* range(x) = int(min(log10(huge(x)), -log10(tiny))  */
 275       mpfr_log10 (a, real_info->huge, GFC_RND_MODE);
 276       mpfr_log10 (b, real_info->tiny, GFC_RND_MODE);
 277       mpfr_neg (b, b, GFC_RND_MODE);
 278
 279       /* a = min(a, b)  */
 280       if (mpfr_cmp (a, b) > 0)
 281         mpfr_set (a, b, GFC_RND_MODE);
 282
 283       mpfr_trunc (a, a);
 284       gfc_mpfr_to_mpz (r, a);
 285       real_info->range = mpz_get_si (r);
 286
 287       /* precision(x) = int((p - 1) * log10(b)) + k  */
 288       mpfr_set_ui (a, real_info->radix, GFC_RND_MODE);
 289       mpfr_log10 (a, a, GFC_RND_MODE);
 290
 291       mpfr_mul_ui (a, a, real_info->digits - 1, GFC_RND_MODE);
 292       mpfr_trunc (a, a);
 293       gfc_mpfr_to_mpz (r, a);
 294       real_info->precision = mpz_get_si (r);
 295
 296       /* If the radix is an integral power of 10, add one to the precision.  */
 297       for (i = 10; i <= real_info->radix; i *= 10)
 298         if (i == real_info->radix)
 299           real_info->precision++;
 300
 301       mpfr_clear (a);
 302       mpfr_clear (b);
 303       mpfr_clear (c);
 304     }
 305
 306   mpz_clear (r);
 307 }
 308
 309
 310 /* Clean up, get rid of numeric constants.  */
 311
 312 void
 313 gfc_arith_done_1 (void)
 314 {
 315   gfc_integer_info *ip;
 316   gfc_real_info *rp;
 317
 318   for (ip = gfc_integer_kinds; ip->kind; ip++)
 319     {
 320       mpz_clear (ip->min_int);
 321       mpz_clear (ip->pedantic_min_int);
 322       mpz_clear (ip->huge);
 323     }
 324
 325   for (rp = gfc_real_kinds; rp->kind; rp++)
 326     {
 327       mpfr_clear (rp->epsilon);
 328       mpfr_clear (rp->huge);
 329       mpfr_clear (rp->tiny);
 330       mpfr_clear (rp->subnormal);
 331     }
 332 }
 333
 334
 335 /* Given an integer and a kind, make sure that the integer lies within
 336    the range of the kind.  Returns ARITH_OK, ARITH_ASYMMETRIC or
 337    ARITH_OVERFLOW.  */
 338
 339 arith
 340 gfc_check_integer_range (mpz_t p, int kind)
 341 {
 342   arith result;
 343   int i;
 344
 345   i = gfc_validate_kind (BT_INTEGER, kind, false);
 346   result = ARITH_OK;
 347
 348   if (pedantic)
 349     {
 350       if (mpz_cmp (p, gfc_integer_kinds[i].pedantic_min_int) < 0)
 351         result = ARITH_ASYMMETRIC;
 352     }
 353
 354   if (mpz_cmp (p, gfc_integer_kinds[i].min_int) < 0
 355       || mpz_cmp (p, gfc_integer_kinds[i].huge) > 0)
 356     result = ARITH_OVERFLOW;
 357
 358   return result;
 359 }
 360
 361
 362 /* Given a real and a kind, make sure that the real lies within the
 363    range of the kind.  Returns ARITH_OK, ARITH_OVERFLOW or
 364    ARITH_UNDERFLOW.  */
 365
 366 static arith
 367 gfc_check_real_range (mpfr_t p, int kind)
 368 {
 369   arith retval;
 370   mpfr_t q;
 371   int i;
 372
 373   i = gfc_validate_kind (BT_REAL, kind, false);
 374
 375   gfc_set_model (p);
 376   mpfr_init (q);
 377   mpfr_abs (q, p, GFC_RND_MODE);
 378
 379   if (mpfr_inf_p (p))
 380     {
 381       if (gfc_option.flag_range_check == 0)
 382         retval = ARITH_OK;
 383       else
 384         retval = ARITH_OVERFLOW;
 385     }
 386   else if (mpfr_nan_p (p))
 387     {
 388       if (gfc_option.flag_range_check == 0)
 389         retval = ARITH_OK;
 390       else
 391         retval = ARITH_NAN;
 392     }
 393   else if (mpfr_sgn (q) == 0)
 394     retval = ARITH_OK;
 395   else if (mpfr_cmp (q, gfc_real_kinds[i].huge) > 0)
 396     {
 397       if (gfc_option.flag_range_check == 0)
 398         retval = ARITH_OK;
 399       else
 400         retval = ARITH_OVERFLOW;
 401     }
 402   else if (mpfr_cmp (q, gfc_real_kinds[i].subnormal) < 0)
 403     {
 404       if (gfc_option.flag_range_check == 0)
 405         retval = ARITH_OK;
 406       else
 407         retval = ARITH_UNDERFLOW;
 408     }
 409   else if (mpfr_cmp (q, gfc_real_kinds[i].tiny) < 0)
 410     {
 411 #if MPFR_VERSION_MAJOR < 2 || (MPFR_VERSION_MAJOR == 2 && MPFR_VERSION_MINOR < 2)
 412       /* MPFR operates on a number with a given precision and enormous
 413         exponential range.  To represent subnormal numbers, the exponent is
 414         allowed to become smaller than emin, but always retains the full
 415         precision.  This code resets unused bits to 0 to alleviate
 416         rounding problems.  Note, a future version of MPFR will have a
 417         mpfr_subnormalize() function, which handles this truncation in a
 418         more efficient and robust way.  */
 419
 420       int j, k;
 421       char *bin, *s;
 422       mp_exp_t e;
 423
 424       bin = mpfr_get_str (NULL, &e, gfc_real_kinds[i].radix, 0, q, GMP_RNDN);
 425       k = gfc_real_kinds[i].digits - (gfc_real_kinds[i].min_exponent - e);
 426       for (j = k; j < gfc_real_kinds[i].digits; j++)
 427         bin[j] = '0';
 428       /* Need space for '0.', bin, 'E', and e */
 429       s = (char *) gfc_getmem (strlen(bin) + 10);
 430       sprintf (s, "0.%sE%d", bin, (int) e);
 431       mpfr_set_str (q, s, gfc_real_kinds[i].radix, GMP_RNDN);
 432
 433       gfc_free (s);
 434       gfc_free (bin);
 435 #else
 436       mp_exp_t emin, emax;
 437
 438       /* Save current values of emin and emax.  */
 439       emin = mpfr_get_emin ();
 440       emax = mpfr_get_emax ();
 441
 442       /* Set emin and emax for the current model number.  */
 443       mpfr_set_emin ((mp_exp_t) gfc_real_kinds[i].min_exponent - 1);
 444       mpfr_set_emax ((mp_exp_t) gfc_real_kinds[i].max_exponent - 1);
 445       mpfr_subnormalize (q, 0, GFC_RND_MODE);
 446
 447       /* Reset emin and emax.  */
 448       mpfr_set_emin (emin);
 449       mpfr_set_emax (emax);
 450 #endif
 451
 452       /* Copy sign if needed.  */
 453       if (mpfr_sgn (p) < 0)
 454         mpfr_neg (p, q, GMP_RNDN);
 455       else
 456         mpfr_set (p, q, GMP_RNDN);
 457
 458       retval = ARITH_OK;
 459     }
 460   else
 461     retval = ARITH_OK;
 462
 463   mpfr_clear (q);
 464
 465   return retval;
 466 }
 467
 468
 469 /* Function to return a constant expression node of a given type and kind.  */
 470
 471 gfc_expr *
 472 gfc_constant_result (bt type, int kind, locus * where)
 473 {
 474   gfc_expr *result;
 475
 476   if (!where)
 477     gfc_internal_error
 478       ("gfc_constant_result(): locus 'where' cannot be NULL");
 479
 480   result = gfc_get_expr ();
 481
 482   result->expr_type = EXPR_CONSTANT;
 483   result->ts.type = type;
 484   result->ts.kind = kind;
 485   result->where = *where;
 486
 487   switch (type)
 488     {
 489     case BT_INTEGER:
 490       mpz_init (result->value.integer);
 491       break;
 492
 493     case BT_REAL:
 494       gfc_set_model_kind (kind);
 495       mpfr_init (result->value.real);
 496       break;
 497
 498     case BT_COMPLEX:
 499       gfc_set_model_kind (kind);
 500       mpfr_init (result->value.complex.r);
 501       mpfr_init (result->value.complex.i);
 502       break;
 503
 504     default:
 505       break;
 506     }
 507
 508   return result;
 509 }
 510
 511
 512 /* Low-level arithmetic functions.  All of these subroutines assume
 513    that all operands are of the same type and return an operand of the
 514    same type.  The other thing about these subroutines is that they
 515    can fail in various ways -- overflow, underflow, division by zero,
 516    zero raised to the zero, etc.  */
 517
 518 static arith
 519 gfc_arith_not (gfc_expr * op1, gfc_expr ** resultp)
 520 {
 521   gfc_expr *result;
 522
 523   result = gfc_constant_result (BT_LOGICAL, op1->ts.kind, &op1->where);
 524   result->value.logical = !op1->value.logical;
 525   *resultp = result;
 526
 527   return ARITH_OK;
 528 }
 529
 530
 531 static arith
 532 gfc_arith_and (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
 533 {
 534   gfc_expr *result;
 535
 536   result = gfc_constant_result (BT_LOGICAL, gfc_kind_max (op1, op2),
 537                                 &op1->where);
 538   result->value.logical = op1->value.logical && op2->value.logical;
 539   *resultp = result;
 540
 541   return ARITH_OK;
 542 }
 543
 544
 545 static arith
 546 gfc_arith_or (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
 547 {
 548   gfc_expr *result;
 549
 550   result = gfc_constant_result (BT_LOGICAL, gfc_kind_max (op1, op2),
 551                                 &op1->where);
 552   result->value.logical = op1->value.logical || op2->value.logical;
 553   *resultp = result;
 554
 555   return ARITH_OK;
 556 }
 557
 558
 559 static arith
 560 gfc_arith_eqv (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
 561 {
 562   gfc_expr *result;
 563
 564   result = gfc_constant_result (BT_LOGICAL, gfc_kind_max (op1, op2),
 565                                 &op1->where);
 566   result->value.logical = op1->value.logical == op2->value.logical;
 567   *resultp = result;
 568
 569   return ARITH_OK;
 570 }
 571
 572
 573 static arith
 574 gfc_arith_neqv (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
 575 {
 576   gfc_expr *result;
 577
 578   result = gfc_constant_result (BT_LOGICAL, gfc_kind_max (op1, op2),
 579                                 &op1->where);
 580   result->value.logical = op1->value.logical != op2->value.logical;
 581   *resultp = result;
 582
 583   return ARITH_OK;
 584 }
 585
 586
 587 /* Make sure a constant numeric expression is within the range for
 588    its type and kind.  Note that there's also a gfc_check_range(),
 589    but that one deals with the intrinsic RANGE function.  */
 590
 591 arith
 592 gfc_range_check (gfc_expr * e)
 593 {
 594   arith rc;
 595
 596   switch (e->ts.type)
 597     {
 598     case BT_INTEGER:
 599       rc = gfc_check_integer_range (e->value.integer, e->ts.kind);
 600       break;
 601
 602     case BT_REAL:
 603       rc = gfc_check_real_range (e->value.real, e->ts.kind);
 604       if (rc == ARITH_UNDERFLOW)
 605         mpfr_set_ui (e->value.real, 0, GFC_RND_MODE);
 606       if (rc == ARITH_OVERFLOW)
 607         mpfr_set_inf (e->value.real, mpfr_sgn (e->value.real));
 608       if (rc == ARITH_NAN)
 609         mpfr_set_nan (e->value.real);
 610       break;
 611
 612     case BT_COMPLEX:
 613       rc = gfc_check_real_range (e->value.complex.r, e->ts.kind);
 614       if (rc == ARITH_UNDERFLOW)
 615         mpfr_set_ui (e->value.complex.r, 0, GFC_RND_MODE);
 616       if (rc == ARITH_OVERFLOW)
 617         mpfr_set_inf (e->value.complex.r, mpfr_sgn (e->value.complex.r));
 618       if (rc == ARITH_NAN)
 619         mpfr_set_nan (e->value.complex.r);
 620
 621       rc = gfc_check_real_range (e->value.complex.i, e->ts.kind);
 622       if (rc == ARITH_UNDERFLOW)
 623         mpfr_set_ui (e->value.complex.i, 0, GFC_RND_MODE);
 624       if (rc == ARITH_OVERFLOW)
 625         mpfr_set_inf (e->value.complex.i, mpfr_sgn (e->value.complex.i));
 626       if (rc == ARITH_NAN)
 627         mpfr_set_nan (e->value.complex.i);
 628       break;
 629
 630     default:
 631       gfc_internal_error ("gfc_range_check(): Bad type");
 632     }
 633
 634   return rc;
 635 }
 636
 637
 638 /* Several of the following routines use the same set of statements to
 639    check the validity of the result.  Encapsulate the checking here.  */
 640
 641 static arith
 642 check_result (arith rc, gfc_expr * x, gfc_expr * r, gfc_expr ** rp)
 643 {
 644   arith val = rc;
 645
 646   if (val == ARITH_UNDERFLOW)
 647     {
 648       if (gfc_option.warn_underflow)
 649         gfc_warning (gfc_arith_error (val), &x->where);
 650       val = ARITH_OK;
 651     }
 652
 653   if (val == ARITH_ASYMMETRIC)
 654     {
 655       gfc_warning (gfc_arith_error (val), &x->where);
 656       val = ARITH_OK;
 657     }
 658
 659   if (val != ARITH_OK)
 660     gfc_free_expr (r);
 661   else
 662     *rp = r;
 663
 664   return val;
 665 }
 666
 667
 668 /* It may seem silly to have a subroutine that actually computes the
 669    unary plus of a constant, but it prevents us from making exceptions
 670    in the code elsewhere.  */
 671
 672 static arith
 673 gfc_arith_uplus (gfc_expr * op1, gfc_expr ** resultp)
 674 {
 675   *resultp = gfc_copy_expr (op1);
 676   return ARITH_OK;
 677 }
 678
 679
 680 static arith
 681 gfc_arith_uminus (gfc_expr * op1, gfc_expr ** resultp)
 682 {
 683   gfc_expr *result;
 684   arith rc;
 685
 686   result = gfc_constant_result (op1->ts.type, op1->ts.kind, &op1->where);
 687
 688   switch (op1->ts.type)
 689     {
 690     case BT_INTEGER:
 691       mpz_neg (result->value.integer, op1->value.integer);
 692       break;
 693
 694     case BT_REAL:
 695       mpfr_neg (result->value.real, op1->value.real, GFC_RND_MODE);
 696       break;
 697
 698     case BT_COMPLEX:
 699       mpfr_neg (result->value.complex.r, op1->value.complex.r, GFC_RND_MODE);
 700       mpfr_neg (result->value.complex.i, op1->value.complex.i, GFC_RND_MODE);
 701       break;
 702
 703     default:
 704       gfc_internal_error ("gfc_arith_uminus(): Bad basic type");
 705     }
 706
 707   rc = gfc_range_check (result);
 708
 709   return check_result (rc, op1, result, resultp);
 710 }
 711
 712
 713 static arith
 714 gfc_arith_plus (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
 715 {
 716   gfc_expr *result;
 717   arith rc;
 718
 719   result = gfc_constant_result (op1->ts.type, op1->ts.kind, &op1->where);
 720
 721   switch (op1->ts.type)
 722     {
 723     case BT_INTEGER:
 724       mpz_add (result->value.integer, op1->value.integer, op2->value.integer);
 725       break;
 726
 727     case BT_REAL:
 728       mpfr_add (result->value.real, op1->value.real, op2->value.real,
 729                GFC_RND_MODE);
 730       break;
 731
 732     case BT_COMPLEX:
 733       mpfr_add (result->value.complex.r, op1->value.complex.r,
 734                op2->value.complex.r, GFC_RND_MODE);
 735
 736       mpfr_add (result->value.complex.i, op1->value.complex.i,
 737                op2->value.complex.i, GFC_RND_MODE);
 738       break;
 739
 740     default:
 741       gfc_internal_error ("gfc_arith_plus(): Bad basic type");
 742     }
 743
 744   rc = gfc_range_check (result);
 745
 746   return check_result (rc, op1, result, resultp);
 747 }
 748
 749
 750 static arith
 751 gfc_arith_minus (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
 752 {
 753   gfc_expr *result;
 754   arith rc;
 755
 756   result = gfc_constant_result (op1->ts.type, op1->ts.kind, &op1->where);
 757
 758   switch (op1->ts.type)
 759     {
 760     case BT_INTEGER:
 761       mpz_sub (result->value.integer, op1->value.integer, op2->value.integer);
 762       break;
 763
 764     case BT_REAL:
 765       mpfr_sub (result->value.real, op1->value.real, op2->value.real,
 766                 GFC_RND_MODE);
 767       break;
 768
 769     case BT_COMPLEX:
 770       mpfr_sub (result->value.complex.r, op1->value.complex.r,
 771                op2->value.complex.r, GFC_RND_MODE);
 772
 773       mpfr_sub (result->value.complex.i, op1->value.complex.i,
 774                op2->value.complex.i, GFC_RND_MODE);
 775       break;
 776
 777     default:
 778       gfc_internal_error ("gfc_arith_minus(): Bad basic type");
 779     }
 780
 781   rc = gfc_range_check (result);
 782
 783   return check_result (rc, op1, result, resultp);
 784 }
 785
 786
 787 static arith
 788 gfc_arith_times (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
 789 {
 790   gfc_expr *result;
 791   mpfr_t x, y;
 792   arith rc;
 793
 794   result = gfc_constant_result (op1->ts.type, op1->ts.kind, &op1->where);
 795
 796   switch (op1->ts.type)
 797     {
 798     case BT_INTEGER:
 799       mpz_mul (result->value.integer, op1->value.integer, op2->value.integer);
 800       break;
 801
 802     case BT_REAL:
 803       mpfr_mul (result->value.real, op1->value.real, op2->value.real,
 804                GFC_RND_MODE);
 805       break;
 806
 807     case BT_COMPLEX:
 808       gfc_set_model (op1->value.complex.r);
 809       mpfr_init (x);
 810       mpfr_init (y);
 811
 812       mpfr_mul (x, op1->value.complex.r, op2->value.complex.r, GFC_RND_MODE);
 813       mpfr_mul (y, op1->value.complex.i, op2->value.complex.i, GFC_RND_MODE);
 814       mpfr_sub (result->value.complex.r, x, y, GFC_RND_MODE);
 815
 816       mpfr_mul (x, op1->value.complex.r, op2->value.complex.i, GFC_RND_MODE);
 817       mpfr_mul (y, op1->value.complex.i, op2->value.complex.r, GFC_RND_MODE);
 818       mpfr_add (result->value.complex.i, x, y, GFC_RND_MODE);
 819
 820       mpfr_clear (x);
 821       mpfr_clear (y);
 822       break;
 823
 824     default:
 825       gfc_internal_error ("gfc_arith_times(): Bad basic type");
 826     }
 827
 828   rc = gfc_range_check (result);
 829
 830   return check_result (rc, op1, result, resultp);
 831 }
 832
 833
 834 static arith
 835 gfc_arith_divide (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
 836 {
 837   gfc_expr *result;
 838   mpfr_t x, y, div;
 839   arith rc;
 840
 841   rc = ARITH_OK;
 842
 843   result = gfc_constant_result (op1->ts.type, op1->ts.kind, &op1->where);
 844
 845   switch (op1->ts.type)
 846     {
 847     case BT_INTEGER:
 848       if (mpz_sgn (op2->value.integer) == 0)
 849         {
 850           rc = ARITH_DIV0;
 851           break;
 852         }
 853
 854       mpz_tdiv_q (result->value.integer, op1->value.integer,
 855                   op2->value.integer);
 856       break;
 857
 858     case BT_REAL:
 859       if (mpfr_sgn (op2->value.real) == 0
 860           && gfc_option.flag_range_check == 1)
 861         {
 862           rc = ARITH_DIV0;
 863           break;
 864         }
 865
 866       mpfr_div (result->value.real, op1->value.real, op2->value.real,
 867                GFC_RND_MODE);
 868       break;
 869
 870     case BT_COMPLEX:
 871       if (mpfr_sgn (op2->value.complex.r) == 0
 872           && mpfr_sgn (op2->value.complex.i) == 0
 873           && gfc_option.flag_range_check == 1)
 874         {
 875           rc = ARITH_DIV0;
 876           break;
 877         }
 878
 879       gfc_set_model (op1->value.complex.r);
 880       mpfr_init (x);
 881       mpfr_init (y);
 882       mpfr_init (div);
 883
 884       mpfr_mul (x, op2->value.complex.r, op2->value.complex.r, GFC_RND_MODE);
 885       mpfr_mul (y, op2->value.complex.i, op2->value.complex.i, GFC_RND_MODE);
 886       mpfr_add (div, x, y, GFC_RND_MODE);
 887
 888       mpfr_mul (x, op1->value.complex.r, op2->value.complex.r, GFC_RND_MODE);
 889       mpfr_mul (y, op1->value.complex.i, op2->value.complex.i, GFC_RND_MODE);
 890       mpfr_add (result->value.complex.r, x, y, GFC_RND_MODE);
 891       mpfr_div (result->value.complex.r, result->value.complex.r, div,
 892                 GFC_RND_MODE);
 893
 894       mpfr_mul (x, op1->value.complex.i, op2->value.complex.r, GFC_RND_MODE);
 895       mpfr_mul (y, op1->value.complex.r, op2->value.complex.i, GFC_RND_MODE);
 896       mpfr_sub (result->value.complex.i, x, y, GFC_RND_MODE);
 897       mpfr_div (result->value.complex.i, result->value.complex.i, div,
 898                 GFC_RND_MODE);
 899
 900       mpfr_clear (x);
 901       mpfr_clear (y);
 902       mpfr_clear (div);
 903       break;
 904
 905     default:
 906       gfc_internal_error ("gfc_arith_divide(): Bad basic type");
 907     }
 908
 909   if (rc == ARITH_OK)
 910     rc = gfc_range_check (result);
 911
 912   return check_result (rc, op1, result, resultp);
 913 }
 914
 915
 916 /* Compute the reciprocal of a complex number (guaranteed nonzero).  */
 917
 918 static void
 919 complex_reciprocal (gfc_expr * op)
 920 {
 921   mpfr_t mod, a, re, im;
 922
 923   gfc_set_model (op->value.complex.r);
 924   mpfr_init (mod);
 925   mpfr_init (a);
 926   mpfr_init (re);
 927   mpfr_init (im);
 928
 929   mpfr_mul (mod, op->value.complex.r, op->value.complex.r, GFC_RND_MODE);
 930   mpfr_mul (a, op->value.complex.i, op->value.complex.i, GFC_RND_MODE);
 931   mpfr_add (mod, mod, a, GFC_RND_MODE);
 932
 933   mpfr_div (re, op->value.complex.r, mod, GFC_RND_MODE);
 934
 935   mpfr_neg (im, op->value.complex.i, GFC_RND_MODE);
 936   mpfr_div (im, im, mod, GFC_RND_MODE);
 937
 938   mpfr_set (op->value.complex.r, re, GFC_RND_MODE);
 939   mpfr_set (op->value.complex.i, im, GFC_RND_MODE);
 940
 941   mpfr_clear (re);
 942   mpfr_clear (im);
 943   mpfr_clear (mod);
 944   mpfr_clear (a);
 945 }
 946
 947
 948 /* Raise a complex number to positive power.  */
 949
 950 static void
 951 complex_pow_ui (gfc_expr * base, int power, gfc_expr * result)
 952 {
 953   mpfr_t re, im, a;
 954
 955   gfc_set_model (base->value.complex.r);
 956   mpfr_init (re);
 957   mpfr_init (im);
 958   mpfr_init (a);
 959
 960   mpfr_set_ui (result->value.complex.r, 1, GFC_RND_MODE);
 961   mpfr_set_ui (result->value.complex.i, 0, GFC_RND_MODE);
 962
 963   for (; power > 0; power--)
 964     {
 965       mpfr_mul (re, base->value.complex.r, result->value.complex.r,
 966                 GFC_RND_MODE);
 967       mpfr_mul (a, base->value.complex.i, result->value.complex.i,
 968                 GFC_RND_MODE);
 969       mpfr_sub (re, re, a, GFC_RND_MODE);
 970
 971       mpfr_mul (im, base->value.complex.r, result->value.complex.i,
 972                 GFC_RND_MODE);
 973       mpfr_mul (a, base->value.complex.i, result->value.complex.r,
 974                 GFC_RND_MODE);
 975       mpfr_add (im, im, a, GFC_RND_MODE);
 976
 977       mpfr_set (result->value.complex.r, re, GFC_RND_MODE);
 978       mpfr_set (result->value.complex.i, im, GFC_RND_MODE);
 979     }
 980
 981   mpfr_clear (re);
 982   mpfr_clear (im);
 983   mpfr_clear (a);
 984 }
 985
 986
 987 /* Raise a number to an integer power.  */
 988
 989 static arith
 990 gfc_arith_power (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
 991 {
 992   int power, apower;
 993   gfc_expr *result;
 994   mpz_t unity_z;
 995   mpfr_t unity_f;
 996   arith rc;
 997
 998   rc = ARITH_OK;
 999
1000   if (gfc_extract_int (op2, &power) != NULL)
1001     gfc_internal_error ("gfc_arith_power(): Bad exponent");
1002
1003   result = gfc_constant_result (op1->ts.type, op1->ts.kind, &op1->where);
1004
1005   if (power == 0)
1006     {
1007       /* Handle something to the zeroth power.  Since we're dealing
1008          with integral exponents, there is no ambiguity in the
1009          limiting procedure used to determine the value of 0**0.  */
1010       switch (op1->ts.type)
1011         {
1012         case BT_INTEGER:
1013           mpz_set_ui (result->value.integer, 1);
1014           break;
1015
1016         case BT_REAL:
1017           mpfr_set_ui (result->value.real, 1, GFC_RND_MODE);
1018           break;
1019
1020         case BT_COMPLEX:
1021           mpfr_set_ui (result->value.complex.r, 1, GFC_RND_MODE);
1022           mpfr_set_ui (result->value.complex.i, 0, GFC_RND_MODE);
1023           break;
1024
1025         default:
1026           gfc_internal_error ("gfc_arith_power(): Bad base");
1027         }
1028     }
1029   else
1030     {
1031       apower = power;
1032       if (power < 0)
1033         apower = -power;
1034
1035       switch (op1->ts.type)
1036         {
1037         case BT_INTEGER:
1038           mpz_pow_ui (result->value.integer, op1->value.integer, apower);
1039
1040           if (power < 0)
1041             {
1042               mpz_init_set_ui (unity_z, 1);
1043               mpz_tdiv_q (result->value.integer, unity_z,
1044                           result->value.integer);
1045               mpz_clear (unity_z);
1046             }
1047           break;
1048
1049         case BT_REAL:
1050           mpfr_pow_ui (result->value.real, op1->value.real, apower,
1051                        GFC_RND_MODE);
1052
1053           if (power < 0)
1054             {
1055               gfc_set_model (op1->value.real);
1056               mpfr_init (unity_f);
1057               mpfr_set_ui (unity_f, 1, GFC_RND_MODE);
1058               mpfr_div (result->value.real, unity_f, result->value.real,
1059                         GFC_RND_MODE);
1060               mpfr_clear (unity_f);
1061             }
1062           break;
1063
1064         case BT_COMPLEX:
1065           complex_pow_ui (op1, apower, result);
1066           if (power < 0)
1067             complex_reciprocal (result);
1068           break;
1069
1070         default:
1071           break;
1072         }
1073     }
1074
1075   if (rc == ARITH_OK)
1076     rc = gfc_range_check (result);
1077
1078   return check_result (rc, op1, result, resultp);
1079 }
1080
1081
1082 /* Concatenate two string constants.  */
1083
1084 static arith
1085 gfc_arith_concat (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
1086 {
1087   gfc_expr *result;
1088   int len;
1089
1090   result = gfc_constant_result (BT_CHARACTER, gfc_default_character_kind,
1091                                 &op1->where);
1092
1093   len = op1->value.character.length + op2->value.character.length;
1094
1095   result->value.character.string = gfc_getmem (len + 1);
1096   result->value.character.length = len;
1097
1098   memcpy (result->value.character.string, op1->value.character.string,
1099           op1->value.character.length);
1100
1101   memcpy (result->value.character.string + op1->value.character.length,
1102           op2->value.character.string, op2->value.character.length);
1103
1104   result->value.character.string[len] = '\0';
1105
1106   *resultp = result;
1107
1108   return ARITH_OK;
1109 }
1110
1111
1112 /* Comparison operators.  Assumes that the two expression nodes
1113    contain two constants of the same type.  */
1114
1115 int
1116 gfc_compare_expr (gfc_expr * op1, gfc_expr * op2)
1117 {
1118   int rc;
1119
1120   switch (op1->ts.type)
1121     {
1122     case BT_INTEGER:
1123       rc = mpz_cmp (op1->value.integer, op2->value.integer);
1124       break;
1125
1126     case BT_REAL:
1127       rc = mpfr_cmp (op1->value.real, op2->value.real);
1128       break;
1129
1130     case BT_CHARACTER:
1131       rc = gfc_compare_string (op1, op2, NULL);
1132       break;
1133
1134     case BT_LOGICAL:
1135       rc = ((!op1->value.logical && op2->value.logical)
1136             || (op1->value.logical && !op2->value.logical));
1137       break;
1138
1139     default:
1140       gfc_internal_error ("gfc_compare_expr(): Bad basic type");
1141     }
1142
1143   return rc;
1144 }
1145
1146
1147 /* Compare a pair of complex numbers.  Naturally, this is only for
1148    equality and nonequality.  */
1149
1150 static int
1151 compare_complex (gfc_expr * op1, gfc_expr * op2)
1152 {
1153   return (mpfr_cmp (op1->value.complex.r, op2->value.complex.r) == 0
1154           && mpfr_cmp (op1->value.complex.i, op2->value.complex.i) == 0);
1155 }
1156
1157
1158 /* Given two constant strings and the inverse collating sequence, compare the
1159    strings.  We return -1 for a < b, 0 for a == b and 1 for a > b.  If the
1160    xcoll_table is NULL, we use the processor's default collating sequence.  */
1161
1162 int
1163 gfc_compare_string (gfc_expr * a, gfc_expr * b, const int * xcoll_table)
1164 {
1165   int len, alen, blen, i, ac, bc;
1166
1167   alen = a->value.character.length;
1168   blen = b->value.character.length;
1169
1170   len = (alen > blen) ? alen : blen;
1171
1172   for (i = 0; i < len; i++)
1173     {
1174       /* We cast to unsigned char because default char, if it is signed,
1175          would lead to ac < 0 for string[i] > 127.  */
1176       ac = (unsigned char) ((i < alen) ? a->value.character.string[i] : ' ');
1177       bc = (unsigned char) ((i < blen) ? b->value.character.string[i] : ' ');
1178
1179       if (xcoll_table != NULL)
1180         {
1181           ac = xcoll_table[ac];
1182           bc = xcoll_table[bc];
1183         }
1184
1185       if (ac < bc)
1186         return -1;
1187       if (ac > bc)
1188         return 1;
1189     }
1190
1191   /* Strings are equal */
1192
1193   return 0;
1194 }
1195
1196
1197 /* Specific comparison subroutines.  */
1198
1199 static arith
1200 gfc_arith_eq (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
1201 {
1202   gfc_expr *result;
1203
1204   result = gfc_constant_result (BT_LOGICAL, gfc_default_logical_kind,
1205                                 &op1->where);
1206   result->value.logical = (op1->ts.type == BT_COMPLEX) ?
1207     compare_complex (op1, op2) : (gfc_compare_expr (op1, op2) == 0);
1208
1209   *resultp = result;
1210   return ARITH_OK;
1211 }
1212
1213
1214 static arith
1215 gfc_arith_ne (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
1216 {
1217   gfc_expr *result;
1218
1219   result = gfc_constant_result (BT_LOGICAL, gfc_default_logical_kind,
1220                                 &op1->where);
1221   result->value.logical = (op1->ts.type == BT_COMPLEX) ?
1222     !compare_complex (op1, op2) : (gfc_compare_expr (op1, op2) != 0);
1223
1224   *resultp = result;
1225   return ARITH_OK;
1226 }
1227
1228
1229 static arith
1230 gfc_arith_gt (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
1231 {
1232   gfc_expr *result;
1233
1234   result = gfc_constant_result (BT_LOGICAL, gfc_default_logical_kind,
1235                                 &op1->where);
1236   result->value.logical = (gfc_compare_expr (op1, op2) > 0);
1237   *resultp = result;
1238
1239   return ARITH_OK;
1240 }
1241
1242
1243 static arith
1244 gfc_arith_ge (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
1245 {
1246   gfc_expr *result;
1247
1248   result = gfc_constant_result (BT_LOGICAL, gfc_default_logical_kind,
1249                                 &op1->where);
1250   result->value.logical = (gfc_compare_expr (op1, op2) >= 0);
1251   *resultp = result;
1252
1253   return ARITH_OK;
1254 }
1255
1256
1257 static arith
1258 gfc_arith_lt (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
1259 {
1260   gfc_expr *result;
1261
1262   result = gfc_constant_result (BT_LOGICAL, gfc_default_logical_kind,
1263                                 &op1->where);
1264   result->value.logical = (gfc_compare_expr (op1, op2) < 0);
1265   *resultp = result;
1266
1267   return ARITH_OK;
1268 }
1269
1270
1271 static arith
1272 gfc_arith_le (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
1273 {
1274   gfc_expr *result;
1275
1276   result = gfc_constant_result (BT_LOGICAL, gfc_default_logical_kind,
1277                                 &op1->where);
1278   result->value.logical = (gfc_compare_expr (op1, op2) <= 0);
1279   *resultp = result;
1280
1281   return ARITH_OK;
1282 }
1283
1284
1285 static arith
1286 reduce_unary (arith (*eval) (gfc_expr *, gfc_expr **), gfc_expr * op,
1287               gfc_expr ** result)
1288 {
1289   gfc_constructor *c, *head;
1290   gfc_expr *r;
1291   arith rc;
1292
1293   if (op->expr_type == EXPR_CONSTANT)
1294     return eval (op, result);
1295
1296   rc = ARITH_OK;
1297   head = gfc_copy_constructor (op->value.constructor);
1298
1299   for (c = head; c; c = c->next)
1300     {
1301       rc = eval (c->expr, &r);
1302       if (rc != ARITH_OK)
1303         break;
1304
1305       gfc_replace_expr (c->expr, r);
1306     }
1307
1308   if (rc != ARITH_OK)
1309     gfc_free_constructor (head);
1310   else
1311     {
1312       r = gfc_get_expr ();
1313       r->expr_type = EXPR_ARRAY;
1314       r->value.constructor = head;
1315       r->shape = gfc_copy_shape (op->shape, op->rank);
1316
1317       r->ts = head->expr->ts;
1318       r->where = op->where;
1319       r->rank = op->rank;
1320
1321       *result = r;
1322     }
1323
1324   return rc;
1325 }
1326
1327
1328 static arith
1329 reduce_binary_ac (arith (*eval) (gfc_expr *, gfc_expr *, gfc_expr **),
1330                   gfc_expr * op1, gfc_expr * op2,
1331                   gfc_expr ** result)
1332 {
1333   gfc_constructor *c, *head;
1334   gfc_expr *r;
1335   arith rc;
1336
1337   head = gfc_copy_constructor (op1->value.constructor);
1338   rc = ARITH_OK;
1339
1340   for (c = head; c; c = c->next)
1341     {
1342       rc = eval (c->expr, op2, &r);
1343       if (rc != ARITH_OK)
1344         break;
1345
1346       gfc_replace_expr (c->expr, r);
1347     }
1348
1349   if (rc != ARITH_OK)
1350     gfc_free_constructor (head);
1351   else
1352     {
1353       r = gfc_get_expr ();
1354       r->expr_type = EXPR_ARRAY;
1355       r->value.constructor = head;
1356       r->shape = gfc_copy_shape (op1->shape, op1->rank);
1357
1358       r->ts = head->expr->ts;
1359       r->where = op1->where;
1360       r->rank = op1->rank;
1361
1362       *result = r;
1363     }
1364
1365   return rc;
1366 }
1367
1368
1369 static arith
1370 reduce_binary_ca (arith (*eval) (gfc_expr *, gfc_expr *, gfc_expr **),
1371                   gfc_expr * op1, gfc_expr * op2,
1372                   gfc_expr ** result)
1373 {
1374   gfc_constructor *c, *head;
1375   gfc_expr *r;
1376   arith rc;
1377
1378   head = gfc_copy_constructor (op2->value.constructor);
1379   rc = ARITH_OK;
1380
1381   for (c = head; c; c = c->next)
1382     {
1383       rc = eval (op1, c->expr, &r);
1384       if (rc != ARITH_OK)
1385         break;
1386
1387       gfc_replace_expr (c->expr, r);
1388     }
1389
1390   if (rc != ARITH_OK)
1391     gfc_free_constructor (head);
1392   else
1393     {
1394       r = gfc_get_expr ();
1395       r->expr_type = EXPR_ARRAY;
1396       r->value.constructor = head;
1397       r->shape = gfc_copy_shape (op2->shape, op2->rank);
1398
1399       r->ts = head->expr->ts;
1400       r->where = op2->where;
1401       r->rank = op2->rank;
1402
1403       *result = r;
1404     }
1405
1406   return rc;
1407 }
1408
1409
1410 static arith
1411 reduce_binary_aa (arith (*eval) (gfc_expr *, gfc_expr *, gfc_expr **),
1412                   gfc_expr * op1, gfc_expr * op2,
1413                   gfc_expr ** result)
1414 {
1415   gfc_constructor *c, *d, *head;
1416   gfc_expr *r;
1417   arith rc;
1418
1419   head = gfc_copy_constructor (op1->value.constructor);
1420
1421   rc = ARITH_OK;
1422   d = op2->value.constructor;
1423
1424   if (gfc_check_conformance ("Elemental binary operation", op1, op2)
1425       != SUCCESS)
1426     rc = ARITH_INCOMMENSURATE;
1427   else
1428     {
1429
1430       for (c = head; c; c = c->next, d = d->next)
1431         {
1432           if (d == NULL)
1433             {
1434               rc = ARITH_INCOMMENSURATE;
1435               break;
1436             }
1437
1438           rc = eval (c->expr, d->expr, &r);
1439           if (rc != ARITH_OK)
1440             break;
1441
1442           gfc_replace_expr (c->expr, r);
1443         }
1444
1445       if (d != NULL)
1446         rc = ARITH_INCOMMENSURATE;
1447     }
1448
1449   if (rc != ARITH_OK)
1450     gfc_free_constructor (head);
1451   else
1452     {
1453       r = gfc_get_expr ();
1454       r->expr_type = EXPR_ARRAY;
1455       r->value.constructor = head;
1456       r->shape = gfc_copy_shape (op1->shape, op1->rank);
1457
1458       r->ts = head->expr->ts;
1459       r->where = op1->where;
1460       r->rank = op1->rank;
1461
1462       *result = r;
1463     }
1464
1465   return rc;
1466 }
1467
1468
1469 static arith
1470 reduce_binary (arith (*eval) (gfc_expr *, gfc_expr *, gfc_expr **),
1471                gfc_expr * op1, gfc_expr * op2,
1472                gfc_expr ** result)
1473 {
1474   if (op1->expr_type == EXPR_CONSTANT && op2->expr_type == EXPR_CONSTANT)
1475     return eval (op1, op2, result);
1476
1477   if (op1->expr_type == EXPR_CONSTANT && op2->expr_type == EXPR_ARRAY)
1478     return reduce_binary_ca (eval, op1, op2, result);
1479
1480   if (op1->expr_type == EXPR_ARRAY && op2->expr_type == EXPR_CONSTANT)
1481     return reduce_binary_ac (eval, op1, op2, result);
1482
1483   return reduce_binary_aa (eval, op1, op2, result);
1484 }
1485
1486
1487 typedef union
1488 {
1489   arith (*f2)(gfc_expr *, gfc_expr **);
1490   arith (*f3)(gfc_expr *, gfc_expr *, gfc_expr **);
1491 }
1492 eval_f;
1493
1494 /* High level arithmetic subroutines.  These subroutines go into
1495    eval_intrinsic(), which can do one of several things to its
1496    operands.  If the operands are incompatible with the intrinsic
1497    operation, we return a node pointing to the operands and hope that
1498    an operator interface is found during resolution.
1499
1500    If the operands are compatible and are constants, then we try doing
1501    the arithmetic.  We also handle the cases where either or both
1502    operands are array constructors.  */
1503
1504 static gfc_expr *
1505 eval_intrinsic (gfc_intrinsic_op operator,
1506                 eval_f eval, gfc_expr * op1, gfc_expr * op2)
1507 {
1508   gfc_expr temp, *result;
1509   int unary;
1510   arith rc;
1511
1512   gfc_clear_ts (&temp.ts);
1513
1514   switch (operator)
1515     {
1516     /* Logical unary  */
1517     case INTRINSIC_NOT:
1518       if (op1->ts.type != BT_LOGICAL)
1519         goto runtime;
1520
1521       temp.ts.type = BT_LOGICAL;
1522       temp.ts.kind = gfc_default_logical_kind;
1523
1524       unary = 1;
1525       break;
1526
1527     /* Logical binary operators  */
1528     case INTRINSIC_OR:
1529     case INTRINSIC_AND:
1530     case INTRINSIC_NEQV:
1531     case INTRINSIC_EQV:
1532       if (op1->ts.type != BT_LOGICAL || op2->ts.type != BT_LOGICAL)
1533         goto runtime;
1534
1535       temp.ts.type = BT_LOGICAL;
1536       temp.ts.kind = gfc_default_logical_kind;
1537
1538       unary = 0;
1539       break;
1540
1541     /* Numeric unary  */
1542     case INTRINSIC_UPLUS:
1543     case INTRINSIC_UMINUS:
1544       if (!gfc_numeric_ts (&op1->ts))
1545         goto runtime;
1546
1547       temp.ts = op1->ts;
1548
1549       unary = 1;
1550       break;
1551
1552     case INTRINSIC_PARENTHESES:
1553       temp.ts = op1->ts;
1554
1555       unary = 1;
1556       break;
1557
1558     /* Additional restrictions for ordering relations.  */
1559     case INTRINSIC_GE:
1560     case INTRINSIC_LT:
1561     case INTRINSIC_LE:
1562     case INTRINSIC_GT:
1563       if (op1->ts.type == BT_COMPLEX || op2->ts.type == BT_COMPLEX)
1564         {
1565           temp.ts.type = BT_LOGICAL;
1566           temp.ts.kind = gfc_default_logical_kind;
1567           goto runtime;
1568         }
1569
1570     /* Fall through  */
1571     case INTRINSIC_EQ:
1572     case INTRINSIC_NE:
1573       if (op1->ts.type == BT_CHARACTER && op2->ts.type == BT_CHARACTER)
1574         {
1575           unary = 0;
1576           temp.ts.type = BT_LOGICAL;
1577           temp.ts.kind = gfc_default_logical_kind;
1578           break;
1579         }
1580
1581     /* Fall through  */
1582     /* Numeric binary  */
1583     case INTRINSIC_PLUS:
1584     case INTRINSIC_MINUS:
1585     case INTRINSIC_TIMES:
1586     case INTRINSIC_DIVIDE:
1587     case INTRINSIC_POWER:
1588       if (!gfc_numeric_ts (&op1->ts) || !gfc_numeric_ts (&op2->ts))
1589         goto runtime;
1590
1591       /* Insert any necessary type conversions to make the operands
1592          compatible.  */
1593
1594       temp.expr_type = EXPR_OP;
1595       gfc_clear_ts (&temp.ts);
1596       temp.value.op.operator = operator;
1597
1598       temp.value.op.op1 = op1;
1599       temp.value.op.op2 = op2;
1600
1601       gfc_type_convert_binary (&temp);
1602
1603       if (operator == INTRINSIC_EQ || operator == INTRINSIC_NE
1604           || operator == INTRINSIC_GE || operator == INTRINSIC_GT
1605           || operator == INTRINSIC_LE || operator == INTRINSIC_LT)
1606         {
1607           temp.ts.type = BT_LOGICAL;
1608           temp.ts.kind = gfc_default_logical_kind;
1609         }
1610
1611       unary = 0;
1612       break;
1613
1614     /* Character binary  */
1615     case INTRINSIC_CONCAT:
1616       if (op1->ts.type != BT_CHARACTER || op2->ts.type != BT_CHARACTER)
1617         goto runtime;
1618
1619       temp.ts.type = BT_CHARACTER;
1620       temp.ts.kind = gfc_default_character_kind;
1621
1622       unary = 0;
1623       break;
1624
1625     case INTRINSIC_USER:
1626       goto runtime;
1627
1628     default:
1629       gfc_internal_error ("eval_intrinsic(): Bad operator");
1630     }
1631
1632   /* Try to combine the operators.  */
1633   if (operator == INTRINSIC_POWER && op2->ts.type != BT_INTEGER)
1634     goto runtime;
1635
1636   if (op1->from_H
1637       || (op1->expr_type != EXPR_CONSTANT
1638           && (op1->expr_type != EXPR_ARRAY
1639               || !gfc_is_constant_expr (op1)
1640               || !gfc_expanded_ac (op1))))
1641     goto runtime;
1642
1643   if (op2 != NULL
1644       && (op2->from_H
1645           || (op2->expr_type != EXPR_CONSTANT
1646               && (op2->expr_type != EXPR_ARRAY
1647               || !gfc_is_constant_expr (op2)
1648               || !gfc_expanded_ac (op2)))))
1649     goto runtime;
1650
1651   if (unary)
1652     rc = reduce_unary (eval.f2, op1, &result);
1653   else
1654     rc = reduce_binary (eval.f3, op1, op2, &result);
1655
1656   if (rc != ARITH_OK)
1657     { /* Something went wrong.  */
1658       gfc_error (gfc_arith_error (rc), &op1->where);
1659       return NULL;
1660     }
1661
1662   gfc_free_expr (op1);
1663   gfc_free_expr (op2);
1664   return result;
1665
1666 runtime:
1667   /* Create a run-time expression.  */
1668   result = gfc_get_expr ();
1669   result->ts = temp.ts;
1670
1671   result->expr_type = EXPR_OP;
1672   result->value.op.operator = operator;
1673
1674   result->value.op.op1 = op1;
1675   result->value.op.op2 = op2;
1676
1677   result->where = op1->where;
1678
1679   return result;
1680 }
1681
1682
1683 /* Modify type of expression for zero size array.  */
1684
1685 static gfc_expr *
1686 eval_type_intrinsic0 (gfc_intrinsic_op operator, gfc_expr * op)
1687 {
1688   if (op == NULL)
1689     gfc_internal_error ("eval_type_intrinsic0(): op NULL");
1690
1691   switch (operator)
1692     {
1693     case INTRINSIC_GE:
1694     case INTRINSIC_LT:
1695     case INTRINSIC_LE:
1696     case INTRINSIC_GT:
1697     case INTRINSIC_EQ:
1698     case INTRINSIC_NE:
1699       op->ts.type = BT_LOGICAL;
1700       op->ts.kind = gfc_default_logical_kind;
1701       break;
1702
1703     default:
1704       break;
1705     }
1706
1707   return op;
1708 }
1709
1710
1711 /* Return nonzero if the expression is a zero size array.  */
1712
1713 static int
1714 gfc_zero_size_array (gfc_expr * e)
1715 {
1716   if (e->expr_type != EXPR_ARRAY)
1717     return 0;
1718
1719   return e->value.constructor == NULL;
1720 }
1721
1722
1723 /* Reduce a binary expression where at least one of the operands
1724    involves a zero-length array.  Returns NULL if neither of the
1725    operands is a zero-length array.  */
1726
1727 static gfc_expr *
1728 reduce_binary0 (gfc_expr * op1, gfc_expr * op2)
1729 {
1730   if (gfc_zero_size_array (op1))
1731     {
1732       gfc_free_expr (op2);
1733       return op1;
1734     }
1735
1736   if (gfc_zero_size_array (op2))
1737     {
1738       gfc_free_expr (op1);
1739       return op2;
1740     }
1741
1742   return NULL;
1743 }
1744
1745
1746 static gfc_expr *
1747 eval_intrinsic_f2 (gfc_intrinsic_op operator,
1748                    arith (*eval) (gfc_expr *, gfc_expr **),
1749                    gfc_expr * op1, gfc_expr * op2)
1750 {
1751   gfc_expr *result;
1752   eval_f f;
1753
1754   if (op2 == NULL)
1755     {
1756       if (gfc_zero_size_array (op1))
1757         return eval_type_intrinsic0 (operator, op1);
1758     }
1759   else
1760     {
1761       result = reduce_binary0 (op1, op2);
1762       if (result != NULL)
1763         return eval_type_intrinsic0 (operator, result);
1764     }
1765
1766   f.f2 = eval;
1767   return eval_intrinsic (operator, f, op1, op2);
1768 }
1769
1770
1771 static gfc_expr *
1772 eval_intrinsic_f3 (gfc_intrinsic_op operator,
1773                    arith (*eval) (gfc_expr *, gfc_expr *, gfc_expr **),
1774                    gfc_expr * op1, gfc_expr * op2)
1775 {
1776   gfc_expr *result;
1777   eval_f f;
1778
1779   result = reduce_binary0 (op1, op2);
1780   if (result != NULL)
1781     return eval_type_intrinsic0(operator, result);
1782
1783   f.f3 = eval;
1784   return eval_intrinsic (operator, f, op1, op2);
1785 }
1786
1787
1788 gfc_expr *
1789 gfc_uplus (gfc_expr * op)
1790 {
1791   return eval_intrinsic_f2 (INTRINSIC_UPLUS, gfc_arith_uplus, op, NULL);
1792 }
1793
1794
1795 gfc_expr *
1796 gfc_uminus (gfc_expr * op)
1797 {
1798   return eval_intrinsic_f2 (INTRINSIC_UMINUS, gfc_arith_uminus, op, NULL);
1799 }
1800
1801
1802 gfc_expr *
1803 gfc_add (gfc_expr * op1, gfc_expr * op2)
1804 {
1805   return eval_intrinsic_f3 (INTRINSIC_PLUS, gfc_arith_plus, op1, op2);
1806 }
1807
1808
1809 gfc_expr *
1810 gfc_subtract (gfc_expr * op1, gfc_expr * op2)
1811 {
1812   return eval_intrinsic_f3 (INTRINSIC_MINUS, gfc_arith_minus, op1, op2);
1813 }
1814
1815
1816 gfc_expr *
1817 gfc_multiply (gfc_expr * op1, gfc_expr * op2)
1818 {
1819   return eval_intrinsic_f3 (INTRINSIC_TIMES, gfc_arith_times, op1, op2);
1820 }
1821
1822
1823 gfc_expr *
1824 gfc_divide (gfc_expr * op1, gfc_expr * op2)
1825 {
1826   return eval_intrinsic_f3 (INTRINSIC_DIVIDE, gfc_arith_divide, op1, op2);
1827 }
1828
1829
1830 gfc_expr *
1831 gfc_power (gfc_expr * op1, gfc_expr * op2)
1832 {
1833   return eval_intrinsic_f3 (INTRINSIC_POWER, gfc_arith_power, op1, op2);
1834 }
1835
1836
1837 gfc_expr *
1838 gfc_concat (gfc_expr * op1, gfc_expr * op2)
1839 {
1840   return eval_intrinsic_f3 (INTRINSIC_CONCAT, gfc_arith_concat, op1, op2);
1841 }
1842
1843
1844 gfc_expr *
1845 gfc_and (gfc_expr * op1, gfc_expr * op2)
1846 {
1847   return eval_intrinsic_f3 (INTRINSIC_AND, gfc_arith_and, op1, op2);
1848 }
1849
1850
1851 gfc_expr *
1852 gfc_or (gfc_expr * op1, gfc_expr * op2)
1853 {
1854   return eval_intrinsic_f3 (INTRINSIC_OR, gfc_arith_or, op1, op2);
1855 }
1856
1857
1858 gfc_expr *
1859 gfc_not (gfc_expr * op1)
1860 {
1861   return eval_intrinsic_f2 (INTRINSIC_NOT, gfc_arith_not, op1, NULL);
1862 }
1863
1864
1865 gfc_expr *
1866 gfc_eqv (gfc_expr * op1, gfc_expr * op2)
1867 {
1868   return eval_intrinsic_f3 (INTRINSIC_EQV, gfc_arith_eqv, op1, op2);
1869 }
1870
1871
1872 gfc_expr *
1873 gfc_neqv (gfc_expr * op1, gfc_expr * op2)
1874 {
1875   return eval_intrinsic_f3 (INTRINSIC_NEQV, gfc_arith_neqv, op1, op2);
1876 }
1877
1878
1879 gfc_expr *
1880 gfc_eq (gfc_expr * op1, gfc_expr * op2)
1881 {
1882   return eval_intrinsic_f3 (INTRINSIC_EQ, gfc_arith_eq, op1, op2);
1883 }
1884
1885
1886 gfc_expr *
1887 gfc_ne (gfc_expr * op1, gfc_expr * op2)
1888 {
1889   return eval_intrinsic_f3 (INTRINSIC_NE, gfc_arith_ne, op1, op2);
1890 }
1891
1892
1893 gfc_expr *
1894 gfc_gt (gfc_expr * op1, gfc_expr * op2)
1895 {
1896   return eval_intrinsic_f3 (INTRINSIC_GT, gfc_arith_gt, op1, op2);
1897 }
1898
1899
1900 gfc_expr *
1901 gfc_ge (gfc_expr * op1, gfc_expr * op2)
1902 {
1903   return eval_intrinsic_f3 (INTRINSIC_GE, gfc_arith_ge, op1, op2);
1904 }
1905
1906
1907 gfc_expr *
1908 gfc_lt (gfc_expr * op1, gfc_expr * op2)
1909 {
1910   return eval_intrinsic_f3 (INTRINSIC_LT, gfc_arith_lt, op1, op2);
1911 }
1912
1913
1914 gfc_expr *
1915 gfc_le (gfc_expr * op1, gfc_expr * op2)
1916 {
1917   return eval_intrinsic_f3 (INTRINSIC_LE, gfc_arith_le, op1, op2);
1918 }
1919
1920
1921 /* Convert an integer string to an expression node.  */
1922
1923 gfc_expr *
1924 gfc_convert_integer (const char * buffer, int kind, int radix, locus * where)
1925 {
1926   gfc_expr *e;
1927   const char *t;
1928
1929   e = gfc_constant_result (BT_INTEGER, kind, where);
1930   /* A leading plus is allowed, but not by mpz_set_str.  */
1931   if (buffer[0] == '+')
1932     t = buffer + 1;
1933   else
1934     t = buffer;
1935   mpz_set_str (e->value.integer, t, radix);
1936
1937   return e;
1938 }
1939
1940
1941 /* Convert a real string to an expression node.  */
1942
1943 gfc_expr *
1944 gfc_convert_real (const char * buffer, int kind, locus * where)
1945 {
1946   gfc_expr *e;
1947
1948   e = gfc_constant_result (BT_REAL, kind, where);
1949   mpfr_set_str (e->value.real, buffer, 10, GFC_RND_MODE);
1950
1951   return e;
1952 }
1953
1954
1955 /* Convert a pair of real, constant expression nodes to a single
1956    complex expression node.  */
1957
1958 gfc_expr *
1959 gfc_convert_complex (gfc_expr * real, gfc_expr * imag, int kind)
1960 {
1961   gfc_expr *e;
1962
1963   e = gfc_constant_result (BT_COMPLEX, kind, &real->where);
1964   mpfr_set (e->value.complex.r, real->value.real, GFC_RND_MODE);
1965   mpfr_set (e->value.complex.i, imag->value.real, GFC_RND_MODE);
1966
1967   return e;
1968 }
1969
1970
1971 /******* Simplification of intrinsic functions with constant arguments *****/
1972
1973
1974 /* Deal with an arithmetic error.  */
1975
1976 static void
1977 arith_error (arith rc, gfc_typespec * from, gfc_typespec * to, locus * where)
1978 {
1979   switch (rc)
1980     {
1981     case ARITH_OK:
1982       gfc_error ("Arithmetic OK converting %s to %s at %L",
1983                  gfc_typename (from), gfc_typename (to), where);
1984       break;
1985     case ARITH_OVERFLOW:
1986       gfc_error ("Arithmetic overflow converting %s to %s at %L",
1987                  gfc_typename (from), gfc_typename (to), where);
1988       break;
1989     case ARITH_UNDERFLOW:
1990       gfc_error ("Arithmetic underflow converting %s to %s at %L",
1991                  gfc_typename (from), gfc_typename (to), where);
1992       break;
1993     case ARITH_NAN:
1994       gfc_error ("Arithmetic NaN converting %s to %s at %L",
1995                  gfc_typename (from), gfc_typename (to), where);
1996       break;
1997     case ARITH_DIV0:
1998       gfc_error ("Division by zero converting %s to %s at %L",
1999                  gfc_typename (from), gfc_typename (to), where);
2000       break;
2001     case ARITH_INCOMMENSURATE:
2002       gfc_error ("Array operands are incommensurate converting %s to %s at %L",
2003                  gfc_typename (from), gfc_typename (to), where);
2004       break;
2005     case ARITH_ASYMMETRIC:
2006       gfc_error ("Integer outside symmetric range implied by Standard Fortran"
2007                  " converting %s to %s at %L",
2008                  gfc_typename (from), gfc_typename (to), where);
2009       break;
2010     default:
2011       gfc_internal_error ("gfc_arith_error(): Bad error code");
2012     }
2013
2014   /* TODO: Do something about the error, ie, throw exception, return
2015      NaN, etc.  */
2016 }
2017
2018
2019 /* Convert integers to integers.  */
2020
2021 gfc_expr *
2022 gfc_int2int (gfc_expr * src, int kind)
2023 {
2024   gfc_expr *result;
2025   arith rc;
2026
2027   result = gfc_constant_result (BT_INTEGER, kind, &src->where);
2028
2029   mpz_set (result->value.integer, src->value.integer);
2030
2031   if ((rc = gfc_check_integer_range (result->value.integer, kind))
2032       != ARITH_OK)
2033     {
2034       if (rc == ARITH_ASYMMETRIC)
2035         {
2036           gfc_warning (gfc_arith_error (rc), &src->where);
2037         }
2038       else
2039         {
2040           arith_error (rc, &src->ts, &result->ts, &src->where);
2041           gfc_free_expr (result);
2042           return NULL;
2043         }
2044     }
2045
2046   return result;
2047 }
2048
2049
2050 /* Convert integers to reals.  */
2051
2052 gfc_expr *
2053 gfc_int2real (gfc_expr * src, int kind)
2054 {
2055   gfc_expr *result;
2056   arith rc;
2057
2058   result = gfc_constant_result (BT_REAL, kind, &src->where);
2059
2060   mpfr_set_z (result->value.real, src->value.integer, GFC_RND_MODE);
2061
2062   if ((rc = gfc_check_real_range (result->value.real, kind)) != ARITH_OK)
2063     {
2064       arith_error (rc, &src->ts, &result->ts, &src->where);
2065       gfc_free_expr (result);
2066       return NULL;
2067     }
2068
2069   return result;
2070 }
2071
2072
2073 /* Convert default integer to default complex.  */
2074
2075 gfc_expr *
2076 gfc_int2complex (gfc_expr * src, int kind)
2077 {
2078   gfc_expr *result;
2079   arith rc;
2080
2081   result = gfc_constant_result (BT_COMPLEX, kind, &src->where);
2082
2083   mpfr_set_z (result->value.complex.r, src->value.integer, GFC_RND_MODE);
2084   mpfr_set_ui (result->value.complex.i, 0, GFC_RND_MODE);
2085
2086   if ((rc = gfc_check_real_range (result->value.complex.r, kind)) != ARITH_OK)
2087     {
2088       arith_error (rc, &src->ts, &result->ts, &src->where);
2089       gfc_free_expr (result);
2090       return NULL;
2091     }
2092
2093   return result;
2094 }
2095
2096
2097 /* Convert default real to default integer.  */
2098
2099 gfc_expr *
2100 gfc_real2int (gfc_expr * src, int kind)
2101 {
2102   gfc_expr *result;
2103   arith rc;
2104
2105   result = gfc_constant_result (BT_INTEGER, kind, &src->where);
2106
2107   gfc_mpfr_to_mpz (result->value.integer, src->value.real);
2108
2109   if ((rc = gfc_check_integer_range (result->value.integer, kind))
2110       != ARITH_OK)
2111     {
2112       arith_error (rc, &src->ts, &result->ts, &src->where);
2113       gfc_free_expr (result);
2114       return NULL;
2115     }
2116
2117   return result;
2118 }
2119
2120
2121 /* Convert real to real.  */
2122
2123 gfc_expr *
2124 gfc_real2real (gfc_expr * src, int kind)
2125 {
2126   gfc_expr *result;
2127   arith rc;
2128
2129   result = gfc_constant_result (BT_REAL, kind, &src->where);
2130
2131   mpfr_set (result->value.real, src->value.real, GFC_RND_MODE);
2132
2133   rc = gfc_check_real_range (result->value.real, kind);
2134
2135   if (rc == ARITH_UNDERFLOW)
2136     {
2137       if (gfc_option.warn_underflow)
2138         gfc_warning (gfc_arith_error (rc), &src->where);
2139       mpfr_set_ui (result->value.real, 0, GFC_RND_MODE);
2140     }
2141   else if (rc != ARITH_OK)
2142     {
2143       arith_error (rc, &src->ts, &result->ts, &src->where);
2144       gfc_free_expr (result);
2145       return NULL;
2146     }
2147
2148   return result;
2149 }
2150
2151
2152 /* Convert real to complex.  */
2153
2154 gfc_expr *
2155 gfc_real2complex (gfc_expr * src, int kind)
2156 {
2157   gfc_expr *result;
2158   arith rc;
2159
2160   result = gfc_constant_result (BT_COMPLEX, kind, &src->where);
2161
2162   mpfr_set (result->value.complex.r, src->value.real, GFC_RND_MODE);
2163   mpfr_set_ui (result->value.complex.i, 0, GFC_RND_MODE);
2164
2165   rc = gfc_check_real_range (result->value.complex.r, kind);
2166
2167   if (rc == ARITH_UNDERFLOW)
2168     {
2169       if (gfc_option.warn_underflow)
2170         gfc_warning (gfc_arith_error (rc), &src->where);
2171       mpfr_set_ui (result->value.complex.r, 0, GFC_RND_MODE);
2172     }
2173   else if (rc != ARITH_OK)
2174     {
2175       arith_error (rc, &src->ts, &result->ts, &src->where);
2176       gfc_free_expr (result);
2177       return NULL;
2178     }
2179
2180   return result;
2181 }
2182
2183
2184 /* Convert complex to integer.  */
2185
2186 gfc_expr *
2187 gfc_complex2int (gfc_expr * src, int kind)
2188 {
2189   gfc_expr *result;
2190   arith rc;
2191
2192   result = gfc_constant_result (BT_INTEGER, kind, &src->where);
2193
2194   gfc_mpfr_to_mpz (result->value.integer, src->value.complex.r);
2195
2196   if ((rc = gfc_check_integer_range (result->value.integer, kind))
2197       != ARITH_OK)
2198     {
2199       arith_error (rc, &src->ts, &result->ts, &src->where);
2200       gfc_free_expr (result);
2201       return NULL;
2202     }
2203
2204   return result;
2205 }
2206
2207
2208 /* Convert complex to real.  */
2209
2210 gfc_expr *
2211 gfc_complex2real (gfc_expr * src, int kind)
2212 {
2213   gfc_expr *result;
2214   arith rc;
2215
2216   result = gfc_constant_result (BT_REAL, kind, &src->where);
2217
2218   mpfr_set (result->value.real, src->value.complex.r, GFC_RND_MODE);
2219
2220   rc = gfc_check_real_range (result->value.real, kind);
2221
2222   if (rc == ARITH_UNDERFLOW)
2223     {
2224       if (gfc_option.warn_underflow)
2225         gfc_warning (gfc_arith_error (rc), &src->where);
2226       mpfr_set_ui (result->value.real, 0, GFC_RND_MODE);
2227     }
2228   if (rc != ARITH_OK)
2229     {
2230       arith_error (rc, &src->ts, &result->ts, &src->where);
2231       gfc_free_expr (result);
2232       return NULL;
2233     }
2234
2235   return result;
2236 }
2237
2238
2239 /* Convert complex to complex.  */
2240
2241 gfc_expr *
2242 gfc_complex2complex (gfc_expr * src, int kind)
2243 {
2244   gfc_expr *result;
2245   arith rc;
2246
2247   result = gfc_constant_result (BT_COMPLEX, kind, &src->where);
2248
2249   mpfr_set (result->value.complex.r, src->value.complex.r, GFC_RND_MODE);
2250   mpfr_set (result->value.complex.i, src->value.complex.i, GFC_RND_MODE);
2251
2252   rc = gfc_check_real_range (result->value.complex.r, kind);
2253
2254   if (rc == ARITH_UNDERFLOW)
2255     {
2256       if (gfc_option.warn_underflow)
2257         gfc_warning (gfc_arith_error (rc), &src->where);
2258       mpfr_set_ui (result->value.complex.r, 0, GFC_RND_MODE);
2259     }
2260   else if (rc != ARITH_OK)
2261     {
2262       arith_error (rc, &src->ts, &result->ts, &src->where);
2263       gfc_free_expr (result);
2264       return NULL;
2265     }
2266
2267   rc = gfc_check_real_range (result->value.complex.i, kind);
2268
2269   if (rc == ARITH_UNDERFLOW)
2270     {
2271       if (gfc_option.warn_underflow)
2272         gfc_warning (gfc_arith_error (rc), &src->where);
2273       mpfr_set_ui (result->value.complex.i, 0, GFC_RND_MODE);
2274     }
2275   else if (rc != ARITH_OK)
2276     {
2277       arith_error (rc, &src->ts, &result->ts, &src->where);
2278       gfc_free_expr (result);
2279       return NULL;
2280     }
2281
2282   return result;
2283 }
2284
2285
2286 /* Logical kind conversion.  */
2287
2288 gfc_expr *
2289 gfc_log2log (gfc_expr * src, int kind)
2290 {
2291   gfc_expr *result;
2292
2293   result = gfc_constant_result (BT_LOGICAL, kind, &src->where);
2294   result->value.logical = src->value.logical;
2295
2296   return result;
2297 }
2298
2299
2300 /* Convert logical to integer.  */
2301
2302 gfc_expr *
2303 gfc_log2int (gfc_expr *src, int kind)
2304 {
2305   gfc_expr *result;
2306
2307   result = gfc_constant_result (BT_INTEGER, kind, &src->where);
2308   mpz_set_si (result->value.integer, src->value.logical);
2309
2310   return result;
2311 }
2312
2313
2314 /* Convert integer to logical.  */
2315
2316 gfc_expr *
2317 gfc_int2log (gfc_expr *src, int kind)
2318 {
2319   gfc_expr *result;
2320
2321   result = gfc_constant_result (BT_LOGICAL, kind, &src->where);
2322   result->value.logical = (mpz_cmp_si (src->value.integer, 0) != 0);
2323
2324   return result;
2325 }
2326
2327
2328 /* Convert Hollerith to integer. The constant will be padded or truncated.  */
2329
2330 gfc_expr *
2331 gfc_hollerith2int (gfc_expr * src, int kind)
2332 {
2333   gfc_expr *result;
2334   int len;
2335
2336   len = src->value.character.length;
2337
2338   result = gfc_get_expr ();
2339   result->expr_type = EXPR_CONSTANT;
2340   result->ts.type = BT_INTEGER;
2341   result->ts.kind = kind;
2342   result->where = src->where;
2343   result->from_H = 1;
2344
2345   if (len > kind)
2346     {
2347       gfc_warning ("The Hollerith constant at %L is too long to convert to %s",
2348                 &src->where, gfc_typename(&result->ts));
2349     }
2350   result->value.character.string = gfc_getmem (kind + 1);
2351   memcpy (result->value.character.string, src->value.character.string,
2352         MIN (kind, len));
2353
2354   if (len < kind)
2355     memset (&result->value.character.string[len], ' ', kind - len);
2356
2357   result->value.character.string[kind] = '\0'; /* For debugger  */
2358   result->value.character.length = kind;
2359
2360   return result;
2361 }
2362
2363
2364 /* Convert Hollerith to real. The constant will be padded or truncated.  */
2365
2366 gfc_expr *
2367 gfc_hollerith2real (gfc_expr * src, int kind)
2368 {
2369   gfc_expr *result;
2370   int len;
2371
2372   len = src->value.character.length;
2373
2374   result = gfc_get_expr ();
2375   result->expr_type = EXPR_CONSTANT;
2376   result->ts.type = BT_REAL;
2377   result->ts.kind = kind;
2378   result->where = src->where;
2379   result->from_H = 1;
2380
2381   if (len > kind)
2382     {
2383       gfc_warning ("The Hollerith constant at %L is too long to convert to %s",
2384                 &src->where, gfc_typename(&result->ts));
2385     }
2386   result->value.character.string = gfc_getmem (kind + 1);
2387   memcpy (result->value.character.string, src->value.character.string,
2388         MIN (kind, len));
2389
2390   if (len < kind)
2391     memset (&result->value.character.string[len], ' ', kind - len);
2392
2393   result->value.character.string[kind] = '\0'; /* For debugger.  */
2394   result->value.character.length = kind;
2395
2396   return result;
2397 }
2398
2399
2400 /* Convert Hollerith to complex. The constant will be padded or truncated.  */
2401
2402 gfc_expr *
2403 gfc_hollerith2complex (gfc_expr * src, int kind)
2404 {
2405   gfc_expr *result;
2406   int len;
2407
2408   len = src->value.character.length;
2409
2410   result = gfc_get_expr ();
2411   result->expr_type = EXPR_CONSTANT;
2412   result->ts.type = BT_COMPLEX;
2413   result->ts.kind = kind;
2414   result->where = src->where;
2415   result->from_H = 1;
2416
2417   kind = kind * 2;
2418
2419   if (len > kind)
2420     {
2421       gfc_warning ("The Hollerith constant at %L is too long to convert to %s",
2422                 &src->where, gfc_typename(&result->ts));
2423     }
2424   result->value.character.string = gfc_getmem (kind + 1);
2425   memcpy (result->value.character.string, src->value.character.string,
2426         MIN (kind, len));
2427
2428   if (len < kind)
2429     memset (&result->value.character.string[len], ' ', kind - len);
2430
2431   result->value.character.string[kind] = '\0'; /* For debugger  */
2432   result->value.character.length = kind;
2433
2434   return result;
2435 }
2436
2437
2438 /* Convert Hollerith to character. */
2439
2440 gfc_expr *
2441 gfc_hollerith2character (gfc_expr * src, int kind)
2442 {
2443   gfc_expr *result;
2444
2445   result = gfc_copy_expr (src);
2446   result->ts.type = BT_CHARACTER;
2447   result->ts.kind = kind;
2448   result->from_H = 1;
2449
2450   return result;
2451 }
2452
2453
2454 /* Convert Hollerith to logical. The constant will be padded or truncated.  */
2455
2456 gfc_expr *
2457 gfc_hollerith2logical (gfc_expr * src, int kind)
2458 {
2459   gfc_expr *result;
2460   int len;
2461
2462   len = src->value.character.length;
2463
2464   result = gfc_get_expr ();
2465   result->expr_type = EXPR_CONSTANT;
2466   result->ts.type = BT_LOGICAL;
2467   result->ts.kind = kind;
2468   result->where = src->where;
2469   result->from_H = 1;
2470
2471   if (len > kind)
2472     {
2473       gfc_warning ("The Hollerith constant at %L is too long to convert to %s",
2474                 &src->where, gfc_typename(&result->ts));
2475     }
2476   result->value.character.string = gfc_getmem (kind + 1);
2477   memcpy (result->value.character.string, src->value.character.string,
2478         MIN (kind, len));
2479
2480   if (len < kind)
2481     memset (&result->value.character.string[len], ' ', kind - len);
2482
2483   result->value.character.string[kind] = '\0'; /* For debugger  */
2484   result->value.character.length = kind;
2485
2486   return result;
2487 }
2488
2489
2490 /* Returns an initializer whose value is one higher than the value of the
2491    LAST_INITIALIZER argument.  If the argument is NULL, the
2492    initializers value will be set to zero.  The initializer's kind
2493    will be set to gfc_c_int_kind.
2494
2495    If -fshort-enums is given, the appropriate kind will be selected
2496    later after all enumerators have been parsed.  A warning is issued
2497    here if an initializer exceeds gfc_c_int_kind.  */
2498
2499 gfc_expr *
2500 gfc_enum_initializer (gfc_expr * last_initializer, locus where)
2501 {
2502   gfc_expr *result;
2503
2504   result = gfc_get_expr ();
2505   result->expr_type = EXPR_CONSTANT;
2506   result->ts.type = BT_INTEGER;
2507   result->ts.kind = gfc_c_int_kind;
2508   result->where = where;
2509
2510   mpz_init (result->value.integer);
2511
2512   if (last_initializer != NULL)
2513     {
2514       mpz_add_ui (result->value.integer, last_initializer->value.integer, 1);
2515       result->where = last_initializer->where;
2516
2517       if (gfc_check_integer_range (result->value.integer,
2518              gfc_c_int_kind) != ARITH_OK)
2519         {
2520           gfc_error ("Enumerator exceeds the C integer type at %C");
2521           return NULL;
2522         }
2523     }
2524   else
2525     {
2526       /* Control comes here, if it's the very first enumerator and no
2527          initializer has been given.  It will be initialized to zero.  */
2528       mpz_set_si (result->value.integer, 0);
2529     }
2530
2531   return result;
2532 }