gcc/fortran/arith.c

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