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 "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.operator = operator;
1602
1603       temp.op1 = op1;
1604       temp.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->operator = operator;
1675
1676   result->op1 = op1;
1677   result->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   const char *t;
1932
1933   e = gfc_constant_result (BT_REAL, kind, where);
1934   /* A leading plus is allowed in Fortran, but not by mpfr_set_str */
1935   if (buffer[0] == '+')
1936     t = buffer + 1;
1937   else
1938     t = buffer;
1939   mpfr_set_str (e->value.real, t, 10, GFC_RND_MODE);
1940
1941   return e;
1942 }
1943
1944
1945 /* Convert a pair of real, constant expression nodes to a single
1946    complex expression node.  */
1947
1948 gfc_expr *
1949 gfc_convert_complex (gfc_expr * real, gfc_expr * imag, int kind)
1950 {
1951   gfc_expr *e;
1952
1953   e = gfc_constant_result (BT_COMPLEX, kind, &real->where);
1954   mpfr_set (e->value.complex.r, real->value.real, GFC_RND_MODE);
1955   mpfr_set (e->value.complex.i, imag->value.real, GFC_RND_MODE);
1956
1957   return e;
1958 }
1959
1960
1961 /******* Simplification of intrinsic functions with constant arguments *****/
1962
1963
1964 /* Deal with an arithmetic error.  */
1965
1966 static void
1967 arith_error (arith rc, gfc_typespec * from, gfc_typespec * to, locus * where)
1968 {
1969   gfc_error ("%s converting %s to %s at %L", gfc_arith_error (rc),
1970              gfc_typename (from), gfc_typename (to), where);
1971
1972   /* TODO: Do something about the error, ie, throw exception, return
1973      NaN, etc.  */
1974 }
1975
1976 /* Convert integers to integers.  */
1977
1978 gfc_expr *
1979 gfc_int2int (gfc_expr * src, int kind)
1980 {
1981   gfc_expr *result;
1982   arith rc;
1983
1984   result = gfc_constant_result (BT_INTEGER, kind, &src->where);
1985
1986   mpz_set (result->value.integer, src->value.integer);
1987
1988   if ((rc = gfc_check_integer_range (result->value.integer, kind))
1989       != ARITH_OK)
1990     {
1991       if (rc == ARITH_ASYMMETRIC)
1992         {
1993           gfc_warning ("%s at %L", gfc_arith_error (rc), &src->where);
1994         }
1995       else
1996         {
1997           arith_error (rc, &src->ts, &result->ts, &src->where);
1998           gfc_free_expr (result);
1999           return NULL;
2000         }
2001     }
2002
2003   return result;
2004 }
2005
2006
2007 /* Convert integers to reals.  */
2008
2009 gfc_expr *
2010 gfc_int2real (gfc_expr * src, int kind)
2011 {
2012   gfc_expr *result;
2013   arith rc;
2014
2015   result = gfc_constant_result (BT_REAL, kind, &src->where);
2016
2017   mpfr_set_z (result->value.real, src->value.integer, GFC_RND_MODE);
2018
2019   if ((rc = gfc_check_real_range (result->value.real, kind)) != ARITH_OK)
2020     {
2021       arith_error (rc, &src->ts, &result->ts, &src->where);
2022       gfc_free_expr (result);
2023       return NULL;
2024     }
2025
2026   return result;
2027 }
2028
2029
2030 /* Convert default integer to default complex.  */
2031
2032 gfc_expr *
2033 gfc_int2complex (gfc_expr * src, int kind)
2034 {
2035   gfc_expr *result;
2036   arith rc;
2037
2038   result = gfc_constant_result (BT_COMPLEX, kind, &src->where);
2039
2040   mpfr_set_z (result->value.complex.r, src->value.integer, GFC_RND_MODE);
2041   mpfr_set_ui (result->value.complex.i, 0, GFC_RND_MODE);
2042
2043   if ((rc = gfc_check_real_range (result->value.complex.r, kind)) != ARITH_OK)
2044     {
2045       arith_error (rc, &src->ts, &result->ts, &src->where);
2046       gfc_free_expr (result);
2047       return NULL;
2048     }
2049
2050   return result;
2051 }
2052
2053
2054 /* Convert default real to default integer.  */
2055
2056 gfc_expr *
2057 gfc_real2int (gfc_expr * src, int kind)
2058 {
2059   gfc_expr *result;
2060   arith rc;
2061
2062   result = gfc_constant_result (BT_INTEGER, kind, &src->where);
2063
2064   gfc_mpfr_to_mpz (result->value.integer, src->value.real);
2065
2066   if ((rc = gfc_check_integer_range (result->value.integer, kind))
2067       != 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 real to real.  */
2079
2080 gfc_expr *
2081 gfc_real2real (gfc_expr * src, int kind)
2082 {
2083   gfc_expr *result;
2084   arith rc;
2085
2086   result = gfc_constant_result (BT_REAL, kind, &src->where);
2087
2088   mpfr_set (result->value.real, src->value.real, GFC_RND_MODE);
2089
2090   rc = gfc_check_real_range (result->value.real, kind);
2091
2092   if (rc == ARITH_UNDERFLOW)
2093     {
2094       if (gfc_option.warn_underflow)
2095         gfc_warning ("%s at %L", gfc_arith_error (rc), &src->where);
2096       mpfr_set_ui (result->value.real, 0, GFC_RND_MODE);
2097     }
2098   else if (rc != ARITH_OK)
2099     {
2100       arith_error (rc, &src->ts, &result->ts, &src->where);
2101       gfc_free_expr (result);
2102       return NULL;
2103     }
2104
2105   return result;
2106 }
2107
2108
2109 /* Convert real to complex.  */
2110
2111 gfc_expr *
2112 gfc_real2complex (gfc_expr * src, int kind)
2113 {
2114   gfc_expr *result;
2115   arith rc;
2116
2117   result = gfc_constant_result (BT_COMPLEX, kind, &src->where);
2118
2119   mpfr_set (result->value.complex.r, src->value.real, GFC_RND_MODE);
2120   mpfr_set_ui (result->value.complex.i, 0, GFC_RND_MODE);
2121
2122   rc = gfc_check_real_range (result->value.complex.r, kind);
2123
2124   if (rc == ARITH_UNDERFLOW)
2125     {
2126       if (gfc_option.warn_underflow)
2127         gfc_warning ("%s at %L", gfc_arith_error (rc), &src->where);
2128       mpfr_set_ui (result->value.complex.r, 0, GFC_RND_MODE);
2129     }
2130   else if (rc != ARITH_OK)
2131     {
2132       arith_error (rc, &src->ts, &result->ts, &src->where);
2133       gfc_free_expr (result);
2134       return NULL;
2135     }
2136
2137   return result;
2138 }
2139
2140
2141 /* Convert complex to integer.  */
2142
2143 gfc_expr *
2144 gfc_complex2int (gfc_expr * src, int kind)
2145 {
2146   gfc_expr *result;
2147   arith rc;
2148
2149   result = gfc_constant_result (BT_INTEGER, kind, &src->where);
2150
2151   gfc_mpfr_to_mpz (result->value.integer, src->value.complex.r);
2152
2153   if ((rc = gfc_check_integer_range (result->value.integer, kind))
2154       != ARITH_OK)
2155     {
2156       arith_error (rc, &src->ts, &result->ts, &src->where);
2157       gfc_free_expr (result);
2158       return NULL;
2159     }
2160
2161   return result;
2162 }
2163
2164
2165 /* Convert complex to real.  */
2166
2167 gfc_expr *
2168 gfc_complex2real (gfc_expr * src, int kind)
2169 {
2170   gfc_expr *result;
2171   arith rc;
2172
2173   result = gfc_constant_result (BT_REAL, kind, &src->where);
2174
2175   mpfr_set (result->value.real, src->value.complex.r, GFC_RND_MODE);
2176
2177   rc = gfc_check_real_range (result->value.real, kind);
2178
2179   if (rc == ARITH_UNDERFLOW)
2180     {
2181       if (gfc_option.warn_underflow)
2182         gfc_warning ("%s at %L", gfc_arith_error (rc), &src->where);
2183       mpfr_set_ui (result->value.real, 0, GFC_RND_MODE);
2184     }
2185   if (rc != ARITH_OK)
2186     {
2187       arith_error (rc, &src->ts, &result->ts, &src->where);
2188       gfc_free_expr (result);
2189       return NULL;
2190     }
2191
2192   return result;
2193 }
2194
2195
2196 /* Convert complex to complex.  */
2197
2198 gfc_expr *
2199 gfc_complex2complex (gfc_expr * src, int kind)
2200 {
2201   gfc_expr *result;
2202   arith rc;
2203
2204   result = gfc_constant_result (BT_COMPLEX, kind, &src->where);
2205
2206   mpfr_set (result->value.complex.r, src->value.complex.r, GFC_RND_MODE);
2207   mpfr_set (result->value.complex.i, src->value.complex.i, GFC_RND_MODE);
2208
2209   rc = gfc_check_real_range (result->value.complex.r, kind);
2210
2211   if (rc == ARITH_UNDERFLOW)
2212     {
2213       if (gfc_option.warn_underflow)
2214         gfc_warning ("%s at %L", gfc_arith_error (rc), &src->where);
2215       mpfr_set_ui (result->value.complex.r, 0, GFC_RND_MODE);
2216     }
2217   else if (rc != ARITH_OK)
2218     {
2219       arith_error (rc, &src->ts, &result->ts, &src->where);
2220       gfc_free_expr (result);
2221       return NULL;
2222     }
2223
2224   rc = gfc_check_real_range (result->value.complex.i, kind);
2225
2226   if (rc == ARITH_UNDERFLOW)
2227     {
2228       if (gfc_option.warn_underflow)
2229         gfc_warning ("%s at %L", gfc_arith_error (rc), &src->where);
2230       mpfr_set_ui (result->value.complex.i, 0, GFC_RND_MODE);
2231     }
2232   else if (rc != ARITH_OK)
2233     {
2234       arith_error (rc, &src->ts, &result->ts, &src->where);
2235       gfc_free_expr (result);
2236       return NULL;
2237     }
2238
2239   return result;
2240 }
2241
2242
2243 /* Logical kind conversion.  */
2244
2245 gfc_expr *
2246 gfc_log2log (gfc_expr * src, int kind)
2247 {
2248   gfc_expr *result;
2249
2250   result = gfc_constant_result (BT_LOGICAL, kind, &src->where);
2251   result->value.logical = src->value.logical;
2252
2253   return result;
2254 }