gcc/fortran/arith.c

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