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. This check "
1985                  "can be disabled with the option -fno-range-check",
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 }