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].subnormal, GFC_RND_MODE);
2297       else
2298         mpfr_sub (result->value.real,
2299                   x->value.real, gfc_real_kinds[k].subnormal, GFC_RND_MODE);
2300     }
2301   else
2302     {
2303       if (sgn < 0)
2304         {
2305           direction = -direction;
2306           mpfr_neg (result->value.real, result->value.real, GFC_RND_MODE);
2307         }
2308
2309       if (direction > 0)
2310         mpfr_add_one_ulp (result->value.real, GFC_RND_MODE);
2311       else
2312         {
2313           /* In this case the exponent can shrink, which makes us skip
2314              over one number because we subtract one ulp with the
2315              larger exponent.  Thus we need to compensate for this.  */
2316           mpfr_init_set (tmp, result->value.real, GFC_RND_MODE);
2317
2318           mpfr_sub_one_ulp (result->value.real, GFC_RND_MODE);
2319           mpfr_add_one_ulp (result->value.real, GFC_RND_MODE);
2320
2321           /* If we're back to where we started, the spacing is one
2322              ulp, and we get the correct result by subtracting.  */
2323           if (mpfr_cmp (tmp, result->value.real) == 0)
2324             mpfr_sub_one_ulp (result->value.real, GFC_RND_MODE);
2325
2326           mpfr_clear (tmp);
2327         }
2328
2329       if (sgn < 0)
2330         mpfr_neg (result->value.real, result->value.real, GFC_RND_MODE);
2331     }
2332
2333   return range_check (result, "NEAREST");
2334 }
2335
2336
2337 static gfc_expr *
2338 simplify_nint (const char *name, gfc_expr * e, gfc_expr * k)
2339 {
2340   gfc_expr *itrunc, *result;
2341   int kind;
2342
2343   kind = get_kind (BT_INTEGER, k, name, gfc_default_integer_kind);
2344   if (kind == -1)
2345     return &gfc_bad_expr;
2346
2347   if (e->expr_type != EXPR_CONSTANT)
2348     return NULL;
2349
2350   result = gfc_constant_result (BT_INTEGER, kind, &e->where);
2351
2352   itrunc = gfc_copy_expr (e);
2353
2354   mpfr_round(itrunc->value.real, e->value.real);
2355
2356   gfc_mpfr_to_mpz (result->value.integer, itrunc->value.real);
2357
2358   gfc_free_expr (itrunc);
2359
2360   return range_check (result, name);
2361 }
2362
2363
2364 gfc_expr *
2365 gfc_simplify_nint (gfc_expr * e, gfc_expr * k)
2366 {
2367   return simplify_nint ("NINT", e, k);
2368 }
2369
2370
2371 gfc_expr *
2372 gfc_simplify_idnint (gfc_expr * e)
2373 {
2374   return simplify_nint ("IDNINT", e, NULL);
2375 }
2376
2377
2378 gfc_expr *
2379 gfc_simplify_not (gfc_expr * e)
2380 {
2381   gfc_expr *result;
2382   int i;
2383
2384   if (e->expr_type != EXPR_CONSTANT)
2385     return NULL;
2386
2387   result = gfc_constant_result (e->ts.type, e->ts.kind, &e->where);
2388
2389   mpz_com (result->value.integer, e->value.integer);
2390
2391   /* Because of how GMP handles numbers, the result must be ANDed with
2392      the max_int mask.  For radices <> 2, this will require change.  */
2393
2394   i = gfc_validate_kind (BT_INTEGER, e->ts.kind, false);
2395
2396   mpz_and (result->value.integer, result->value.integer,
2397            gfc_integer_kinds[i].max_int);
2398
2399   return range_check (result, "NOT");
2400 }
2401
2402
2403 gfc_expr *
2404 gfc_simplify_null (gfc_expr * mold)
2405 {
2406   gfc_expr *result;
2407
2408   result = gfc_get_expr ();
2409   result->expr_type = EXPR_NULL;
2410
2411   if (mold == NULL)
2412     result->ts.type = BT_UNKNOWN;
2413   else
2414     {
2415       result->ts = mold->ts;
2416       result->where = mold->where;
2417     }
2418
2419   return result;
2420 }
2421
2422
2423 gfc_expr *
2424 gfc_simplify_precision (gfc_expr * e)
2425 {
2426   gfc_expr *result;
2427   int i;
2428
2429   i = gfc_validate_kind (e->ts.type, e->ts.kind, false);
2430
2431   result = gfc_int_expr (gfc_real_kinds[i].precision);
2432   result->where = e->where;
2433
2434   return result;
2435 }
2436
2437
2438 gfc_expr *
2439 gfc_simplify_radix (gfc_expr * e)
2440 {
2441   gfc_expr *result;
2442   int i;
2443
2444   i = gfc_validate_kind (e->ts.type, e->ts.kind, false);
2445   switch (e->ts.type)
2446     {
2447     case BT_INTEGER:
2448       i = gfc_integer_kinds[i].radix;
2449       break;
2450
2451     case BT_REAL:
2452       i = gfc_real_kinds[i].radix;
2453       break;
2454
2455     default:
2456       gcc_unreachable ();
2457     }
2458
2459   result = gfc_int_expr (i);
2460   result->where = e->where;
2461
2462   return result;
2463 }
2464
2465
2466 gfc_expr *
2467 gfc_simplify_range (gfc_expr * e)
2468 {
2469   gfc_expr *result;
2470   int i;
2471   long j;
2472
2473   i = gfc_validate_kind (e->ts.type, e->ts.kind, false);
2474
2475   switch (e->ts.type)
2476     {
2477     case BT_INTEGER:
2478       j = gfc_integer_kinds[i].range;
2479       break;
2480
2481     case BT_REAL:
2482     case BT_COMPLEX:
2483       j = gfc_real_kinds[i].range;
2484       break;
2485
2486     default:
2487       gcc_unreachable ();
2488     }
2489
2490   result = gfc_int_expr (j);
2491   result->where = e->where;
2492
2493   return result;
2494 }
2495
2496
2497 gfc_expr *
2498 gfc_simplify_real (gfc_expr * e, gfc_expr * k)
2499 {
2500   gfc_expr *result;
2501   int kind;
2502
2503   if (e->ts.type == BT_COMPLEX)
2504     kind = get_kind (BT_REAL, k, "REAL", e->ts.kind);
2505   else
2506     kind = get_kind (BT_REAL, k, "REAL", gfc_default_real_kind);
2507
2508   if (kind == -1)
2509     return &gfc_bad_expr;
2510
2511   if (e->expr_type != EXPR_CONSTANT)
2512     return NULL;
2513
2514   switch (e->ts.type)
2515     {
2516     case BT_INTEGER:
2517       result = gfc_int2real (e, kind);
2518       break;
2519
2520     case BT_REAL:
2521       result = gfc_real2real (e, kind);
2522       break;
2523
2524     case BT_COMPLEX:
2525       result = gfc_complex2real (e, kind);
2526       break;
2527
2528     default:
2529       gfc_internal_error ("bad type in REAL");
2530       /* Not reached */
2531     }
2532
2533   return range_check (result, "REAL");
2534 }
2535
2536 gfc_expr *
2537 gfc_simplify_repeat (gfc_expr * e, gfc_expr * n)
2538 {
2539   gfc_expr *result;
2540   int i, j, len, ncopies, nlen;
2541
2542   if (e->expr_type != EXPR_CONSTANT || n->expr_type != EXPR_CONSTANT)
2543     return NULL;
2544
2545   if (n != NULL && (gfc_extract_int (n, &ncopies) != NULL || ncopies < 0))
2546     {
2547       gfc_error ("Invalid second argument of REPEAT at %L", &n->where);
2548       return &gfc_bad_expr;
2549     }
2550
2551   len = e->value.character.length;
2552   nlen = ncopies * len;
2553
2554   result = gfc_constant_result (BT_CHARACTER, e->ts.kind, &e->where);
2555
2556   if (ncopies == 0)
2557     {
2558       result->value.character.string = gfc_getmem (1);
2559       result->value.character.length = 0;
2560       result->value.character.string[0] = '\0';
2561       return result;
2562     }
2563
2564   result->value.character.length = nlen;
2565   result->value.character.string = gfc_getmem (nlen + 1);
2566
2567   for (i = 0; i < ncopies; i++)
2568     for (j = 0; j < len; j++)
2569       result->value.character.string[j + i * len] =
2570         e->value.character.string[j];
2571
2572   result->value.character.string[nlen] = '\0';  /* For debugger */
2573   return result;
2574 }
2575
2576
2577 /* This one is a bear, but mainly has to do with shuffling elements.  */
2578
2579 gfc_expr *
2580 gfc_simplify_reshape (gfc_expr * source, gfc_expr * shape_exp,
2581                       gfc_expr * pad, gfc_expr * order_exp)
2582 {
2583
2584   int order[GFC_MAX_DIMENSIONS], shape[GFC_MAX_DIMENSIONS];
2585   int i, rank, npad, x[GFC_MAX_DIMENSIONS];
2586   gfc_constructor *head, *tail;
2587   mpz_t index, size;
2588   unsigned long j;
2589   size_t nsource;
2590   gfc_expr *e;
2591
2592   /* Unpack the shape array.  */
2593   if (source->expr_type != EXPR_ARRAY || !gfc_is_constant_expr (source))
2594     return NULL;
2595
2596   if (shape_exp->expr_type != EXPR_ARRAY || !gfc_is_constant_expr (shape_exp))
2597     return NULL;
2598
2599   if (pad != NULL
2600       && (pad->expr_type != EXPR_ARRAY
2601           || !gfc_is_constant_expr (pad)))
2602     return NULL;
2603
2604   if (order_exp != NULL
2605       && (order_exp->expr_type != EXPR_ARRAY
2606           || !gfc_is_constant_expr (order_exp)))
2607     return NULL;
2608
2609   mpz_init (index);
2610   rank = 0;
2611   head = tail = NULL;
2612
2613   for (;;)
2614     {
2615       e = gfc_get_array_element (shape_exp, rank);
2616       if (e == NULL)
2617         break;
2618
2619       if (gfc_extract_int (e, &shape[rank]) != NULL)
2620         {
2621           gfc_error ("Integer too large in shape specification at %L",
2622                      &e->where);
2623           gfc_free_expr (e);
2624           goto bad_reshape;
2625         }
2626
2627       gfc_free_expr (e);
2628
2629       if (rank >= GFC_MAX_DIMENSIONS)
2630         {
2631           gfc_error ("Too many dimensions in shape specification for RESHAPE "
2632                      "at %L", &e->where);
2633
2634           goto bad_reshape;
2635         }
2636
2637       if (shape[rank] < 0)
2638         {
2639           gfc_error ("Shape specification at %L cannot be negative",
2640                      &e->where);
2641           goto bad_reshape;
2642         }
2643
2644       rank++;
2645     }
2646
2647   if (rank == 0)
2648     {
2649       gfc_error ("Shape specification at %L cannot be the null array",
2650                  &shape_exp->where);
2651       goto bad_reshape;
2652     }
2653
2654   /* Now unpack the order array if present.  */
2655   if (order_exp == NULL)
2656     {
2657       for (i = 0; i < rank; i++)
2658         order[i] = i;
2659
2660     }
2661   else
2662     {
2663
2664       for (i = 0; i < rank; i++)
2665         x[i] = 0;
2666
2667       for (i = 0; i < rank; i++)
2668         {
2669           e = gfc_get_array_element (order_exp, i);
2670           if (e == NULL)
2671             {
2672               gfc_error
2673                 ("ORDER parameter of RESHAPE at %L is not the same size "
2674                  "as SHAPE parameter", &order_exp->where);
2675               goto bad_reshape;
2676             }
2677
2678           if (gfc_extract_int (e, &order[i]) != NULL)
2679             {
2680               gfc_error ("Error in ORDER parameter of RESHAPE at %L",
2681                          &e->where);
2682               gfc_free_expr (e);
2683               goto bad_reshape;
2684             }
2685
2686           gfc_free_expr (e);
2687
2688           if (order[i] < 1 || order[i] > rank)
2689             {
2690               gfc_error ("ORDER parameter of RESHAPE at %L is out of range",
2691                          &e->where);
2692               goto bad_reshape;
2693             }
2694
2695           order[i]--;
2696
2697           if (x[order[i]])
2698             {
2699               gfc_error ("Invalid permutation in ORDER parameter at %L",
2700                          &e->where);
2701               goto bad_reshape;
2702             }
2703
2704           x[order[i]] = 1;
2705         }
2706     }
2707
2708   /* Count the elements in the source and padding arrays.  */
2709
2710   npad = 0;
2711   if (pad != NULL)
2712     {
2713       gfc_array_size (pad, &size);
2714       npad = mpz_get_ui (size);
2715       mpz_clear (size);
2716     }
2717
2718   gfc_array_size (source, &size);
2719   nsource = mpz_get_ui (size);
2720   mpz_clear (size);
2721
2722   /* If it weren't for that pesky permutation we could just loop
2723      through the source and round out any shortage with pad elements.
2724      But no, someone just had to have the compiler do something the
2725      user should be doing.  */
2726
2727   for (i = 0; i < rank; i++)
2728     x[i] = 0;
2729
2730   for (;;)
2731     {
2732       /* Figure out which element to extract.  */
2733       mpz_set_ui (index, 0);
2734
2735       for (i = rank - 1; i >= 0; i--)
2736         {
2737           mpz_add_ui (index, index, x[order[i]]);
2738           if (i != 0)
2739             mpz_mul_ui (index, index, shape[order[i - 1]]);
2740         }
2741
2742       if (mpz_cmp_ui (index, INT_MAX) > 0)
2743         gfc_internal_error ("Reshaped array too large at %L", &e->where);
2744
2745       j = mpz_get_ui (index);
2746
2747       if (j < nsource)
2748         e = gfc_get_array_element (source, j);
2749       else
2750         {
2751           j = j - nsource;
2752
2753           if (npad == 0)
2754             {
2755               gfc_error
2756                 ("PAD parameter required for short SOURCE parameter at %L",
2757                  &source->where);
2758               goto bad_reshape;
2759             }
2760
2761           j = j % npad;
2762           e = gfc_get_array_element (pad, j);
2763         }
2764
2765       if (head == NULL)
2766         head = tail = gfc_get_constructor ();
2767       else
2768         {
2769           tail->next = gfc_get_constructor ();
2770           tail = tail->next;
2771         }
2772
2773       if (e == NULL)
2774         goto bad_reshape;
2775
2776       tail->where = e->where;
2777       tail->expr = e;
2778
2779       /* Calculate the next element.  */
2780       i = 0;
2781
2782 inc:
2783       if (++x[i] < shape[i])
2784         continue;
2785       x[i++] = 0;
2786       if (i < rank)
2787         goto inc;
2788
2789       break;
2790     }
2791
2792   mpz_clear (index);
2793
2794   e = gfc_get_expr ();
2795   e->where = source->where;
2796   e->expr_type = EXPR_ARRAY;
2797   e->value.constructor = head;
2798   e->shape = gfc_get_shape (rank);
2799
2800   for (i = 0; i < rank; i++)
2801     mpz_init_set_ui (e->shape[i], shape[i]);
2802
2803   e->ts = head->expr->ts;
2804   e->rank = rank;
2805
2806   return e;
2807
2808 bad_reshape:
2809   gfc_free_constructor (head);
2810   mpz_clear (index);
2811   return &gfc_bad_expr;
2812 }
2813
2814
2815 gfc_expr *
2816 gfc_simplify_rrspacing (gfc_expr * x)
2817 {
2818   gfc_expr *result;
2819   mpfr_t absv, log2, exp, frac, pow2;
2820   int i, p;
2821
2822   if (x->expr_type != EXPR_CONSTANT)
2823     return NULL;
2824
2825   i = gfc_validate_kind (x->ts.type, x->ts.kind, false);
2826
2827   result = gfc_constant_result (BT_REAL, x->ts.kind, &x->where);
2828
2829   p = gfc_real_kinds[i].digits;
2830
2831   gfc_set_model_kind (x->ts.kind);
2832
2833   if (mpfr_sgn (x->value.real) == 0)
2834     {
2835       mpfr_ui_div (result->value.real, 1, gfc_real_kinds[i].tiny, GFC_RND_MODE);
2836       return result;
2837     }
2838
2839   mpfr_init (log2);
2840   mpfr_init (absv);
2841   mpfr_init (frac);
2842   mpfr_init (pow2);
2843
2844   mpfr_abs (absv, x->value.real, GFC_RND_MODE);
2845   mpfr_log2 (log2, absv, GFC_RND_MODE);
2846
2847   mpfr_trunc (log2, log2);
2848   mpfr_add_ui (exp, log2, 1, GFC_RND_MODE);
2849
2850   mpfr_ui_pow (pow2, 2, exp, GFC_RND_MODE);
2851   mpfr_div (frac, absv, pow2, GFC_RND_MODE);
2852
2853   mpfr_mul_2exp (result->value.real, frac, (unsigned long)p, GFC_RND_MODE);
2854
2855   mpfr_clear (log2);
2856   mpfr_clear (absv);
2857   mpfr_clear (frac);
2858   mpfr_clear (pow2);
2859
2860   return range_check (result, "RRSPACING");
2861 }
2862
2863
2864 gfc_expr *
2865 gfc_simplify_scale (gfc_expr * x, gfc_expr * i)
2866 {
2867   int k, neg_flag, power, exp_range;
2868   mpfr_t scale, radix;
2869   gfc_expr *result;
2870
2871   if (x->expr_type != EXPR_CONSTANT || i->expr_type != EXPR_CONSTANT)
2872     return NULL;
2873
2874   result = gfc_constant_result (BT_REAL, x->ts.kind, &x->where);
2875
2876   if (mpfr_sgn (x->value.real) == 0)
2877     {
2878       mpfr_set_ui (result->value.real, 0, GFC_RND_MODE);
2879       return result;
2880     }
2881
2882   k = gfc_validate_kind (BT_REAL, x->ts.kind, false);
2883
2884   exp_range = gfc_real_kinds[k].max_exponent - gfc_real_kinds[k].min_exponent;
2885
2886   /* This check filters out values of i that would overflow an int.  */
2887   if (mpz_cmp_si (i->value.integer, exp_range + 2) > 0
2888       || mpz_cmp_si (i->value.integer, -exp_range - 2) < 0)
2889     {
2890       gfc_error ("Result of SCALE overflows its kind at %L", &result->where);
2891       return &gfc_bad_expr;
2892     }
2893
2894   /* Compute scale = radix ** power.  */
2895   power = mpz_get_si (i->value.integer);
2896
2897   if (power >= 0)
2898     neg_flag = 0;
2899   else
2900     {
2901       neg_flag = 1;
2902       power = -power;
2903     }
2904
2905   gfc_set_model_kind (x->ts.kind);
2906   mpfr_init (scale);
2907   mpfr_init (radix);
2908   mpfr_set_ui (radix, gfc_real_kinds[k].radix, GFC_RND_MODE);
2909   mpfr_pow_ui (scale, radix, power, GFC_RND_MODE);
2910
2911   if (neg_flag)
2912     mpfr_div (result->value.real, x->value.real, scale, GFC_RND_MODE);
2913   else
2914     mpfr_mul (result->value.real, x->value.real, scale, GFC_RND_MODE);
2915
2916   mpfr_clear (scale);
2917   mpfr_clear (radix);
2918
2919   return range_check (result, "SCALE");
2920 }
2921
2922
2923 gfc_expr *
2924 gfc_simplify_scan (gfc_expr * e, gfc_expr * c, gfc_expr * b)
2925 {
2926   gfc_expr *result;
2927   int back;
2928   size_t i;
2929   size_t indx, len, lenc;
2930
2931   if (e->expr_type != EXPR_CONSTANT || c->expr_type != EXPR_CONSTANT)
2932     return NULL;
2933
2934   if (b != NULL && b->value.logical != 0)
2935     back = 1;
2936   else
2937     back = 0;
2938
2939   result = gfc_constant_result (BT_INTEGER, gfc_default_integer_kind,
2940                                 &e->where);
2941
2942   len = e->value.character.length;
2943   lenc = c->value.character.length;
2944
2945   if (len == 0 || lenc == 0)
2946     {
2947       indx = 0;
2948     }
2949   else
2950     {
2951       if (back == 0)
2952         {
2953           indx =
2954             strcspn (e->value.character.string, c->value.character.string) + 1;
2955           if (indx > len)
2956             indx = 0;
2957         }
2958       else
2959         {
2960           i = 0;
2961           for (indx = len; indx > 0; indx--)
2962             {
2963               for (i = 0; i < lenc; i++)
2964                 {
2965                   if (c->value.character.string[i]
2966                         == e->value.character.string[indx - 1])
2967                     break;
2968                 }
2969               if (i < lenc)
2970                 break;
2971             }
2972         }
2973     }
2974   mpz_set_ui (result->value.integer, indx);
2975   return range_check (result, "SCAN");
2976 }
2977
2978
2979 gfc_expr *
2980 gfc_simplify_selected_int_kind (gfc_expr * e)
2981 {
2982   int i, kind, range;
2983   gfc_expr *result;
2984
2985   if (e->expr_type != EXPR_CONSTANT || gfc_extract_int (e, &range) != NULL)
2986     return NULL;
2987
2988   kind = INT_MAX;
2989
2990   for (i = 0; gfc_integer_kinds[i].kind != 0; i++)
2991     if (gfc_integer_kinds[i].range >= range
2992         && gfc_integer_kinds[i].kind < kind)
2993       kind = gfc_integer_kinds[i].kind;
2994
2995   if (kind == INT_MAX)
2996     kind = -1;
2997
2998   result = gfc_int_expr (kind);
2999   result->where = e->where;
3000
3001   return result;
3002 }
3003
3004
3005 gfc_expr *
3006 gfc_simplify_selected_real_kind (gfc_expr * p, gfc_expr * q)
3007 {
3008   int range, precision, i, kind, found_precision, found_range;
3009   gfc_expr *result;
3010
3011   if (p == NULL)
3012     precision = 0;
3013   else
3014     {
3015       if (p->expr_type != EXPR_CONSTANT
3016           || gfc_extract_int (p, &precision) != NULL)
3017         return NULL;
3018     }
3019
3020   if (q == NULL)
3021     range = 0;
3022   else
3023     {
3024       if (q->expr_type != EXPR_CONSTANT
3025           || gfc_extract_int (q, &range) != NULL)
3026         return NULL;
3027     }
3028
3029   kind = INT_MAX;
3030   found_precision = 0;
3031   found_range = 0;
3032
3033   for (i = 0; gfc_real_kinds[i].kind != 0; i++)
3034     {
3035       if (gfc_real_kinds[i].precision >= precision)
3036         found_precision = 1;
3037
3038       if (gfc_real_kinds[i].range >= range)
3039         found_range = 1;
3040
3041       if (gfc_real_kinds[i].precision >= precision
3042           && gfc_real_kinds[i].range >= range && gfc_real_kinds[i].kind < kind)
3043         kind = gfc_real_kinds[i].kind;
3044     }
3045
3046   if (kind == INT_MAX)
3047     {
3048       kind = 0;
3049
3050       if (!found_precision)
3051         kind = -1;
3052       if (!found_range)
3053         kind -= 2;
3054     }
3055
3056   result = gfc_int_expr (kind);
3057   result->where = (p != NULL) ? p->where : q->where;
3058
3059   return result;
3060 }
3061
3062
3063 gfc_expr *
3064 gfc_simplify_set_exponent (gfc_expr * x, gfc_expr * i)
3065 {
3066   gfc_expr *result;
3067   mpfr_t exp, absv, log2, pow2, frac;
3068   unsigned long exp2;
3069
3070   if (x->expr_type != EXPR_CONSTANT || i->expr_type != EXPR_CONSTANT)
3071     return NULL;
3072
3073   result = gfc_constant_result (BT_REAL, x->ts.kind, &x->where);
3074
3075   gfc_set_model_kind (x->ts.kind);
3076
3077   if (mpfr_sgn (x->value.real) == 0)
3078     {
3079       mpfr_set_ui (result->value.real, 0, GFC_RND_MODE);
3080       return result;
3081     }
3082
3083   mpfr_init (absv);
3084   mpfr_init (log2);
3085   mpfr_init (exp);
3086   mpfr_init (pow2);
3087   mpfr_init (frac);
3088
3089   mpfr_abs (absv, x->value.real, GFC_RND_MODE);
3090   mpfr_log2 (log2, absv, GFC_RND_MODE);
3091
3092   mpfr_trunc (log2, log2);
3093   mpfr_add_ui (exp, log2, 1, GFC_RND_MODE);
3094
3095   /* Old exponent value, and fraction.  */
3096   mpfr_ui_pow (pow2, 2, exp, GFC_RND_MODE);
3097
3098   mpfr_div (frac, absv, pow2, GFC_RND_MODE);
3099
3100   /* New exponent.  */
3101   exp2 = (unsigned long) mpz_get_d (i->value.integer);
3102   mpfr_mul_2exp (result->value.real, frac, exp2, GFC_RND_MODE);
3103
3104   mpfr_clear (absv);
3105   mpfr_clear (log2);
3106   mpfr_clear (pow2);
3107   mpfr_clear (frac);
3108
3109   return range_check (result, "SET_EXPONENT");
3110 }
3111
3112
3113 gfc_expr *
3114 gfc_simplify_shape (gfc_expr * source)
3115 {
3116   mpz_t shape[GFC_MAX_DIMENSIONS];
3117   gfc_expr *result, *e, *f;
3118   gfc_array_ref *ar;
3119   int n;
3120   try t;
3121
3122   if (source->rank == 0 || source->expr_type != EXPR_VARIABLE)
3123     return NULL;
3124
3125   result = gfc_start_constructor (BT_INTEGER, gfc_default_integer_kind,
3126                                   &source->where);
3127
3128   ar = gfc_find_array_ref (source);
3129
3130   t = gfc_array_ref_shape (ar, shape);
3131
3132   for (n = 0; n < source->rank; n++)
3133     {
3134       e = gfc_constant_result (BT_INTEGER, gfc_default_integer_kind,
3135                                &source->where);
3136
3137       if (t == SUCCESS)
3138         {
3139           mpz_set (e->value.integer, shape[n]);
3140           mpz_clear (shape[n]);
3141         }
3142       else
3143         {
3144           mpz_set_ui (e->value.integer, n + 1);
3145
3146           f = gfc_simplify_size (source, e);
3147           gfc_free_expr (e);
3148           if (f == NULL)
3149             {
3150               gfc_free_expr (result);
3151               return NULL;
3152             }
3153           else
3154             {
3155               e = f;
3156             }
3157         }
3158
3159       gfc_append_constructor (result, e);
3160     }
3161
3162   return result;
3163 }
3164
3165
3166 gfc_expr *
3167 gfc_simplify_size (gfc_expr * array, gfc_expr * dim)
3168 {
3169   mpz_t size;
3170   gfc_expr *result;
3171   int d;
3172
3173   if (dim == NULL)
3174     {
3175       if (gfc_array_size (array, &size) == FAILURE)
3176         return NULL;
3177     }
3178   else
3179     {
3180       if (dim->expr_type != EXPR_CONSTANT)
3181         return NULL;
3182
3183       d = mpz_get_ui (dim->value.integer) - 1;
3184       if (gfc_array_dimen_size (array, d, &size) == FAILURE)
3185         return NULL;
3186     }
3187
3188   result = gfc_constant_result (BT_INTEGER, gfc_default_integer_kind,
3189                                 &array->where);
3190
3191   mpz_set (result->value.integer, size);
3192
3193   return result;
3194 }
3195
3196
3197 gfc_expr *
3198 gfc_simplify_sign (gfc_expr * x, gfc_expr * y)
3199 {
3200   gfc_expr *result;
3201
3202   if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
3203     return NULL;
3204
3205   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
3206
3207   switch (x->ts.type)
3208     {
3209     case BT_INTEGER:
3210       mpz_abs (result->value.integer, x->value.integer);
3211       if (mpz_sgn (y->value.integer) < 0)
3212         mpz_neg (result->value.integer, result->value.integer);
3213
3214       break;
3215
3216     case BT_REAL:
3217       /* TODO: Handle -0.0 and +0.0 correctly on machines that support
3218          it.  */
3219       mpfr_abs (result->value.real, x->value.real, GFC_RND_MODE);
3220       if (mpfr_sgn (y->value.real) < 0)
3221         mpfr_neg (result->value.real, result->value.real, GFC_RND_MODE);
3222
3223       break;
3224
3225     default:
3226       gfc_internal_error ("Bad type in gfc_simplify_sign");
3227     }
3228
3229   return result;
3230 }
3231
3232
3233 gfc_expr *
3234 gfc_simplify_sin (gfc_expr * x)
3235 {
3236   gfc_expr *result;
3237   mpfr_t xp, xq;
3238
3239   if (x->expr_type != EXPR_CONSTANT)
3240     return NULL;
3241
3242   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
3243
3244   switch (x->ts.type)
3245     {
3246     case BT_REAL:
3247       mpfr_sin (result->value.real, x->value.real, GFC_RND_MODE);
3248       break;
3249
3250     case BT_COMPLEX:
3251       gfc_set_model (x->value.real);
3252       mpfr_init (xp);
3253       mpfr_init (xq);
3254
3255       mpfr_sin  (xp, x->value.complex.r, GFC_RND_MODE);
3256       mpfr_cosh (xq, x->value.complex.i, GFC_RND_MODE);
3257       mpfr_mul (result->value.complex.r, xp, xq, GFC_RND_MODE);
3258
3259       mpfr_cos  (xp, x->value.complex.r, GFC_RND_MODE);
3260       mpfr_sinh (xq, x->value.complex.i, GFC_RND_MODE);
3261       mpfr_mul (result->value.complex.i, xp, xq, GFC_RND_MODE);
3262
3263       mpfr_clear (xp);
3264       mpfr_clear (xq);
3265       break;
3266
3267     default:
3268       gfc_internal_error ("in gfc_simplify_sin(): Bad type");
3269     }
3270
3271   return range_check (result, "SIN");
3272 }
3273
3274
3275 gfc_expr *
3276 gfc_simplify_sinh (gfc_expr * x)
3277 {
3278   gfc_expr *result;
3279
3280   if (x->expr_type != EXPR_CONSTANT)
3281     return NULL;
3282
3283   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
3284
3285   mpfr_sinh(result->value.real, x->value.real, GFC_RND_MODE);
3286
3287   return range_check (result, "SINH");
3288 }
3289
3290
3291 /* The argument is always a double precision real that is converted to
3292    single precision.  TODO: Rounding!  */
3293
3294 gfc_expr *
3295 gfc_simplify_sngl (gfc_expr * a)
3296 {
3297   gfc_expr *result;
3298
3299   if (a->expr_type != EXPR_CONSTANT)
3300     return NULL;
3301
3302   result = gfc_real2real (a, gfc_default_real_kind);
3303   return range_check (result, "SNGL");
3304 }
3305
3306
3307 gfc_expr *
3308 gfc_simplify_spacing (gfc_expr * x)
3309 {
3310   gfc_expr *result;
3311   mpfr_t absv, log2;
3312   long diff;
3313   int i, p;
3314
3315   if (x->expr_type != EXPR_CONSTANT)
3316     return NULL;
3317
3318   i = gfc_validate_kind (x->ts.type, x->ts.kind, false);
3319
3320   p = gfc_real_kinds[i].digits;
3321
3322   result = gfc_constant_result (BT_REAL, x->ts.kind, &x->where);
3323
3324   gfc_set_model_kind (x->ts.kind);
3325
3326   if (mpfr_sgn (x->value.real) == 0)
3327     {
3328       mpfr_set (result->value.real, gfc_real_kinds[i].tiny, GFC_RND_MODE);
3329       return result;
3330     }
3331
3332   mpfr_init (log2);
3333   mpfr_init (absv);
3334
3335   mpfr_abs (absv, x->value.real, GFC_RND_MODE);
3336   mpfr_log2 (log2, absv, GFC_RND_MODE);
3337   mpfr_trunc (log2, log2);
3338
3339   mpfr_add_ui (log2, log2, 1, GFC_RND_MODE);
3340
3341   /* FIXME: We should be using mpfr_get_si here, but this function is
3342      not available with the version of mpfr distributed with gmp (as of
3343      2004-09-17). Replace once mpfr has been imported into the gcc cvs
3344      tree.  */
3345   diff = (long)mpfr_get_d (log2, GFC_RND_MODE) - (long)p;
3346   mpfr_set_ui (result->value.real, 1, GFC_RND_MODE);
3347   mpfr_mul_2si (result->value.real, result->value.real, diff, GFC_RND_MODE);
3348
3349   mpfr_clear (log2);
3350   mpfr_clear (absv);
3351
3352   if (mpfr_cmp (result->value.real, gfc_real_kinds[i].tiny) < 0)
3353     mpfr_set (result->value.real, gfc_real_kinds[i].tiny, GFC_RND_MODE);
3354
3355   return range_check (result, "SPACING");
3356 }
3357
3358
3359 gfc_expr *
3360 gfc_simplify_sqrt (gfc_expr * e)
3361 {
3362   gfc_expr *result;
3363   mpfr_t ac, ad, s, t, w;
3364
3365   if (e->expr_type != EXPR_CONSTANT)
3366     return NULL;
3367
3368   result = gfc_constant_result (e->ts.type, e->ts.kind, &e->where);
3369
3370   switch (e->ts.type)
3371     {
3372     case BT_REAL:
3373       if (mpfr_cmp_si (e->value.real, 0) < 0)
3374         goto negative_arg;
3375       mpfr_sqrt (result->value.real, e->value.real, GFC_RND_MODE);
3376
3377       break;
3378
3379     case BT_COMPLEX:
3380       /* Formula taken from Numerical Recipes to avoid over- and
3381          underflow.  */
3382
3383       gfc_set_model (e->value.real);
3384       mpfr_init (ac);
3385       mpfr_init (ad);
3386       mpfr_init (s);
3387       mpfr_init (t);
3388       mpfr_init (w);
3389
3390       if (mpfr_cmp_ui (e->value.complex.r, 0) == 0
3391           && mpfr_cmp_ui (e->value.complex.i, 0) == 0)
3392         {
3393
3394           mpfr_set_ui (result->value.complex.r, 0, GFC_RND_MODE);
3395           mpfr_set_ui (result->value.complex.i, 0, GFC_RND_MODE);
3396           break;
3397         }
3398
3399       mpfr_abs (ac, e->value.complex.r, GFC_RND_MODE);
3400       mpfr_abs (ad, e->value.complex.i, GFC_RND_MODE);
3401
3402       if (mpfr_cmp (ac, ad) >= 0)
3403         {
3404           mpfr_div (t, e->value.complex.i, e->value.complex.r, GFC_RND_MODE);
3405           mpfr_mul (t, t, t, GFC_RND_MODE);
3406           mpfr_add_ui (t, t, 1, GFC_RND_MODE);
3407           mpfr_sqrt (t, t, GFC_RND_MODE);
3408           mpfr_add_ui (t, t, 1, GFC_RND_MODE);
3409           mpfr_div_ui (t, t, 2, GFC_RND_MODE);
3410           mpfr_sqrt (t, t, GFC_RND_MODE);
3411           mpfr_sqrt (s, ac, GFC_RND_MODE);
3412           mpfr_mul (w, s, t, GFC_RND_MODE);
3413         }
3414       else
3415         {
3416           mpfr_div (s, e->value.complex.r, e->value.complex.i, GFC_RND_MODE);
3417           mpfr_mul (t, s, s, GFC_RND_MODE);
3418           mpfr_add_ui (t, t, 1, GFC_RND_MODE);
3419           mpfr_sqrt (t, t, GFC_RND_MODE);
3420           mpfr_abs (s, s, GFC_RND_MODE);
3421           mpfr_add (t, t, s, GFC_RND_MODE);
3422           mpfr_div_ui (t, t, 2, GFC_RND_MODE);
3423           mpfr_sqrt (t, t, GFC_RND_MODE);
3424           mpfr_sqrt (s, ad, GFC_RND_MODE);
3425           mpfr_mul (w, s, t, GFC_RND_MODE);
3426         }
3427
3428       if (mpfr_cmp_ui (w, 0) != 0 && mpfr_cmp_ui (e->value.complex.r, 0) >= 0)
3429         {
3430           mpfr_mul_ui (t, w, 2, GFC_RND_MODE);
3431           mpfr_div (result->value.complex.i, e->value.complex.i, t, GFC_RND_MODE);
3432           mpfr_set (result->value.complex.r, w, GFC_RND_MODE);
3433         }
3434       else if (mpfr_cmp_ui (w, 0) != 0
3435                && mpfr_cmp_ui (e->value.complex.r, 0) < 0
3436                && mpfr_cmp_ui (e->value.complex.i, 0) >= 0)
3437         {
3438           mpfr_mul_ui (t, w, 2, GFC_RND_MODE);
3439           mpfr_div (result->value.complex.r, e->value.complex.i, t, GFC_RND_MODE);
3440           mpfr_set (result->value.complex.i, w, GFC_RND_MODE);
3441         }
3442       else if (mpfr_cmp_ui (w, 0) != 0
3443                && mpfr_cmp_ui (e->value.complex.r, 0) < 0
3444                && mpfr_cmp_ui (e->value.complex.i, 0) < 0)
3445         {
3446           mpfr_mul_ui (t, w, 2, GFC_RND_MODE);
3447           mpfr_div (result->value.complex.r, ad, t, GFC_RND_MODE);
3448           mpfr_neg (w, w, GFC_RND_MODE);
3449           mpfr_set (result->value.complex.i, w, GFC_RND_MODE);
3450         }
3451       else
3452         gfc_internal_error ("invalid complex argument of SQRT at %L",
3453                             &e->where);
3454
3455       mpfr_clear (s);
3456       mpfr_clear (t);
3457       mpfr_clear (ac);
3458       mpfr_clear (ad);
3459       mpfr_clear (w);
3460
3461       break;
3462
3463     default:
3464       gfc_internal_error ("invalid argument of SQRT at %L", &e->where);
3465     }
3466
3467   return range_check (result, "SQRT");
3468
3469 negative_arg:
3470   gfc_free_expr (result);
3471   gfc_error ("Argument of SQRT at %L has a negative value", &e->where);
3472   return &gfc_bad_expr;
3473 }
3474
3475
3476 gfc_expr *
3477 gfc_simplify_tan (gfc_expr * x)
3478 {
3479   int i;
3480   gfc_expr *result;
3481
3482   if (x->expr_type != EXPR_CONSTANT)
3483     return NULL;
3484
3485   i = gfc_validate_kind (BT_REAL, x->ts.kind, false);
3486
3487   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
3488
3489   mpfr_tan (result->value.real, x->value.real, GFC_RND_MODE);
3490
3491   return range_check (result, "TAN");
3492 }
3493
3494
3495 gfc_expr *
3496 gfc_simplify_tanh (gfc_expr * x)
3497 {
3498   gfc_expr *result;
3499
3500   if (x->expr_type != EXPR_CONSTANT)
3501     return NULL;
3502
3503   result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
3504
3505   mpfr_tanh (result->value.real, x->value.real, GFC_RND_MODE);
3506
3507   return range_check (result, "TANH");
3508
3509 }
3510
3511
3512 gfc_expr *
3513 gfc_simplify_tiny (gfc_expr * e)
3514 {
3515   gfc_expr *result;
3516   int i;
3517
3518   i = gfc_validate_kind (BT_REAL, e->ts.kind, false);
3519
3520   result = gfc_constant_result (BT_REAL, e->ts.kind, &e->where);
3521   mpfr_set (result->value.real, gfc_real_kinds[i].tiny, GFC_RND_MODE);
3522
3523   return result;
3524 }
3525
3526
3527 gfc_expr *
3528 gfc_simplify_trim (gfc_expr * e)
3529 {
3530   gfc_expr *result;
3531   int count, i, len, lentrim;
3532
3533   if (e->expr_type != EXPR_CONSTANT)
3534     return NULL;
3535
3536   len = e->value.character.length;
3537
3538   result = gfc_constant_result (BT_CHARACTER, e->ts.kind, &e->where);
3539
3540   for (count = 0, i = 1; i <= len; ++i)
3541     {
3542       if (e->value.character.string[len - i] == ' ')
3543         count++;
3544       else
3545         break;
3546     }
3547
3548   lentrim = len - count;
3549
3550   result->value.character.length = lentrim;
3551   result->value.character.string = gfc_getmem (lentrim + 1);
3552
3553   for (i = 0; i < lentrim; i++)
3554     result->value.character.string[i] = e->value.character.string[i];
3555
3556   result->value.character.string[lentrim] = '\0';       /* For debugger */
3557
3558   return result;
3559 }
3560
3561
3562 gfc_expr *
3563 gfc_simplify_ubound (gfc_expr * array, gfc_expr * dim)
3564 {
3565   return simplify_bound (array, dim, 1);
3566 }
3567
3568
3569 gfc_expr *
3570 gfc_simplify_verify (gfc_expr * s, gfc_expr * set, gfc_expr * b)
3571 {
3572   gfc_expr *result;
3573   int back;
3574   size_t index, len, lenset;
3575   size_t i;
3576
3577   if (s->expr_type != EXPR_CONSTANT || set->expr_type != EXPR_CONSTANT)
3578     return NULL;
3579
3580   if (b != NULL && b->value.logical != 0)
3581     back = 1;
3582   else
3583     back = 0;
3584
3585   result = gfc_constant_result (BT_INTEGER, gfc_default_integer_kind,
3586                                 &s->where);
3587
3588   len = s->value.character.length;
3589   lenset = set->value.character.length;
3590
3591   if (len == 0)
3592     {
3593       mpz_set_ui (result->value.integer, 0);
3594       return result;
3595     }
3596
3597   if (back == 0)
3598     {
3599       if (lenset == 0)
3600         {
3601           mpz_set_ui (result->value.integer, len);
3602           return result;
3603         }
3604
3605       index =
3606         strspn (s->value.character.string, set->value.character.string) + 1;
3607       if (index > len)
3608         index = 0;
3609
3610     }
3611   else
3612     {
3613       if (lenset == 0)
3614         {
3615           mpz_set_ui (result->value.integer, 1);
3616           return result;
3617         }
3618       for (index = len; index > 0; index --)
3619         {
3620           for (i = 0; i < lenset; i++)
3621             {
3622               if (s->value.character.string[index - 1]
3623                     == set->value.character.string[i])
3624                 break;
3625             }
3626           if (i == lenset)
3627             break;
3628         }
3629     }
3630
3631   mpz_set_ui (result->value.integer, index);
3632   return result;
3633 }
3634
3635 /****************** Constant simplification *****************/
3636
3637 /* Master function to convert one constant to another.  While this is
3638    used as a simplification function, it requires the destination type
3639    and kind information which is supplied by a special case in
3640    do_simplify().  */
3641
3642 gfc_expr *
3643 gfc_convert_constant (gfc_expr * e, bt type, int kind)
3644 {
3645   gfc_expr *g, *result, *(*f) (gfc_expr *, int);
3646   gfc_constructor *head, *c, *tail = NULL;
3647
3648   switch (e->ts.type)
3649     {
3650     case BT_INTEGER:
3651       switch (type)
3652         {
3653         case BT_INTEGER:
3654           f = gfc_int2int;
3655           break;
3656         case BT_REAL:
3657           f = gfc_int2real;
3658           break;
3659         case BT_COMPLEX:
3660           f = gfc_int2complex;
3661           break;
3662         default:
3663           goto oops;
3664         }
3665       break;
3666
3667     case BT_REAL:
3668       switch (type)
3669         {
3670         case BT_INTEGER:
3671           f = gfc_real2int;
3672           break;
3673         case BT_REAL:
3674           f = gfc_real2real;
3675           break;
3676         case BT_COMPLEX:
3677           f = gfc_real2complex;
3678           break;
3679         default:
3680           goto oops;
3681         }
3682       break;
3683
3684     case BT_COMPLEX:
3685       switch (type)
3686         {
3687         case BT_INTEGER:
3688           f = gfc_complex2int;
3689           break;
3690         case BT_REAL:
3691           f = gfc_complex2real;
3692           break;
3693         case BT_COMPLEX:
3694           f = gfc_complex2complex;
3695           break;
3696
3697         default:
3698           goto oops;
3699         }
3700       break;
3701
3702     case BT_LOGICAL:
3703       if (type != BT_LOGICAL)
3704         goto oops;
3705       f = gfc_log2log;
3706       break;
3707
3708     default:
3709     oops:
3710       gfc_internal_error ("gfc_convert_constant(): Unexpected type");
3711     }
3712
3713   result = NULL;
3714
3715   switch (e->expr_type)
3716     {
3717     case EXPR_CONSTANT:
3718       result = f (e, kind);
3719       if (result == NULL)
3720         return &gfc_bad_expr;
3721       break;
3722
3723     case EXPR_ARRAY:
3724       if (!gfc_is_constant_expr (e))
3725         break;
3726
3727       head = NULL;
3728
3729       for (c = e->value.constructor; c; c = c->next)
3730         {
3731           if (head == NULL)
3732             head = tail = gfc_get_constructor ();
3733           else
3734             {
3735               tail->next = gfc_get_constructor ();
3736               tail = tail->next;
3737             }
3738
3739           tail->where = c->where;
3740
3741           if (c->iterator == NULL)
3742             tail->expr = f (c->expr, kind);
3743           else
3744             {
3745               g = gfc_convert_constant (c->expr, type, kind);
3746               if (g == &gfc_bad_expr)
3747                 return g;
3748               tail->expr = g;
3749             }
3750
3751           if (tail->expr == NULL)
3752             {
3753               gfc_free_constructor (head);
3754               return NULL;
3755             }
3756         }
3757
3758       result = gfc_get_expr ();
3759       result->ts.type = type;
3760       result->ts.kind = kind;
3761       result->expr_type = EXPR_ARRAY;
3762       result->value.constructor = head;
3763       result->shape = gfc_copy_shape (e->shape, e->rank);
3764       result->where = e->where;
3765       result->rank = e->rank;
3766       break;
3767
3768     default:
3769       break;
3770     }
3771
3772   return result;
3773 }
3774
3775
3776 /****************** Helper functions ***********************/
3777
3778 /* Given a collating table, create the inverse table.  */
3779
3780 static void
3781 invert_table (const int *table, int *xtable)
3782 {
3783   int i;
3784
3785   for (i = 0; i < 256; i++)
3786     xtable[i] = 0;
3787
3788   for (i = 0; i < 256; i++)
3789     xtable[table[i]] = i;
3790 }
3791
3792
3793 void
3794 gfc_simplify_init_1 (void)
3795 {
3796
3797   invert_table (ascii_table, xascii_table);
3798 }