gcc/fortran/arith.c

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