gcc/fortran/arith.c

   1 /* Compiler arithmetic
   2    Copyright (C) 2000, 2001, 2002, 2003, 2004 Free Software Foundation,
   3    Inc.
   4    Contributed by Andy Vaught
   5
   6 This file is part of GCC.
   7
   8 GCC is free software; you can redistribute it and/or modify it under
   9 the terms of the GNU General Public License as published by the Free
  10 Software Foundation; either version 2, or (at your option) any later
  11 version.
  12
  13 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
  14 WARRANTY; without even the implied warranty of MERCHANTABILITY or
  15 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  16 for more details.
  17
  18 You should have received a copy of the GNU General Public License
  19 along with GCC; see the file COPYING.  If not, write to the Free
  20 Software Foundation, 59 Temple Place - Suite 330, Boston, MA
  21 02111-1307, USA.  */
  22
  23 /* Since target arithmetic must be done on the host, there has to
  24    be some way of evaluating arithmetic expressions as the host
  25    would evaluate them.  We use the GNU MP library to do arithmetic,
  26    and this file provides the interface.  */
  27
  28 #include "config.h"
  29
  30 #include <string.h>
  31
  32 #include "gfortran.h"
  33 #include "arith.h"
  34
  35 mpf_t pi, half_pi, two_pi, e;
  36
  37 /* The gfc_(integer|real)_kinds[] structures have everything the front
  38    end needs to know about integers and real numbers on the target.
  39    Other entries of the structure are calculated from these values.
  40    The first entry is the default kind, the second entry of the real
  41    structure is the default double kind.  */
  42
  43 #define MPZ_NULL {{0,0,0}}
  44 #define MPF_NULL {{0,0,0,0}}
  45
  46 #define DEF_GFC_INTEGER_KIND(KIND,RADIX,DIGITS,BIT_SIZE)                \
  47         {KIND, RADIX, DIGITS, BIT_SIZE, 0, MPZ_NULL, MPZ_NULL, MPZ_NULL}
  48
  49 #define DEF_GFC_LOGICAL_KIND(KIND,BIT_SIZE)                             \
  50         {KIND, BIT_SIZE}
  51
  52 #define DEF_GFC_REAL_KIND(KIND,RADIX,DIGITS,MIN_EXP, MAX_EXP)           \
  53         {KIND, RADIX, DIGITS, MIN_EXP, MAX_EXP,                         \
  54          0, 0, MPF_NULL, MPF_NULL, MPF_NULL}
  55
  56 gfc_integer_info gfc_integer_kinds[] = {
  57   DEF_GFC_INTEGER_KIND (4, 2, 31, 32),
  58   DEF_GFC_INTEGER_KIND (8, 2, 63, 64),
  59   DEF_GFC_INTEGER_KIND (2, 2, 15, 16),
  60   DEF_GFC_INTEGER_KIND (1, 2,  7,  8),
  61   DEF_GFC_INTEGER_KIND (0, 0,  0,  0)
  62 };
  63
  64 gfc_logical_info gfc_logical_kinds[] = {
  65   DEF_GFC_LOGICAL_KIND (4, 32),
  66   DEF_GFC_LOGICAL_KIND (8, 64),
  67   DEF_GFC_LOGICAL_KIND (2, 16),
  68   DEF_GFC_LOGICAL_KIND (1,  8),
  69   DEF_GFC_LOGICAL_KIND (0,  0)
  70 };
  71
  72 gfc_real_info gfc_real_kinds[] = {
  73   DEF_GFC_REAL_KIND (4, 2, 24,  -125,  128),
  74   DEF_GFC_REAL_KIND (8, 2, 53, -1021, 1024),
  75   DEF_GFC_REAL_KIND (0, 0,  0,     0,    0)
  76 };
  77
  78
  79 /* The integer kind to use for array indices.  This will be set to the
  80    proper value based on target information from the backend.  */
  81
  82 int gfc_index_integer_kind;
  83
  84
  85 /* Compute the natural log of arg.
  86
  87    We first get the argument into the range 0.5 to 1.5 by successive
  88    multiplications or divisions by e.  Then we use the series:
  89
  90      ln(x) = (x-1) - (x-1)^/2 + (x-1)^3/3 - (x-1)^4/4 + ...
  91
  92    Because we are expanding in powers of (x-1), and 0.5 < x < 1.5, we
  93    have -0.5 < (x-1) < 0.5.  Ignoring the harmonic term, this means
  94    that each term is at most 1/(2^i), meaning one bit is gained per
  95    iteration.
  96
  97    Not very efficient, but it doesn't have to be.  */
  98
  99 void
 100 natural_logarithm (mpf_t * arg, mpf_t * result)
 101 {
 102   mpf_t x, xp, t, log;
 103   int i, p;
 104
 105   mpf_init_set (x, *arg);
 106   mpf_init (t);
 107
 108   p = 0;
 109
 110   mpf_set_str (t, "0.5", 10);
 111   while (mpf_cmp (x, t) < 0)
 112     {
 113       mpf_mul (x, x, e);
 114       p--;
 115     }
 116
 117   mpf_set_str (t, "1.5", 10);
 118   while (mpf_cmp (x, t) > 0)
 119     {
 120       mpf_div (x, x, e);
 121       p++;
 122     }
 123
 124   mpf_sub_ui (x, x, 1);
 125   mpf_init_set_ui (log, 0);
 126   mpf_init_set_ui (xp, 1);
 127
 128   for (i = 1; i < GFC_REAL_BITS; i++)
 129     {
 130       mpf_mul (xp, xp, x);
 131       mpf_div_ui (t, xp, i);
 132
 133       if (i % 2 == 0)
 134         mpf_sub (log, log, t);
 135       else
 136         mpf_add (log, log, t);
 137     }
 138
 139   /* Add in the log (e^p) = p */
 140
 141   if (p < 0)
 142     mpf_sub_ui (log, log, -p);
 143   else
 144     mpf_add_ui (log, log, p);
 145
 146   mpf_clear (x);
 147   mpf_clear (xp);
 148   mpf_clear (t);
 149
 150   mpf_set (*result, log);
 151   mpf_clear (log);
 152 }
 153
 154
 155 /* Calculate the common logarithm of arg.  We use the natural
 156    logaritm of arg and of 10:
 157
 158    log10(arg) = log(arg)/log(10)  */
 159
 160 void
 161 common_logarithm (mpf_t * arg, mpf_t * result)
 162 {
 163   mpf_t i10, log10;
 164
 165   natural_logarithm (arg, result);
 166
 167   mpf_init_set_ui (i10, 10);
 168   mpf_init (log10);
 169   natural_logarithm (&i10, &log10);
 170
 171   mpf_div (*result, *result, log10);
 172   mpf_clear (i10);
 173   mpf_clear (log10);
 174 }
 175
 176 /* Calculate exp(arg).
 177
 178    We use a reduction of the form
 179
 180      x = Nln2 + r
 181
 182    Then we obtain exp(r) from the McLaurin series.
 183    exp(x) is then recovered from the identity
 184
 185      exp(x) = 2^N*exp(r).  */
 186
 187 void
 188 exponential (mpf_t * arg, mpf_t * result)
 189 {
 190   mpf_t two, ln2, power, q, r, num, denom, term, x, xp;
 191   int i;
 192   long n;
 193   unsigned long p, mp;
 194
 195
 196   mpf_init_set (x, *arg);
 197
 198   if (mpf_cmp_ui (x, 0) == 0)
 199     {
 200       mpf_set_ui (*result, 1);
 201     }
 202   else if (mpf_cmp_ui (x, 1) == 0)
 203     {
 204       mpf_set (*result, e);
 205     }
 206   else
 207     {
 208       mpf_init_set_ui (two, 2);
 209       mpf_init (ln2);
 210       mpf_init (q);
 211       mpf_init (r);
 212       mpf_init (power);
 213       mpf_init (term);
 214
 215       natural_logarithm (&two, &ln2);
 216
 217       mpf_div (q, x, ln2);
 218       mpf_floor (power, q);
 219       mpf_mul (q, power, ln2);
 220       mpf_sub (r, x, q);
 221
 222       mpf_init_set_ui (xp, 1);
 223       mpf_init_set_ui (num, 1);
 224       mpf_init_set_ui (denom, 1);
 225
 226       for (i = 1; i <= GFC_REAL_BITS + 10; i++)
 227         {
 228           mpf_mul (num, num, r);
 229           mpf_mul_ui (denom, denom, i);
 230           mpf_div (term, num, denom);
 231           mpf_add (xp, xp, term);
 232         }
 233
 234       /* Reconstruction step */
 235       n = (long) mpf_get_d (power);
 236
 237       if (n > 0)
 238         {
 239           p = (unsigned int) n;
 240           mpf_mul_2exp (*result, xp, p);
 241         }
 242       else
 243         {
 244           mp = (unsigned int) (-n);
 245           mpf_div_2exp (*result, xp, mp);
 246         }
 247
 248       mpf_clear (two);
 249       mpf_clear (ln2);
 250       mpf_clear (q);
 251       mpf_clear (r);
 252       mpf_clear (power);
 253       mpf_clear (num);
 254       mpf_clear (denom);
 255       mpf_clear (term);
 256       mpf_clear (xp);
 257     }
 258
 259   mpf_clear (x);
 260 }
 261
 262
 263 /* Calculate sin(arg).
 264
 265    We use a reduction of the form
 266
 267      x= N*2pi + r
 268
 269    Then we obtain sin(r) from the McLaurin series.  */
 270
 271 void
 272 sine (mpf_t * arg, mpf_t * result)
 273 {
 274   mpf_t factor, q, r, num, denom, term, x, xp;
 275   int i, sign;
 276
 277   mpf_init_set (x, *arg);
 278
 279   /* Special case (we do not treat multiples of pi due to roundoff issues) */
 280   if (mpf_cmp_ui (x, 0) == 0)
 281     {
 282       mpf_set_ui (*result, 0);
 283     }
 284   else
 285     {
 286       mpf_init (q);
 287       mpf_init (r);
 288       mpf_init (factor);
 289       mpf_init (term);
 290
 291       mpf_div (q, x, two_pi);
 292       mpf_floor (factor, q);
 293       mpf_mul (q, factor, two_pi);
 294       mpf_sub (r, x, q);
 295
 296       mpf_init_set_ui (xp, 0);
 297       mpf_init_set_ui (num, 1);
 298       mpf_init_set_ui (denom, 1);
 299
 300       sign = -1;
 301       for (i = 1; i < GFC_REAL_BITS + 10; i++)
 302         {
 303           mpf_mul (num, num, r);
 304           mpf_mul_ui (denom, denom, i);
 305           if (i % 2 == 0)
 306             continue;
 307
 308           sign = -sign;
 309           mpf_div (term, num, denom);
 310           if (sign > 0)
 311             mpf_add (xp, xp, term);
 312           else
 313             mpf_sub (xp, xp, term);
 314         }
 315
 316       mpf_set (*result, xp);
 317
 318       mpf_clear (q);
 319       mpf_clear (r);
 320       mpf_clear (factor);
 321       mpf_clear (num);
 322       mpf_clear (denom);
 323       mpf_clear (term);
 324       mpf_clear (xp);
 325     }
 326
 327   mpf_clear (x);
 328 }
 329
 330
 331 /* Calculate cos(arg).
 332
 333    Similar to sine.  */
 334
 335 void
 336 cosine (mpf_t * arg, mpf_t * result)
 337 {
 338   mpf_t factor, q, r, num, denom, term, x, xp;
 339   int i, sign;
 340
 341   mpf_init_set (x, *arg);
 342
 343   /* Special case (we do not treat multiples of pi due to roundoff issues) */
 344   if (mpf_cmp_ui (x, 0) == 0)
 345     {
 346       mpf_set_ui (*result, 1);
 347     }
 348   else
 349     {
 350       mpf_init (q);
 351       mpf_init (r);
 352       mpf_init (factor);
 353       mpf_init (term);
 354
 355       mpf_div (q, x, two_pi);
 356       mpf_floor (factor, q);
 357       mpf_mul (q, factor, two_pi);
 358       mpf_sub (r, x, q);
 359
 360       mpf_init_set_ui (xp, 1);
 361       mpf_init_set_ui (num, 1);
 362       mpf_init_set_ui (denom, 1);
 363
 364       sign = 1;
 365       for (i = 1; i < GFC_REAL_BITS + 10; i++)
 366         {
 367           mpf_mul (num, num, r);
 368           mpf_mul_ui (denom, denom, i);
 369           if (i % 2 != 0)
 370             continue;
 371
 372           sign = -sign;
 373           mpf_div (term, num, denom);
 374           if (sign > 0)
 375             mpf_add (xp, xp, term);
 376           else
 377             mpf_sub (xp, xp, term);
 378         }
 379       mpf_set (*result, xp);
 380
 381       mpf_clear (q);
 382       mpf_clear (r);
 383       mpf_clear (factor);
 384       mpf_clear (num);
 385       mpf_clear (denom);
 386       mpf_clear (term);
 387       mpf_clear (xp);
 388     }
 389
 390   mpf_clear (x);
 391 }
 392
 393
 394 /* Calculate atan(arg).
 395
 396    Similar to sine but requires special handling for x near 1.  */
 397
 398 void
 399 arctangent (mpf_t * arg, mpf_t * result)
 400 {
 401   mpf_t absval, convgu, convgl, num, term, x, xp;
 402   int i, sign;
 403
 404   mpf_init_set (x, *arg);
 405
 406   /* Special cases */
 407   if (mpf_cmp_ui (x, 0) == 0)
 408     {
 409       mpf_set_ui (*result, 0);
 410     }
 411   else if (mpf_cmp_ui (x, 1) == 0)
 412     {
 413       mpf_init (num);
 414       mpf_div_ui (num, half_pi, 2);
 415       mpf_set (*result, num);
 416       mpf_clear (num);
 417     }
 418   else if (mpf_cmp_si (x, -1) == 0)
 419     {
 420       mpf_init (num);
 421       mpf_div_ui (num, half_pi, 2);
 422       mpf_neg (*result, num);
 423       mpf_clear (num);
 424     }
 425   else
 426     {                           /* General cases */
 427
 428       mpf_init (absval);
 429       mpf_abs (absval, x);
 430
 431       mpf_init_set_d (convgu, 1.5);
 432       mpf_init_set_d (convgl, 0.5);
 433       mpf_init_set_ui (num, 1);
 434       mpf_init (term);
 435
 436       if (mpf_cmp (absval, convgl) < 0)
 437         {
 438           mpf_init_set_ui (xp, 0);
 439           sign = -1;
 440           for (i = 1; i < GFC_REAL_BITS + 10; i++)
 441             {
 442               mpf_mul (num, num, absval);
 443               if (i % 2 == 0)
 444                 continue;
 445
 446               sign = -sign;
 447               mpf_div_ui (term, num, i);
 448               if (sign > 0)
 449                 mpf_add (xp, xp, term);
 450               else
 451                 mpf_sub (xp, xp, term);
 452             }
 453         }
 454       else if (mpf_cmp (absval, convgu) >= 0)
 455         {
 456           mpf_init_set (xp, half_pi);
 457           sign = 1;
 458           for (i = 1; i < GFC_REAL_BITS + 10; i++)
 459             {
 460               mpf_div (num, num, absval);
 461               if (i % 2 == 0)
 462                 continue;
 463
 464               sign = -sign;
 465               mpf_div_ui (term, num, i);
 466               if (sign > 0)
 467                 mpf_add (xp, xp, term);
 468               else
 469                 mpf_sub (xp, xp, term);
 470             }
 471         }
 472       else
 473         {
 474           mpf_init_set_ui (xp, 0);
 475
 476           mpf_sub_ui (num, absval, 1);
 477           mpf_add_ui (term, absval, 1);
 478           mpf_div (absval, num, term);
 479
 480           mpf_set_ui (num, 1);
 481
 482           sign = -1;
 483           for (i = 1; i < GFC_REAL_BITS + 10; i++)
 484             {
 485               mpf_mul (num, num, absval);
 486               if (i % 2 == 0)
 487                 continue;
 488               sign = -sign;
 489               mpf_div_ui (term, num, i);
 490               if (sign > 0)
 491                 mpf_add (xp, xp, term);
 492               else
 493                 mpf_sub (xp, xp, term);
 494             }
 495
 496           mpf_div_ui (term, half_pi, 2);
 497           mpf_add (xp, term, xp);
 498         }
 499
 500       /* This makes sure to preserve the identity arctan(-x) = -arctan(x)
 501          and improves accuracy to boot.  */
 502
 503       if (mpf_cmp_ui (x, 0) > 0)
 504         mpf_set (*result, xp);
 505       else
 506         mpf_neg (*result, xp);
 507
 508       mpf_clear (absval);
 509       mpf_clear (convgl);
 510       mpf_clear (convgu);
 511       mpf_clear (num);
 512       mpf_clear (term);
 513       mpf_clear (xp);
 514     }
 515   mpf_clear (x);
 516 }
 517
 518
 519 /* Calculate atan2 (y, x)
 520
 521 atan2(y, x) = atan(y/x)                         if x > 0,
 522               sign(y)*(pi - atan(|y/x|))        if x < 0,
 523               0                                 if x = 0 && y == 0,
 524               sign(y)*pi/2                      if x = 0 && y != 0.
 525 */
 526
 527 void
 528 arctangent2 (mpf_t * y, mpf_t * x, mpf_t * result)
 529 {
 530   mpf_t t;
 531
 532   mpf_init (t);
 533
 534   switch (mpf_sgn (*x))
 535     {
 536     case 1:
 537       mpf_div (t, *y, *x);
 538       arctangent (&t, result);
 539       break;
 540     case -1:
 541       mpf_div (t, *y, *x);
 542       mpf_abs (t, t);
 543       arctangent (&t, &t);
 544       mpf_sub (*result, pi, t);
 545       if (mpf_sgn (*y) == -1)
 546         mpf_neg (*result, *result);
 547       break;
 548     case 0:
 549       if (mpf_sgn (*y) == 0)
 550         mpf_set_ui (*result, 0);
 551       else
 552         {
 553           mpf_set (*result, half_pi);
 554           if (mpf_sgn (*y) == -1)
 555             mpf_neg (*result, *result);
 556         }
 557        break;
 558     }
 559   mpf_clear (t);
 560 }
 561
 562 /* Calculate cosh(arg). */
 563
 564 void
 565 hypercos (mpf_t * arg, mpf_t * result)
 566 {
 567   mpf_t neg, term1, term2, x, xp;
 568
 569   mpf_init_set (x, *arg);
 570
 571   mpf_init (neg);
 572   mpf_init (term1);
 573   mpf_init (term2);
 574   mpf_init (xp);
 575
 576   mpf_neg (neg, x);
 577
 578   exponential (&x, &term1);
 579   exponential (&neg, &term2);
 580
 581   mpf_add (xp, term1, term2);
 582   mpf_div_ui (*result, xp, 2);
 583
 584   mpf_clear (neg);
 585   mpf_clear (term1);
 586   mpf_clear (term2);
 587   mpf_clear (x);
 588   mpf_clear (xp);
 589 }
 590
 591
 592 /* Calculate sinh(arg). */
 593
 594 void
 595 hypersine (mpf_t * arg, mpf_t * result)
 596 {
 597   mpf_t neg, term1, term2, x, xp;
 598
 599   mpf_init_set (x, *arg);
 600
 601   mpf_init (neg);
 602   mpf_init (term1);
 603   mpf_init (term2);
 604   mpf_init (xp);
 605
 606   mpf_neg (neg, x);
 607
 608   exponential (&x, &term1);
 609   exponential (&neg, &term2);
 610
 611   mpf_sub (xp, term1, term2);
 612   mpf_div_ui (*result, xp, 2);
 613
 614   mpf_clear (neg);
 615   mpf_clear (term1);
 616   mpf_clear (term2);
 617   mpf_clear (x);
 618   mpf_clear (xp);
 619 }
 620
 621
 622 /* Given an arithmetic error code, return a pointer to a string that
 623    explains the error.  */
 624
 625 static const char *
 626 gfc_arith_error (arith code)
 627 {
 628   const char *p;
 629
 630   switch (code)
 631     {
 632     case ARITH_OK:
 633       p = "Arithmetic OK";
 634       break;
 635     case ARITH_OVERFLOW:
 636       p = "Arithmetic overflow";
 637       break;
 638     case ARITH_UNDERFLOW:
 639       p = "Arithmetic underflow";
 640       break;
 641     case ARITH_DIV0:
 642       p = "Division by zero";
 643       break;
 644     case ARITH_0TO0:
 645       p = "Indeterminate form 0 ** 0";
 646       break;
 647     case ARITH_INCOMMENSURATE:
 648       p = "Array operands are incommensurate";
 649       break;
 650     default:
 651       gfc_internal_error ("gfc_arith_error(): Bad error code");
 652     }
 653
 654   return p;
 655 }
 656
 657
 658 /* Get things ready to do math.  */
 659
 660 void
 661 gfc_arith_init_1 (void)
 662 {
 663   gfc_integer_info *int_info;
 664   gfc_real_info *real_info;
 665   mpf_t a, b;
 666   mpz_t r;
 667   int i, n, limit;
 668
 669   /* Set the default precision for GMP computations.  */
 670   mpf_set_default_prec (GFC_REAL_BITS + 30);
 671
 672   /* Calculate e, needed by the natural_logarithm() subroutine.  */
 673   mpf_init (b);
 674   mpf_init_set_ui (e, 0);
 675   mpf_init_set_ui (a, 1);
 676
 677   for (i = 1; i < 100; i++)
 678     {
 679       mpf_add (e, e, a);
 680       mpf_div_ui (a, a, i);     /* 1/(i!) */
 681     }
 682
 683   /* Calculate pi, 2pi, pi/2, and -pi/2, needed for trigonometric
 684      functions.
 685
 686      We use the Bailey, Borwein and Plouffe formula:
 687
 688        pi = \sum{n=0}^\infty (1/16)^n [4/(8n+1) - 2/(8n+4) - 1/(8n+5) - 1/(8n+6)]
 689
 690      which gives about four bits per iteration.  */
 691
 692   mpf_init_set_ui (pi, 0);
 693
 694   mpf_init (two_pi);
 695   mpf_init (half_pi);
 696
 697   limit = (GFC_REAL_BITS / 4) + 10;     /* (1/16)^n gives 4 bits per iteration */
 698
 699   for (n = 0; n < limit; n++)
 700     {
 701       mpf_set_ui (b, 4);
 702       mpf_div_ui (b, b, 8 * n + 1);     /* 4/(8n+1) */
 703
 704       mpf_set_ui (a, 2);
 705       mpf_div_ui (a, a, 8 * n + 4);     /* 2/(8n+4) */
 706       mpf_sub (b, b, a);
 707
 708       mpf_set_ui (a, 1);
 709       mpf_div_ui (a, a, 8 * n + 5);     /* 1/(8n+5) */
 710       mpf_sub (b, b, a);
 711
 712       mpf_set_ui (a, 1);
 713       mpf_div_ui (a, a, 8 * n + 6);     /* 1/(8n+6) */
 714       mpf_sub (b, b, a);
 715
 716       mpf_set_ui (a, 16);
 717       mpf_pow_ui (a, a, n);     /* 16^n */
 718
 719       mpf_div (b, b, a);
 720
 721       mpf_add (pi, pi, b);
 722     }
 723
 724   mpf_mul_ui (two_pi, pi, 2);
 725   mpf_div_ui (half_pi, pi, 2);
 726
 727   /* Convert the minimum/maximum values for each kind into their
 728      GNU MP representation.  */
 729   mpz_init (r);
 730
 731   for (int_info = gfc_integer_kinds; int_info->kind != 0; int_info++)
 732     {
 733       /* Huge */
 734       mpz_set_ui (r, int_info->radix);
 735       mpz_pow_ui (r, r, int_info->digits);
 736
 737       mpz_init (int_info->huge);
 738       mpz_sub_ui (int_info->huge, r, 1);
 739
 740       /* These are the numbers that are actually representable by the
 741          target.  For bases other than two, this needs to be changed.  */
 742       if (int_info->radix != 2)
 743         gfc_internal_error ("Fix min_int, max_int calculation");
 744
 745       mpz_init (int_info->min_int);
 746       mpz_neg (int_info->min_int, int_info->huge);
 747       /* No -1 here, because the representation is symmetric.  */
 748
 749       mpz_init (int_info->max_int);
 750       mpz_add (int_info->max_int, int_info->huge, int_info->huge);
 751       mpz_add_ui (int_info->max_int, int_info->max_int, 1);
 752
 753       /* Range */
 754       mpf_set_z (a, int_info->huge);
 755       common_logarithm (&a, &a);
 756       mpf_trunc (a, a);
 757       mpz_set_f (r, a);
 758       int_info->range = mpz_get_si (r);
 759     }
 760
 761   /*  mpf_set_default_prec(GFC_REAL_BITS); */
 762   for (real_info = gfc_real_kinds; real_info->kind != 0; real_info++)
 763     {
 764       /* Huge */
 765       mpf_set_ui (a, real_info->radix);
 766       mpf_set_ui (b, real_info->radix);
 767
 768       mpf_pow_ui (a, a, real_info->max_exponent);
 769       mpf_pow_ui (b, b, real_info->max_exponent - real_info->digits);
 770
 771       mpf_init (real_info->huge);
 772       mpf_sub (real_info->huge, a, b);
 773
 774       /* Tiny */
 775       mpf_set_ui (b, real_info->radix);
 776       mpf_pow_ui (b, b, 1 - real_info->min_exponent);
 777
 778       mpf_init (real_info->tiny);
 779       mpf_ui_div (real_info->tiny, 1, b);
 780
 781       /* Epsilon */
 782       mpf_set_ui (b, real_info->radix);
 783       mpf_pow_ui (b, b, real_info->digits - 1);
 784
 785       mpf_init (real_info->epsilon);
 786       mpf_ui_div (real_info->epsilon, 1, b);
 787
 788       /* Range */
 789       common_logarithm (&real_info->huge, &a);
 790       common_logarithm (&real_info->tiny, &b);
 791       mpf_neg (b, b);
 792
 793       if (mpf_cmp (a, b) > 0)
 794         mpf_set (a, b);         /* a = min(a, b) */
 795
 796       mpf_trunc (a, a);
 797       mpz_set_f (r, a);
 798       real_info->range = mpz_get_si (r);
 799
 800       /* Precision */
 801       mpf_set_ui (a, real_info->radix);
 802       common_logarithm (&a, &a);
 803
 804       mpf_mul_ui (a, a, real_info->digits - 1);
 805       mpf_trunc (a, a);
 806       mpz_set_f (r, a);
 807       real_info->precision = mpz_get_si (r);
 808
 809       /* If the radix is an integral power of 10, add one to the
 810          precision.  */
 811       for (i = 10; i <= real_info->radix; i *= 10)
 812         if (i == real_info->radix)
 813           real_info->precision++;
 814     }
 815
 816   mpz_clear (r);
 817   mpf_clear (a);
 818   mpf_clear (b);
 819 }
 820
 821
 822 /* Clean up, get rid of numeric constants.  */
 823
 824 void
 825 gfc_arith_done_1 (void)
 826 {
 827   gfc_integer_info *ip;
 828   gfc_real_info *rp;
 829
 830   mpf_clear (e);
 831
 832   mpf_clear (pi);
 833   mpf_clear (half_pi);
 834   mpf_clear (two_pi);
 835
 836   for (ip = gfc_integer_kinds; ip->kind; ip++)
 837     {
 838       mpz_clear (ip->min_int);
 839       mpz_clear (ip->max_int);
 840       mpz_clear (ip->huge);
 841     }
 842
 843   for (rp = gfc_real_kinds; rp->kind; rp++)
 844     {
 845       mpf_clear (rp->epsilon);
 846       mpf_clear (rp->huge);
 847       mpf_clear (rp->tiny);
 848     }
 849 }
 850
 851
 852 /* Return default kinds.  */
 853
 854 int
 855 gfc_default_integer_kind (void)
 856 {
 857   return gfc_integer_kinds[gfc_option.i8 ? 1 : 0].kind;
 858 }
 859
 860 int
 861 gfc_default_real_kind (void)
 862 {
 863   return gfc_real_kinds[gfc_option.r8 ? 1 : 0].kind;
 864 }
 865
 866 int
 867 gfc_default_double_kind (void)
 868 {
 869   return gfc_real_kinds[1].kind;
 870 }
 871
 872 int
 873 gfc_default_character_kind (void)
 874 {
 875   return 1;
 876 }
 877
 878 int
 879 gfc_default_logical_kind (void)
 880 {
 881   return gfc_logical_kinds[gfc_option.i8 ? 1 : 0].kind;
 882 }
 883
 884 int
 885 gfc_default_complex_kind (void)
 886 {
 887   return gfc_default_real_kind ();
 888 }
 889
 890
 891 /* Make sure that a valid kind is present.  Returns an index into the
 892    gfc_integer_kinds array, -1 if the kind is not present.  */
 893
 894 static int
 895 validate_integer (int kind)
 896 {
 897   int i;
 898
 899   for (i = 0;; i++)
 900     {
 901       if (gfc_integer_kinds[i].kind == 0)
 902         {
 903           i = -1;
 904           break;
 905         }
 906       if (gfc_integer_kinds[i].kind == kind)
 907         break;
 908     }
 909
 910   return i;
 911 }
 912
 913
 914 static int
 915 validate_real (int kind)
 916 {
 917   int i;
 918
 919   for (i = 0;; i++)
 920     {
 921       if (gfc_real_kinds[i].kind == 0)
 922         {
 923           i = -1;
 924           break;
 925         }
 926       if (gfc_real_kinds[i].kind == kind)
 927         break;
 928     }
 929
 930   return i;
 931 }
 932
 933
 934 static int
 935 validate_logical (int kind)
 936 {
 937   int i;
 938
 939   for (i = 0;; i++)
 940     {
 941       if (gfc_logical_kinds[i].kind == 0)
 942         {
 943           i = -1;
 944           break;
 945         }
 946       if (gfc_logical_kinds[i].kind == kind)
 947         break;
 948     }
 949
 950   return i;
 951 }
 952
 953
 954 static int
 955 validate_character (int kind)
 956 {
 957
 958   if (kind == gfc_default_character_kind ())
 959     return 0;
 960   return -1;
 961 }
 962
 963
 964 /* Validate a kind given a basic type.  The return value is the same
 965    for the child functions, with -1 indicating nonexistence of the
 966    type.  */
 967
 968 int
 969 gfc_validate_kind (bt type, int kind)
 970 {
 971   int rc;
 972
 973   switch (type)
 974     {
 975     case BT_REAL:               /* Fall through */
 976     case BT_COMPLEX:
 977       rc = validate_real (kind);
 978       break;
 979     case BT_INTEGER:
 980       rc = validate_integer (kind);
 981       break;
 982     case BT_LOGICAL:
 983       rc = validate_logical (kind);
 984       break;
 985     case BT_CHARACTER:
 986       rc = validate_character (kind);
 987       break;
 988
 989     default:
 990       gfc_internal_error ("gfc_validate_kind(): Got bad type");
 991     }
 992
 993   return rc;
 994 }
 995
 996
 997 /* Given an integer and a kind, make sure that the integer lies within
 998    the range of the kind.  Returns ARITH_OK or ARITH_OVERFLOW.  */
 999
1000 static arith
1001 gfc_check_integer_range (mpz_t p, int kind)
1002 {
1003   arith result;
1004   int i;
1005
1006   i = validate_integer (kind);
1007   if (i == -1)
1008     gfc_internal_error ("gfc_check_integer_range(): Bad kind");
1009
1010   result = ARITH_OK;
1011
1012   if (mpz_cmp (p, gfc_integer_kinds[i].min_int) < 0
1013       || mpz_cmp (p, gfc_integer_kinds[i].max_int) > 0)
1014     result = ARITH_OVERFLOW;
1015
1016   return result;
1017 }
1018
1019
1020 /* Given a real and a kind, make sure that the real lies within the
1021    range of the kind.  Returns ARITH_OK, ARITH_OVERFLOW or
1022    ARITH_UNDERFLOW.  */
1023
1024 static arith
1025 gfc_check_real_range (mpf_t p, int kind)
1026 {
1027   arith retval;
1028   mpf_t q;
1029   int i;
1030
1031   mpf_init (q);
1032   mpf_abs (q, p);
1033
1034   i = validate_real (kind);
1035   if (i == -1)
1036     gfc_internal_error ("gfc_check_real_range(): Bad kind");
1037
1038   retval = ARITH_OK;
1039   if (mpf_sgn (q) == 0)
1040     goto done;
1041
1042   if (mpf_cmp (q, gfc_real_kinds[i].huge) == 1)
1043     {
1044       retval = ARITH_OVERFLOW;
1045       goto done;
1046     }
1047
1048   if (mpf_cmp (q, gfc_real_kinds[i].tiny) == -1)
1049     retval = ARITH_UNDERFLOW;
1050
1051 done:
1052   mpf_clear (q);
1053
1054   return retval;
1055 }
1056
1057
1058 /* Function to return a constant expression node of a given type and
1059    kind.  */
1060
1061 gfc_expr *
1062 gfc_constant_result (bt type, int kind, locus * where)
1063 {
1064   gfc_expr *result;
1065
1066   if (!where)
1067     gfc_internal_error
1068       ("gfc_constant_result(): locus 'where' cannot be NULL");
1069
1070   result = gfc_get_expr ();
1071
1072   result->expr_type = EXPR_CONSTANT;
1073   result->ts.type = type;
1074   result->ts.kind = kind;
1075   result->where = *where;
1076
1077   switch (type)
1078     {
1079     case BT_INTEGER:
1080       mpz_init (result->value.integer);
1081       break;
1082
1083     case BT_REAL:
1084       mpf_init (result->value.real);
1085       break;
1086
1087     case BT_COMPLEX:
1088       mpf_init (result->value.complex.r);
1089       mpf_init (result->value.complex.i);
1090       break;
1091
1092     default:
1093       break;
1094     }
1095
1096   return result;
1097 }
1098
1099
1100 /* Low-level arithmetic functions.  All of these subroutines assume
1101    that all operands are of the same type and return an operand of the
1102    same type.  The other thing about these subroutines is that they
1103    can fail in various ways -- overflow, underflow, division by zero,
1104    zero raised to the zero, etc.  */
1105
1106 static arith
1107 gfc_arith_not (gfc_expr * op1, gfc_expr ** resultp)
1108 {
1109   gfc_expr *result;
1110
1111   result = gfc_constant_result (BT_LOGICAL, op1->ts.kind, &op1->where);
1112   result->value.logical = !op1->value.logical;
1113   *resultp = result;
1114
1115   return ARITH_OK;
1116 }
1117
1118
1119 static arith
1120 gfc_arith_and (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
1121 {
1122   gfc_expr *result;
1123
1124   result = gfc_constant_result (BT_LOGICAL, gfc_kind_max (op1, op2),
1125                                 &op1->where);
1126   result->value.logical = op1->value.logical && op2->value.logical;
1127   *resultp = result;
1128
1129   return ARITH_OK;
1130 }
1131
1132
1133 static arith
1134 gfc_arith_or (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
1135 {
1136   gfc_expr *result;
1137
1138   result = gfc_constant_result (BT_LOGICAL, gfc_kind_max (op1, op2),
1139                                 &op1->where);
1140   result->value.logical = op1->value.logical || op2->value.logical;
1141   *resultp = result;
1142
1143   return ARITH_OK;
1144 }
1145
1146
1147 static arith
1148 gfc_arith_eqv (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
1149 {
1150   gfc_expr *result;
1151
1152   result = gfc_constant_result (BT_LOGICAL, gfc_kind_max (op1, op2),
1153                                 &op1->where);
1154   result->value.logical = op1->value.logical == op2->value.logical;
1155   *resultp = result;
1156
1157   return ARITH_OK;
1158 }
1159
1160
1161 static arith
1162 gfc_arith_neqv (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
1163 {
1164   gfc_expr *result;
1165
1166   result = gfc_constant_result (BT_LOGICAL, gfc_kind_max (op1, op2),
1167                                 &op1->where);
1168   result->value.logical = op1->value.logical != op2->value.logical;
1169   *resultp = result;
1170
1171   return ARITH_OK;
1172 }
1173
1174
1175 /* Make sure a constant numeric expression is within the range for
1176    it's type and kind.  Note that there's also a gfc_check_range(),
1177    but that one deals with the intrinsic RANGE function.  */
1178
1179 arith
1180 gfc_range_check (gfc_expr * e)
1181 {
1182   arith rc;
1183
1184   switch (e->ts.type)
1185     {
1186     case BT_INTEGER:
1187       rc = gfc_check_integer_range (e->value.integer, e->ts.kind);
1188       break;
1189
1190     case BT_REAL:
1191       rc = gfc_check_real_range (e->value.real, e->ts.kind);
1192       break;
1193
1194     case BT_COMPLEX:
1195       rc = gfc_check_real_range (e->value.complex.r, e->ts.kind);
1196       if (rc != ARITH_OK)
1197         rc = gfc_check_real_range (e->value.complex.i, e->ts.kind);
1198
1199       break;
1200
1201     default:
1202       gfc_internal_error ("gfc_range_check(): Bad type");
1203     }
1204
1205   return rc;
1206 }
1207
1208
1209 /* It may seem silly to have a subroutine that actually computes the
1210    unary plus of a constant, but it prevents us from making exceptions
1211    in the code elsewhere.  */
1212
1213 static arith
1214 gfc_arith_uplus (gfc_expr * op1, gfc_expr ** resultp)
1215 {
1216
1217   *resultp = gfc_copy_expr (op1);
1218   return ARITH_OK;
1219 }
1220
1221
1222 static arith
1223 gfc_arith_uminus (gfc_expr * op1, gfc_expr ** resultp)
1224 {
1225   gfc_expr *result;
1226   arith rc;
1227
1228   result = gfc_constant_result (op1->ts.type, op1->ts.kind, &op1->where);
1229
1230   switch (op1->ts.type)
1231     {
1232     case BT_INTEGER:
1233       mpz_neg (result->value.integer, op1->value.integer);
1234       break;
1235
1236     case BT_REAL:
1237       mpf_neg (result->value.real, op1->value.real);
1238       break;
1239
1240     case BT_COMPLEX:
1241       mpf_neg (result->value.complex.r, op1->value.complex.r);
1242       mpf_neg (result->value.complex.i, op1->value.complex.i);
1243       break;
1244
1245     default:
1246       gfc_internal_error ("gfc_arith_uminus(): Bad basic type");
1247     }
1248
1249   rc = gfc_range_check (result);
1250
1251   if (rc != ARITH_OK)
1252     gfc_free_expr (result);
1253   else
1254     *resultp = result;
1255
1256   return rc;
1257 }
1258
1259
1260 static arith
1261 gfc_arith_plus (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
1262 {
1263   gfc_expr *result;
1264   arith rc;
1265
1266   result = gfc_constant_result (op1->ts.type, op1->ts.kind, &op1->where);
1267
1268   switch (op1->ts.type)
1269     {
1270     case BT_INTEGER:
1271       mpz_add (result->value.integer, op1->value.integer, op2->value.integer);
1272       break;
1273
1274     case BT_REAL:
1275       mpf_add (result->value.real, op1->value.real, op2->value.real);
1276       break;
1277
1278     case BT_COMPLEX:
1279       mpf_add (result->value.complex.r, op1->value.complex.r,
1280                op2->value.complex.r);
1281
1282       mpf_add (result->value.complex.i, op1->value.complex.i,
1283                op2->value.complex.i);
1284       break;
1285
1286     default:
1287       gfc_internal_error ("gfc_arith_plus(): Bad basic type");
1288     }
1289
1290   rc = gfc_range_check (result);
1291
1292   if (rc != ARITH_OK)
1293     gfc_free_expr (result);
1294   else
1295     *resultp = result;
1296
1297   return rc;
1298 }
1299
1300
1301 static arith
1302 gfc_arith_minus (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
1303 {
1304   gfc_expr *result;
1305   arith rc;
1306
1307   result = gfc_constant_result (op1->ts.type, op1->ts.kind, &op1->where);
1308
1309   switch (op1->ts.type)
1310     {
1311     case BT_INTEGER:
1312       mpz_sub (result->value.integer, op1->value.integer, op2->value.integer);
1313       break;
1314
1315     case BT_REAL:
1316       mpf_sub (result->value.real, op1->value.real, op2->value.real);
1317       break;
1318
1319     case BT_COMPLEX:
1320       mpf_sub (result->value.complex.r, op1->value.complex.r,
1321                op2->value.complex.r);
1322
1323       mpf_sub (result->value.complex.i, op1->value.complex.i,
1324                op2->value.complex.i);
1325
1326       break;
1327
1328     default:
1329       gfc_internal_error ("gfc_arith_minus(): Bad basic type");
1330     }
1331
1332   rc = gfc_range_check (result);
1333
1334   if (rc != ARITH_OK)
1335     gfc_free_expr (result);
1336   else
1337     *resultp = result;
1338
1339   return rc;
1340 }
1341
1342
1343 static arith
1344 gfc_arith_times (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
1345 {
1346   gfc_expr *result;
1347   mpf_t x, y;
1348   arith rc;
1349
1350   result = gfc_constant_result (op1->ts.type, op1->ts.kind, &op1->where);
1351
1352   switch (op1->ts.type)
1353     {
1354     case BT_INTEGER:
1355       mpz_mul (result->value.integer, op1->value.integer, op2->value.integer);
1356       break;
1357
1358     case BT_REAL:
1359       mpf_mul (result->value.real, op1->value.real, op2->value.real);
1360       break;
1361
1362     case BT_COMPLEX:
1363       mpf_init (x);
1364       mpf_init (y);
1365
1366       mpf_mul (x, op1->value.complex.r, op2->value.complex.r);
1367       mpf_mul (y, op1->value.complex.i, op2->value.complex.i);
1368       mpf_sub (result->value.complex.r, x, y);
1369
1370       mpf_mul (x, op1->value.complex.r, op2->value.complex.i);
1371       mpf_mul (y, op1->value.complex.i, op2->value.complex.r);
1372       mpf_add (result->value.complex.i, x, y);
1373
1374       mpf_clear (x);
1375       mpf_clear (y);
1376
1377       break;
1378
1379     default:
1380       gfc_internal_error ("gfc_arith_times(): Bad basic type");
1381     }
1382
1383   rc = gfc_range_check (result);
1384
1385   if (rc != ARITH_OK)
1386     gfc_free_expr (result);
1387   else
1388     *resultp = result;
1389
1390   return rc;
1391 }
1392
1393
1394 static arith
1395 gfc_arith_divide (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
1396 {
1397   gfc_expr *result;
1398   mpf_t x, y, div;
1399   arith rc;
1400
1401   rc = ARITH_OK;
1402
1403   result = gfc_constant_result (op1->ts.type, op1->ts.kind, &op1->where);
1404
1405   switch (op1->ts.type)
1406     {
1407     case BT_INTEGER:
1408       if (mpz_sgn (op2->value.integer) == 0)
1409         {
1410           rc = ARITH_DIV0;
1411           break;
1412         }
1413
1414       mpz_tdiv_q (result->value.integer, op1->value.integer,
1415                   op2->value.integer);
1416       break;
1417
1418     case BT_REAL:
1419       if (mpf_sgn (op2->value.real) == 0)
1420         {
1421           rc = ARITH_DIV0;
1422           break;
1423         }
1424
1425       mpf_div (result->value.real, op1->value.real, op2->value.real);
1426       break;
1427
1428     case BT_COMPLEX:
1429       if (mpf_sgn (op2->value.complex.r) == 0
1430           && mpf_sgn (op2->value.complex.i) == 0)
1431         {
1432           rc = ARITH_DIV0;
1433           break;
1434         }
1435
1436       mpf_init (x);
1437       mpf_init (y);
1438       mpf_init (div);
1439
1440       mpf_mul (x, op2->value.complex.r, op2->value.complex.r);
1441       mpf_mul (y, op2->value.complex.i, op2->value.complex.i);
1442       mpf_add (div, x, y);
1443
1444       mpf_mul (x, op1->value.complex.r, op2->value.complex.r);
1445       mpf_mul (y, op1->value.complex.i, op2->value.complex.i);
1446       mpf_add (result->value.complex.r, x, y);
1447       mpf_div (result->value.complex.r, result->value.complex.r, div);
1448
1449       mpf_mul (x, op1->value.complex.i, op2->value.complex.r);
1450       mpf_mul (y, op1->value.complex.r, op2->value.complex.i);
1451       mpf_sub (result->value.complex.i, x, y);
1452       mpf_div (result->value.complex.i, result->value.complex.i, div);
1453
1454       mpf_clear (x);
1455       mpf_clear (y);
1456       mpf_clear (div);
1457
1458       break;
1459
1460     default:
1461       gfc_internal_error ("gfc_arith_divide(): Bad basic type");
1462     }
1463
1464   if (rc == ARITH_OK)
1465     rc = gfc_range_check (result);
1466
1467   if (rc != ARITH_OK)
1468     gfc_free_expr (result);
1469   else
1470     *resultp = result;
1471
1472   return rc;
1473 }
1474
1475
1476 /* Compute the reciprocal of a complex number (guaranteed nonzero).  */
1477
1478 static void
1479 complex_reciprocal (gfc_expr * op)
1480 {
1481   mpf_t mod, a, result_r, result_i;
1482
1483   mpf_init (mod);
1484   mpf_init (a);
1485
1486   mpf_mul (mod, op->value.complex.r, op->value.complex.r);
1487   mpf_mul (a, op->value.complex.i, op->value.complex.i);
1488   mpf_add (mod, mod, a);
1489
1490   mpf_init (result_r);
1491   mpf_div (result_r, op->value.complex.r, mod);
1492
1493   mpf_init (result_i);
1494   mpf_neg (result_i, op->value.complex.i);
1495   mpf_div (result_i, result_i, mod);
1496
1497   mpf_set (op->value.complex.r, result_r);
1498   mpf_set (op->value.complex.i, result_i);
1499
1500   mpf_clear (result_r);
1501   mpf_clear (result_i);
1502
1503   mpf_clear (mod);
1504   mpf_clear (a);
1505 }
1506
1507
1508 /* Raise a complex number to positive power.  */
1509
1510 static void
1511 complex_pow_ui (gfc_expr * base, int power, gfc_expr * result)
1512 {
1513   mpf_t temp_r, temp_i, a;
1514
1515   mpf_set_ui (result->value.complex.r, 1);
1516   mpf_set_ui (result->value.complex.i, 0);
1517
1518   mpf_init (temp_r);
1519   mpf_init (temp_i);
1520   mpf_init (a);
1521
1522   for (; power > 0; power--)
1523     {
1524       mpf_mul (temp_r, base->value.complex.r, result->value.complex.r);
1525       mpf_mul (a, base->value.complex.i, result->value.complex.i);
1526       mpf_sub (temp_r, temp_r, a);
1527
1528       mpf_mul (temp_i, base->value.complex.r, result->value.complex.i);
1529       mpf_mul (a, base->value.complex.i, result->value.complex.r);
1530       mpf_add (temp_i, temp_i, a);
1531
1532       mpf_set (result->value.complex.r, temp_r);
1533       mpf_set (result->value.complex.i, temp_i);
1534     }
1535
1536   mpf_clear (temp_r);
1537   mpf_clear (temp_i);
1538   mpf_clear (a);
1539 }
1540
1541
1542 /* Raise a number to an integer power.  */
1543
1544 static arith
1545 gfc_arith_power (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
1546 {
1547   int power, apower;
1548   gfc_expr *result;
1549   mpz_t unity_z;
1550   mpf_t unity_f;
1551   arith rc;
1552
1553   rc = ARITH_OK;
1554
1555   if (gfc_extract_int (op2, &power) != NULL)
1556     gfc_internal_error ("gfc_arith_power(): Bad exponent");
1557
1558   result = gfc_constant_result (op1->ts.type, op1->ts.kind, &op1->where);
1559
1560   if (power == 0)
1561     {                           /* Handle something to the zeroth power */
1562       switch (op1->ts.type)
1563         {
1564         case BT_INTEGER:
1565           if (mpz_sgn (op1->value.integer) == 0)
1566             rc = ARITH_0TO0;
1567           else
1568             mpz_set_ui (result->value.integer, 1);
1569
1570           break;
1571
1572         case BT_REAL:
1573           if (mpf_sgn (op1->value.real) == 0)
1574             rc = ARITH_0TO0;
1575           else
1576             mpf_set_ui (result->value.real, 1);
1577
1578           break;
1579
1580         case BT_COMPLEX:
1581           if (mpf_sgn (op1->value.complex.r) == 0
1582               && mpf_sgn (op1->value.complex.i) == 0)
1583             rc = ARITH_0TO0;
1584           else
1585             {
1586               mpf_set_ui (result->value.complex.r, 1);
1587               mpf_set_ui (result->value.complex.r, 0);
1588             }
1589
1590           break;
1591
1592         default:
1593           gfc_internal_error ("gfc_arith_power(): Bad base");
1594         }
1595     }
1596
1597   if (power != 0)
1598     {
1599       apower = power;
1600       if (power < 0)
1601         apower = -power;
1602
1603       switch (op1->ts.type)
1604         {
1605         case BT_INTEGER:
1606           mpz_pow_ui (result->value.integer, op1->value.integer, apower);
1607
1608           if (power < 0)
1609             {
1610               mpz_init_set_ui (unity_z, 1);
1611               mpz_tdiv_q (result->value.integer, unity_z,
1612                           result->value.integer);
1613               mpz_clear (unity_z);
1614             }
1615
1616           break;
1617
1618         case BT_REAL:
1619           mpf_pow_ui (result->value.real, op1->value.real, apower);
1620
1621           if (power < 0)
1622             {
1623               mpf_init_set_ui (unity_f, 1);
1624               mpf_div (result->value.real, unity_f, result->value.real);
1625               mpf_clear (unity_f);
1626             }
1627
1628           break;
1629
1630         case BT_COMPLEX:
1631           complex_pow_ui (op1, apower, result);
1632           if (power < 0)
1633             complex_reciprocal (result);
1634
1635           break;
1636
1637         default:
1638           break;
1639         }
1640     }
1641
1642   if (rc == ARITH_OK)
1643     rc = gfc_range_check (result);
1644
1645   if (rc != ARITH_OK)
1646     gfc_free_expr (result);
1647   else
1648     *resultp = result;
1649
1650   return rc;
1651 }
1652
1653
1654 /* Concatenate two string constants.  */
1655
1656 static arith
1657 gfc_arith_concat (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
1658 {
1659   gfc_expr *result;
1660   int len;
1661
1662   result = gfc_constant_result (BT_CHARACTER, gfc_default_character_kind (),
1663                                 &op1->where);
1664
1665   len = op1->value.character.length + op2->value.character.length;
1666
1667   result->value.character.string = gfc_getmem (len + 1);
1668   result->value.character.length = len;
1669
1670   memcpy (result->value.character.string, op1->value.character.string,
1671           op1->value.character.length);
1672
1673   memcpy (result->value.character.string + op1->value.character.length,
1674           op2->value.character.string, op2->value.character.length);
1675
1676   result->value.character.string[len] = '\0';
1677
1678   *resultp = result;
1679
1680   return ARITH_OK;
1681 }
1682
1683
1684 /* Comparison operators.  Assumes that the two expression nodes
1685    contain two constants of the same type.  */
1686
1687 int
1688 gfc_compare_expr (gfc_expr * op1, gfc_expr * op2)
1689 {
1690   int rc;
1691
1692   switch (op1->ts.type)
1693     {
1694     case BT_INTEGER:
1695       rc = mpz_cmp (op1->value.integer, op2->value.integer);
1696       break;
1697
1698     case BT_REAL:
1699       rc = mpf_cmp (op1->value.real, op2->value.real);
1700       break;
1701
1702     case BT_CHARACTER:
1703       rc = gfc_compare_string (op1, op2, NULL);
1704       break;
1705
1706     case BT_LOGICAL:
1707       rc = ((!op1->value.logical && op2->value.logical)
1708             || (op1->value.logical && !op2->value.logical));
1709       break;
1710
1711     default:
1712       gfc_internal_error ("gfc_compare_expr(): Bad basic type");
1713     }
1714
1715   return rc;
1716 }
1717
1718
1719 /* Compare a pair of complex numbers.  Naturally, this is only for
1720    equality/nonequality.  */
1721
1722 static int
1723 compare_complex (gfc_expr * op1, gfc_expr * op2)
1724 {
1725
1726   return (mpf_cmp (op1->value.complex.r, op2->value.complex.r) == 0
1727           && mpf_cmp (op1->value.complex.i, op2->value.complex.i) == 0);
1728 }
1729
1730
1731 /* Given two constant strings and the inverse collating sequence,
1732    compare the strings.  We return -1 for a<b, 0 for a==b and 1 for
1733    a>b.  If the xcoll_table is NULL, we use the processor's default
1734    collating sequence.  */
1735
1736 int
1737 gfc_compare_string (gfc_expr * a, gfc_expr * b, const int *xcoll_table)
1738 {
1739   int len, alen, blen, i, ac, bc;
1740
1741   alen = a->value.character.length;
1742   blen = b->value.character.length;
1743
1744   len = (alen > blen) ? alen : blen;
1745
1746   for (i = 0; i < len; i++)
1747     {
1748       ac = (i < alen) ? a->value.character.string[i] : ' ';
1749       bc = (i < blen) ? b->value.character.string[i] : ' ';
1750
1751       if (xcoll_table != NULL)
1752         {
1753           ac = xcoll_table[ac];
1754           bc = xcoll_table[bc];
1755         }
1756
1757       if (ac < bc)
1758         return -1;
1759       if (ac > bc)
1760         return 1;
1761     }
1762
1763   /* Strings are equal */
1764
1765   return 0;
1766 }
1767
1768
1769 /* Specific comparison subroutines.  */
1770
1771 static arith
1772 gfc_arith_eq (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
1773 {
1774   gfc_expr *result;
1775
1776   result = gfc_constant_result (BT_LOGICAL, gfc_default_logical_kind (),
1777                                 &op1->where);
1778   result->value.logical = (op1->ts.type == BT_COMPLEX) ?
1779     compare_complex (op1, op2) : (gfc_compare_expr (op1, op2) == 0);
1780
1781   *resultp = result;
1782   return ARITH_OK;
1783 }
1784
1785
1786 static arith
1787 gfc_arith_ne (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
1788 {
1789   gfc_expr *result;
1790
1791   result = gfc_constant_result (BT_LOGICAL, gfc_default_logical_kind (),
1792                                 &op1->where);
1793   result->value.logical = (op1->ts.type == BT_COMPLEX) ?
1794     !compare_complex (op1, op2) : (gfc_compare_expr (op1, op2) != 0);
1795
1796   *resultp = result;
1797   return ARITH_OK;
1798 }
1799
1800
1801 static arith
1802 gfc_arith_gt (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
1803 {
1804   gfc_expr *result;
1805
1806   result = gfc_constant_result (BT_LOGICAL, gfc_default_logical_kind (),
1807                                 &op1->where);
1808   result->value.logical = (gfc_compare_expr (op1, op2) > 0);
1809   *resultp = result;
1810
1811   return ARITH_OK;
1812 }
1813
1814
1815 static arith
1816 gfc_arith_ge (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
1817 {
1818   gfc_expr *result;
1819
1820   result = gfc_constant_result (BT_LOGICAL, gfc_default_logical_kind (),
1821                                 &op1->where);
1822   result->value.logical = (gfc_compare_expr (op1, op2) >= 0);
1823   *resultp = result;
1824
1825   return ARITH_OK;
1826 }
1827
1828
1829 static arith
1830 gfc_arith_lt (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
1831 {
1832   gfc_expr *result;
1833
1834   result = gfc_constant_result (BT_LOGICAL, gfc_default_logical_kind (),
1835                                 &op1->where);
1836   result->value.logical = (gfc_compare_expr (op1, op2) < 0);
1837   *resultp = result;
1838
1839   return ARITH_OK;
1840 }
1841
1842
1843 static arith
1844 gfc_arith_le (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
1845 {
1846   gfc_expr *result;
1847
1848   result = gfc_constant_result (BT_LOGICAL, gfc_default_logical_kind (),
1849                                 &op1->where);
1850   result->value.logical = (gfc_compare_expr (op1, op2) <= 0);
1851   *resultp = result;
1852
1853   return ARITH_OK;
1854 }
1855
1856
1857 static arith
1858 reduce_unary (arith (*eval) (gfc_expr *, gfc_expr **), gfc_expr * op,
1859               gfc_expr ** result)
1860 {
1861   gfc_constructor *c, *head;
1862   gfc_expr *r;
1863   arith rc;
1864
1865   if (op->expr_type == EXPR_CONSTANT)
1866     return eval (op, result);
1867
1868   rc = ARITH_OK;
1869   head = gfc_copy_constructor (op->value.constructor);
1870
1871   for (c = head; c; c = c->next)
1872     {
1873       rc = eval (c->expr, &r);
1874       if (rc != ARITH_OK)
1875         break;
1876
1877       gfc_replace_expr (c->expr, r);
1878     }
1879
1880   if (rc != ARITH_OK)
1881     gfc_free_constructor (head);
1882   else
1883     {
1884       r = gfc_get_expr ();
1885       r->expr_type = EXPR_ARRAY;
1886       r->value.constructor = head;
1887       r->shape = gfc_copy_shape (op->shape, op->rank);
1888
1889       r->ts = head->expr->ts;
1890       r->where = op->where;
1891       r->rank = op->rank;
1892
1893       *result = r;
1894     }
1895
1896   return rc;
1897 }
1898
1899
1900 static arith
1901 reduce_binary_ac (arith (*eval) (gfc_expr *, gfc_expr *, gfc_expr **),
1902                   gfc_expr * op1, gfc_expr * op2,
1903                   gfc_expr ** result)
1904 {
1905   gfc_constructor *c, *head;
1906   gfc_expr *r;
1907   arith rc;
1908
1909   head = gfc_copy_constructor (op1->value.constructor);
1910   rc = ARITH_OK;
1911
1912   for (c = head; c; c = c->next)
1913     {
1914       rc = eval (c->expr, op2, &r);
1915       if (rc != ARITH_OK)
1916         break;
1917
1918       gfc_replace_expr (c->expr, r);
1919     }
1920
1921   if (rc != ARITH_OK)
1922     gfc_free_constructor (head);
1923   else
1924     {
1925       r = gfc_get_expr ();
1926       r->expr_type = EXPR_ARRAY;
1927       r->value.constructor = head;
1928       r->shape = gfc_copy_shape (op1->shape, op1->rank);
1929
1930       r->ts = head->expr->ts;
1931       r->where = op1->where;
1932       r->rank = op1->rank;
1933
1934       *result = r;
1935     }
1936
1937   return rc;
1938 }
1939
1940
1941 static arith
1942 reduce_binary_ca (arith (*eval) (gfc_expr *, gfc_expr *, gfc_expr **),
1943                   gfc_expr * op1, gfc_expr * op2,
1944                   gfc_expr ** result)
1945 {
1946   gfc_constructor *c, *head;
1947   gfc_expr *r;
1948   arith rc;
1949
1950   head = gfc_copy_constructor (op2->value.constructor);
1951   rc = ARITH_OK;
1952
1953   for (c = head; c; c = c->next)
1954     {
1955       rc = eval (op1, c->expr, &r);
1956       if (rc != ARITH_OK)
1957         break;
1958
1959       gfc_replace_expr (c->expr, r);
1960     }
1961
1962   if (rc != ARITH_OK)
1963     gfc_free_constructor (head);
1964   else
1965     {
1966       r = gfc_get_expr ();
1967       r->expr_type = EXPR_ARRAY;
1968       r->value.constructor = head;
1969       r->shape = gfc_copy_shape (op2->shape, op2->rank);
1970
1971       r->ts = head->expr->ts;
1972       r->where = op2->where;
1973       r->rank = op2->rank;
1974
1975       *result = r;
1976     }
1977
1978   return rc;
1979 }
1980
1981
1982 static arith
1983 reduce_binary_aa (arith (*eval) (gfc_expr *, gfc_expr *, gfc_expr **),
1984                   gfc_expr * op1, gfc_expr * op2,
1985                   gfc_expr ** result)
1986 {
1987   gfc_constructor *c, *d, *head;
1988   gfc_expr *r;
1989   arith rc;
1990
1991   head = gfc_copy_constructor (op1->value.constructor);
1992
1993   rc = ARITH_OK;
1994   d = op2->value.constructor;
1995
1996   if (gfc_check_conformance ("Elemental binary operation", op1, op2)
1997       != SUCCESS)
1998     rc = ARITH_INCOMMENSURATE;
1999   else
2000     {
2001
2002       for (c = head; c; c = c->next, d = d->next)
2003         {
2004           if (d == NULL)
2005             {
2006               rc = ARITH_INCOMMENSURATE;
2007               break;
2008             }
2009
2010           rc = eval (c->expr, d->expr, &r);
2011           if (rc != ARITH_OK)
2012             break;
2013
2014           gfc_replace_expr (c->expr, r);
2015         }
2016
2017       if (d != NULL)
2018         rc = ARITH_INCOMMENSURATE;
2019     }
2020
2021   if (rc != ARITH_OK)
2022     gfc_free_constructor (head);
2023   else
2024     {
2025       r = gfc_get_expr ();
2026       r->expr_type = EXPR_ARRAY;
2027       r->value.constructor = head;
2028       r->shape = gfc_copy_shape (op1->shape, op1->rank);
2029
2030       r->ts = head->expr->ts;
2031       r->where = op1->where;
2032       r->rank = op1->rank;
2033
2034       *result = r;
2035     }
2036
2037   return rc;
2038 }
2039
2040
2041 static arith
2042 reduce_binary (arith (*eval) (gfc_expr *, gfc_expr *, gfc_expr **),
2043                gfc_expr * op1, gfc_expr * op2,
2044                gfc_expr ** result)
2045 {
2046
2047   if (op1->expr_type == EXPR_CONSTANT && op2->expr_type == EXPR_CONSTANT)
2048     return eval (op1, op2, result);
2049
2050   if (op1->expr_type == EXPR_CONSTANT && op2->expr_type == EXPR_ARRAY)
2051     return reduce_binary_ca (eval, op1, op2, result);
2052
2053   if (op1->expr_type == EXPR_ARRAY && op2->expr_type == EXPR_CONSTANT)
2054     return reduce_binary_ac (eval, op1, op2, result);
2055
2056   return reduce_binary_aa (eval, op1, op2, result);
2057 }
2058
2059
2060 typedef union
2061 {
2062   arith (*f2)(gfc_expr *, gfc_expr **);
2063   arith (*f3)(gfc_expr *, gfc_expr *, gfc_expr **);
2064 }
2065 eval_f;
2066
2067 /* High level arithmetic subroutines.  These subroutines go into
2068    eval_intrinsic(), which can do one of several things to its
2069    operands.  If the operands are incompatible with the intrinsic
2070    operation, we return a node pointing to the operands and hope that
2071    an operator interface is found during resolution.
2072
2073    If the operands are compatible and are constants, then we try doing
2074    the arithmetic.  We also handle the cases where either or both
2075    operands are array constructors.  */
2076
2077 static gfc_expr *
2078 eval_intrinsic (gfc_intrinsic_op operator,
2079                 eval_f eval, gfc_expr * op1, gfc_expr * op2)
2080 {
2081   gfc_expr temp, *result;
2082   int unary;
2083   arith rc;
2084
2085   gfc_clear_ts (&temp.ts);
2086
2087   switch (operator)
2088     {
2089     case INTRINSIC_NOT: /* Logical unary */
2090       if (op1->ts.type != BT_LOGICAL)
2091         goto runtime;
2092
2093       temp.ts.type = BT_LOGICAL;
2094       temp.ts.kind = gfc_default_logical_kind ();
2095
2096       unary = 1;
2097       break;
2098
2099       /* Logical binary operators */
2100     case INTRINSIC_OR:
2101     case INTRINSIC_AND:
2102     case INTRINSIC_NEQV:
2103     case INTRINSIC_EQV:
2104       if (op1->ts.type != BT_LOGICAL || op2->ts.type != BT_LOGICAL)
2105         goto runtime;
2106
2107       temp.ts.type = BT_LOGICAL;
2108       temp.ts.kind = gfc_default_logical_kind ();
2109
2110       unary = 0;
2111       break;
2112
2113     case INTRINSIC_UPLUS:
2114     case INTRINSIC_UMINUS:      /* Numeric unary */
2115       if (!gfc_numeric_ts (&op1->ts))
2116         goto runtime;
2117
2118       temp.ts = op1->ts;
2119
2120       unary = 1;
2121       break;
2122
2123     case INTRINSIC_GE:
2124     case INTRINSIC_LT:          /* Additional restrictions  */
2125     case INTRINSIC_LE:          /* for ordering relations.  */
2126     case INTRINSIC_GT:
2127       if (op1->ts.type == BT_COMPLEX || op2->ts.type == BT_COMPLEX)
2128         {
2129           temp.ts.type = BT_LOGICAL;
2130           temp.ts.kind = gfc_default_logical_kind();
2131           goto runtime;
2132         }
2133
2134       /* else fall through */
2135
2136     case INTRINSIC_EQ:
2137     case INTRINSIC_NE:
2138       if (op1->ts.type == BT_CHARACTER && op2->ts.type == BT_CHARACTER)
2139         {
2140           unary = 0;
2141           temp.ts.type = BT_LOGICAL;
2142           temp.ts.kind = gfc_default_logical_kind();
2143           break;
2144         }
2145
2146       /* else fall through */
2147
2148     case INTRINSIC_PLUS:
2149     case INTRINSIC_MINUS:
2150     case INTRINSIC_TIMES:
2151     case INTRINSIC_DIVIDE:
2152     case INTRINSIC_POWER:       /* Numeric binary */
2153       if (!gfc_numeric_ts (&op1->ts) || !gfc_numeric_ts (&op2->ts))
2154         goto runtime;
2155
2156       /* Insert any necessary type conversions to make the operands compatible.  */
2157
2158       temp.expr_type = EXPR_OP;
2159       gfc_clear_ts (&temp.ts);
2160       temp.operator = operator;
2161
2162       temp.op1 = op1;
2163       temp.op2 = op2;
2164
2165       gfc_type_convert_binary (&temp);
2166
2167       if (operator == INTRINSIC_EQ || operator == INTRINSIC_NE
2168           || operator == INTRINSIC_GE || operator == INTRINSIC_GT
2169           || operator == INTRINSIC_LE || operator == INTRINSIC_LT)
2170         {
2171           temp.ts.type = BT_LOGICAL;
2172           temp.ts.kind = gfc_default_logical_kind ();
2173         }
2174
2175       unary = 0;
2176       break;
2177
2178     case INTRINSIC_CONCAT:      /* Character binary */
2179       if (op1->ts.type != BT_CHARACTER || op2->ts.type != BT_CHARACTER)
2180         goto runtime;
2181
2182       temp.ts.type = BT_CHARACTER;
2183       temp.ts.kind = gfc_default_character_kind ();
2184
2185       unary = 0;
2186       break;
2187
2188     case INTRINSIC_USER:
2189       goto runtime;
2190
2191     default:
2192       gfc_internal_error ("eval_intrinsic(): Bad operator");
2193     }
2194
2195   /* Try to combine the operators.  */
2196   if (operator == INTRINSIC_POWER && op2->ts.type != BT_INTEGER)
2197     goto runtime;
2198
2199   if (op1->expr_type != EXPR_CONSTANT
2200       && (op1->expr_type != EXPR_ARRAY
2201           || !gfc_is_constant_expr (op1)
2202           || !gfc_expanded_ac (op1)))
2203     goto runtime;
2204
2205   if (op2 != NULL
2206       && op2->expr_type != EXPR_CONSTANT
2207       && (op2->expr_type != EXPR_ARRAY
2208           || !gfc_is_constant_expr (op2)
2209           || !gfc_expanded_ac (op2)))
2210     goto runtime;
2211
2212   if (unary)
2213     rc = reduce_unary (eval.f2, op1, &result);
2214   else
2215     rc = reduce_binary (eval.f3, op1, op2, &result);
2216
2217   if (rc != ARITH_OK)
2218     {                           /* Something went wrong */
2219       gfc_error ("%s at %L", gfc_arith_error (rc), &op1->where);
2220       return NULL;
2221     }
2222
2223   gfc_free_expr (op1);
2224   gfc_free_expr (op2);
2225   return result;
2226
2227 runtime:
2228   /* Create a run-time expression */
2229   result = gfc_get_expr ();
2230   result->ts = temp.ts;
2231
2232   result->expr_type = EXPR_OP;
2233   result->operator = operator;
2234
2235   result->op1 = op1;
2236   result->op2 = op2;
2237
2238   result->where = op1->where;
2239
2240   return result;
2241 }
2242
2243
2244 /* Modify type of expression for zero size array.  */
2245 static gfc_expr *
2246 eval_type_intrinsic0 (gfc_intrinsic_op operator, gfc_expr *op)
2247 {
2248   if (op == NULL)
2249     gfc_internal_error("eval_type_intrinsic0(): op NULL");
2250
2251   switch(operator)
2252     {
2253     case INTRINSIC_GE:
2254     case INTRINSIC_LT:
2255     case INTRINSIC_LE:
2256     case INTRINSIC_GT:
2257     case INTRINSIC_EQ:
2258     case INTRINSIC_NE:
2259       op->ts.type = BT_LOGICAL;
2260       op->ts.kind = gfc_default_logical_kind();
2261       break;
2262
2263     default:
2264       break;
2265     }
2266
2267   return op;
2268 }
2269
2270
2271 /* Return nonzero if the expression is a zero size array.  */
2272
2273 static int
2274 gfc_zero_size_array (gfc_expr * e)
2275 {
2276
2277   if (e->expr_type != EXPR_ARRAY)
2278     return 0;
2279
2280   return e->value.constructor == NULL;
2281 }
2282
2283
2284 /* Reduce a binary expression where at least one of the operands
2285    involves a zero-length array.  Returns NULL if neither of the
2286    operands is a zero-length array.  */
2287
2288 static gfc_expr *
2289 reduce_binary0 (gfc_expr * op1, gfc_expr * op2)
2290 {
2291
2292   if (gfc_zero_size_array (op1))
2293     {
2294       gfc_free_expr (op2);
2295       return op1;
2296     }
2297
2298   if (gfc_zero_size_array (op2))
2299     {
2300       gfc_free_expr (op1);
2301       return op2;
2302     }
2303
2304   return NULL;
2305 }
2306
2307
2308 static gfc_expr *
2309 eval_intrinsic_f2 (gfc_intrinsic_op operator,
2310                    arith (*eval) (gfc_expr *, gfc_expr **),
2311                    gfc_expr * op1, gfc_expr * op2)
2312 {
2313   gfc_expr *result;
2314   eval_f f;
2315
2316   if (op2 == NULL)
2317     {
2318       if (gfc_zero_size_array (op1))
2319         return eval_type_intrinsic0(operator, op1);
2320     }
2321   else
2322     {
2323       result = reduce_binary0 (op1, op2);
2324       if (result != NULL)
2325         return eval_type_intrinsic0(operator, result);
2326     }
2327
2328   f.f2 = eval;
2329   return eval_intrinsic (operator, f, op1, op2);
2330 }
2331
2332
2333 static gfc_expr *
2334 eval_intrinsic_f3 (gfc_intrinsic_op operator,
2335                    arith (*eval) (gfc_expr *, gfc_expr *, gfc_expr **),
2336                    gfc_expr * op1, gfc_expr * op2)
2337 {
2338   gfc_expr *result;
2339   eval_f f;
2340
2341   result = reduce_binary0 (op1, op2);
2342   if (result != NULL)
2343     return eval_type_intrinsic0(operator, result);
2344
2345   f.f3 = eval;
2346   return eval_intrinsic (operator, f, op1, op2);
2347 }
2348
2349
2350
2351 gfc_expr *
2352 gfc_uplus (gfc_expr * op)
2353 {
2354   return eval_intrinsic_f2 (INTRINSIC_UPLUS, gfc_arith_uplus, op, NULL);
2355 }
2356
2357 gfc_expr *
2358 gfc_uminus (gfc_expr * op)
2359 {
2360   return eval_intrinsic_f2 (INTRINSIC_UMINUS, gfc_arith_uminus, op, NULL);
2361 }
2362
2363 gfc_expr *
2364 gfc_add (gfc_expr * op1, gfc_expr * op2)
2365 {
2366   return eval_intrinsic_f3 (INTRINSIC_PLUS, gfc_arith_plus, op1, op2);
2367 }
2368
2369 gfc_expr *
2370 gfc_subtract (gfc_expr * op1, gfc_expr * op2)
2371 {
2372   return eval_intrinsic_f3 (INTRINSIC_MINUS, gfc_arith_minus, op1, op2);
2373 }
2374
2375 gfc_expr *
2376 gfc_multiply (gfc_expr * op1, gfc_expr * op2)
2377 {
2378   return eval_intrinsic_f3 (INTRINSIC_TIMES, gfc_arith_times, op1, op2);
2379 }
2380
2381 gfc_expr *
2382 gfc_divide (gfc_expr * op1, gfc_expr * op2)
2383 {
2384   return eval_intrinsic_f3 (INTRINSIC_DIVIDE, gfc_arith_divide, op1, op2);
2385 }
2386
2387 gfc_expr *
2388 gfc_power (gfc_expr * op1, gfc_expr * op2)
2389 {
2390   return eval_intrinsic_f3 (INTRINSIC_POWER, gfc_arith_power, op1, op2);
2391 }
2392
2393 gfc_expr *
2394 gfc_concat (gfc_expr * op1, gfc_expr * op2)
2395 {
2396   return eval_intrinsic_f3 (INTRINSIC_CONCAT, gfc_arith_concat, op1, op2);
2397 }
2398
2399 gfc_expr *
2400 gfc_and (gfc_expr * op1, gfc_expr * op2)
2401 {
2402   return eval_intrinsic_f3 (INTRINSIC_AND, gfc_arith_and, op1, op2);
2403 }
2404
2405 gfc_expr *
2406 gfc_or (gfc_expr * op1, gfc_expr * op2)
2407 {
2408   return eval_intrinsic_f3 (INTRINSIC_OR, gfc_arith_or, op1, op2);
2409 }
2410
2411 gfc_expr *
2412 gfc_not (gfc_expr * op1)
2413 {
2414   return eval_intrinsic_f2 (INTRINSIC_NOT, gfc_arith_not, op1, NULL);
2415 }
2416
2417 gfc_expr *
2418 gfc_eqv (gfc_expr * op1, gfc_expr * op2)
2419 {
2420   return eval_intrinsic_f3 (INTRINSIC_EQV, gfc_arith_eqv, op1, op2);
2421 }
2422
2423 gfc_expr *
2424 gfc_neqv (gfc_expr * op1, gfc_expr * op2)
2425 {
2426   return eval_intrinsic_f3 (INTRINSIC_NEQV, gfc_arith_neqv, op1, op2);
2427 }
2428
2429 gfc_expr *
2430 gfc_eq (gfc_expr * op1, gfc_expr * op2)
2431 {
2432   return eval_intrinsic_f3 (INTRINSIC_EQ, gfc_arith_eq, op1, op2);
2433 }
2434
2435 gfc_expr *
2436 gfc_ne (gfc_expr * op1, gfc_expr * op2)
2437 {
2438   return eval_intrinsic_f3 (INTRINSIC_NE, gfc_arith_ne, op1, op2);
2439 }
2440
2441 gfc_expr *
2442 gfc_gt (gfc_expr * op1, gfc_expr * op2)
2443 {
2444   return eval_intrinsic_f3 (INTRINSIC_GT, gfc_arith_gt, op1, op2);
2445 }
2446
2447 gfc_expr *
2448 gfc_ge (gfc_expr * op1, gfc_expr * op2)
2449 {
2450   return eval_intrinsic_f3 (INTRINSIC_GE, gfc_arith_ge, op1, op2);
2451 }
2452
2453 gfc_expr *
2454 gfc_lt (gfc_expr * op1, gfc_expr * op2)
2455 {
2456   return eval_intrinsic_f3 (INTRINSIC_LT, gfc_arith_lt, op1, op2);
2457 }
2458
2459 gfc_expr *
2460 gfc_le (gfc_expr * op1, gfc_expr * op2)
2461 {
2462   return eval_intrinsic_f3 (INTRINSIC_LE, gfc_arith_le, op1, op2);
2463 }
2464
2465
2466 /* Convert an integer string to an expression node.  */
2467
2468 gfc_expr *
2469 gfc_convert_integer (const char *buffer, int kind, int radix, locus * where)
2470 {
2471   gfc_expr *e;
2472   const char *t;
2473
2474   e = gfc_constant_result (BT_INTEGER, kind, where);
2475   /* a leading plus is allowed, but not by mpz_set_str */
2476   if (buffer[0] == '+')
2477     t = buffer + 1;
2478   else
2479     t = buffer;
2480   mpz_set_str (e->value.integer, t, radix);
2481
2482   return e;
2483 }
2484
2485
2486 /* Convert a real string to an expression node.  */
2487
2488 gfc_expr *
2489 gfc_convert_real (const char *buffer, int kind, locus * where)
2490 {
2491   gfc_expr *e;
2492   const char *t;
2493
2494   e = gfc_constant_result (BT_REAL, kind, where);
2495   /* a leading plus is allowed, but not by mpf_set_str */
2496   if (buffer[0] == '+')
2497     t = buffer + 1;
2498   else
2499     t = buffer;
2500   mpf_set_str (e->value.real, t, 10);
2501
2502   return e;
2503 }
2504
2505
2506 /* Convert a pair of real, constant expression nodes to a single
2507    complex expression node.  */
2508
2509 gfc_expr *
2510 gfc_convert_complex (gfc_expr * real, gfc_expr * imag, int kind)
2511 {
2512   gfc_expr *e;
2513
2514   e = gfc_constant_result (BT_COMPLEX, kind, &real->where);
2515   mpf_set (e->value.complex.r, real->value.real);
2516   mpf_set (e->value.complex.i, imag->value.real);
2517
2518   return e;
2519 }
2520
2521
2522 /******* Simplification of intrinsic functions with constant arguments *****/
2523
2524
2525 /* Deal with an arithmetic error.  */
2526
2527 static void
2528 arith_error (arith rc, gfc_typespec * from, gfc_typespec * to, locus * where)
2529 {
2530
2531   gfc_error ("%s converting %s to %s at %L", gfc_arith_error (rc),
2532              gfc_typename (from), gfc_typename (to), where);
2533
2534   /* TODO: Do something about the error, ie underflow rounds to 0,
2535      throw exception, return NaN, etc.  */
2536 }
2537
2538 /* Convert integers to integers.  */
2539
2540 gfc_expr *
2541 gfc_int2int (gfc_expr * src, int kind)
2542 {
2543   gfc_expr *result;
2544   arith rc;
2545
2546   result = gfc_constant_result (BT_INTEGER, kind, &src->where);
2547
2548   mpz_set (result->value.integer, src->value.integer);
2549
2550   if ((rc = gfc_check_integer_range (result->value.integer, kind))
2551       != ARITH_OK)
2552     {
2553       arith_error (rc, &src->ts, &result->ts, &src->where);
2554       gfc_free_expr (result);
2555       return NULL;
2556     }
2557
2558   return result;
2559 }
2560
2561
2562 /* Convert integers to reals.  */
2563
2564 gfc_expr *
2565 gfc_int2real (gfc_expr * src, int kind)
2566 {
2567   gfc_expr *result;
2568   arith rc;
2569
2570   result = gfc_constant_result (BT_REAL, kind, &src->where);
2571
2572   mpf_set_z (result->value.real, src->value.integer);
2573
2574   if ((rc = gfc_check_real_range (result->value.real, kind)) != ARITH_OK)
2575     {
2576       arith_error (rc, &src->ts, &result->ts, &src->where);
2577       gfc_free_expr (result);
2578       return NULL;
2579     }
2580
2581   return result;
2582 }
2583
2584
2585 /* Convert default integer to default complex.  */
2586
2587 gfc_expr *
2588 gfc_int2complex (gfc_expr * src, int kind)
2589 {
2590   gfc_expr *result;
2591   arith rc;
2592
2593   result = gfc_constant_result (BT_COMPLEX, kind, &src->where);
2594
2595   mpf_set_z (result->value.complex.r, src->value.integer);
2596   mpf_set_ui (result->value.complex.i, 0);
2597
2598   if ((rc = gfc_check_real_range (result->value.complex.i, kind)) != ARITH_OK)
2599     {
2600       arith_error (rc, &src->ts, &result->ts, &src->where);
2601       gfc_free_expr (result);
2602       return NULL;
2603     }
2604
2605   return result;
2606 }
2607
2608
2609 /* Convert default real to default integer.  */
2610
2611 gfc_expr *
2612 gfc_real2int (gfc_expr * src, int kind)
2613 {
2614   gfc_expr *result;
2615   arith rc;
2616
2617   result = gfc_constant_result (BT_INTEGER, kind, &src->where);
2618
2619   mpz_set_f (result->value.integer, src->value.real);
2620
2621   if ((rc = gfc_check_integer_range (result->value.integer, kind))
2622       != ARITH_OK)
2623     {
2624       arith_error (rc, &src->ts, &result->ts, &src->where);
2625       gfc_free_expr (result);
2626       return NULL;
2627     }
2628
2629   return result;
2630 }
2631
2632
2633 /* Convert real to real.  */
2634
2635 gfc_expr *
2636 gfc_real2real (gfc_expr * src, int kind)
2637 {
2638   gfc_expr *result;
2639   arith rc;
2640
2641   result = gfc_constant_result (BT_REAL, kind, &src->where);
2642
2643   mpf_set (result->value.real, src->value.real);
2644
2645   if ((rc = gfc_check_real_range (result->value.real, kind)) != ARITH_OK)
2646     {
2647       arith_error (rc, &src->ts, &result->ts, &src->where);
2648       gfc_free_expr (result);
2649       return NULL;
2650     }
2651
2652   return result;
2653 }
2654
2655
2656 /* Convert real to complex.  */
2657
2658 gfc_expr *
2659 gfc_real2complex (gfc_expr * src, int kind)
2660 {
2661   gfc_expr *result;
2662   arith rc;
2663
2664   result = gfc_constant_result (BT_COMPLEX, kind, &src->where);
2665
2666   mpf_set (result->value.complex.r, src->value.real);
2667   mpf_set_ui (result->value.complex.i, 0);
2668
2669   if ((rc = gfc_check_real_range (result->value.complex.i, kind)) != ARITH_OK)
2670     {
2671       arith_error (rc, &src->ts, &result->ts, &src->where);
2672       gfc_free_expr (result);
2673       return NULL;
2674     }
2675
2676   return result;
2677 }
2678
2679
2680 /* Convert complex to integer.  */
2681
2682 gfc_expr *
2683 gfc_complex2int (gfc_expr * src, int kind)
2684 {
2685   gfc_expr *result;
2686   arith rc;
2687
2688   result = gfc_constant_result (BT_INTEGER, kind, &src->where);
2689
2690   mpz_set_f (result->value.integer, src->value.complex.r);
2691
2692   if ((rc = gfc_check_integer_range (result->value.integer, kind))
2693       != ARITH_OK)
2694     {
2695       arith_error (rc, &src->ts, &result->ts, &src->where);
2696       gfc_free_expr (result);
2697       return NULL;
2698     }
2699
2700   return result;
2701 }
2702
2703
2704 /* Convert complex to real.  */
2705
2706 gfc_expr *
2707 gfc_complex2real (gfc_expr * src, int kind)
2708 {
2709   gfc_expr *result;
2710   arith rc;
2711
2712   result = gfc_constant_result (BT_REAL, kind, &src->where);
2713
2714   mpf_set (result->value.real, src->value.complex.r);
2715
2716   if ((rc = gfc_check_real_range (result->value.real, kind)) != ARITH_OK)
2717     {
2718       arith_error (rc, &src->ts, &result->ts, &src->where);
2719       gfc_free_expr (result);
2720       return NULL;
2721     }
2722
2723   return result;
2724 }
2725
2726
2727 /* Convert complex to complex.  */
2728
2729 gfc_expr *
2730 gfc_complex2complex (gfc_expr * src, int kind)
2731 {
2732   gfc_expr *result;
2733   arith rc;
2734
2735   result = gfc_constant_result (BT_COMPLEX, kind, &src->where);
2736
2737   mpf_set (result->value.complex.r, src->value.complex.r);
2738   mpf_set (result->value.complex.i, src->value.complex.i);
2739
2740   if ((rc = gfc_check_real_range (result->value.complex.r, kind)) != ARITH_OK
2741       || (rc =
2742           gfc_check_real_range (result->value.complex.i, kind)) != ARITH_OK)
2743     {
2744       arith_error (rc, &src->ts, &result->ts, &src->where);
2745       gfc_free_expr (result);
2746       return NULL;
2747     }
2748
2749   return result;
2750 }
2751
2752
2753 /* Logical kind conversion.  */
2754
2755 gfc_expr *
2756 gfc_log2log (gfc_expr * src, int kind)
2757 {
2758   gfc_expr *result;
2759
2760   result = gfc_constant_result (BT_LOGICAL, kind, &src->where);
2761   result->value.logical = src->value.logical;
2762
2763   return result;
2764 }