gcc/fortran/arith.c

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