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