gcc/fortran/arith.c

   1 /* Compiler arithmetic
   2    Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007
   3    Free Software Foundation, Inc.
   4    Contributed by Andy Vaught
   5
   6 This file is part of GCC.
   7
   8 GCC is free software; you can redistribute it and/or modify it under
   9 the terms of the GNU General Public License as published by the Free
  10 Software Foundation; either version 2, or (at your option) any later
  11 version.
  12
  13 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
  14 WARRANTY; without even the implied warranty of MERCHANTABILITY or
  15 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  16 for more details.
  17
  18 You should have received a copy of the GNU General Public License
  19 along with GCC; see the file COPYING.  If not, write to the Free
  20 Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA
  21 02110-1301, USA.  */
  22
  23 /* Since target arithmetic must be done on the host, there has to
  24    be some way of evaluating arithmetic expressions as the host
  25    would evaluate them.  We use the GNU MP library and the MPFR
  26    library to do arithmetic, and this file provides the interface.  */
  27
  28 #include "config.h"
  29 #include "system.h"
  30 #include "flags.h"
  31 #include "gfortran.h"
  32 #include "arith.h"
  33
  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_LT:
1542     case INTRINSIC_LE:
1543     case INTRINSIC_GT:
1544       if (op1->ts.type == BT_COMPLEX || op2->ts.type == BT_COMPLEX)
1545         {
1546           temp.ts.type = BT_LOGICAL;
1547           temp.ts.kind = gfc_default_logical_kind;
1548           goto runtime;
1549         }
1550
1551     /* Fall through  */
1552     case INTRINSIC_EQ:
1553     case INTRINSIC_NE:
1554       if (op1->ts.type == BT_CHARACTER && op2->ts.type == BT_CHARACTER)
1555         {
1556           unary = 0;
1557           temp.ts.type = BT_LOGICAL;
1558           temp.ts.kind = gfc_default_logical_kind;
1559           break;
1560         }
1561
1562     /* Fall through  */
1563     /* Numeric binary  */
1564     case INTRINSIC_PLUS:
1565     case INTRINSIC_MINUS:
1566     case INTRINSIC_TIMES:
1567     case INTRINSIC_DIVIDE:
1568     case INTRINSIC_POWER:
1569       if (!gfc_numeric_ts (&op1->ts) || !gfc_numeric_ts (&op2->ts))
1570         goto runtime;
1571
1572       /* Insert any necessary type conversions to make the operands
1573          compatible.  */
1574
1575       temp.expr_type = EXPR_OP;
1576       gfc_clear_ts (&temp.ts);
1577       temp.value.op.operator = operator;
1578
1579       temp.value.op.op1 = op1;
1580       temp.value.op.op2 = op2;
1581
1582       gfc_type_convert_binary (&temp);
1583
1584       if (operator == INTRINSIC_EQ || operator == INTRINSIC_NE
1585           || operator == INTRINSIC_GE || operator == INTRINSIC_GT
1586           || operator == INTRINSIC_LE || operator == INTRINSIC_LT)
1587         {
1588           temp.ts.type = BT_LOGICAL;
1589           temp.ts.kind = gfc_default_logical_kind;
1590         }
1591
1592       unary = 0;
1593       break;
1594
1595     /* Character binary  */
1596     case INTRINSIC_CONCAT:
1597       if (op1->ts.type != BT_CHARACTER || op2->ts.type != BT_CHARACTER)
1598         goto runtime;
1599
1600       temp.ts.type = BT_CHARACTER;
1601       temp.ts.kind = gfc_default_character_kind;
1602       unary = 0;
1603       break;
1604
1605     case INTRINSIC_USER:
1606       goto runtime;
1607
1608     default:
1609       gfc_internal_error ("eval_intrinsic(): Bad operator");
1610     }
1611
1612   /* Try to combine the operators.  */
1613   if (operator == INTRINSIC_POWER && op2->ts.type != BT_INTEGER)
1614     goto runtime;
1615
1616   if (op1->from_H
1617       || (op1->expr_type != EXPR_CONSTANT
1618           && (op1->expr_type != EXPR_ARRAY
1619               || !gfc_is_constant_expr (op1) || !gfc_expanded_ac (op1))))
1620     goto runtime;
1621
1622   if (op2 != NULL
1623       && (op2->from_H
1624           || (op2->expr_type != EXPR_CONSTANT
1625               && (op2->expr_type != EXPR_ARRAY
1626                   || !gfc_is_constant_expr (op2) || !gfc_expanded_ac (op2)))))
1627     goto runtime;
1628
1629   if (unary)
1630     rc = reduce_unary (eval.f2, op1, &result);
1631   else
1632     rc = reduce_binary (eval.f3, op1, op2, &result);
1633
1634   if (rc != ARITH_OK)
1635     { /* Something went wrong.  */
1636       gfc_error (gfc_arith_error (rc), &op1->where);
1637       return NULL;
1638     }
1639
1640   gfc_free_expr (op1);
1641   gfc_free_expr (op2);
1642   return result;
1643
1644 runtime:
1645   /* Create a run-time expression.  */
1646   result = gfc_get_expr ();
1647   result->ts = temp.ts;
1648
1649   result->expr_type = EXPR_OP;
1650   result->value.op.operator = operator;
1651
1652   result->value.op.op1 = op1;
1653   result->value.op.op2 = op2;
1654
1655   result->where = op1->where;
1656
1657   return result;
1658 }
1659
1660
1661 /* Modify type of expression for zero size array.  */
1662
1663 static gfc_expr *
1664 eval_type_intrinsic0 (gfc_intrinsic_op operator, gfc_expr *op)
1665 {
1666   if (op == NULL)
1667     gfc_internal_error ("eval_type_intrinsic0(): op NULL");
1668
1669   switch (operator)
1670     {
1671     case INTRINSIC_GE:
1672     case INTRINSIC_LT:
1673     case INTRINSIC_LE:
1674     case INTRINSIC_GT:
1675     case INTRINSIC_EQ:
1676     case INTRINSIC_NE:
1677       op->ts.type = BT_LOGICAL;
1678       op->ts.kind = gfc_default_logical_kind;
1679       break;
1680
1681     default:
1682       break;
1683     }
1684
1685   return op;
1686 }
1687
1688
1689 /* Return nonzero if the expression is a zero size array.  */
1690
1691 static int
1692 gfc_zero_size_array (gfc_expr *e)
1693 {
1694   if (e->expr_type != EXPR_ARRAY)
1695     return 0;
1696
1697   return e->value.constructor == NULL;
1698 }
1699
1700
1701 /* Reduce a binary expression where at least one of the operands
1702    involves a zero-length array.  Returns NULL if neither of the
1703    operands is a zero-length array.  */
1704
1705 static gfc_expr *
1706 reduce_binary0 (gfc_expr *op1, gfc_expr *op2)
1707 {
1708   if (gfc_zero_size_array (op1))
1709     {
1710       gfc_free_expr (op2);
1711       return op1;
1712     }
1713
1714   if (gfc_zero_size_array (op2))
1715     {
1716       gfc_free_expr (op1);
1717       return op2;
1718     }
1719
1720   return NULL;
1721 }
1722
1723
1724 static gfc_expr *
1725 eval_intrinsic_f2 (gfc_intrinsic_op operator,
1726                    arith (*eval) (gfc_expr *, gfc_expr **),
1727                    gfc_expr *op1, gfc_expr *op2)
1728 {
1729   gfc_expr *result;
1730   eval_f f;
1731
1732   if (op2 == NULL)
1733     {
1734       if (gfc_zero_size_array (op1))
1735         return eval_type_intrinsic0 (operator, op1);
1736     }
1737   else
1738     {
1739       result = reduce_binary0 (op1, op2);
1740       if (result != NULL)
1741         return eval_type_intrinsic0 (operator, result);
1742     }
1743
1744   f.f2 = eval;
1745   return eval_intrinsic (operator, f, op1, op2);
1746 }
1747
1748
1749 static gfc_expr *
1750 eval_intrinsic_f3 (gfc_intrinsic_op operator,
1751                    arith (*eval) (gfc_expr *, gfc_expr *, gfc_expr **),
1752                    gfc_expr *op1, gfc_expr *op2)
1753 {
1754   gfc_expr *result;
1755   eval_f f;
1756
1757   result = reduce_binary0 (op1, op2);
1758   if (result != NULL)
1759     return eval_type_intrinsic0(operator, result);
1760
1761   f.f3 = eval;
1762   return eval_intrinsic (operator, f, op1, op2);
1763 }
1764
1765
1766 gfc_expr *
1767 gfc_parentheses (gfc_expr *op)
1768 {
1769   return eval_intrinsic_f2 (INTRINSIC_PARENTHESES, gfc_arith_identity,
1770                             op, NULL);
1771 }
1772
1773 gfc_expr *
1774 gfc_uplus (gfc_expr *op)
1775 {
1776   return eval_intrinsic_f2 (INTRINSIC_UPLUS, gfc_arith_identity, op, NULL);
1777 }
1778
1779
1780 gfc_expr *
1781 gfc_uminus (gfc_expr *op)
1782 {
1783   return eval_intrinsic_f2 (INTRINSIC_UMINUS, gfc_arith_uminus, op, NULL);
1784 }
1785
1786
1787 gfc_expr *
1788 gfc_add (gfc_expr *op1, gfc_expr *op2)
1789 {
1790   return eval_intrinsic_f3 (INTRINSIC_PLUS, gfc_arith_plus, op1, op2);
1791 }
1792
1793
1794 gfc_expr *
1795 gfc_subtract (gfc_expr *op1, gfc_expr *op2)
1796 {
1797   return eval_intrinsic_f3 (INTRINSIC_MINUS, gfc_arith_minus, op1, op2);
1798 }
1799
1800
1801 gfc_expr *
1802 gfc_multiply (gfc_expr *op1, gfc_expr *op2)
1803 {
1804   return eval_intrinsic_f3 (INTRINSIC_TIMES, gfc_arith_times, op1, op2);
1805 }
1806
1807
1808 gfc_expr *
1809 gfc_divide (gfc_expr *op1, gfc_expr *op2)
1810 {
1811   return eval_intrinsic_f3 (INTRINSIC_DIVIDE, gfc_arith_divide, op1, op2);
1812 }
1813
1814
1815 gfc_expr *
1816 gfc_power (gfc_expr *op1, gfc_expr *op2)
1817 {
1818   return eval_intrinsic_f3 (INTRINSIC_POWER, gfc_arith_power, op1, op2);
1819 }
1820
1821
1822 gfc_expr *
1823 gfc_concat (gfc_expr *op1, gfc_expr *op2)
1824 {
1825   return eval_intrinsic_f3 (INTRINSIC_CONCAT, gfc_arith_concat, op1, op2);
1826 }
1827
1828
1829 gfc_expr *
1830 gfc_and (gfc_expr *op1, gfc_expr *op2)
1831 {
1832   return eval_intrinsic_f3 (INTRINSIC_AND, gfc_arith_and, op1, op2);
1833 }
1834
1835
1836 gfc_expr *
1837 gfc_or (gfc_expr *op1, gfc_expr *op2)
1838 {
1839   return eval_intrinsic_f3 (INTRINSIC_OR, gfc_arith_or, op1, op2);
1840 }
1841
1842
1843 gfc_expr *
1844 gfc_not (gfc_expr *op1)
1845 {
1846   return eval_intrinsic_f2 (INTRINSIC_NOT, gfc_arith_not, op1, NULL);
1847 }
1848
1849
1850 gfc_expr *
1851 gfc_eqv (gfc_expr *op1, gfc_expr *op2)
1852 {
1853   return eval_intrinsic_f3 (INTRINSIC_EQV, gfc_arith_eqv, op1, op2);
1854 }
1855
1856
1857 gfc_expr *
1858 gfc_neqv (gfc_expr *op1, gfc_expr *op2)
1859 {
1860   return eval_intrinsic_f3 (INTRINSIC_NEQV, gfc_arith_neqv, op1, op2);
1861 }
1862
1863
1864 gfc_expr *
1865 gfc_eq (gfc_expr *op1, gfc_expr *op2)
1866 {
1867   return eval_intrinsic_f3 (INTRINSIC_EQ, gfc_arith_eq, op1, op2);
1868 }
1869
1870
1871 gfc_expr *
1872 gfc_ne (gfc_expr *op1, gfc_expr *op2)
1873 {
1874   return eval_intrinsic_f3 (INTRINSIC_NE, gfc_arith_ne, op1, op2);
1875 }
1876
1877
1878 gfc_expr *
1879 gfc_gt (gfc_expr *op1, gfc_expr *op2)
1880 {
1881   return eval_intrinsic_f3 (INTRINSIC_GT, gfc_arith_gt, op1, op2);
1882 }
1883
1884
1885 gfc_expr *
1886 gfc_ge (gfc_expr *op1, gfc_expr *op2)
1887 {
1888   return eval_intrinsic_f3 (INTRINSIC_GE, gfc_arith_ge, op1, op2);
1889 }
1890
1891
1892 gfc_expr *
1893 gfc_lt (gfc_expr *op1, gfc_expr *op2)
1894 {
1895   return eval_intrinsic_f3 (INTRINSIC_LT, gfc_arith_lt, op1, op2);
1896 }
1897
1898
1899 gfc_expr *
1900 gfc_le (gfc_expr *op1, gfc_expr *op2)
1901 {
1902   return eval_intrinsic_f3 (INTRINSIC_LE, gfc_arith_le, op1, op2);
1903 }
1904
1905
1906 /* Convert an integer string to an expression node.  */
1907
1908 gfc_expr *
1909 gfc_convert_integer (const char *buffer, int kind, int radix, locus *where)
1910 {
1911   gfc_expr *e;
1912   const char *t;
1913
1914   e = gfc_constant_result (BT_INTEGER, kind, where);
1915   /* A leading plus is allowed, but not by mpz_set_str.  */
1916   if (buffer[0] == '+')
1917     t = buffer + 1;
1918   else
1919     t = buffer;
1920   mpz_set_str (e->value.integer, t, radix);
1921
1922   return e;
1923 }
1924
1925
1926 /* Convert a real string to an expression node.  */
1927
1928 gfc_expr *
1929 gfc_convert_real (const char *buffer, int kind, locus *where)
1930 {
1931   gfc_expr *e;
1932
1933   e = gfc_constant_result (BT_REAL, kind, where);
1934   mpfr_set_str (e->value.real, buffer, 10, GFC_RND_MODE);
1935
1936   return e;
1937 }
1938
1939
1940 /* Convert a pair of real, constant expression nodes to a single
1941    complex expression node.  */
1942
1943 gfc_expr *
1944 gfc_convert_complex (gfc_expr *real, gfc_expr *imag, int kind)
1945 {
1946   gfc_expr *e;
1947
1948   e = gfc_constant_result (BT_COMPLEX, kind, &real->where);
1949   mpfr_set (e->value.complex.r, real->value.real, GFC_RND_MODE);
1950   mpfr_set (e->value.complex.i, imag->value.real, GFC_RND_MODE);
1951
1952   return e;
1953 }
1954
1955
1956 /******* Simplification of intrinsic functions with constant arguments *****/
1957
1958
1959 /* Deal with an arithmetic error.  */
1960
1961 static void
1962 arith_error (arith rc, gfc_typespec *from, gfc_typespec *to, locus *where)
1963 {
1964   switch (rc)
1965     {
1966     case ARITH_OK:
1967       gfc_error ("Arithmetic OK converting %s to %s at %L",
1968                  gfc_typename (from), gfc_typename (to), where);
1969       break;
1970     case ARITH_OVERFLOW:
1971       gfc_error ("Arithmetic overflow converting %s to %s at %L",
1972                  gfc_typename (from), gfc_typename (to), where);
1973       break;
1974     case ARITH_UNDERFLOW:
1975       gfc_error ("Arithmetic underflow converting %s to %s at %L",
1976                  gfc_typename (from), gfc_typename (to), where);
1977       break;
1978     case ARITH_NAN:
1979       gfc_error ("Arithmetic NaN converting %s to %s at %L",
1980                  gfc_typename (from), gfc_typename (to), where);
1981       break;
1982     case ARITH_DIV0:
1983       gfc_error ("Division by zero converting %s to %s at %L",
1984                  gfc_typename (from), gfc_typename (to), where);
1985       break;
1986     case ARITH_INCOMMENSURATE:
1987       gfc_error ("Array operands are incommensurate converting %s to %s at %L",
1988                  gfc_typename (from), gfc_typename (to), where);
1989       break;
1990     case ARITH_ASYMMETRIC:
1991       gfc_error ("Integer outside symmetric range implied by Standard Fortran"
1992                  " converting %s to %s at %L",
1993                  gfc_typename (from), gfc_typename (to), where);
1994       break;
1995     default:
1996       gfc_internal_error ("gfc_arith_error(): Bad error code");
1997     }
1998
1999   /* TODO: Do something about the error, ie, throw exception, return
2000      NaN, etc.  */
2001 }
2002
2003
2004 /* Convert integers to integers.  */
2005
2006 gfc_expr *
2007 gfc_int2int (gfc_expr *src, int kind)
2008 {
2009   gfc_expr *result;
2010   arith rc;
2011
2012   result = gfc_constant_result (BT_INTEGER, kind, &src->where);
2013
2014   mpz_set (result->value.integer, src->value.integer);
2015
2016   if ((rc = gfc_check_integer_range (result->value.integer, kind)) != ARITH_OK)
2017     {
2018       if (rc == ARITH_ASYMMETRIC)
2019         {
2020           gfc_warning (gfc_arith_error (rc), &src->where);
2021         }
2022       else
2023         {
2024           arith_error (rc, &src->ts, &result->ts, &src->where);
2025           gfc_free_expr (result);
2026           return NULL;
2027         }
2028     }
2029
2030   return result;
2031 }
2032
2033
2034 /* Convert integers to reals.  */
2035
2036 gfc_expr *
2037 gfc_int2real (gfc_expr *src, int kind)
2038 {
2039   gfc_expr *result;
2040   arith rc;
2041
2042   result = gfc_constant_result (BT_REAL, kind, &src->where);
2043
2044   mpfr_set_z (result->value.real, src->value.integer, GFC_RND_MODE);
2045
2046   if ((rc = gfc_check_real_range (result->value.real, kind)) != ARITH_OK)
2047     {
2048       arith_error (rc, &src->ts, &result->ts, &src->where);
2049       gfc_free_expr (result);
2050       return NULL;
2051     }
2052
2053   return result;
2054 }
2055
2056
2057 /* Convert default integer to default complex.  */
2058
2059 gfc_expr *
2060 gfc_int2complex (gfc_expr *src, int kind)
2061 {
2062   gfc_expr *result;
2063   arith rc;
2064
2065   result = gfc_constant_result (BT_COMPLEX, kind, &src->where);
2066
2067   mpfr_set_z (result->value.complex.r, src->value.integer, GFC_RND_MODE);
2068   mpfr_set_ui (result->value.complex.i, 0, GFC_RND_MODE);
2069
2070   if ((rc = gfc_check_real_range (result->value.complex.r, kind)) != ARITH_OK)
2071     {
2072       arith_error (rc, &src->ts, &result->ts, &src->where);
2073       gfc_free_expr (result);
2074       return NULL;
2075     }
2076
2077   return result;
2078 }
2079
2080
2081 /* Convert default real to default integer.  */
2082
2083 gfc_expr *
2084 gfc_real2int (gfc_expr *src, int kind)
2085 {
2086   gfc_expr *result;
2087   arith rc;
2088
2089   result = gfc_constant_result (BT_INTEGER, kind, &src->where);
2090
2091   gfc_mpfr_to_mpz (result->value.integer, src->value.real);
2092
2093   if ((rc = gfc_check_integer_range (result->value.integer, kind)) != ARITH_OK)
2094     {
2095       arith_error (rc, &src->ts, &result->ts, &src->where);
2096       gfc_free_expr (result);
2097       return NULL;
2098     }
2099
2100   return result;
2101 }
2102
2103
2104 /* Convert real to real.  */
2105
2106 gfc_expr *
2107 gfc_real2real (gfc_expr *src, int kind)
2108 {
2109   gfc_expr *result;
2110   arith rc;
2111
2112   result = gfc_constant_result (BT_REAL, kind, &src->where);
2113
2114   mpfr_set (result->value.real, src->value.real, GFC_RND_MODE);
2115
2116   rc = gfc_check_real_range (result->value.real, kind);
2117
2118   if (rc == ARITH_UNDERFLOW)
2119     {
2120       if (gfc_option.warn_underflow)
2121         gfc_warning (gfc_arith_error (rc), &src->where);
2122       mpfr_set_ui (result->value.real, 0, GFC_RND_MODE);
2123     }
2124   else if (rc != ARITH_OK)
2125     {
2126       arith_error (rc, &src->ts, &result->ts, &src->where);
2127       gfc_free_expr (result);
2128       return NULL;
2129     }
2130
2131   return result;
2132 }
2133
2134
2135 /* Convert real to complex.  */
2136
2137 gfc_expr *
2138 gfc_real2complex (gfc_expr *src, int kind)
2139 {
2140   gfc_expr *result;
2141   arith rc;
2142
2143   result = gfc_constant_result (BT_COMPLEX, kind, &src->where);
2144
2145   mpfr_set (result->value.complex.r, src->value.real, GFC_RND_MODE);
2146   mpfr_set_ui (result->value.complex.i, 0, GFC_RND_MODE);
2147
2148   rc = gfc_check_real_range (result->value.complex.r, kind);
2149
2150   if (rc == ARITH_UNDERFLOW)
2151     {
2152       if (gfc_option.warn_underflow)
2153         gfc_warning (gfc_arith_error (rc), &src->where);
2154       mpfr_set_ui (result->value.complex.r, 0, GFC_RND_MODE);
2155     }
2156   else if (rc != ARITH_OK)
2157     {
2158       arith_error (rc, &src->ts, &result->ts, &src->where);
2159       gfc_free_expr (result);
2160       return NULL;
2161     }
2162
2163   return result;
2164 }
2165
2166
2167 /* Convert complex to integer.  */
2168
2169 gfc_expr *
2170 gfc_complex2int (gfc_expr *src, int kind)
2171 {
2172   gfc_expr *result;
2173   arith rc;
2174
2175   result = gfc_constant_result (BT_INTEGER, kind, &src->where);
2176
2177   gfc_mpfr_to_mpz (result->value.integer, src->value.complex.r);
2178
2179   if ((rc = gfc_check_integer_range (result->value.integer, kind)) != ARITH_OK)
2180     {
2181       arith_error (rc, &src->ts, &result->ts, &src->where);
2182       gfc_free_expr (result);
2183       return NULL;
2184     }
2185
2186   return result;
2187 }
2188
2189
2190 /* Convert complex to real.  */
2191
2192 gfc_expr *
2193 gfc_complex2real (gfc_expr *src, int kind)
2194 {
2195   gfc_expr *result;
2196   arith rc;
2197
2198   result = gfc_constant_result (BT_REAL, kind, &src->where);
2199
2200   mpfr_set (result->value.real, src->value.complex.r, GFC_RND_MODE);
2201
2202   rc = gfc_check_real_range (result->value.real, kind);
2203
2204   if (rc == ARITH_UNDERFLOW)
2205     {
2206       if (gfc_option.warn_underflow)
2207         gfc_warning (gfc_arith_error (rc), &src->where);
2208       mpfr_set_ui (result->value.real, 0, GFC_RND_MODE);
2209     }
2210   if (rc != ARITH_OK)
2211     {
2212       arith_error (rc, &src->ts, &result->ts, &src->where);
2213       gfc_free_expr (result);
2214       return NULL;
2215     }
2216
2217   return result;
2218 }
2219
2220
2221 /* Convert complex to complex.  */
2222
2223 gfc_expr *
2224 gfc_complex2complex (gfc_expr *src, int kind)
2225 {
2226   gfc_expr *result;
2227   arith rc;
2228
2229   result = gfc_constant_result (BT_COMPLEX, kind, &src->where);
2230
2231   mpfr_set (result->value.complex.r, src->value.complex.r, GFC_RND_MODE);
2232   mpfr_set (result->value.complex.i, src->value.complex.i, GFC_RND_MODE);
2233
2234   rc = gfc_check_real_range (result->value.complex.r, kind);
2235
2236   if (rc == ARITH_UNDERFLOW)
2237     {
2238       if (gfc_option.warn_underflow)
2239         gfc_warning (gfc_arith_error (rc), &src->where);
2240       mpfr_set_ui (result->value.complex.r, 0, GFC_RND_MODE);
2241     }
2242   else if (rc != ARITH_OK)
2243     {
2244       arith_error (rc, &src->ts, &result->ts, &src->where);
2245       gfc_free_expr (result);
2246       return NULL;
2247     }
2248
2249   rc = gfc_check_real_range (result->value.complex.i, kind);
2250
2251   if (rc == ARITH_UNDERFLOW)
2252     {
2253       if (gfc_option.warn_underflow)
2254         gfc_warning (gfc_arith_error (rc), &src->where);
2255       mpfr_set_ui (result->value.complex.i, 0, GFC_RND_MODE);
2256     }
2257   else if (rc != ARITH_OK)
2258     {
2259       arith_error (rc, &src->ts, &result->ts, &src->where);
2260       gfc_free_expr (result);
2261       return NULL;
2262     }
2263
2264   return result;
2265 }
2266
2267
2268 /* Logical kind conversion.  */
2269
2270 gfc_expr *
2271 gfc_log2log (gfc_expr *src, int kind)
2272 {
2273   gfc_expr *result;
2274
2275   result = gfc_constant_result (BT_LOGICAL, kind, &src->where);
2276   result->value.logical = src->value.logical;
2277
2278   return result;
2279 }
2280
2281
2282 /* Convert logical to integer.  */
2283
2284 gfc_expr *
2285 gfc_log2int (gfc_expr *src, int kind)
2286 {
2287   gfc_expr *result;
2288
2289   result = gfc_constant_result (BT_INTEGER, kind, &src->where);
2290   mpz_set_si (result->value.integer, src->value.logical);
2291
2292   return result;
2293 }
2294
2295
2296 /* Convert integer to logical.  */
2297
2298 gfc_expr *
2299 gfc_int2log (gfc_expr *src, int kind)
2300 {
2301   gfc_expr *result;
2302
2303   result = gfc_constant_result (BT_LOGICAL, kind, &src->where);
2304   result->value.logical = (mpz_cmp_si (src->value.integer, 0) != 0);
2305
2306   return result;
2307 }
2308
2309
2310 /* Convert Hollerith to integer. The constant will be padded or truncated.  */
2311
2312 gfc_expr *
2313 gfc_hollerith2int (gfc_expr *src, int kind)
2314 {
2315   gfc_expr *result;
2316   int len;
2317
2318   len = src->value.character.length;
2319
2320   result = gfc_get_expr ();
2321   result->expr_type = EXPR_CONSTANT;
2322   result->ts.type = BT_INTEGER;
2323   result->ts.kind = kind;
2324   result->where = src->where;
2325   result->from_H = 1;
2326
2327   if (len > kind)
2328     {
2329       gfc_warning ("The Hollerith constant at %L is too long to convert to %s",
2330                    &src->where, gfc_typename(&result->ts));
2331     }
2332   result->value.character.string = gfc_getmem (kind + 1);
2333   memcpy (result->value.character.string, src->value.character.string,
2334         MIN (kind, len));
2335
2336   if (len < kind)
2337     memset (&result->value.character.string[len], ' ', kind - len);
2338
2339   result->value.character.string[kind] = '\0'; /* For debugger  */
2340   result->value.character.length = kind;
2341
2342   return result;
2343 }
2344
2345
2346 /* Convert Hollerith to real. The constant will be padded or truncated.  */
2347
2348 gfc_expr *
2349 gfc_hollerith2real (gfc_expr *src, int kind)
2350 {
2351   gfc_expr *result;
2352   int len;
2353
2354   len = src->value.character.length;
2355
2356   result = gfc_get_expr ();
2357   result->expr_type = EXPR_CONSTANT;
2358   result->ts.type = BT_REAL;
2359   result->ts.kind = kind;
2360   result->where = src->where;
2361   result->from_H = 1;
2362
2363   if (len > kind)
2364     {
2365       gfc_warning ("The Hollerith constant at %L is too long to convert to %s",
2366                    &src->where, gfc_typename(&result->ts));
2367     }
2368   result->value.character.string = gfc_getmem (kind + 1);
2369   memcpy (result->value.character.string, src->value.character.string,
2370         MIN (kind, len));
2371
2372   if (len < kind)
2373     memset (&result->value.character.string[len], ' ', kind - len);
2374
2375   result->value.character.string[kind] = '\0'; /* For debugger.  */
2376   result->value.character.length = kind;
2377
2378   return result;
2379 }
2380
2381
2382 /* Convert Hollerith to complex. The constant will be padded or truncated.  */
2383
2384 gfc_expr *
2385 gfc_hollerith2complex (gfc_expr *src, int kind)
2386 {
2387   gfc_expr *result;
2388   int len;
2389
2390   len = src->value.character.length;
2391
2392   result = gfc_get_expr ();
2393   result->expr_type = EXPR_CONSTANT;
2394   result->ts.type = BT_COMPLEX;
2395   result->ts.kind = kind;
2396   result->where = src->where;
2397   result->from_H = 1;
2398
2399   kind = kind * 2;
2400
2401   if (len > kind)
2402     {
2403       gfc_warning ("The Hollerith constant at %L is too long to convert to %s",
2404                    &src->where, gfc_typename(&result->ts));
2405     }
2406   result->value.character.string = gfc_getmem (kind + 1);
2407   memcpy (result->value.character.string, src->value.character.string,
2408           MIN (kind, len));
2409
2410   if (len < kind)
2411     memset (&result->value.character.string[len], ' ', kind - len);
2412
2413   result->value.character.string[kind] = '\0'; /* For debugger  */
2414   result->value.character.length = kind;
2415
2416   return result;
2417 }
2418
2419
2420 /* Convert Hollerith to character. */
2421
2422 gfc_expr *
2423 gfc_hollerith2character (gfc_expr *src, int kind)
2424 {
2425   gfc_expr *result;
2426
2427   result = gfc_copy_expr (src);
2428   result->ts.type = BT_CHARACTER;
2429   result->ts.kind = kind;
2430   result->from_H = 1;
2431
2432   return result;
2433 }
2434
2435
2436 /* Convert Hollerith to logical. The constant will be padded or truncated.  */
2437
2438 gfc_expr *
2439 gfc_hollerith2logical (gfc_expr *src, int kind)
2440 {
2441   gfc_expr *result;
2442   int len;
2443
2444   len = src->value.character.length;
2445
2446   result = gfc_get_expr ();
2447   result->expr_type = EXPR_CONSTANT;
2448   result->ts.type = BT_LOGICAL;
2449   result->ts.kind = kind;
2450   result->where = src->where;
2451   result->from_H = 1;
2452
2453   if (len > kind)
2454     {
2455       gfc_warning ("The Hollerith constant at %L is too long to convert to %s",
2456                    &src->where, gfc_typename(&result->ts));
2457     }
2458   result->value.character.string = gfc_getmem (kind + 1);
2459   memcpy (result->value.character.string, src->value.character.string,
2460         MIN (kind, len));
2461
2462   if (len < kind)
2463     memset (&result->value.character.string[len], ' ', kind - len);
2464
2465   result->value.character.string[kind] = '\0'; /* For debugger  */
2466   result->value.character.length = kind;
2467
2468   return result;
2469 }
2470
2471
2472 /* Returns an initializer whose value is one higher than the value of the
2473    LAST_INITIALIZER argument.  If the argument is NULL, the
2474    initializers value will be set to zero.  The initializer's kind
2475    will be set to gfc_c_int_kind.
2476
2477    If -fshort-enums is given, the appropriate kind will be selected
2478    later after all enumerators have been parsed.  A warning is issued
2479    here if an initializer exceeds gfc_c_int_kind.  */
2480
2481 gfc_expr *
2482 gfc_enum_initializer (gfc_expr *last_initializer, locus where)
2483 {
2484   gfc_expr *result;
2485
2486   result = gfc_get_expr ();
2487   result->expr_type = EXPR_CONSTANT;
2488   result->ts.type = BT_INTEGER;
2489   result->ts.kind = gfc_c_int_kind;
2490   result->where = where;
2491
2492   mpz_init (result->value.integer);
2493
2494   if (last_initializer != NULL)
2495     {
2496       mpz_add_ui (result->value.integer, last_initializer->value.integer, 1);
2497       result->where = last_initializer->where;
2498
2499       if (gfc_check_integer_range (result->value.integer,
2500              gfc_c_int_kind) != ARITH_OK)
2501         {
2502           gfc_error ("Enumerator exceeds the C integer type at %C");
2503           return NULL;
2504         }
2505     }
2506   else
2507     {
2508       /* Control comes here, if it's the very first enumerator and no
2509          initializer has been given.  It will be initialized to zero.  */
2510       mpz_set_si (result->value.integer, 0);
2511     }
2512
2513   return result;
2514 }