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