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