gcc/fortran/arith.c

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