gcc/fortran/arith.c

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