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