gcc/fortran/arith.c

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