gcc/fortran/arith.c

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