gcc/fortran/arith.c

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