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