gcc/fortran/arith.c

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