gcc/fortran/arith.c

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