gcc/fortran/simplify.c

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