gcc/fortran/simplify.c

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