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