gcc/fortran/simplify.c

   1 /* Simplify intrinsic functions at compile-time.
   2    Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006 Free Software
   3    Foundation, Inc.
   4    Contributed by Andy Vaught & Katherine Holcomb
   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 #include "config.h"
  24 #include "system.h"
  25 #include "flags.h"
  26 #include "gfortran.h"
  27 #include "arith.h"
  28 #include "intrinsic.h"
  29
  30 gfc_expr gfc_bad_expr;
  31
  32
  33 /* Note that 'simplification' is not just transforming expressions.
  34    For functions that are not simplified at compile time, range
  35    checking is done if possible.
  36
  37    The return convention is that each simplification function returns:
  38
  39      A new expression node corresponding to the simplified arguments.
  40      The original arguments are destroyed by the caller, and must not
  41      be a part of the new expression.
  42
  43      NULL pointer indicating that no simplification was possible and
  44      the original expression should remain intact.  If the
  45      simplification function sets the type and/or the function name
  46      via the pointer gfc_simple_expression, then this type is
  47      retained.
  48
  49      An expression pointer to gfc_bad_expr (a static placeholder)
  50      indicating that some error has prevented simplification.  For
  51      example, sqrt(-1.0).  The error is generated within the function
  52      and should be propagated upwards
  53
  54    By the time a simplification function gets control, it has been
  55    decided that the function call is really supposed to be the
  56    intrinsic.  No type checking is strictly necessary, since only
  57    valid types will be passed on.  On the other hand, a simplification
  58    subroutine may have to look at the type of an argument as part of
  59    its processing.
  60
  61    Array arguments are never passed to these subroutines.
  62
  63    The functions in this file don't have much comment with them, but
  64    everything is reasonably straight-forward.  The Standard, chapter 13
  65    is the best comment you'll find for this file anyway.  */
  66
  67 /* Static table for converting non-ascii character sets to ascii.
  68    The xascii_table[] is the inverse table.  */
  69
  70 static int ascii_table[256] = {
  71   '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0',
  72   '\b', '\t', '\n', '\v', '\0', '\r', '\0', '\0',
  73   '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0',
  74   '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0',
  75   ' ', '!', '"', '#', '$', '%', '&', '\'',
  76   '(', ')', '*', '+', ',', '-', '.', '/',
  77   '0', '1', '2', '3', '4', '5', '6', '7',
  78   '8', '9', ':', ';', '<', '=', '>', '?',
  79   '@', 'A', 'B', 'C', 'D', 'E', 'F', 'G',
  80   'H', 'I', 'J', 'K', 'L', 'M', 'N', 'O',
  81   'P', 'Q', 'R', 'S', 'T', 'U', 'V', 'W',
  82   'X', 'Y', 'Z', '[', '\\', ']', '^', '_',
  83   '`', 'a', 'b', 'c', 'd', 'e', 'f', 'g',
  84   'h', 'i', 'j', 'k', 'l', 'm', 'n', 'o',
  85   'p', 'q', 'r', 's', 't', 'u', 'v', 'w',
  86   'x', 'y', 'z', '{', '|', '}', '~', '\?'
  87 };
  88
  89 static int xascii_table[256];
  90
  91
  92 /* Range checks an expression node.  If all goes well, returns the
  93    node, otherwise returns &gfc_bad_expr and frees the node.  */
  94
  95 static gfc_expr *
  96 range_check (gfc_expr * result, const char *name)
  97 {
  98
  99   switch (gfc_range_check (result))
 100     {
 101       case ARITH_OK:
 102         return result;
 103
 104       case ARITH_OVERFLOW:
 105         gfc_error ("Result of %s overflows its kind at %L", name, &result->where);
 106         break;
 107
 108       case ARITH_UNDERFLOW:
 109         gfc_error ("Result of %s underflows its kind at %L", name, &result->where);
 110         break;
 111
 112       case ARITH_NAN:
 113         gfc_error ("Result of %s is NaN at %L", name, &result->where);
 114         break;
 115
 116       default:
 117         gfc_error ("Result of %s gives range error for its kind at %L", name, &result->where);
 118         break;
 119     }
 120
 121   gfc_free_expr (result);
 122   return &gfc_bad_expr;
 123 }
 124
 125
 126 /* A helper function that gets an optional and possibly missing
 127    kind parameter.  Returns the kind, -1 if something went wrong.  */
 128
 129 static int
 130 get_kind (bt type, gfc_expr * k, const char *name, int default_kind)
 131 {
 132   int kind;
 133
 134   if (k == NULL)
 135     return default_kind;
 136
 137   if (k->expr_type != EXPR_CONSTANT)
 138     {
 139       gfc_error ("KIND parameter of %s at %L must be an initialization "
 140                  "expression", name, &k->where);
 141
 142       return -1;
 143     }
 144
 145   if (gfc_extract_int (k, &kind) != NULL
 146       || gfc_validate_kind (type, kind, true) < 0)
 147     {
 148
 149       gfc_error ("Invalid KIND parameter of %s at %L", name, &k->where);
 150       return -1;
 151     }
 152
 153   return kind;
 154 }
 155
 156
 157 /* Checks if X, which is assumed to represent a two's complement
 158    integer of binary width BITSIZE, has the signbit set.  If so, makes
 159    X the corresponding negative number.  */
 160
 161 static void
 162 twos_complement (mpz_t x, int bitsize)
 163 {
 164   mpz_t mask;
 165
 166   if (mpz_tstbit (x, bitsize - 1) == 1)
 167     {
 168       mpz_init_set_ui(mask, 1);
 169       mpz_mul_2exp(mask, mask, bitsize);
 170       mpz_sub_ui(mask, mask, 1);
 171
 172       /* We negate the number by hand, zeroing the high bits, that is
 173         make it the corresponding positive number, and then have it
 174         negated by GMP, giving the correct representation of the
 175         negative number.  */
 176       mpz_com (x, x);
 177       mpz_add_ui (x, x, 1);
 178       mpz_and (x, x, mask);
 179
 180       mpz_neg (x, x);
 181
 182       mpz_clear (mask);
 183     }
 184 }
 185
 186
 187 /********************** Simplification functions *****************************/
 188
 189 gfc_expr *
 190 gfc_simplify_abs (gfc_expr * e)
 191 {
 192   gfc_expr *result;
 193
 194   if (e->expr_type != EXPR_CONSTANT)
 195     return NULL;
 196
 197   switch (e->ts.type)
 198     {
 199     case BT_INTEGER:
 200       result = gfc_constant_result (BT_INTEGER, e->ts.kind, &e->where);
 201
 202       mpz_abs (result->value.integer, e->value.integer);
 203
 204       result = range_check (result, "IABS");
 205       break;
 206
 207     case BT_REAL:
 208       result = gfc_constant_result (BT_REAL, e->ts.kind, &e->where);
 209
 210       mpfr_abs (result->value.real, e->value.real, GFC_RND_MODE);
 211
 212       result = range_check (result, "ABS");
 213       break;
 214
 215     case BT_COMPLEX:
 216       result = gfc_constant_result (BT_REAL, e->ts.kind, &e->where);
 217
 218       gfc_set_model_kind (e->ts.kind);
 219
 220       mpfr_hypot (result->value.real, e->value.complex.r,
 221                   e->value.complex.i, GFC_RND_MODE);
 222       result = range_check (result, "CABS");
 223       break;
 224
 225     default:
 226       gfc_internal_error ("gfc_simplify_abs(): Bad type");
 227     }
 228
 229   return result;
 230 }
 231
 232
 233 gfc_expr *
 234 gfc_simplify_achar (gfc_expr * e)
 235 {
 236   gfc_expr *result;
 237   int index;
 238
 239   if (e->expr_type != EXPR_CONSTANT)
 240     return NULL;
 241
 242   /* We cannot assume that the native character set is ASCII in this
 243      function.  */
 244   if (gfc_extract_int (e, &index) != NULL || index < 0 || index > 127)
 245     {
 246       gfc_error ("Extended ASCII not implemented: argument of ACHAR at %L "
 247                  "must be between 0 and 127", &e->where);
 248       return &gfc_bad_expr;
 249     }
 250
 251   result = gfc_constant_result (BT_CHARACTER, gfc_default_character_kind,
 252                                 &e->where);
 253
 254   result->value.character.string = gfc_getmem (2);
 255
 256   result->value.character.length = 1;
 257   result->value.character.string[0] = ascii_table[index];
 258   result->value.character.string[1] = '\0';     /* For debugger */
 259   return result;
 260 }
 261
 262
 263 gfc_expr *
 264 gfc_simplify_acos (gfc_expr * x)
 265 {
 266   gfc_expr *result;
 267
 268   if (x->expr_type != EXPR_CONSTANT)
 269     return NULL;
 270
 271   if (mpfr_cmp_si (x->value.real, 1) > 0 || mpfr_cmp_si (x->value.real, -1) < 0)
 272     {
 273       gfc_error ("Argument of ACOS at %L must be between -1 and 1",
 274                  &x->where);
 275       return &gfc_bad_expr;
 276     }
 277
 278   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
 279
 280   mpfr_acos (result->value.real, x->value.real, GFC_RND_MODE);
 281
 282   return range_check (result, "ACOS");
 283 }
 284
 285 gfc_expr *
 286 gfc_simplify_acosh (gfc_expr * x)
 287 {
 288   gfc_expr *result;
 289
 290   if (x->expr_type != EXPR_CONSTANT)
 291     return NULL;
 292
 293   if (mpfr_cmp_si (x->value.real, 1) < 0)
 294     {
 295       gfc_error ("Argument of ACOSH at %L must not be less than 1",
 296                  &x->where);
 297       return &gfc_bad_expr;
 298     }
 299
 300   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
 301
 302   mpfr_acosh (result->value.real, x->value.real, GFC_RND_MODE);
 303
 304   return range_check (result, "ACOSH");
 305 }
 306
 307 gfc_expr *
 308 gfc_simplify_adjustl (gfc_expr * e)
 309 {
 310   gfc_expr *result;
 311   int count, i, len;
 312   char ch;
 313
 314   if (e->expr_type != EXPR_CONSTANT)
 315     return NULL;
 316
 317   len = e->value.character.length;
 318
 319   result = gfc_constant_result (BT_CHARACTER, e->ts.kind, &e->where);
 320
 321   result->value.character.length = len;
 322   result->value.character.string = gfc_getmem (len + 1);
 323
 324   for (count = 0, i = 0; i < len; ++i)
 325     {
 326       ch = e->value.character.string[i];
 327       if (ch != ' ')
 328         break;
 329       ++count;
 330     }
 331
 332   for (i = 0; i < len - count; ++i)
 333     {
 334       result->value.character.string[i] =
 335         e->value.character.string[count + i];
 336     }
 337
 338   for (i = len - count; i < len; ++i)
 339     {
 340       result->value.character.string[i] = ' ';
 341     }
 342
 343   result->value.character.string[len] = '\0';   /* For debugger */
 344
 345   return result;
 346 }
 347
 348
 349 gfc_expr *
 350 gfc_simplify_adjustr (gfc_expr * e)
 351 {
 352   gfc_expr *result;
 353   int count, i, len;
 354   char ch;
 355
 356   if (e->expr_type != EXPR_CONSTANT)
 357     return NULL;
 358
 359   len = e->value.character.length;
 360
 361   result = gfc_constant_result (BT_CHARACTER, e->ts.kind, &e->where);
 362
 363   result->value.character.length = len;
 364   result->value.character.string = gfc_getmem (len + 1);
 365
 366   for (count = 0, i = len - 1; i >= 0; --i)
 367     {
 368       ch = e->value.character.string[i];
 369       if (ch != ' ')
 370         break;
 371       ++count;
 372     }
 373
 374   for (i = 0; i < count; ++i)
 375     {
 376       result->value.character.string[i] = ' ';
 377     }
 378
 379   for (i = count; i < len; ++i)
 380     {
 381       result->value.character.string[i] =
 382         e->value.character.string[i - count];
 383     }
 384
 385   result->value.character.string[len] = '\0';   /* For debugger */
 386
 387   return result;
 388 }
 389
 390
 391 gfc_expr *
 392 gfc_simplify_aimag (gfc_expr * e)
 393 {
 394
 395   gfc_expr *result;
 396
 397   if (e->expr_type != EXPR_CONSTANT)
 398     return NULL;
 399
 400   result = gfc_constant_result (BT_REAL, e->ts.kind, &e->where);
 401   mpfr_set (result->value.real, e->value.complex.i, GFC_RND_MODE);
 402
 403   return range_check (result, "AIMAG");
 404 }
 405
 406
 407 gfc_expr *
 408 gfc_simplify_aint (gfc_expr * e, gfc_expr * k)
 409 {
 410   gfc_expr *rtrunc, *result;
 411   int kind;
 412
 413   kind = get_kind (BT_REAL, k, "AINT", e->ts.kind);
 414   if (kind == -1)
 415     return &gfc_bad_expr;
 416
 417   if (e->expr_type != EXPR_CONSTANT)
 418     return NULL;
 419
 420   rtrunc = gfc_copy_expr (e);
 421
 422   mpfr_trunc (rtrunc->value.real, e->value.real);
 423
 424   result = gfc_real2real (rtrunc, kind);
 425   gfc_free_expr (rtrunc);
 426
 427   return range_check (result, "AINT");
 428 }
 429
 430
 431 gfc_expr *
 432 gfc_simplify_dint (gfc_expr * e)
 433 {
 434   gfc_expr *rtrunc, *result;
 435
 436   if (e->expr_type != EXPR_CONSTANT)
 437     return NULL;
 438
 439   rtrunc = gfc_copy_expr (e);
 440
 441   mpfr_trunc (rtrunc->value.real, e->value.real);
 442
 443   result = gfc_real2real (rtrunc, gfc_default_double_kind);
 444   gfc_free_expr (rtrunc);
 445
 446   return range_check (result, "DINT");
 447 }
 448
 449
 450 gfc_expr *
 451 gfc_simplify_anint (gfc_expr * e, gfc_expr * k)
 452 {
 453   gfc_expr *result;
 454   int kind;
 455
 456   kind = get_kind (BT_REAL, k, "ANINT", e->ts.kind);
 457   if (kind == -1)
 458     return &gfc_bad_expr;
 459
 460   if (e->expr_type != EXPR_CONSTANT)
 461     return NULL;
 462
 463   result = gfc_constant_result (e->ts.type, kind, &e->where);
 464
 465   mpfr_round (result->value.real, e->value.real);
 466
 467   return range_check (result, "ANINT");
 468 }
 469
 470
 471 gfc_expr *
 472 gfc_simplify_and (gfc_expr * x, gfc_expr * y)
 473 {
 474   gfc_expr *result;
 475   int kind;
 476
 477   if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
 478     return NULL;
 479
 480   kind = x->ts.kind > y->ts.kind ? x->ts.kind : y->ts.kind;
 481   if (x->ts.type == BT_INTEGER)
 482     {
 483       result = gfc_constant_result (BT_INTEGER, kind, &x->where);
 484       mpz_and (result->value.integer, x->value.integer, y->value.integer);
 485     }
 486   else /* BT_LOGICAL */
 487     {
 488       result = gfc_constant_result (BT_LOGICAL, kind, &x->where);
 489       result->value.logical = x->value.logical && y->value.logical;
 490     }
 491
 492   return range_check (result, "AND");
 493 }
 494
 495
 496 gfc_expr *
 497 gfc_simplify_dnint (gfc_expr * e)
 498 {
 499   gfc_expr *result;
 500
 501   if (e->expr_type != EXPR_CONSTANT)
 502     return NULL;
 503
 504   result = gfc_constant_result (BT_REAL, gfc_default_double_kind, &e->where);
 505
 506   mpfr_round (result->value.real, e->value.real);
 507
 508   return range_check (result, "DNINT");
 509 }
 510
 511
 512 gfc_expr *
 513 gfc_simplify_asin (gfc_expr * x)
 514 {
 515   gfc_expr *result;
 516
 517   if (x->expr_type != EXPR_CONSTANT)
 518     return NULL;
 519
 520   if (mpfr_cmp_si (x->value.real, 1) > 0 || mpfr_cmp_si (x->value.real, -1) < 0)
 521     {
 522       gfc_error ("Argument of ASIN at %L must be between -1 and 1",
 523                  &x->where);
 524       return &gfc_bad_expr;
 525     }
 526
 527   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
 528
 529   mpfr_asin(result->value.real, x->value.real, GFC_RND_MODE);
 530
 531   return range_check (result, "ASIN");
 532 }
 533
 534
 535 gfc_expr *
 536 gfc_simplify_asinh (gfc_expr * x)
 537 {
 538   gfc_expr *result;
 539
 540   if (x->expr_type != EXPR_CONSTANT)
 541     return NULL;
 542
 543   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
 544
 545   mpfr_asinh(result->value.real, x->value.real, GFC_RND_MODE);
 546
 547   return range_check (result, "ASINH");
 548 }
 549
 550
 551 gfc_expr *
 552 gfc_simplify_atan (gfc_expr * x)
 553 {
 554   gfc_expr *result;
 555
 556   if (x->expr_type != EXPR_CONSTANT)
 557     return NULL;
 558
 559   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
 560
 561   mpfr_atan(result->value.real, x->value.real, GFC_RND_MODE);
 562
 563   return range_check (result, "ATAN");
 564 }
 565
 566
 567 gfc_expr *
 568 gfc_simplify_atanh (gfc_expr * x)
 569 {
 570   gfc_expr *result;
 571
 572   if (x->expr_type != EXPR_CONSTANT)
 573     return NULL;
 574
 575   if (mpfr_cmp_si (x->value.real, 1) >= 0 ||
 576       mpfr_cmp_si (x->value.real, -1) <= 0)
 577     {
 578       gfc_error ("Argument of ATANH at %L must be inside the range -1 to 1",
 579                  &x->where);
 580       return &gfc_bad_expr;
 581     }
 582
 583   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
 584
 585   mpfr_atanh(result->value.real, x->value.real, GFC_RND_MODE);
 586
 587   return range_check (result, "ATANH");
 588 }
 589
 590
 591 gfc_expr *
 592 gfc_simplify_atan2 (gfc_expr * y, gfc_expr * x)
 593 {
 594   gfc_expr *result;
 595
 596   if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
 597     return NULL;
 598
 599   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
 600
 601   if (mpfr_sgn (y->value.real) == 0 && mpfr_sgn (x->value.real) == 0)
 602     {
 603       gfc_error
 604         ("If first argument of ATAN2 %L is zero, then the second argument "
 605           "must not be zero", &x->where);
 606       gfc_free_expr (result);
 607       return &gfc_bad_expr;
 608     }
 609
 610   arctangent2 (y->value.real, x->value.real, result->value.real);
 611
 612   return range_check (result, "ATAN2");
 613 }
 614
 615
 616 gfc_expr *
 617 gfc_simplify_bit_size (gfc_expr * e)
 618 {
 619   gfc_expr *result;
 620   int i;
 621
 622   i = gfc_validate_kind (e->ts.type, e->ts.kind, false);
 623   result = gfc_constant_result (BT_INTEGER, e->ts.kind, &e->where);
 624   mpz_set_ui (result->value.integer, gfc_integer_kinds[i].bit_size);
 625
 626   return result;
 627 }
 628
 629
 630 gfc_expr *
 631 gfc_simplify_btest (gfc_expr * e, gfc_expr * bit)
 632 {
 633   int b;
 634
 635   if (e->expr_type != EXPR_CONSTANT || bit->expr_type != EXPR_CONSTANT)
 636     return NULL;
 637
 638   if (gfc_extract_int (bit, &b) != NULL || b < 0)
 639     return gfc_logical_expr (0, &e->where);
 640
 641   return gfc_logical_expr (mpz_tstbit (e->value.integer, b), &e->where);
 642 }
 643
 644
 645 gfc_expr *
 646 gfc_simplify_ceiling (gfc_expr * e, gfc_expr * k)
 647 {
 648   gfc_expr *ceil, *result;
 649   int kind;
 650
 651   kind = get_kind (BT_INTEGER, k, "CEILING", gfc_default_integer_kind);
 652   if (kind == -1)
 653     return &gfc_bad_expr;
 654
 655   if (e->expr_type != EXPR_CONSTANT)
 656     return NULL;
 657
 658   result = gfc_constant_result (BT_INTEGER, kind, &e->where);
 659
 660   ceil = gfc_copy_expr (e);
 661
 662   mpfr_ceil (ceil->value.real, e->value.real);
 663   gfc_mpfr_to_mpz(result->value.integer, ceil->value.real);
 664
 665   gfc_free_expr (ceil);
 666
 667   return range_check (result, "CEILING");
 668 }
 669
 670
 671 gfc_expr *
 672 gfc_simplify_char (gfc_expr * e, gfc_expr * k)
 673 {
 674   gfc_expr *result;
 675   int c, kind;
 676
 677   kind = get_kind (BT_CHARACTER, k, "CHAR", gfc_default_character_kind);
 678   if (kind == -1)
 679     return &gfc_bad_expr;
 680
 681   if (e->expr_type != EXPR_CONSTANT)
 682     return NULL;
 683
 684   if (gfc_extract_int (e, &c) != NULL || c < 0 || c > UCHAR_MAX)
 685     {
 686       gfc_error ("Bad character in CHAR function at %L", &e->where);
 687       return &gfc_bad_expr;
 688     }
 689
 690   result = gfc_constant_result (BT_CHARACTER, kind, &e->where);
 691
 692   result->value.character.length = 1;
 693   result->value.character.string = gfc_getmem (2);
 694
 695   result->value.character.string[0] = c;
 696   result->value.character.string[1] = '\0';     /* For debugger */
 697
 698   return result;
 699 }
 700
 701
 702 /* Common subroutine for simplifying CMPLX and DCMPLX.  */
 703
 704 static gfc_expr *
 705 simplify_cmplx (const char *name, gfc_expr * x, gfc_expr * y, int kind)
 706 {
 707   gfc_expr *result;
 708
 709   result = gfc_constant_result (BT_COMPLEX, kind, &x->where);
 710
 711   mpfr_set_ui (result->value.complex.i, 0, GFC_RND_MODE);
 712
 713   switch (x->ts.type)
 714     {
 715     case BT_INTEGER:
 716       mpfr_set_z (result->value.complex.r, x->value.integer, GFC_RND_MODE);
 717       break;
 718
 719     case BT_REAL:
 720       mpfr_set (result->value.complex.r, x->value.real, GFC_RND_MODE);
 721       break;
 722
 723     case BT_COMPLEX:
 724       mpfr_set (result->value.complex.r, x->value.complex.r, GFC_RND_MODE);
 725       mpfr_set (result->value.complex.i, x->value.complex.i, GFC_RND_MODE);
 726       break;
 727
 728     default:
 729       gfc_internal_error ("gfc_simplify_dcmplx(): Bad type (x)");
 730     }
 731
 732   if (y != NULL)
 733     {
 734       switch (y->ts.type)
 735         {
 736         case BT_INTEGER:
 737           mpfr_set_z (result->value.complex.i, y->value.integer, GFC_RND_MODE);
 738           break;
 739
 740         case BT_REAL:
 741           mpfr_set (result->value.complex.i, y->value.real, GFC_RND_MODE);
 742           break;
 743
 744         default:
 745           gfc_internal_error ("gfc_simplify_dcmplx(): Bad type (y)");
 746         }
 747     }
 748
 749   return range_check (result, name);
 750 }
 751
 752
 753 gfc_expr *
 754 gfc_simplify_cmplx (gfc_expr * x, gfc_expr * y, gfc_expr * k)
 755 {
 756   int kind;
 757
 758   if (x->expr_type != EXPR_CONSTANT
 759       || (y != NULL && y->expr_type != EXPR_CONSTANT))
 760     return NULL;
 761
 762   kind = get_kind (BT_REAL, k, "CMPLX", gfc_default_real_kind);
 763   if (kind == -1)
 764     return &gfc_bad_expr;
 765
 766   return simplify_cmplx ("CMPLX", x, y, kind);
 767 }
 768
 769
 770 gfc_expr *
 771 gfc_simplify_complex (gfc_expr * x, gfc_expr * y)
 772 {
 773   int kind;
 774
 775   if (x->expr_type != EXPR_CONSTANT
 776       || (y != NULL && y->expr_type != EXPR_CONSTANT))
 777     return NULL;
 778
 779   if (x->ts.type == BT_INTEGER)
 780     {
 781       if (y->ts.type == BT_INTEGER)
 782         kind = gfc_default_real_kind;
 783       else
 784         kind = y->ts.kind;
 785     }
 786   else
 787     {
 788       if (y->ts.type == BT_REAL)
 789         kind = (x->ts.kind > y->ts.kind) ? x->ts.kind : y->ts.kind;
 790       else
 791         kind = x->ts.kind;
 792     }
 793
 794   return simplify_cmplx ("COMPLEX", x, y, kind);
 795 }
 796
 797
 798 gfc_expr *
 799 gfc_simplify_conjg (gfc_expr * e)
 800 {
 801   gfc_expr *result;
 802
 803   if (e->expr_type != EXPR_CONSTANT)
 804     return NULL;
 805
 806   result = gfc_copy_expr (e);
 807   mpfr_neg (result->value.complex.i, result->value.complex.i, GFC_RND_MODE);
 808
 809   return range_check (result, "CONJG");
 810 }
 811
 812
 813 gfc_expr *
 814 gfc_simplify_cos (gfc_expr * x)
 815 {
 816   gfc_expr *result;
 817   mpfr_t xp, xq;
 818
 819   if (x->expr_type != EXPR_CONSTANT)
 820     return NULL;
 821
 822   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
 823
 824   switch (x->ts.type)
 825     {
 826     case BT_REAL:
 827       mpfr_cos (result->value.real, x->value.real, GFC_RND_MODE);
 828       break;
 829     case BT_COMPLEX:
 830       gfc_set_model_kind (x->ts.kind);
 831       mpfr_init (xp);
 832       mpfr_init (xq);
 833
 834       mpfr_cos  (xp, x->value.complex.r, GFC_RND_MODE);
 835       mpfr_cosh (xq, x->value.complex.i, GFC_RND_MODE);
 836       mpfr_mul(result->value.complex.r, xp, xq, GFC_RND_MODE);
 837
 838       mpfr_sin  (xp, x->value.complex.r, GFC_RND_MODE);
 839       mpfr_sinh (xq, x->value.complex.i, GFC_RND_MODE);
 840       mpfr_mul (xp, xp, xq, GFC_RND_MODE);
 841       mpfr_neg (result->value.complex.i, xp, GFC_RND_MODE );
 842
 843       mpfr_clear (xp);
 844       mpfr_clear (xq);
 845       break;
 846     default:
 847       gfc_internal_error ("in gfc_simplify_cos(): Bad type");
 848     }
 849
 850   return range_check (result, "COS");
 851
 852 }
 853
 854
 855 gfc_expr *
 856 gfc_simplify_cosh (gfc_expr * x)
 857 {
 858   gfc_expr *result;
 859
 860   if (x->expr_type != EXPR_CONSTANT)
 861     return NULL;
 862
 863   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
 864
 865   mpfr_cosh (result->value.real, x->value.real, GFC_RND_MODE);
 866
 867   return range_check (result, "COSH");
 868 }
 869
 870
 871 gfc_expr *
 872 gfc_simplify_dcmplx (gfc_expr * x, gfc_expr * y)
 873 {
 874
 875   if (x->expr_type != EXPR_CONSTANT
 876       || (y != NULL && y->expr_type != EXPR_CONSTANT))
 877     return NULL;
 878
 879   return simplify_cmplx ("DCMPLX", x, y, gfc_default_double_kind);
 880 }
 881
 882
 883 gfc_expr *
 884 gfc_simplify_dble (gfc_expr * e)
 885 {
 886   gfc_expr *result;
 887
 888   if (e->expr_type != EXPR_CONSTANT)
 889     return NULL;
 890
 891   switch (e->ts.type)
 892     {
 893     case BT_INTEGER:
 894       result = gfc_int2real (e, gfc_default_double_kind);
 895       break;
 896
 897     case BT_REAL:
 898       result = gfc_real2real (e, gfc_default_double_kind);
 899       break;
 900
 901     case BT_COMPLEX:
 902       result = gfc_complex2real (e, gfc_default_double_kind);
 903       break;
 904
 905     default:
 906       gfc_internal_error ("gfc_simplify_dble(): bad type at %L", &e->where);
 907     }
 908
 909   return range_check (result, "DBLE");
 910 }
 911
 912
 913 gfc_expr *
 914 gfc_simplify_digits (gfc_expr * x)
 915 {
 916   int i, digits;
 917
 918   i = gfc_validate_kind (x->ts.type, x->ts.kind, false);
 919   switch (x->ts.type)
 920     {
 921     case BT_INTEGER:
 922       digits = gfc_integer_kinds[i].digits;
 923       break;
 924
 925     case BT_REAL:
 926     case BT_COMPLEX:
 927       digits = gfc_real_kinds[i].digits;
 928       break;
 929
 930     default:
 931       gcc_unreachable ();
 932     }
 933
 934   return gfc_int_expr (digits);
 935 }
 936
 937
 938 gfc_expr *
 939 gfc_simplify_dim (gfc_expr * x, gfc_expr * y)
 940 {
 941   gfc_expr *result;
 942   int kind;
 943
 944   if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
 945     return NULL;
 946
 947   kind = x->ts.kind > y->ts.kind ? x->ts.kind : y->ts.kind;
 948   result = gfc_constant_result (x->ts.type, kind, &x->where);
 949
 950   switch (x->ts.type)
 951     {
 952     case BT_INTEGER:
 953       if (mpz_cmp (x->value.integer, y->value.integer) > 0)
 954         mpz_sub (result->value.integer, x->value.integer, y->value.integer);
 955       else
 956         mpz_set_ui (result->value.integer, 0);
 957
 958       break;
 959
 960     case BT_REAL:
 961       if (mpfr_cmp (x->value.real, y->value.real) > 0)
 962         mpfr_sub (result->value.real, x->value.real, y->value.real, GFC_RND_MODE);
 963       else
 964         mpfr_set_ui (result->value.real, 0, GFC_RND_MODE);
 965
 966       break;
 967
 968     default:
 969       gfc_internal_error ("gfc_simplify_dim(): Bad type");
 970     }
 971
 972   return range_check (result, "DIM");
 973 }
 974
 975
 976 gfc_expr *
 977 gfc_simplify_dprod (gfc_expr * x, gfc_expr * y)
 978 {
 979   gfc_expr *a1, *a2, *result;
 980
 981   if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
 982     return NULL;
 983
 984   result =
 985     gfc_constant_result (BT_REAL, gfc_default_double_kind, &x->where);
 986
 987   a1 = gfc_real2real (x, gfc_default_double_kind);
 988   a2 = gfc_real2real (y, gfc_default_double_kind);
 989
 990   mpfr_mul (result->value.real, a1->value.real, a2->value.real, GFC_RND_MODE);
 991
 992   gfc_free_expr (a1);
 993   gfc_free_expr (a2);
 994
 995   return range_check (result, "DPROD");
 996 }
 997
 998
 999 gfc_expr *
1000 gfc_simplify_epsilon (gfc_expr * e)
1001 {
1002   gfc_expr *result;
1003   int i;
1004
1005   i = gfc_validate_kind (e->ts.type, e->ts.kind, false);
1006
1007   result = gfc_constant_result (BT_REAL, e->ts.kind, &e->where);
1008
1009   mpfr_set (result->value.real, gfc_real_kinds[i].epsilon, GFC_RND_MODE);
1010
1011   return range_check (result, "EPSILON");
1012 }
1013
1014
1015 gfc_expr *
1016 gfc_simplify_exp (gfc_expr * x)
1017 {
1018   gfc_expr *result;
1019   mpfr_t xp, xq;
1020
1021   if (x->expr_type != EXPR_CONSTANT)
1022     return NULL;
1023
1024   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
1025
1026   switch (x->ts.type)
1027     {
1028     case BT_REAL:
1029       mpfr_exp(result->value.real, x->value.real, GFC_RND_MODE);
1030       break;
1031
1032     case BT_COMPLEX:
1033       gfc_set_model_kind (x->ts.kind);
1034       mpfr_init (xp);
1035       mpfr_init (xq);
1036       mpfr_exp (xq, x->value.complex.r, GFC_RND_MODE);
1037       mpfr_cos (xp, x->value.complex.i, GFC_RND_MODE);
1038       mpfr_mul (result->value.complex.r, xq, xp, GFC_RND_MODE);
1039       mpfr_sin (xp, x->value.complex.i, GFC_RND_MODE);
1040       mpfr_mul (result->value.complex.i, xq, xp, GFC_RND_MODE);
1041       mpfr_clear (xp);
1042       mpfr_clear (xq);
1043       break;
1044
1045     default:
1046       gfc_internal_error ("in gfc_simplify_exp(): Bad type");
1047     }
1048
1049   return range_check (result, "EXP");
1050 }
1051
1052 /* FIXME:  MPFR should be able to do this better */
1053 gfc_expr *
1054 gfc_simplify_exponent (gfc_expr * x)
1055 {
1056   int i;
1057   mpfr_t tmp;
1058   gfc_expr *result;
1059
1060   if (x->expr_type != EXPR_CONSTANT)
1061     return NULL;
1062
1063   result = gfc_constant_result (BT_INTEGER, gfc_default_integer_kind,
1064                                 &x->where);
1065
1066   gfc_set_model (x->value.real);
1067
1068   if (mpfr_sgn (x->value.real) == 0)
1069     {
1070       mpz_set_ui (result->value.integer, 0);
1071       return result;
1072     }
1073
1074   mpfr_init (tmp);
1075
1076   mpfr_abs (tmp, x->value.real, GFC_RND_MODE);
1077   mpfr_log2 (tmp, tmp, GFC_RND_MODE);
1078
1079   gfc_mpfr_to_mpz (result->value.integer, tmp);
1080
1081   /* The model number for tiny(x) is b**(emin - 1) where b is the base and emin
1082      is the smallest exponent value.  So, we need to add 1 if x is tiny(x).  */
1083   i = gfc_validate_kind (x->ts.type, x->ts.kind, false);
1084   if (mpfr_cmp (x->value.real, gfc_real_kinds[i].tiny) == 0)
1085     mpz_add_ui (result->value.integer,result->value.integer, 1);
1086
1087   mpfr_clear (tmp);
1088
1089   return range_check (result, "EXPONENT");
1090 }
1091
1092
1093 gfc_expr *
1094 gfc_simplify_float (gfc_expr * a)
1095 {
1096   gfc_expr *result;
1097
1098   if (a->expr_type != EXPR_CONSTANT)
1099     return NULL;
1100
1101   result = gfc_int2real (a, gfc_default_real_kind);
1102   return range_check (result, "FLOAT");
1103 }
1104
1105
1106 gfc_expr *
1107 gfc_simplify_floor (gfc_expr * e, gfc_expr * k)
1108 {
1109   gfc_expr *result;
1110   mpfr_t floor;
1111   int kind;
1112
1113   kind = get_kind (BT_INTEGER, k, "FLOOR", gfc_default_integer_kind);
1114   if (kind == -1)
1115     gfc_internal_error ("gfc_simplify_floor(): Bad kind");
1116
1117   if (e->expr_type != EXPR_CONSTANT)
1118     return NULL;
1119
1120   result = gfc_constant_result (BT_INTEGER, kind, &e->where);
1121
1122   gfc_set_model_kind (kind);
1123   mpfr_init (floor);
1124   mpfr_floor (floor, e->value.real);
1125
1126   gfc_mpfr_to_mpz (result->value.integer, floor);
1127
1128   mpfr_clear (floor);
1129
1130   return range_check (result, "FLOOR");
1131 }
1132
1133
1134 gfc_expr *
1135 gfc_simplify_fraction (gfc_expr * x)
1136 {
1137   gfc_expr *result;
1138   mpfr_t absv, exp, pow2;
1139
1140   if (x->expr_type != EXPR_CONSTANT)
1141     return NULL;
1142
1143   result = gfc_constant_result (BT_REAL, x->ts.kind, &x->where);
1144
1145   gfc_set_model_kind (x->ts.kind);
1146
1147   if (mpfr_sgn (x->value.real) == 0)
1148     {
1149       mpfr_set_ui (result->value.real, 0, GFC_RND_MODE);
1150       return result;
1151     }
1152
1153   mpfr_init (exp);
1154   mpfr_init (absv);
1155   mpfr_init (pow2);
1156
1157   mpfr_abs (absv, x->value.real, GFC_RND_MODE);
1158   mpfr_log2 (exp, absv, GFC_RND_MODE);
1159
1160   mpfr_trunc (exp, exp);
1161   mpfr_add_ui (exp, exp, 1, GFC_RND_MODE);
1162
1163   mpfr_ui_pow (pow2, 2, exp, GFC_RND_MODE);
1164
1165   mpfr_div (result->value.real, absv, pow2, GFC_RND_MODE);
1166
1167   mpfr_clear (exp);
1168   mpfr_clear (absv);
1169   mpfr_clear (pow2);
1170
1171   return range_check (result, "FRACTION");
1172 }
1173
1174
1175 gfc_expr *
1176 gfc_simplify_huge (gfc_expr * e)
1177 {
1178   gfc_expr *result;
1179   int i;
1180
1181   i = gfc_validate_kind (e->ts.type, e->ts.kind, false);
1182
1183   result = gfc_constant_result (e->ts.type, e->ts.kind, &e->where);
1184
1185   switch (e->ts.type)
1186     {
1187     case BT_INTEGER:
1188       mpz_set (result->value.integer, gfc_integer_kinds[i].huge);
1189       break;
1190
1191     case BT_REAL:
1192       mpfr_set (result->value.real, gfc_real_kinds[i].huge, GFC_RND_MODE);
1193       break;
1194
1195     default:
1196       gcc_unreachable ();
1197     }
1198
1199   return result;
1200 }
1201
1202
1203 gfc_expr *
1204 gfc_simplify_iachar (gfc_expr * e)
1205 {
1206   gfc_expr *result;
1207   int index;
1208
1209   if (e->expr_type != EXPR_CONSTANT)
1210     return NULL;
1211
1212   if (e->value.character.length != 1)
1213     {
1214       gfc_error ("Argument of IACHAR at %L must be of length one", &e->where);
1215       return &gfc_bad_expr;
1216     }
1217
1218   index = xascii_table[(int) e->value.character.string[0] & 0xFF];
1219
1220   result = gfc_int_expr (index);
1221   result->where = e->where;
1222
1223   return range_check (result, "IACHAR");
1224 }
1225
1226
1227 gfc_expr *
1228 gfc_simplify_iand (gfc_expr * x, gfc_expr * y)
1229 {
1230   gfc_expr *result;
1231
1232   if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
1233     return NULL;
1234
1235   result = gfc_constant_result (BT_INTEGER, x->ts.kind, &x->where);
1236
1237   mpz_and (result->value.integer, x->value.integer, y->value.integer);
1238
1239   return range_check (result, "IAND");
1240 }
1241
1242
1243 gfc_expr *
1244 gfc_simplify_ibclr (gfc_expr * x, gfc_expr * y)
1245 {
1246   gfc_expr *result;
1247   int k, pos;
1248
1249   if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
1250     return NULL;
1251
1252   if (gfc_extract_int (y, &pos) != NULL || pos < 0)
1253     {
1254       gfc_error ("Invalid second argument of IBCLR at %L", &y->where);
1255       return &gfc_bad_expr;
1256     }
1257
1258   k = gfc_validate_kind (x->ts.type, x->ts.kind, false);
1259
1260   if (pos > gfc_integer_kinds[k].bit_size)
1261     {
1262       gfc_error ("Second argument of IBCLR exceeds bit size at %L",
1263                  &y->where);
1264       return &gfc_bad_expr;
1265     }
1266
1267   result = gfc_copy_expr (x);
1268
1269   mpz_clrbit (result->value.integer, pos);
1270   return range_check (result, "IBCLR");
1271 }
1272
1273
1274 gfc_expr *
1275 gfc_simplify_ibits (gfc_expr * x, gfc_expr * y, gfc_expr * z)
1276 {
1277   gfc_expr *result;
1278   int pos, len;
1279   int i, k, bitsize;
1280   int *bits;
1281
1282   if (x->expr_type != EXPR_CONSTANT
1283       || y->expr_type != EXPR_CONSTANT
1284       || z->expr_type != EXPR_CONSTANT)
1285     return NULL;
1286
1287   if (gfc_extract_int (y, &pos) != NULL || pos < 0)
1288     {
1289       gfc_error ("Invalid second argument of IBITS at %L", &y->where);
1290       return &gfc_bad_expr;
1291     }
1292
1293   if (gfc_extract_int (z, &len) != NULL || len < 0)
1294     {
1295       gfc_error ("Invalid third argument of IBITS at %L", &z->where);
1296       return &gfc_bad_expr;
1297     }
1298
1299   k = gfc_validate_kind (BT_INTEGER, x->ts.kind, false);
1300
1301   bitsize = gfc_integer_kinds[k].bit_size;
1302
1303   if (pos + len > bitsize)
1304     {
1305       gfc_error
1306         ("Sum of second and third arguments of IBITS exceeds bit size "
1307          "at %L", &y->where);
1308       return &gfc_bad_expr;
1309     }
1310
1311   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
1312
1313   bits = gfc_getmem (bitsize * sizeof (int));
1314
1315   for (i = 0; i < bitsize; i++)
1316     bits[i] = 0;
1317
1318   for (i = 0; i < len; i++)
1319     bits[i] = mpz_tstbit (x->value.integer, i + pos);
1320
1321   for (i = 0; i < bitsize; i++)
1322     {
1323       if (bits[i] == 0)
1324         {
1325           mpz_clrbit (result->value.integer, i);
1326         }
1327       else if (bits[i] == 1)
1328         {
1329           mpz_setbit (result->value.integer, i);
1330         }
1331       else
1332         {
1333           gfc_internal_error ("IBITS: Bad bit");
1334         }
1335     }
1336
1337   gfc_free (bits);
1338
1339   return range_check (result, "IBITS");
1340 }
1341
1342
1343 gfc_expr *
1344 gfc_simplify_ibset (gfc_expr * x, gfc_expr * y)
1345 {
1346   gfc_expr *result;
1347   int k, pos;
1348
1349   if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
1350     return NULL;
1351
1352   if (gfc_extract_int (y, &pos) != NULL || pos < 0)
1353     {
1354       gfc_error ("Invalid second argument of IBSET at %L", &y->where);
1355       return &gfc_bad_expr;
1356     }
1357
1358   k = gfc_validate_kind (x->ts.type, x->ts.kind, false);
1359
1360   if (pos > gfc_integer_kinds[k].bit_size)
1361     {
1362       gfc_error ("Second argument of IBSET exceeds bit size at %L",
1363                  &y->where);
1364       return &gfc_bad_expr;
1365     }
1366
1367   result = gfc_copy_expr (x);
1368
1369   mpz_setbit (result->value.integer, pos);
1370
1371   twos_complement (result->value.integer, gfc_integer_kinds[k].bit_size);
1372
1373   return range_check (result, "IBSET");
1374 }
1375
1376
1377 gfc_expr *
1378 gfc_simplify_ichar (gfc_expr * e)
1379 {
1380   gfc_expr *result;
1381   int index;
1382
1383   if (e->expr_type != EXPR_CONSTANT)
1384     return NULL;
1385
1386   if (e->value.character.length != 1)
1387     {
1388       gfc_error ("Argument of ICHAR at %L must be of length one", &e->where);
1389       return &gfc_bad_expr;
1390     }
1391
1392   index = (unsigned char) e->value.character.string[0];
1393
1394   if (index < 0 || index > UCHAR_MAX)
1395     {
1396       gfc_error ("Argument of ICHAR at %L out of range of this processor",
1397                  &e->where);
1398       return &gfc_bad_expr;
1399     }
1400
1401   result = gfc_int_expr (index);
1402   result->where = e->where;
1403   return range_check (result, "ICHAR");
1404 }
1405
1406
1407 gfc_expr *
1408 gfc_simplify_ieor (gfc_expr * x, gfc_expr * y)
1409 {
1410   gfc_expr *result;
1411
1412   if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
1413     return NULL;
1414
1415   result = gfc_constant_result (BT_INTEGER, x->ts.kind, &x->where);
1416
1417   mpz_xor (result->value.integer, x->value.integer, y->value.integer);
1418
1419   return range_check (result, "IEOR");
1420 }
1421
1422
1423 gfc_expr *
1424 gfc_simplify_index (gfc_expr * x, gfc_expr * y, gfc_expr * b)
1425 {
1426   gfc_expr *result;
1427   int back, len, lensub;
1428   int i, j, k, count, index = 0, start;
1429
1430   if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
1431     return NULL;
1432
1433   if (b != NULL && b->value.logical != 0)
1434     back = 1;
1435   else
1436     back = 0;
1437
1438   result = gfc_constant_result (BT_INTEGER, gfc_default_integer_kind,
1439                                 &x->where);
1440
1441   len = x->value.character.length;
1442   lensub = y->value.character.length;
1443
1444   if (len < lensub)
1445     {
1446       mpz_set_si (result->value.integer, 0);
1447       return result;
1448     }
1449
1450   if (back == 0)
1451     {
1452
1453       if (lensub == 0)
1454         {
1455           mpz_set_si (result->value.integer, 1);
1456           return result;
1457         }
1458       else if (lensub == 1)
1459         {
1460           for (i = 0; i < len; i++)
1461             {
1462               for (j = 0; j < lensub; j++)
1463                 {
1464                   if (y->value.character.string[j] ==
1465                       x->value.character.string[i])
1466                     {
1467                       index = i + 1;
1468                       goto done;
1469                     }
1470                 }
1471             }
1472         }
1473       else
1474         {
1475           for (i = 0; i < len; i++)
1476             {
1477               for (j = 0; j < lensub; j++)
1478                 {
1479                   if (y->value.character.string[j] ==
1480                       x->value.character.string[i])
1481                     {
1482                       start = i;
1483                       count = 0;
1484
1485                       for (k = 0; k < lensub; k++)
1486                         {
1487                           if (y->value.character.string[k] ==
1488                               x->value.character.string[k + start])
1489                             count++;
1490                         }
1491
1492                       if (count == lensub)
1493                         {
1494                           index = start + 1;
1495                           goto done;
1496                         }
1497                     }
1498                 }
1499             }
1500         }
1501
1502     }
1503   else
1504     {
1505
1506       if (lensub == 0)
1507         {
1508           mpz_set_si (result->value.integer, len + 1);
1509           return result;
1510         }
1511       else if (lensub == 1)
1512         {
1513           for (i = 0; i < len; i++)
1514             {
1515               for (j = 0; j < lensub; j++)
1516                 {
1517                   if (y->value.character.string[j] ==
1518                       x->value.character.string[len - i])
1519                     {
1520                       index = len - i + 1;
1521                       goto done;
1522                     }
1523                 }
1524             }
1525         }
1526       else
1527         {
1528           for (i = 0; i < len; i++)
1529             {
1530               for (j = 0; j < lensub; j++)
1531                 {
1532                   if (y->value.character.string[j] ==
1533                       x->value.character.string[len - i])
1534                     {
1535                       start = len - i;
1536                       if (start <= len - lensub)
1537                         {
1538                           count = 0;
1539                           for (k = 0; k < lensub; k++)
1540                             if (y->value.character.string[k] ==
1541                                 x->value.character.string[k + start])
1542                               count++;
1543
1544                           if (count == lensub)
1545                             {
1546                               index = start + 1;
1547                               goto done;
1548                             }
1549                         }
1550                       else
1551                         {
1552                           continue;
1553                         }
1554                     }
1555                 }
1556             }
1557         }
1558     }
1559
1560 done:
1561   mpz_set_si (result->value.integer, index);
1562   return range_check (result, "INDEX");
1563 }
1564
1565
1566 gfc_expr *
1567 gfc_simplify_int (gfc_expr * e, gfc_expr * k)
1568 {
1569   gfc_expr *rpart, *rtrunc, *result;
1570   int kind;
1571
1572   kind = get_kind (BT_INTEGER, k, "INT", gfc_default_integer_kind);
1573   if (kind == -1)
1574     return &gfc_bad_expr;
1575
1576   if (e->expr_type != EXPR_CONSTANT)
1577     return NULL;
1578
1579   result = gfc_constant_result (BT_INTEGER, kind, &e->where);
1580
1581   switch (e->ts.type)
1582     {
1583     case BT_INTEGER:
1584       mpz_set (result->value.integer, e->value.integer);
1585       break;
1586
1587     case BT_REAL:
1588       rtrunc = gfc_copy_expr (e);
1589       mpfr_trunc (rtrunc->value.real, e->value.real);
1590       gfc_mpfr_to_mpz (result->value.integer, rtrunc->value.real);
1591       gfc_free_expr (rtrunc);
1592       break;
1593
1594     case BT_COMPLEX:
1595       rpart = gfc_complex2real (e, kind);
1596       rtrunc = gfc_copy_expr (rpart);
1597       mpfr_trunc (rtrunc->value.real, rpart->value.real);
1598       gfc_mpfr_to_mpz (result->value.integer, rtrunc->value.real);
1599       gfc_free_expr (rpart);
1600       gfc_free_expr (rtrunc);
1601       break;
1602
1603     default:
1604       gfc_error ("Argument of INT at %L is not a valid type", &e->where);
1605       gfc_free_expr (result);
1606       return &gfc_bad_expr;
1607     }
1608
1609   return range_check (result, "INT");
1610 }
1611
1612
1613 gfc_expr *
1614 gfc_simplify_ifix (gfc_expr * e)
1615 {
1616   gfc_expr *rtrunc, *result;
1617
1618   if (e->expr_type != EXPR_CONSTANT)
1619     return NULL;
1620
1621   result = gfc_constant_result (BT_INTEGER, gfc_default_integer_kind,
1622                                 &e->where);
1623
1624   rtrunc = gfc_copy_expr (e);
1625
1626   mpfr_trunc (rtrunc->value.real, e->value.real);
1627   gfc_mpfr_to_mpz (result->value.integer, rtrunc->value.real);
1628
1629   gfc_free_expr (rtrunc);
1630   return range_check (result, "IFIX");
1631 }
1632
1633
1634 gfc_expr *
1635 gfc_simplify_idint (gfc_expr * e)
1636 {
1637   gfc_expr *rtrunc, *result;
1638
1639   if (e->expr_type != EXPR_CONSTANT)
1640     return NULL;
1641
1642   result = gfc_constant_result (BT_INTEGER, gfc_default_integer_kind,
1643                                 &e->where);
1644
1645   rtrunc = gfc_copy_expr (e);
1646
1647   mpfr_trunc (rtrunc->value.real, e->value.real);
1648   gfc_mpfr_to_mpz (result->value.integer, rtrunc->value.real);
1649
1650   gfc_free_expr (rtrunc);
1651   return range_check (result, "IDINT");
1652 }
1653
1654
1655 gfc_expr *
1656 gfc_simplify_ior (gfc_expr * x, gfc_expr * y)
1657 {
1658   gfc_expr *result;
1659
1660   if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
1661     return NULL;
1662
1663   result = gfc_constant_result (BT_INTEGER, x->ts.kind, &x->where);
1664
1665   mpz_ior (result->value.integer, x->value.integer, y->value.integer);
1666   return range_check (result, "IOR");
1667 }
1668
1669
1670 gfc_expr *
1671 gfc_simplify_ishft (gfc_expr * e, gfc_expr * s)
1672 {
1673   gfc_expr *result;
1674   int shift, ashift, isize, k, *bits, i;
1675
1676   if (e->expr_type != EXPR_CONSTANT || s->expr_type != EXPR_CONSTANT)
1677     return NULL;
1678
1679   if (gfc_extract_int (s, &shift) != NULL)
1680     {
1681       gfc_error ("Invalid second argument of ISHFT at %L", &s->where);
1682       return &gfc_bad_expr;
1683     }
1684
1685   k = gfc_validate_kind (BT_INTEGER, e->ts.kind, false);
1686
1687   isize = gfc_integer_kinds[k].bit_size;
1688
1689   if (shift >= 0)
1690     ashift = shift;
1691   else
1692     ashift = -shift;
1693
1694   if (ashift > isize)
1695     {
1696       gfc_error
1697         ("Magnitude of second argument of ISHFT exceeds bit size at %L",
1698          &s->where);
1699       return &gfc_bad_expr;
1700     }
1701
1702   result = gfc_constant_result (e->ts.type, e->ts.kind, &e->where);
1703
1704   if (shift == 0)
1705     {
1706       mpz_set (result->value.integer, e->value.integer);
1707       return range_check (result, "ISHFT");
1708     }
1709
1710   bits = gfc_getmem (isize * sizeof (int));
1711
1712   for (i = 0; i < isize; i++)
1713     bits[i] = mpz_tstbit (e->value.integer, i);
1714
1715   if (shift > 0)
1716     {
1717       for (i = 0; i < shift; i++)
1718         mpz_clrbit (result->value.integer, i);
1719
1720       for (i = 0; i < isize - shift; i++)
1721         {
1722           if (bits[i] == 0)
1723             mpz_clrbit (result->value.integer, i + shift);
1724           else
1725             mpz_setbit (result->value.integer, i + shift);
1726         }
1727     }
1728   else
1729     {
1730       for (i = isize - 1; i >= isize - ashift; i--)
1731         mpz_clrbit (result->value.integer, i);
1732
1733       for (i = isize - 1; i >= ashift; i--)
1734         {
1735           if (bits[i] == 0)
1736             mpz_clrbit (result->value.integer, i - ashift);
1737           else
1738             mpz_setbit (result->value.integer, i - ashift);
1739         }
1740     }
1741
1742   twos_complement (result->value.integer, isize);
1743
1744   gfc_free (bits);
1745   return result;
1746 }
1747
1748
1749 gfc_expr *
1750 gfc_simplify_ishftc (gfc_expr * e, gfc_expr * s, gfc_expr * sz)
1751 {
1752   gfc_expr *result;
1753   int shift, ashift, isize, delta, k;
1754   int i, *bits;
1755
1756   if (e->expr_type != EXPR_CONSTANT || s->expr_type != EXPR_CONSTANT)
1757     return NULL;
1758
1759   if (gfc_extract_int (s, &shift) != NULL)
1760     {
1761       gfc_error ("Invalid second argument of ISHFTC at %L", &s->where);
1762       return &gfc_bad_expr;
1763     }
1764
1765   k = gfc_validate_kind (e->ts.type, e->ts.kind, false);
1766
1767   if (sz != NULL)
1768     {
1769       if (gfc_extract_int (sz, &isize) != NULL || isize < 0)
1770         {
1771           gfc_error ("Invalid third argument of ISHFTC at %L", &sz->where);
1772           return &gfc_bad_expr;
1773         }
1774     }
1775   else
1776     isize = gfc_integer_kinds[k].bit_size;
1777
1778   if (shift >= 0)
1779     ashift = shift;
1780   else
1781     ashift = -shift;
1782
1783   if (ashift > isize)
1784     {
1785       gfc_error
1786         ("Magnitude of second argument of ISHFTC exceeds third argument "
1787          "at %L", &s->where);
1788       return &gfc_bad_expr;
1789     }
1790
1791   result = gfc_constant_result (e->ts.type, e->ts.kind, &e->where);
1792
1793   if (shift == 0)
1794     {
1795       mpz_set (result->value.integer, e->value.integer);
1796       return result;
1797     }
1798
1799   bits = gfc_getmem (isize * sizeof (int));
1800
1801   for (i = 0; i < isize; i++)
1802     bits[i] = mpz_tstbit (e->value.integer, i);
1803
1804   delta = isize - ashift;
1805
1806   if (shift > 0)
1807     {
1808       for (i = 0; i < delta; i++)
1809         {
1810           if (bits[i] == 0)
1811             mpz_clrbit (result->value.integer, i + shift);
1812           else
1813             mpz_setbit (result->value.integer, i + shift);
1814         }
1815
1816       for (i = delta; i < isize; i++)
1817         {
1818           if (bits[i] == 0)
1819             mpz_clrbit (result->value.integer, i - delta);
1820           else
1821             mpz_setbit (result->value.integer, i - delta);
1822         }
1823     }
1824   else
1825     {
1826       for (i = 0; i < ashift; i++)
1827         {
1828           if (bits[i] == 0)
1829             mpz_clrbit (result->value.integer, i + delta);
1830           else
1831             mpz_setbit (result->value.integer, i + delta);
1832         }
1833
1834       for (i = ashift; i < isize; i++)
1835         {
1836           if (bits[i] == 0)
1837             mpz_clrbit (result->value.integer, i + shift);
1838           else
1839             mpz_setbit (result->value.integer, i + shift);
1840         }
1841     }
1842
1843   twos_complement (result->value.integer, isize);
1844
1845   gfc_free (bits);
1846   return result;
1847 }
1848
1849
1850 gfc_expr *
1851 gfc_simplify_kind (gfc_expr * e)
1852 {
1853
1854   if (e->ts.type == BT_DERIVED)
1855     {
1856       gfc_error ("Argument of KIND at %L is a DERIVED type", &e->where);
1857       return &gfc_bad_expr;
1858     }
1859
1860   return gfc_int_expr (e->ts.kind);
1861 }
1862
1863
1864 static gfc_expr *
1865 simplify_bound (gfc_expr * array, gfc_expr * dim, int upper)
1866 {
1867   gfc_ref *ref;
1868   gfc_array_spec *as;
1869   gfc_expr *e;
1870   int d;
1871
1872   if (array->expr_type != EXPR_VARIABLE)
1873     return NULL;
1874
1875   if (dim == NULL)
1876     /* TODO: Simplify constant multi-dimensional bounds.  */
1877     return NULL;
1878
1879   if (dim->expr_type != EXPR_CONSTANT)
1880     return NULL;
1881
1882   /* Follow any component references.  */
1883   as = array->symtree->n.sym->as;
1884   for (ref = array->ref; ref; ref = ref->next)
1885     {
1886       switch (ref->type)
1887         {
1888         case REF_ARRAY:
1889           switch (ref->u.ar.type)
1890             {
1891             case AR_ELEMENT:
1892               as = NULL;
1893               continue;
1894
1895             case AR_FULL:
1896               /* We're done because 'as' has already been set in the
1897                  previous iteration.  */
1898               goto done;
1899
1900             case AR_SECTION:
1901             case AR_UNKNOWN:
1902               return NULL;
1903             }
1904
1905           gcc_unreachable ();
1906
1907         case REF_COMPONENT:
1908           as = ref->u.c.component->as;
1909           continue;
1910
1911         case REF_SUBSTRING:
1912           continue;
1913         }
1914     }
1915
1916   gcc_unreachable ();
1917
1918  done:
1919   if (as->type == AS_DEFERRED || as->type == AS_ASSUMED_SHAPE)
1920     return NULL;
1921
1922   d = mpz_get_si (dim->value.integer);
1923
1924   if (d < 1 || d > as->rank
1925       || (d == as->rank && as->type == AS_ASSUMED_SIZE && upper))
1926     {
1927       gfc_error ("DIM argument at %L is out of bounds", &dim->where);
1928       return &gfc_bad_expr;
1929     }
1930
1931   e = upper ? as->upper[d-1] : as->lower[d-1];
1932
1933   if (e->expr_type != EXPR_CONSTANT)
1934     return NULL;
1935
1936   return gfc_copy_expr (e);
1937 }
1938
1939
1940 gfc_expr *
1941 gfc_simplify_lbound (gfc_expr * array, gfc_expr * dim)
1942 {
1943   return simplify_bound (array, dim, 0);
1944 }
1945
1946
1947 gfc_expr *
1948 gfc_simplify_len (gfc_expr * e)
1949 {
1950   gfc_expr *result;
1951
1952   if (e->expr_type == EXPR_CONSTANT)
1953     {
1954       result = gfc_constant_result (BT_INTEGER, gfc_default_integer_kind,
1955                                     &e->where);
1956       mpz_set_si (result->value.integer, e->value.character.length);
1957       return range_check (result, "LEN");
1958     }
1959
1960   if (e->ts.cl != NULL && e->ts.cl->length != NULL
1961       && e->ts.cl->length->expr_type == EXPR_CONSTANT)
1962     {
1963       result = gfc_constant_result (BT_INTEGER, gfc_default_integer_kind,
1964                                     &e->where);
1965       mpz_set (result->value.integer, e->ts.cl->length->value.integer);
1966       return range_check (result, "LEN");
1967     }
1968
1969   return NULL;
1970 }
1971
1972
1973 gfc_expr *
1974 gfc_simplify_len_trim (gfc_expr * e)
1975 {
1976   gfc_expr *result;
1977   int count, len, lentrim, i;
1978
1979   if (e->expr_type != EXPR_CONSTANT)
1980     return NULL;
1981
1982   result = gfc_constant_result (BT_INTEGER, gfc_default_integer_kind,
1983                                 &e->where);
1984
1985   len = e->value.character.length;
1986
1987   for (count = 0, i = 1; i <= len; i++)
1988     if (e->value.character.string[len - i] == ' ')
1989       count++;
1990     else
1991       break;
1992
1993   lentrim = len - count;
1994
1995   mpz_set_si (result->value.integer, lentrim);
1996   return range_check (result, "LEN_TRIM");
1997 }
1998
1999
2000 gfc_expr *
2001 gfc_simplify_lge (gfc_expr * a, gfc_expr * b)
2002 {
2003
2004   if (a->expr_type != EXPR_CONSTANT || b->expr_type != EXPR_CONSTANT)
2005     return NULL;
2006
2007   return gfc_logical_expr (gfc_compare_string (a, b, xascii_table) >= 0,
2008                            &a->where);
2009 }
2010
2011
2012 gfc_expr *
2013 gfc_simplify_lgt (gfc_expr * a, gfc_expr * b)
2014 {
2015
2016   if (a->expr_type != EXPR_CONSTANT || b->expr_type != EXPR_CONSTANT)
2017     return NULL;
2018
2019   return gfc_logical_expr (gfc_compare_string (a, b, xascii_table) > 0,
2020                            &a->where);
2021 }
2022
2023
2024 gfc_expr *
2025 gfc_simplify_lle (gfc_expr * a, gfc_expr * b)
2026 {
2027
2028   if (a->expr_type != EXPR_CONSTANT || b->expr_type != EXPR_CONSTANT)
2029     return NULL;
2030
2031   return gfc_logical_expr (gfc_compare_string (a, b, xascii_table) <= 0,
2032                            &a->where);
2033 }
2034
2035
2036 gfc_expr *
2037 gfc_simplify_llt (gfc_expr * a, gfc_expr * b)
2038 {
2039
2040   if (a->expr_type != EXPR_CONSTANT || b->expr_type != EXPR_CONSTANT)
2041     return NULL;
2042
2043   return gfc_logical_expr (gfc_compare_string (a, b, xascii_table) < 0,
2044                            &a->where);
2045 }
2046
2047
2048 gfc_expr *
2049 gfc_simplify_log (gfc_expr * x)
2050 {
2051   gfc_expr *result;
2052   mpfr_t xr, xi;
2053
2054   if (x->expr_type != EXPR_CONSTANT)
2055     return NULL;
2056
2057   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
2058
2059   gfc_set_model_kind (x->ts.kind);
2060
2061   switch (x->ts.type)
2062     {
2063     case BT_REAL:
2064       if (mpfr_sgn (x->value.real) <= 0)
2065         {
2066           gfc_error
2067             ("Argument of LOG at %L cannot be less than or equal to zero",
2068              &x->where);
2069           gfc_free_expr (result);
2070           return &gfc_bad_expr;
2071         }
2072
2073       mpfr_log(result->value.real, x->value.real, GFC_RND_MODE);
2074       break;
2075
2076     case BT_COMPLEX:
2077       if ((mpfr_sgn (x->value.complex.r) == 0)
2078           && (mpfr_sgn (x->value.complex.i) == 0))
2079         {
2080           gfc_error ("Complex argument of LOG at %L cannot be zero",
2081                      &x->where);
2082           gfc_free_expr (result);
2083           return &gfc_bad_expr;
2084         }
2085
2086       mpfr_init (xr);
2087       mpfr_init (xi);
2088
2089       arctangent2 (x->value.complex.i, x->value.complex.r,
2090                    result->value.complex.i);
2091
2092       mpfr_mul (xr, x->value.complex.r, x->value.complex.r, GFC_RND_MODE);
2093       mpfr_mul (xi, x->value.complex.i, x->value.complex.i, GFC_RND_MODE);
2094       mpfr_add (xr, xr, xi, GFC_RND_MODE);
2095       mpfr_sqrt (xr, xr, GFC_RND_MODE);
2096       mpfr_log (result->value.complex.r, xr, GFC_RND_MODE);
2097
2098       mpfr_clear (xr);
2099       mpfr_clear (xi);
2100
2101       break;
2102
2103     default:
2104       gfc_internal_error ("gfc_simplify_log: bad type");
2105     }
2106
2107   return range_check (result, "LOG");
2108 }
2109
2110
2111 gfc_expr *
2112 gfc_simplify_log10 (gfc_expr * x)
2113 {
2114   gfc_expr *result;
2115
2116   if (x->expr_type != EXPR_CONSTANT)
2117     return NULL;
2118
2119   gfc_set_model_kind (x->ts.kind);
2120
2121   if (mpfr_sgn (x->value.real) <= 0)
2122     {
2123       gfc_error
2124         ("Argument of LOG10 at %L cannot be less than or equal to zero",
2125          &x->where);
2126       return &gfc_bad_expr;
2127     }
2128
2129   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
2130
2131   mpfr_log10 (result->value.real, x->value.real, GFC_RND_MODE);
2132
2133   return range_check (result, "LOG10");
2134 }
2135
2136
2137 gfc_expr *
2138 gfc_simplify_logical (gfc_expr * e, gfc_expr * k)
2139 {
2140   gfc_expr *result;
2141   int kind;
2142
2143   kind = get_kind (BT_LOGICAL, k, "LOGICAL", gfc_default_logical_kind);
2144   if (kind < 0)
2145     return &gfc_bad_expr;
2146
2147   if (e->expr_type != EXPR_CONSTANT)
2148     return NULL;
2149
2150   result = gfc_constant_result (BT_LOGICAL, kind, &e->where);
2151
2152   result->value.logical = e->value.logical;
2153
2154   return result;
2155 }
2156
2157
2158 /* This function is special since MAX() can take any number of
2159    arguments.  The simplified expression is a rewritten version of the
2160    argument list containing at most one constant element.  Other
2161    constant elements are deleted.  Because the argument list has
2162    already been checked, this function always succeeds.  sign is 1 for
2163    MAX(), -1 for MIN().  */
2164
2165 static gfc_expr *
2166 simplify_min_max (gfc_expr * expr, int sign)
2167 {
2168   gfc_actual_arglist *arg, *last, *extremum;
2169   gfc_intrinsic_sym * specific;
2170
2171   last = NULL;
2172   extremum = NULL;
2173   specific = expr->value.function.isym;
2174
2175   arg = expr->value.function.actual;
2176
2177   for (; arg; last = arg, arg = arg->next)
2178     {
2179       if (arg->expr->expr_type != EXPR_CONSTANT)
2180         continue;
2181
2182       if (extremum == NULL)
2183         {
2184           extremum = arg;
2185           continue;
2186         }
2187
2188       switch (arg->expr->ts.type)
2189         {
2190         case BT_INTEGER:
2191           if (mpz_cmp (arg->expr->value.integer,
2192                        extremum->expr->value.integer) * sign > 0)
2193             mpz_set (extremum->expr->value.integer, arg->expr->value.integer);
2194
2195           break;
2196
2197         case BT_REAL:
2198           if (mpfr_cmp (arg->expr->value.real, extremum->expr->value.real) *
2199               sign > 0)
2200             mpfr_set (extremum->expr->value.real, arg->expr->value.real,
2201                       GFC_RND_MODE);
2202
2203           break;
2204
2205         default:
2206           gfc_internal_error ("gfc_simplify_max(): Bad type in arglist");
2207         }
2208
2209       /* Delete the extra constant argument.  */
2210       if (last == NULL)
2211         expr->value.function.actual = arg->next;
2212       else
2213         last->next = arg->next;
2214
2215       arg->next = NULL;
2216       gfc_free_actual_arglist (arg);
2217       arg = last;
2218     }
2219
2220   /* If there is one value left, replace the function call with the
2221      expression.  */
2222   if (expr->value.function.actual->next != NULL)
2223     return NULL;
2224
2225   /* Convert to the correct type and kind.  */
2226   if (expr->ts.type != BT_UNKNOWN)
2227     return gfc_convert_constant (expr->value.function.actual->expr,
2228         expr->ts.type, expr->ts.kind);
2229
2230   if (specific->ts.type != BT_UNKNOWN)
2231     return gfc_convert_constant (expr->value.function.actual->expr,
2232         specific->ts.type, specific->ts.kind);
2233
2234   return gfc_copy_expr (expr->value.function.actual->expr);
2235 }
2236
2237
2238 gfc_expr *
2239 gfc_simplify_min (gfc_expr * e)
2240 {
2241   return simplify_min_max (e, -1);
2242 }
2243
2244
2245 gfc_expr *
2246 gfc_simplify_max (gfc_expr * e)
2247 {
2248   return simplify_min_max (e, 1);
2249 }
2250
2251
2252 gfc_expr *
2253 gfc_simplify_maxexponent (gfc_expr * x)
2254 {
2255   gfc_expr *result;
2256   int i;
2257
2258   i = gfc_validate_kind (BT_REAL, x->ts.kind, false);
2259
2260   result = gfc_int_expr (gfc_real_kinds[i].max_exponent);
2261   result->where = x->where;
2262
2263   return result;
2264 }
2265
2266
2267 gfc_expr *
2268 gfc_simplify_minexponent (gfc_expr * x)
2269 {
2270   gfc_expr *result;
2271   int i;
2272
2273   i = gfc_validate_kind (BT_REAL, x->ts.kind, false);
2274
2275   result = gfc_int_expr (gfc_real_kinds[i].min_exponent);
2276   result->where = x->where;
2277
2278   return result;
2279 }
2280
2281
2282 gfc_expr *
2283 gfc_simplify_mod (gfc_expr * a, gfc_expr * p)
2284 {
2285   gfc_expr *result;
2286   mpfr_t quot, iquot, term;
2287   int kind;
2288
2289   if (a->expr_type != EXPR_CONSTANT || p->expr_type != EXPR_CONSTANT)
2290     return NULL;
2291
2292   kind = a->ts.kind > p->ts.kind ? a->ts.kind : p->ts.kind;
2293   result = gfc_constant_result (a->ts.type, kind, &a->where);
2294
2295   switch (a->ts.type)
2296     {
2297     case BT_INTEGER:
2298       if (mpz_cmp_ui (p->value.integer, 0) == 0)
2299         {
2300           /* Result is processor-dependent.  */
2301           gfc_error ("Second argument MOD at %L is zero", &a->where);
2302           gfc_free_expr (result);
2303           return &gfc_bad_expr;
2304         }
2305       mpz_tdiv_r (result->value.integer, a->value.integer, p->value.integer);
2306       break;
2307
2308     case BT_REAL:
2309       if (mpfr_cmp_ui (p->value.real, 0) == 0)
2310         {
2311           /* Result is processor-dependent.  */
2312           gfc_error ("Second argument of MOD at %L is zero", &p->where);
2313           gfc_free_expr (result);
2314           return &gfc_bad_expr;
2315         }
2316
2317       gfc_set_model_kind (kind);
2318       mpfr_init (quot);
2319       mpfr_init (iquot);
2320       mpfr_init (term);
2321
2322       mpfr_div (quot, a->value.real, p->value.real, GFC_RND_MODE);
2323       mpfr_trunc (iquot, quot);
2324       mpfr_mul (term, iquot, p->value.real, GFC_RND_MODE);
2325       mpfr_sub (result->value.real, a->value.real, term, GFC_RND_MODE);
2326
2327       mpfr_clear (quot);
2328       mpfr_clear (iquot);
2329       mpfr_clear (term);
2330       break;
2331
2332     default:
2333       gfc_internal_error ("gfc_simplify_mod(): Bad arguments");
2334     }
2335
2336   return range_check (result, "MOD");
2337 }
2338
2339
2340 gfc_expr *
2341 gfc_simplify_modulo (gfc_expr * a, gfc_expr * p)
2342 {
2343   gfc_expr *result;
2344   mpfr_t quot, iquot, term;
2345   int kind;
2346
2347   if (a->expr_type != EXPR_CONSTANT || p->expr_type != EXPR_CONSTANT)
2348     return NULL;
2349
2350   kind = a->ts.kind > p->ts.kind ? a->ts.kind : p->ts.kind;
2351   result = gfc_constant_result (a->ts.type, kind, &a->where);
2352
2353   switch (a->ts.type)
2354     {
2355     case BT_INTEGER:
2356       if (mpz_cmp_ui (p->value.integer, 0) == 0)
2357         {
2358           /* Result is processor-dependent. This processor just opts
2359              to not handle it at all.  */
2360           gfc_error ("Second argument of MODULO at %L is zero", &a->where);
2361           gfc_free_expr (result);
2362           return &gfc_bad_expr;
2363         }
2364       mpz_fdiv_r (result->value.integer, a->value.integer, p->value.integer);
2365
2366       break;
2367
2368     case BT_REAL:
2369       if (mpfr_cmp_ui (p->value.real, 0) == 0)
2370         {
2371           /* Result is processor-dependent.  */
2372           gfc_error ("Second argument of MODULO at %L is zero", &p->where);
2373           gfc_free_expr (result);
2374           return &gfc_bad_expr;
2375         }
2376
2377       gfc_set_model_kind (kind);
2378       mpfr_init (quot);
2379       mpfr_init (iquot);
2380       mpfr_init (term);
2381
2382       mpfr_div (quot, a->value.real, p->value.real, GFC_RND_MODE);
2383       mpfr_floor (iquot, quot);
2384       mpfr_mul (term, iquot, p->value.real, GFC_RND_MODE);
2385       mpfr_sub (result->value.real, a->value.real, term, GFC_RND_MODE);
2386
2387       mpfr_clear (quot);
2388       mpfr_clear (iquot);
2389       mpfr_clear (term);
2390       break;
2391
2392     default:
2393       gfc_internal_error ("gfc_simplify_modulo(): Bad arguments");
2394     }
2395
2396   return range_check (result, "MODULO");
2397 }
2398
2399
2400 /* Exists for the sole purpose of consistency with other intrinsics.  */
2401 gfc_expr *
2402 gfc_simplify_mvbits (gfc_expr * f  ATTRIBUTE_UNUSED,
2403                      gfc_expr * fp ATTRIBUTE_UNUSED,
2404                      gfc_expr * l  ATTRIBUTE_UNUSED,
2405                      gfc_expr * to ATTRIBUTE_UNUSED,
2406                      gfc_expr * tp ATTRIBUTE_UNUSED)
2407 {
2408   return NULL;
2409 }
2410
2411
2412 gfc_expr *
2413 gfc_simplify_nearest (gfc_expr * x, gfc_expr * s)
2414 {
2415   gfc_expr *result;
2416   mpfr_t tmp;
2417   int direction, sgn;
2418
2419   if (x->expr_type != EXPR_CONSTANT || s->expr_type != EXPR_CONSTANT)
2420     return NULL;
2421
2422   gfc_set_model_kind (x->ts.kind);
2423   result = gfc_copy_expr (x);
2424
2425   direction = mpfr_sgn (s->value.real);
2426
2427   if (direction == 0)
2428     {
2429       gfc_error ("Second argument of NEAREST at %L may not be zero",
2430                  &s->where);
2431       gfc_free (result);
2432       return &gfc_bad_expr;
2433     }
2434
2435   /* TODO: Use mpfr_nextabove and mpfr_nextbelow once we move to a
2436      newer version of mpfr.  */
2437
2438   sgn = mpfr_sgn (x->value.real);
2439
2440   if (sgn == 0)
2441     {
2442       int k = gfc_validate_kind (BT_REAL, x->ts.kind, 0);
2443
2444       if (direction > 0)
2445         mpfr_add (result->value.real,
2446                   x->value.real, gfc_real_kinds[k].subnormal, GFC_RND_MODE);
2447       else
2448         mpfr_sub (result->value.real,
2449                   x->value.real, gfc_real_kinds[k].subnormal, GFC_RND_MODE);
2450     }
2451   else
2452     {
2453       if (sgn < 0)
2454         {
2455           direction = -direction;
2456           mpfr_neg (result->value.real, result->value.real, GFC_RND_MODE);
2457         }
2458
2459       if (direction > 0)
2460         mpfr_add_one_ulp (result->value.real, GFC_RND_MODE);
2461       else
2462         {
2463           /* In this case the exponent can shrink, which makes us skip
2464              over one number because we subtract one ulp with the
2465              larger exponent.  Thus we need to compensate for this.  */
2466           mpfr_init_set (tmp, result->value.real, GFC_RND_MODE);
2467
2468           mpfr_sub_one_ulp (result->value.real, GFC_RND_MODE);
2469           mpfr_add_one_ulp (result->value.real, GFC_RND_MODE);
2470
2471           /* If we're back to where we started, the spacing is one
2472              ulp, and we get the correct result by subtracting.  */
2473           if (mpfr_cmp (tmp, result->value.real) == 0)
2474             mpfr_sub_one_ulp (result->value.real, GFC_RND_MODE);
2475
2476           mpfr_clear (tmp);
2477         }
2478
2479       if (sgn < 0)
2480         mpfr_neg (result->value.real, result->value.real, GFC_RND_MODE);
2481     }
2482
2483   return range_check (result, "NEAREST");
2484 }
2485
2486
2487 static gfc_expr *
2488 simplify_nint (const char *name, gfc_expr * e, gfc_expr * k)
2489 {
2490   gfc_expr *itrunc, *result;
2491   int kind;
2492
2493   kind = get_kind (BT_INTEGER, k, name, gfc_default_integer_kind);
2494   if (kind == -1)
2495     return &gfc_bad_expr;
2496
2497   if (e->expr_type != EXPR_CONSTANT)
2498     return NULL;
2499
2500   result = gfc_constant_result (BT_INTEGER, kind, &e->where);
2501
2502   itrunc = gfc_copy_expr (e);
2503
2504   mpfr_round(itrunc->value.real, e->value.real);
2505
2506   gfc_mpfr_to_mpz (result->value.integer, itrunc->value.real);
2507
2508   gfc_free_expr (itrunc);
2509
2510   return range_check (result, name);
2511 }
2512
2513
2514 gfc_expr *
2515 gfc_simplify_nint (gfc_expr * e, gfc_expr * k)
2516 {
2517   return simplify_nint ("NINT", e, k);
2518 }
2519
2520
2521 gfc_expr *
2522 gfc_simplify_idnint (gfc_expr * e)
2523 {
2524   return simplify_nint ("IDNINT", e, NULL);
2525 }
2526
2527
2528 gfc_expr *
2529 gfc_simplify_not (gfc_expr * e)
2530 {
2531   gfc_expr *result;
2532   int i;
2533
2534   if (e->expr_type != EXPR_CONSTANT)
2535     return NULL;
2536
2537   result = gfc_constant_result (e->ts.type, e->ts.kind, &e->where);
2538
2539   mpz_com (result->value.integer, e->value.integer);
2540
2541   /* Because of how GMP handles numbers, the result must be ANDed with
2542      the max_int mask.  For radices <> 2, this will require change.  */
2543
2544   i = gfc_validate_kind (BT_INTEGER, e->ts.kind, false);
2545
2546   mpz_and (result->value.integer, result->value.integer,
2547            gfc_integer_kinds[i].max_int);
2548
2549   twos_complement (result->value.integer, gfc_integer_kinds[i].bit_size);
2550
2551   return range_check (result, "NOT");
2552 }
2553
2554
2555 gfc_expr *
2556 gfc_simplify_null (gfc_expr * mold)
2557 {
2558   gfc_expr *result;
2559
2560   if (mold == NULL)
2561     {
2562       result = gfc_get_expr ();
2563       result->ts.type = BT_UNKNOWN;
2564     }
2565   else
2566     result = gfc_copy_expr (mold);
2567   result->expr_type = EXPR_NULL;
2568
2569   return result;
2570 }
2571
2572
2573 gfc_expr *
2574 gfc_simplify_or (gfc_expr * x, gfc_expr * y)
2575 {
2576   gfc_expr *result;
2577   int kind;
2578
2579   if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
2580     return NULL;
2581
2582   kind = x->ts.kind > y->ts.kind ? x->ts.kind : y->ts.kind;
2583   if (x->ts.type == BT_INTEGER)
2584     {
2585       result = gfc_constant_result (BT_INTEGER, kind, &x->where);
2586       mpz_ior (result->value.integer, x->value.integer, y->value.integer);
2587     }
2588   else /* BT_LOGICAL */
2589     {
2590       result = gfc_constant_result (BT_LOGICAL, kind, &x->where);
2591       result->value.logical = x->value.logical || y->value.logical;
2592     }
2593
2594   return range_check (result, "OR");
2595 }
2596
2597
2598 gfc_expr *
2599 gfc_simplify_precision (gfc_expr * e)
2600 {
2601   gfc_expr *result;
2602   int i;
2603
2604   i = gfc_validate_kind (e->ts.type, e->ts.kind, false);
2605
2606   result = gfc_int_expr (gfc_real_kinds[i].precision);
2607   result->where = e->where;
2608
2609   return result;
2610 }
2611
2612
2613 gfc_expr *
2614 gfc_simplify_radix (gfc_expr * e)
2615 {
2616   gfc_expr *result;
2617   int i;
2618
2619   i = gfc_validate_kind (e->ts.type, e->ts.kind, false);
2620   switch (e->ts.type)
2621     {
2622     case BT_INTEGER:
2623       i = gfc_integer_kinds[i].radix;
2624       break;
2625
2626     case BT_REAL:
2627       i = gfc_real_kinds[i].radix;
2628       break;
2629
2630     default:
2631       gcc_unreachable ();
2632     }
2633
2634   result = gfc_int_expr (i);
2635   result->where = e->where;
2636
2637   return result;
2638 }
2639
2640
2641 gfc_expr *
2642 gfc_simplify_range (gfc_expr * e)
2643 {
2644   gfc_expr *result;
2645   int i;
2646   long j;
2647
2648   i = gfc_validate_kind (e->ts.type, e->ts.kind, false);
2649
2650   switch (e->ts.type)
2651     {
2652     case BT_INTEGER:
2653       j = gfc_integer_kinds[i].range;
2654       break;
2655
2656     case BT_REAL:
2657     case BT_COMPLEX:
2658       j = gfc_real_kinds[i].range;
2659       break;
2660
2661     default:
2662       gcc_unreachable ();
2663     }
2664
2665   result = gfc_int_expr (j);
2666   result->where = e->where;
2667
2668   return result;
2669 }
2670
2671
2672 gfc_expr *
2673 gfc_simplify_real (gfc_expr * e, gfc_expr * k)
2674 {
2675   gfc_expr *result;
2676   int kind;
2677
2678   if (e->ts.type == BT_COMPLEX)
2679     kind = get_kind (BT_REAL, k, "REAL", e->ts.kind);
2680   else
2681     kind = get_kind (BT_REAL, k, "REAL", gfc_default_real_kind);
2682
2683   if (kind == -1)
2684     return &gfc_bad_expr;
2685
2686   if (e->expr_type != EXPR_CONSTANT)
2687     return NULL;
2688
2689   switch (e->ts.type)
2690     {
2691     case BT_INTEGER:
2692       result = gfc_int2real (e, kind);
2693       break;
2694
2695     case BT_REAL:
2696       result = gfc_real2real (e, kind);
2697       break;
2698
2699     case BT_COMPLEX:
2700       result = gfc_complex2real (e, kind);
2701       break;
2702
2703     default:
2704       gfc_internal_error ("bad type in REAL");
2705       /* Not reached */
2706     }
2707
2708   return range_check (result, "REAL");
2709 }
2710
2711
2712 gfc_expr *
2713 gfc_simplify_realpart (gfc_expr * e)
2714 {
2715   gfc_expr *result;
2716
2717   if (e->expr_type != EXPR_CONSTANT)
2718     return NULL;
2719
2720   result = gfc_constant_result (BT_REAL, e->ts.kind, &e->where);
2721   mpfr_set (result->value.real, e->value.complex.r, GFC_RND_MODE);
2722
2723   return range_check (result, "REALPART");
2724 }
2725
2726 gfc_expr *
2727 gfc_simplify_repeat (gfc_expr * e, gfc_expr * n)
2728 {
2729   gfc_expr *result;
2730   int i, j, len, ncopies, nlen;
2731
2732   if (e->expr_type != EXPR_CONSTANT || n->expr_type != EXPR_CONSTANT)
2733     return NULL;
2734
2735   if (n != NULL && (gfc_extract_int (n, &ncopies) != NULL || ncopies < 0))
2736     {
2737       gfc_error ("Invalid second argument of REPEAT at %L", &n->where);
2738       return &gfc_bad_expr;
2739     }
2740
2741   len = e->value.character.length;
2742   nlen = ncopies * len;
2743
2744   result = gfc_constant_result (BT_CHARACTER, e->ts.kind, &e->where);
2745
2746   if (ncopies == 0)
2747     {
2748       result->value.character.string = gfc_getmem (1);
2749       result->value.character.length = 0;
2750       result->value.character.string[0] = '\0';
2751       return result;
2752     }
2753
2754   result->value.character.length = nlen;
2755   result->value.character.string = gfc_getmem (nlen + 1);
2756
2757   for (i = 0; i < ncopies; i++)
2758     for (j = 0; j < len; j++)
2759       result->value.character.string[j + i * len] =
2760         e->value.character.string[j];
2761
2762   result->value.character.string[nlen] = '\0';  /* For debugger */
2763   return result;
2764 }
2765
2766
2767 /* This one is a bear, but mainly has to do with shuffling elements.  */
2768
2769 gfc_expr *
2770 gfc_simplify_reshape (gfc_expr * source, gfc_expr * shape_exp,
2771                       gfc_expr * pad, gfc_expr * order_exp)
2772 {
2773
2774   int order[GFC_MAX_DIMENSIONS], shape[GFC_MAX_DIMENSIONS];
2775   int i, rank, npad, x[GFC_MAX_DIMENSIONS];
2776   gfc_constructor *head, *tail;
2777   mpz_t index, size;
2778   unsigned long j;
2779   size_t nsource;
2780   gfc_expr *e;
2781
2782   /* Unpack the shape array.  */
2783   if (source->expr_type != EXPR_ARRAY || !gfc_is_constant_expr (source))
2784     return NULL;
2785
2786   if (shape_exp->expr_type != EXPR_ARRAY || !gfc_is_constant_expr (shape_exp))
2787     return NULL;
2788
2789   if (pad != NULL
2790       && (pad->expr_type != EXPR_ARRAY
2791           || !gfc_is_constant_expr (pad)))
2792     return NULL;
2793
2794   if (order_exp != NULL
2795       && (order_exp->expr_type != EXPR_ARRAY
2796           || !gfc_is_constant_expr (order_exp)))
2797     return NULL;
2798
2799   mpz_init (index);
2800   rank = 0;
2801   head = tail = NULL;
2802
2803   for (;;)
2804     {
2805       e = gfc_get_array_element (shape_exp, rank);
2806       if (e == NULL)
2807         break;
2808
2809       if (gfc_extract_int (e, &shape[rank]) != NULL)
2810         {
2811           gfc_error ("Integer too large in shape specification at %L",
2812                      &e->where);
2813           gfc_free_expr (e);
2814           goto bad_reshape;
2815         }
2816
2817       gfc_free_expr (e);
2818
2819       if (rank >= GFC_MAX_DIMENSIONS)
2820         {
2821           gfc_error ("Too many dimensions in shape specification for RESHAPE "
2822                      "at %L", &e->where);
2823
2824           goto bad_reshape;
2825         }
2826
2827       if (shape[rank] < 0)
2828         {
2829           gfc_error ("Shape specification at %L cannot be negative",
2830                      &e->where);
2831           goto bad_reshape;
2832         }
2833
2834       rank++;
2835     }
2836
2837   if (rank == 0)
2838     {
2839       gfc_error ("Shape specification at %L cannot be the null array",
2840                  &shape_exp->where);
2841       goto bad_reshape;
2842     }
2843
2844   /* Now unpack the order array if present.  */
2845   if (order_exp == NULL)
2846     {
2847       for (i = 0; i < rank; i++)
2848         order[i] = i;
2849
2850     }
2851   else
2852     {
2853
2854       for (i = 0; i < rank; i++)
2855         x[i] = 0;
2856
2857       for (i = 0; i < rank; i++)
2858         {
2859           e = gfc_get_array_element (order_exp, i);
2860           if (e == NULL)
2861             {
2862               gfc_error
2863                 ("ORDER parameter of RESHAPE at %L is not the same size "
2864                  "as SHAPE parameter", &order_exp->where);
2865               goto bad_reshape;
2866             }
2867
2868           if (gfc_extract_int (e, &order[i]) != NULL)
2869             {
2870               gfc_error ("Error in ORDER parameter of RESHAPE at %L",
2871                          &e->where);
2872               gfc_free_expr (e);
2873               goto bad_reshape;
2874             }
2875
2876           gfc_free_expr (e);
2877
2878           if (order[i] < 1 || order[i] > rank)
2879             {
2880               gfc_error ("ORDER parameter of RESHAPE at %L is out of range",
2881                          &e->where);
2882               goto bad_reshape;
2883             }
2884
2885           order[i]--;
2886
2887           if (x[order[i]])
2888             {
2889               gfc_error ("Invalid permutation in ORDER parameter at %L",
2890                          &e->where);
2891               goto bad_reshape;
2892             }
2893
2894           x[order[i]] = 1;
2895         }
2896     }
2897
2898   /* Count the elements in the source and padding arrays.  */
2899
2900   npad = 0;
2901   if (pad != NULL)
2902     {
2903       gfc_array_size (pad, &size);
2904       npad = mpz_get_ui (size);
2905       mpz_clear (size);
2906     }
2907
2908   gfc_array_size (source, &size);
2909   nsource = mpz_get_ui (size);
2910   mpz_clear (size);
2911
2912   /* If it weren't for that pesky permutation we could just loop
2913      through the source and round out any shortage with pad elements.
2914      But no, someone just had to have the compiler do something the
2915      user should be doing.  */
2916
2917   for (i = 0; i < rank; i++)
2918     x[i] = 0;
2919
2920   for (;;)
2921     {
2922       /* Figure out which element to extract.  */
2923       mpz_set_ui (index, 0);
2924
2925       for (i = rank - 1; i >= 0; i--)
2926         {
2927           mpz_add_ui (index, index, x[order[i]]);
2928           if (i != 0)
2929             mpz_mul_ui (index, index, shape[order[i - 1]]);
2930         }
2931
2932       if (mpz_cmp_ui (index, INT_MAX) > 0)
2933         gfc_internal_error ("Reshaped array too large at %L", &e->where);
2934
2935       j = mpz_get_ui (index);
2936
2937       if (j < nsource)
2938         e = gfc_get_array_element (source, j);
2939       else
2940         {
2941           j = j - nsource;
2942
2943           if (npad == 0)
2944             {
2945               gfc_error
2946                 ("PAD parameter required for short SOURCE parameter at %L",
2947                  &source->where);
2948               goto bad_reshape;
2949             }
2950
2951           j = j % npad;
2952           e = gfc_get_array_element (pad, j);
2953         }
2954
2955       if (head == NULL)
2956         head = tail = gfc_get_constructor ();
2957       else
2958         {
2959           tail->next = gfc_get_constructor ();
2960           tail = tail->next;
2961         }
2962
2963       if (e == NULL)
2964         goto bad_reshape;
2965
2966       tail->where = e->where;
2967       tail->expr = e;
2968
2969       /* Calculate the next element.  */
2970       i = 0;
2971
2972 inc:
2973       if (++x[i] < shape[i])
2974         continue;
2975       x[i++] = 0;
2976       if (i < rank)
2977         goto inc;
2978
2979       break;
2980     }
2981
2982   mpz_clear (index);
2983
2984   e = gfc_get_expr ();
2985   e->where = source->where;
2986   e->expr_type = EXPR_ARRAY;
2987   e->value.constructor = head;
2988   e->shape = gfc_get_shape (rank);
2989
2990   for (i = 0; i < rank; i++)
2991     mpz_init_set_ui (e->shape[i], shape[i]);
2992
2993   e->ts = source->ts;
2994   e->rank = rank;
2995
2996   return e;
2997
2998 bad_reshape:
2999   gfc_free_constructor (head);
3000   mpz_clear (index);
3001   return &gfc_bad_expr;
3002 }
3003
3004
3005 gfc_expr *
3006 gfc_simplify_rrspacing (gfc_expr * x)
3007 {
3008   gfc_expr *result;
3009   mpfr_t absv, log2, exp, frac, pow2;
3010   int i, p;
3011
3012   if (x->expr_type != EXPR_CONSTANT)
3013     return NULL;
3014
3015   i = gfc_validate_kind (x->ts.type, x->ts.kind, false);
3016
3017   result = gfc_constant_result (BT_REAL, x->ts.kind, &x->where);
3018
3019   p = gfc_real_kinds[i].digits;
3020
3021   gfc_set_model_kind (x->ts.kind);
3022
3023   if (mpfr_sgn (x->value.real) == 0)
3024     {
3025       mpfr_ui_div (result->value.real, 1, gfc_real_kinds[i].tiny, GFC_RND_MODE);
3026       return result;
3027     }
3028
3029   mpfr_init (log2);
3030   mpfr_init (absv);
3031   mpfr_init (frac);
3032   mpfr_init (pow2);
3033   mpfr_init (exp);
3034
3035   mpfr_abs (absv, x->value.real, GFC_RND_MODE);
3036   mpfr_log2 (log2, absv, GFC_RND_MODE);
3037
3038   mpfr_trunc (log2, log2);
3039   mpfr_add_ui (exp, log2, 1, GFC_RND_MODE);
3040
3041   mpfr_ui_pow (pow2, 2, exp, GFC_RND_MODE);
3042   mpfr_div (frac, absv, pow2, GFC_RND_MODE);
3043
3044   mpfr_mul_2exp (result->value.real, frac, (unsigned long)p, GFC_RND_MODE);
3045
3046   mpfr_clear (log2);
3047   mpfr_clear (absv);
3048   mpfr_clear (frac);
3049   mpfr_clear (pow2);
3050   mpfr_clear (exp);
3051
3052   return range_check (result, "RRSPACING");
3053 }
3054
3055
3056 gfc_expr *
3057 gfc_simplify_scale (gfc_expr * x, gfc_expr * i)
3058 {
3059   int k, neg_flag, power, exp_range;
3060   mpfr_t scale, radix;
3061   gfc_expr *result;
3062
3063   if (x->expr_type != EXPR_CONSTANT || i->expr_type != EXPR_CONSTANT)
3064     return NULL;
3065
3066   result = gfc_constant_result (BT_REAL, x->ts.kind, &x->where);
3067
3068   if (mpfr_sgn (x->value.real) == 0)
3069     {
3070       mpfr_set_ui (result->value.real, 0, GFC_RND_MODE);
3071       return result;
3072     }
3073
3074   k = gfc_validate_kind (BT_REAL, x->ts.kind, false);
3075
3076   exp_range = gfc_real_kinds[k].max_exponent - gfc_real_kinds[k].min_exponent;
3077
3078   /* This check filters out values of i that would overflow an int.  */
3079   if (mpz_cmp_si (i->value.integer, exp_range + 2) > 0
3080       || mpz_cmp_si (i->value.integer, -exp_range - 2) < 0)
3081     {
3082       gfc_error ("Result of SCALE overflows its kind at %L", &result->where);
3083       return &gfc_bad_expr;
3084     }
3085
3086   /* Compute scale = radix ** power.  */
3087   power = mpz_get_si (i->value.integer);
3088
3089   if (power >= 0)
3090     neg_flag = 0;
3091   else
3092     {
3093       neg_flag = 1;
3094       power = -power;
3095     }
3096
3097   gfc_set_model_kind (x->ts.kind);
3098   mpfr_init (scale);
3099   mpfr_init (radix);
3100   mpfr_set_ui (radix, gfc_real_kinds[k].radix, GFC_RND_MODE);
3101   mpfr_pow_ui (scale, radix, power, GFC_RND_MODE);
3102
3103   if (neg_flag)
3104     mpfr_div (result->value.real, x->value.real, scale, GFC_RND_MODE);
3105   else
3106     mpfr_mul (result->value.real, x->value.real, scale, GFC_RND_MODE);
3107
3108   mpfr_clear (scale);
3109   mpfr_clear (radix);
3110
3111   return range_check (result, "SCALE");
3112 }
3113
3114
3115 gfc_expr *
3116 gfc_simplify_scan (gfc_expr * e, gfc_expr * c, gfc_expr * b)
3117 {
3118   gfc_expr *result;
3119   int back;
3120   size_t i;
3121   size_t indx, len, lenc;
3122
3123   if (e->expr_type != EXPR_CONSTANT || c->expr_type != EXPR_CONSTANT)
3124     return NULL;
3125
3126   if (b != NULL && b->value.logical != 0)
3127     back = 1;
3128   else
3129     back = 0;
3130
3131   result = gfc_constant_result (BT_INTEGER, gfc_default_integer_kind,
3132                                 &e->where);
3133
3134   len = e->value.character.length;
3135   lenc = c->value.character.length;
3136
3137   if (len == 0 || lenc == 0)
3138     {
3139       indx = 0;
3140     }
3141   else
3142     {
3143       if (back == 0)
3144         {
3145           indx =
3146             strcspn (e->value.character.string, c->value.character.string) + 1;
3147           if (indx > len)
3148             indx = 0;
3149         }
3150       else
3151         {
3152           i = 0;
3153           for (indx = len; indx > 0; indx--)
3154             {
3155               for (i = 0; i < lenc; i++)
3156                 {
3157                   if (c->value.character.string[i]
3158                         == e->value.character.string[indx - 1])
3159                     break;
3160                 }
3161               if (i < lenc)
3162                 break;
3163             }
3164         }
3165     }
3166   mpz_set_ui (result->value.integer, indx);
3167   return range_check (result, "SCAN");
3168 }
3169
3170
3171 gfc_expr *
3172 gfc_simplify_selected_int_kind (gfc_expr * e)
3173 {
3174   int i, kind, range;
3175   gfc_expr *result;
3176
3177   if (e->expr_type != EXPR_CONSTANT || gfc_extract_int (e, &range) != NULL)
3178     return NULL;
3179
3180   kind = INT_MAX;
3181
3182   for (i = 0; gfc_integer_kinds[i].kind != 0; i++)
3183     if (gfc_integer_kinds[i].range >= range
3184         && gfc_integer_kinds[i].kind < kind)
3185       kind = gfc_integer_kinds[i].kind;
3186
3187   if (kind == INT_MAX)
3188     kind = -1;
3189
3190   result = gfc_int_expr (kind);
3191   result->where = e->where;
3192
3193   return result;
3194 }
3195
3196
3197 gfc_expr *
3198 gfc_simplify_selected_real_kind (gfc_expr * p, gfc_expr * q)
3199 {
3200   int range, precision, i, kind, found_precision, found_range;
3201   gfc_expr *result;
3202
3203   if (p == NULL)
3204     precision = 0;
3205   else
3206     {
3207       if (p->expr_type != EXPR_CONSTANT
3208           || gfc_extract_int (p, &precision) != NULL)
3209         return NULL;
3210     }
3211
3212   if (q == NULL)
3213     range = 0;
3214   else
3215     {
3216       if (q->expr_type != EXPR_CONSTANT
3217           || gfc_extract_int (q, &range) != NULL)
3218         return NULL;
3219     }
3220
3221   kind = INT_MAX;
3222   found_precision = 0;
3223   found_range = 0;
3224
3225   for (i = 0; gfc_real_kinds[i].kind != 0; i++)
3226     {
3227       if (gfc_real_kinds[i].precision >= precision)
3228         found_precision = 1;
3229
3230       if (gfc_real_kinds[i].range >= range)
3231         found_range = 1;
3232
3233       if (gfc_real_kinds[i].precision >= precision
3234           && gfc_real_kinds[i].range >= range && gfc_real_kinds[i].kind < kind)
3235         kind = gfc_real_kinds[i].kind;
3236     }
3237
3238   if (kind == INT_MAX)
3239     {
3240       kind = 0;
3241
3242       if (!found_precision)
3243         kind = -1;
3244       if (!found_range)
3245         kind -= 2;
3246     }
3247
3248   result = gfc_int_expr (kind);
3249   result->where = (p != NULL) ? p->where : q->where;
3250
3251   return result;
3252 }
3253
3254
3255 gfc_expr *
3256 gfc_simplify_set_exponent (gfc_expr * x, gfc_expr * i)
3257 {
3258   gfc_expr *result;
3259   mpfr_t exp, absv, log2, pow2, frac;
3260   unsigned long exp2;
3261
3262   if (x->expr_type != EXPR_CONSTANT || i->expr_type != EXPR_CONSTANT)
3263     return NULL;
3264
3265   result = gfc_constant_result (BT_REAL, x->ts.kind, &x->where);
3266
3267   gfc_set_model_kind (x->ts.kind);
3268
3269   if (mpfr_sgn (x->value.real) == 0)
3270     {
3271       mpfr_set_ui (result->value.real, 0, GFC_RND_MODE);
3272       return result;
3273     }
3274
3275   mpfr_init (absv);
3276   mpfr_init (log2);
3277   mpfr_init (exp);
3278   mpfr_init (pow2);
3279   mpfr_init (frac);
3280
3281   mpfr_abs (absv, x->value.real, GFC_RND_MODE);
3282   mpfr_log2 (log2, absv, GFC_RND_MODE);
3283
3284   mpfr_trunc (log2, log2);
3285   mpfr_add_ui (exp, log2, 1, GFC_RND_MODE);
3286
3287   /* Old exponent value, and fraction.  */
3288   mpfr_ui_pow (pow2, 2, exp, GFC_RND_MODE);
3289
3290   mpfr_div (frac, absv, pow2, GFC_RND_MODE);
3291
3292   /* New exponent.  */
3293   exp2 = (unsigned long) mpz_get_d (i->value.integer);
3294   mpfr_mul_2exp (result->value.real, frac, exp2, GFC_RND_MODE);
3295
3296   mpfr_clear (absv);
3297   mpfr_clear (log2);
3298   mpfr_clear (pow2);
3299   mpfr_clear (frac);
3300
3301   return range_check (result, "SET_EXPONENT");
3302 }
3303
3304
3305 gfc_expr *
3306 gfc_simplify_shape (gfc_expr * source)
3307 {
3308   mpz_t shape[GFC_MAX_DIMENSIONS];
3309   gfc_expr *result, *e, *f;
3310   gfc_array_ref *ar;
3311   int n;
3312   try t;
3313
3314   if (source->rank == 0 || source->expr_type != EXPR_VARIABLE)
3315     return NULL;
3316
3317   result = gfc_start_constructor (BT_INTEGER, gfc_default_integer_kind,
3318                                   &source->where);
3319
3320   ar = gfc_find_array_ref (source);
3321
3322   t = gfc_array_ref_shape (ar, shape);
3323
3324   for (n = 0; n < source->rank; n++)
3325     {
3326       e = gfc_constant_result (BT_INTEGER, gfc_default_integer_kind,
3327                                &source->where);
3328
3329       if (t == SUCCESS)
3330         {
3331           mpz_set (e->value.integer, shape[n]);
3332           mpz_clear (shape[n]);
3333         }
3334       else
3335         {
3336           mpz_set_ui (e->value.integer, n + 1);
3337
3338           f = gfc_simplify_size (source, e);
3339           gfc_free_expr (e);
3340           if (f == NULL)
3341             {
3342               gfc_free_expr (result);
3343               return NULL;
3344             }
3345           else
3346             {
3347               e = f;
3348             }
3349         }
3350
3351       gfc_append_constructor (result, e);
3352     }
3353
3354   return result;
3355 }
3356
3357
3358 gfc_expr *
3359 gfc_simplify_size (gfc_expr * array, gfc_expr * dim)
3360 {
3361   mpz_t size;
3362   gfc_expr *result;
3363   int d;
3364
3365   if (dim == NULL)
3366     {
3367       if (gfc_array_size (array, &size) == FAILURE)
3368         return NULL;
3369     }
3370   else
3371     {
3372       if (dim->expr_type != EXPR_CONSTANT)
3373         return NULL;
3374
3375       d = mpz_get_ui (dim->value.integer) - 1;
3376       if (gfc_array_dimen_size (array, d, &size) == FAILURE)
3377         return NULL;
3378     }
3379
3380   result = gfc_constant_result (BT_INTEGER, gfc_default_integer_kind,
3381                                 &array->where);
3382
3383   mpz_set (result->value.integer, size);
3384
3385   return result;
3386 }
3387
3388
3389 gfc_expr *
3390 gfc_simplify_sign (gfc_expr * x, gfc_expr * y)
3391 {
3392   gfc_expr *result;
3393
3394   if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
3395     return NULL;
3396
3397   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
3398
3399   switch (x->ts.type)
3400     {
3401     case BT_INTEGER:
3402       mpz_abs (result->value.integer, x->value.integer);
3403       if (mpz_sgn (y->value.integer) < 0)
3404         mpz_neg (result->value.integer, result->value.integer);
3405
3406       break;
3407
3408     case BT_REAL:
3409       /* TODO: Handle -0.0 and +0.0 correctly on machines that support
3410          it.  */
3411       mpfr_abs (result->value.real, x->value.real, GFC_RND_MODE);
3412       if (mpfr_sgn (y->value.real) < 0)
3413         mpfr_neg (result->value.real, result->value.real, GFC_RND_MODE);
3414
3415       break;
3416
3417     default:
3418       gfc_internal_error ("Bad type in gfc_simplify_sign");
3419     }
3420
3421   return result;
3422 }
3423
3424
3425 gfc_expr *
3426 gfc_simplify_sin (gfc_expr * x)
3427 {
3428   gfc_expr *result;
3429   mpfr_t xp, xq;
3430
3431   if (x->expr_type != EXPR_CONSTANT)
3432     return NULL;
3433
3434   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
3435
3436   switch (x->ts.type)
3437     {
3438     case BT_REAL:
3439       mpfr_sin (result->value.real, x->value.real, GFC_RND_MODE);
3440       break;
3441
3442     case BT_COMPLEX:
3443       gfc_set_model (x->value.real);
3444       mpfr_init (xp);
3445       mpfr_init (xq);
3446
3447       mpfr_sin  (xp, x->value.complex.r, GFC_RND_MODE);
3448       mpfr_cosh (xq, x->value.complex.i, GFC_RND_MODE);
3449       mpfr_mul (result->value.complex.r, xp, xq, GFC_RND_MODE);
3450
3451       mpfr_cos  (xp, x->value.complex.r, GFC_RND_MODE);
3452       mpfr_sinh (xq, x->value.complex.i, GFC_RND_MODE);
3453       mpfr_mul (result->value.complex.i, xp, xq, GFC_RND_MODE);
3454
3455       mpfr_clear (xp);
3456       mpfr_clear (xq);
3457       break;
3458
3459     default:
3460       gfc_internal_error ("in gfc_simplify_sin(): Bad type");
3461     }
3462
3463   return range_check (result, "SIN");
3464 }
3465
3466
3467 gfc_expr *
3468 gfc_simplify_sinh (gfc_expr * x)
3469 {
3470   gfc_expr *result;
3471
3472   if (x->expr_type != EXPR_CONSTANT)
3473     return NULL;
3474
3475   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
3476
3477   mpfr_sinh(result->value.real, x->value.real, GFC_RND_MODE);
3478
3479   return range_check (result, "SINH");
3480 }
3481
3482
3483 /* The argument is always a double precision real that is converted to
3484    single precision.  TODO: Rounding!  */
3485
3486 gfc_expr *
3487 gfc_simplify_sngl (gfc_expr * a)
3488 {
3489   gfc_expr *result;
3490
3491   if (a->expr_type != EXPR_CONSTANT)
3492     return NULL;
3493
3494   result = gfc_real2real (a, gfc_default_real_kind);
3495   return range_check (result, "SNGL");
3496 }
3497
3498
3499 gfc_expr *
3500 gfc_simplify_spacing (gfc_expr * x)
3501 {
3502   gfc_expr *result;
3503   mpfr_t absv, log2;
3504   long diff;
3505   int i, p;
3506
3507   if (x->expr_type != EXPR_CONSTANT)
3508     return NULL;
3509
3510   i = gfc_validate_kind (x->ts.type, x->ts.kind, false);
3511
3512   p = gfc_real_kinds[i].digits;
3513
3514   result = gfc_constant_result (BT_REAL, x->ts.kind, &x->where);
3515
3516   gfc_set_model_kind (x->ts.kind);
3517
3518   if (mpfr_sgn (x->value.real) == 0)
3519     {
3520       mpfr_set (result->value.real, gfc_real_kinds[i].tiny, GFC_RND_MODE);
3521       return result;
3522     }
3523
3524   mpfr_init (log2);
3525   mpfr_init (absv);
3526
3527   mpfr_abs (absv, x->value.real, GFC_RND_MODE);
3528   mpfr_log2 (log2, absv, GFC_RND_MODE);
3529   mpfr_trunc (log2, log2);
3530
3531   mpfr_add_ui (log2, log2, 1, GFC_RND_MODE);
3532
3533   /* FIXME: We should be using mpfr_get_si here, but this function is
3534      not available with the version of mpfr distributed with gmp (as of
3535      2004-09-17). Replace once mpfr has been imported into the gcc cvs
3536      tree.  */
3537   diff = (long)mpfr_get_d (log2, GFC_RND_MODE) - (long)p;
3538   mpfr_set_ui (result->value.real, 1, GFC_RND_MODE);
3539   mpfr_mul_2si (result->value.real, result->value.real, diff, GFC_RND_MODE);
3540
3541   mpfr_clear (log2);
3542   mpfr_clear (absv);
3543
3544   if (mpfr_cmp (result->value.real, gfc_real_kinds[i].tiny) < 0)
3545     mpfr_set (result->value.real, gfc_real_kinds[i].tiny, GFC_RND_MODE);
3546
3547   return range_check (result, "SPACING");
3548 }
3549
3550
3551 gfc_expr *
3552 gfc_simplify_sqrt (gfc_expr * e)
3553 {
3554   gfc_expr *result;
3555   mpfr_t ac, ad, s, t, w;
3556
3557   if (e->expr_type != EXPR_CONSTANT)
3558     return NULL;
3559
3560   result = gfc_constant_result (e->ts.type, e->ts.kind, &e->where);
3561
3562   switch (e->ts.type)
3563     {
3564     case BT_REAL:
3565       if (mpfr_cmp_si (e->value.real, 0) < 0)
3566         goto negative_arg;
3567       mpfr_sqrt (result->value.real, e->value.real, GFC_RND_MODE);
3568
3569       break;
3570
3571     case BT_COMPLEX:
3572       /* Formula taken from Numerical Recipes to avoid over- and
3573          underflow.  */
3574
3575       gfc_set_model (e->value.real);
3576       mpfr_init (ac);
3577       mpfr_init (ad);
3578       mpfr_init (s);
3579       mpfr_init (t);
3580       mpfr_init (w);
3581
3582       if (mpfr_cmp_ui (e->value.complex.r, 0) == 0
3583           && mpfr_cmp_ui (e->value.complex.i, 0) == 0)
3584         {
3585
3586           mpfr_set_ui (result->value.complex.r, 0, GFC_RND_MODE);
3587           mpfr_set_ui (result->value.complex.i, 0, GFC_RND_MODE);
3588           break;
3589         }
3590
3591       mpfr_abs (ac, e->value.complex.r, GFC_RND_MODE);
3592       mpfr_abs (ad, e->value.complex.i, GFC_RND_MODE);
3593
3594       if (mpfr_cmp (ac, ad) >= 0)
3595         {
3596           mpfr_div (t, e->value.complex.i, e->value.complex.r, GFC_RND_MODE);
3597           mpfr_mul (t, t, t, GFC_RND_MODE);
3598           mpfr_add_ui (t, t, 1, GFC_RND_MODE);
3599           mpfr_sqrt (t, t, GFC_RND_MODE);
3600           mpfr_add_ui (t, t, 1, GFC_RND_MODE);
3601           mpfr_div_ui (t, t, 2, GFC_RND_MODE);
3602           mpfr_sqrt (t, t, GFC_RND_MODE);
3603           mpfr_sqrt (s, ac, GFC_RND_MODE);
3604           mpfr_mul (w, s, t, GFC_RND_MODE);
3605         }
3606       else
3607         {
3608           mpfr_div (s, e->value.complex.r, e->value.complex.i, GFC_RND_MODE);
3609           mpfr_mul (t, s, s, GFC_RND_MODE);
3610           mpfr_add_ui (t, t, 1, GFC_RND_MODE);
3611           mpfr_sqrt (t, t, GFC_RND_MODE);
3612           mpfr_abs (s, s, GFC_RND_MODE);
3613           mpfr_add (t, t, s, GFC_RND_MODE);
3614           mpfr_div_ui (t, t, 2, GFC_RND_MODE);
3615           mpfr_sqrt (t, t, GFC_RND_MODE);
3616           mpfr_sqrt (s, ad, GFC_RND_MODE);
3617           mpfr_mul (w, s, t, GFC_RND_MODE);
3618         }
3619
3620       if (mpfr_cmp_ui (w, 0) != 0 && mpfr_cmp_ui (e->value.complex.r, 0) >= 0)
3621         {
3622           mpfr_mul_ui (t, w, 2, GFC_RND_MODE);
3623           mpfr_div (result->value.complex.i, e->value.complex.i, t, GFC_RND_MODE);
3624           mpfr_set (result->value.complex.r, w, GFC_RND_MODE);
3625         }
3626       else if (mpfr_cmp_ui (w, 0) != 0
3627                && mpfr_cmp_ui (e->value.complex.r, 0) < 0
3628                && mpfr_cmp_ui (e->value.complex.i, 0) >= 0)
3629         {
3630           mpfr_mul_ui (t, w, 2, GFC_RND_MODE);
3631           mpfr_div (result->value.complex.r, e->value.complex.i, t, GFC_RND_MODE);
3632           mpfr_set (result->value.complex.i, w, GFC_RND_MODE);
3633         }
3634       else if (mpfr_cmp_ui (w, 0) != 0
3635                && mpfr_cmp_ui (e->value.complex.r, 0) < 0
3636                && mpfr_cmp_ui (e->value.complex.i, 0) < 0)
3637         {
3638           mpfr_mul_ui (t, w, 2, GFC_RND_MODE);
3639           mpfr_div (result->value.complex.r, ad, t, GFC_RND_MODE);
3640           mpfr_neg (w, w, GFC_RND_MODE);
3641           mpfr_set (result->value.complex.i, w, GFC_RND_MODE);
3642         }
3643       else
3644         gfc_internal_error ("invalid complex argument of SQRT at %L",
3645                             &e->where);
3646
3647       mpfr_clear (s);
3648       mpfr_clear (t);
3649       mpfr_clear (ac);
3650       mpfr_clear (ad);
3651       mpfr_clear (w);
3652
3653       break;
3654
3655     default:
3656       gfc_internal_error ("invalid argument of SQRT at %L", &e->where);
3657     }
3658
3659   return range_check (result, "SQRT");
3660
3661 negative_arg:
3662   gfc_free_expr (result);
3663   gfc_error ("Argument of SQRT at %L has a negative value", &e->where);
3664   return &gfc_bad_expr;
3665 }
3666
3667
3668 gfc_expr *
3669 gfc_simplify_tan (gfc_expr * x)
3670 {
3671   int i;
3672   gfc_expr *result;
3673
3674   if (x->expr_type != EXPR_CONSTANT)
3675     return NULL;
3676
3677   i = gfc_validate_kind (BT_REAL, x->ts.kind, false);
3678
3679   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
3680
3681   mpfr_tan (result->value.real, x->value.real, GFC_RND_MODE);
3682
3683   return range_check (result, "TAN");
3684 }
3685
3686
3687 gfc_expr *
3688 gfc_simplify_tanh (gfc_expr * x)
3689 {
3690   gfc_expr *result;
3691
3692   if (x->expr_type != EXPR_CONSTANT)
3693     return NULL;
3694
3695   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
3696
3697   mpfr_tanh (result->value.real, x->value.real, GFC_RND_MODE);
3698
3699   return range_check (result, "TANH");
3700
3701 }
3702
3703
3704 gfc_expr *
3705 gfc_simplify_tiny (gfc_expr * e)
3706 {
3707   gfc_expr *result;
3708   int i;
3709
3710   i = gfc_validate_kind (BT_REAL, e->ts.kind, false);
3711
3712   result = gfc_constant_result (BT_REAL, e->ts.kind, &e->where);
3713   mpfr_set (result->value.real, gfc_real_kinds[i].tiny, GFC_RND_MODE);
3714
3715   return result;
3716 }
3717
3718
3719 gfc_expr *
3720 gfc_simplify_transfer (gfc_expr * source, gfc_expr *mold, gfc_expr * size)
3721 {
3722
3723   /* Reference mold and size to suppress warning.  */
3724   if (gfc_init_expr && (mold || size))
3725     gfc_error ("TRANSFER intrinsic not implemented for initialization at %L",
3726                &source->where);
3727
3728   return NULL;
3729 }
3730
3731
3732 gfc_expr *
3733 gfc_simplify_trim (gfc_expr * e)
3734 {
3735   gfc_expr *result;
3736   int count, i, len, lentrim;
3737
3738   if (e->expr_type != EXPR_CONSTANT)
3739     return NULL;
3740
3741   len = e->value.character.length;
3742
3743   result = gfc_constant_result (BT_CHARACTER, e->ts.kind, &e->where);
3744
3745   for (count = 0, i = 1; i <= len; ++i)
3746     {
3747       if (e->value.character.string[len - i] == ' ')
3748         count++;
3749       else
3750         break;
3751     }
3752
3753   lentrim = len - count;
3754
3755   result->value.character.length = lentrim;
3756   result->value.character.string = gfc_getmem (lentrim + 1);
3757