+2004-08-06 Steven G. Kargl <kargls@comcast.net>
+
+ * arith.c: Add #define for model numbers. Remove global GMP variables.
+ (natural_logarithm,common_logarithm,exponential,sine,
+ cosine,arctangent,hypercos,hypersine ): Remove.
+ (gfc_mpfr_to_mpz,gfc_set_model_kind,gfc_set_model): New functions.
+ (arctangent2,gfc_arith_init_1,gfc_arith_done_1
+ gfc_check_real_range, gfc_constant_result, gfc_range_check,
+ gfc_arith_uminus,gfc_arith_plus, gfc_arith_minus, gfc_arith_times,
+ gfc_arith_divide,complex_reciprocal,complex_pow_ui,
+ gfc_arith_power,gfc_compare_expr,compare_complex,gfc_convert_real,
+ gfc_convert_complex,gfc_int2real,gfc_int2complex,
+ gfc_real2int,gfc_real2real,gfc_real2complex,
+ gfc_complex2int,gfc_complex2real,gfc_complex2complex): Convert GMP
+ to MPFR, use new functions.
+ * arith.h: Remove extern global variables.
+ (natural_logarithm,common_logarithm,exponential, sine, cosine,
+ arctangent,hypercos,hypersine): Remove prototypes.
+ (arctangent2): Update prototype from GMP to MPFR.
+ (gfc_mpfr_to_mpz, gfc_set_model_kind,gfc_set_model): Add prototypes.
+ * dump-parse-tree.c (gfc_show_expr): Convert GMP to MPFR.
+ * expr.c (free_expr0,gfc_copy_expr): Convert GMP to MPFR.
+ * gfortran.h (GFC_REAL_BITS): Remove.
+ (arith): Add ARITH_NAN.
+ Include mpfr.h. Define GFC_RND_MODE.
+ Rename GCC_GFORTRAN_H GFC_GFC_H.
+ (gfc_expr): Convert GMP to MPFR.
+ * module.c: Add arith.h, correct type in comment.
+ (mio_gmp_real): Convert GMP to MPFR.
+ (mio_expr): Use gfc_set_model_kind().
+ * primary.c: Update copyright date with 2004.
+ (match_real_constant,match_const_complex_part): Convert GMP to MPFR.
+ * simplify.c: Remove global GMP variables
+ (gfc_simplify_abs,gfc_simplify_acos,gfc_simplify_aimag,
+ gfc_simplify_aint,gfc_simplify_dint,gfc_simplify_anint,
+ gfc_simplify_dnint,gfc_simplify_asin,gfc_simplify_atan,
+ gfc_simplify_atan2,gfc_simplify_ceiling,simplify_cmplx,
+ gfc_simplify_conjg,gfc_simplify_cos,gfc_simplify_cosh,
+ gfc_simplify_dim,gfc_simplify_dprod,gfc_simplify_epsilon,
+ gfc_simplify_exp,gfc_simplify_exponent,gfc_simplify_floor,
+ gfc_simplify_fraction,gfc_simplify_huge,gfc_simplify_int,
+ gfc_simplify_ifix,gfc_simplify_idint,gfc_simplify_log,
+ gfc_simplify_log10,simplify_min_max,gfc_simplify_mod,
+ gfc_simplify_modulo,gfc_simplify_nearest,simplify_nint,
+ gfc_simplify_rrspacing,gfc_simplify_scale,
+ gfc_simplify_set_exponent,gfc_simplify_sign,gfc_simplify_sin,
+ gfc_simplify_sinh,gfc_simplify_spacing,gfc_simplify_sqrt,
+ gfc_simplify_tan,gfc_simplify_tanh,gfc_simplify_tiny,
+ gfc_simplify_init_1,gfc_simplify_done_1): Convert GMP to MPFR.
+ Use new functions.
+ * trans-const.c (gfc_conv_mpfr_to_tree): Rename from
+ gfc_conv_mpf_to_tree. Convert it to use MPFR
+ (gfc_conv_constant_to_tree): Use it.
+ * trans-const.h: Update prototype for gfc_conv_mpfr_to_tree().
+ * trans-intrinsic.c: Add arith.h, remove gmp.h
+ (gfc_conv_intrinsic_aint,gfc_conv_intrinsic_mod): Convert GMP to MPFR.
+
2004-08-06 Victor Leikehman <lei@il.ibm.com>
Paul Brook <paul@codesourcery.com>
#include "gfortran.h"
#include "arith.h"
-mpf_t pi, half_pi, two_pi, e;
-
/* The gfc_(integer|real)_kinds[] structures have everything the front
end needs to know about integers and real numbers on the target.
Other entries of the structure are calculated from these values.
DEF_GFC_LOGICAL_KIND (0, 0)
};
+
+/* IEEE-754 uses 1.xEe representation whereas the fortran standard
+ uses 0.xEe representation. Hence the exponents below are biased
+ by one. */
+
+#define GFC_SP_KIND 4
+#define GFC_SP_PREC 24 /* p = 24, IEEE-754 */
+#define GFC_SP_EMIN -125 /* emin = -126, IEEE-754 */
+#define GFC_SP_EMAX 128 /* emin = 127, IEEE-754 */
+
+/* Double precision model numbers. */
+#define GFC_DP_KIND 8
+#define GFC_DP_PREC 53 /* p = 53, IEEE-754 */
+#define GFC_DP_EMIN -1021 /* emin = -1022, IEEE-754 */
+#define GFC_DP_EMAX 1024 /* emin = 1023, IEEE-754 */
+
+/* Quad precision model numbers. Not used. */
+#define GFC_QP_KIND 16
+#define GFC_QP_PREC 113 /* p = 113, IEEE-754 */
+#define GFC_QP_EMIN -16381 /* emin = -16382, IEEE-754 */
+#define GFC_QP_EMAX 16384 /* emin = 16383, IEEE-754 */
+
gfc_real_info gfc_real_kinds[] = {
- DEF_GFC_REAL_KIND (4, 2, 24, -125, 128),
- DEF_GFC_REAL_KIND (8, 2, 53, -1021, 1024),
+ DEF_GFC_REAL_KIND (GFC_SP_KIND, 2, GFC_SP_PREC, GFC_SP_EMIN, GFC_SP_EMAX),
+ DEF_GFC_REAL_KIND (GFC_DP_KIND, 2, GFC_DP_PREC, GFC_DP_EMIN, GFC_DP_EMAX),
DEF_GFC_REAL_KIND (0, 0, 0, 0, 0)
};
int gfc_index_integer_kind;
-/* Compute the natural log of arg.
-
- We first get the argument into the range 0.5 to 1.5 by successive
- multiplications or divisions by e. Then we use the series:
-
- ln(x) = (x-1) - (x-1)^2/2 + (x-1)^3/3 - (x-1)^4/4 + ...
-
- Because we are expanding in powers of (x-1), and 0.5 < x < 1.5, we
- have -0.5 < (x-1) < 0.5. Ignoring the harmonic term, this means
- that each term is at most 1/(2^i), meaning one bit is gained per
- iteration.
-
- Not very efficient, but it doesn't have to be. */
+/* MPFR does not have a direct replacement for mpz_set_f() from GMP.
+ It's easily implemented with a few calls though. */
void
-natural_logarithm (mpf_t * arg, mpf_t * result)
+gfc_mpfr_to_mpz(mpz_t z, mpfr_t x)
{
- mpf_t x, xp, t, log;
- int i, p;
-
- mpf_init_set (x, *arg);
- mpf_init (t);
-
- p = 0;
-
- mpf_set_str (t, "0.5", 10);
- while (mpf_cmp (x, t) < 0)
- {
- mpf_mul (x, x, e);
- p--;
- }
-
- mpf_set_str (t, "1.5", 10);
- while (mpf_cmp (x, t) > 0)
- {
- mpf_div (x, x, e);
- p++;
- }
-
- mpf_sub_ui (x, x, 1);
- mpf_init_set_ui (log, 0);
- mpf_init_set_ui (xp, 1);
-
- for (i = 1; i < GFC_REAL_BITS; i++)
- {
- mpf_mul (xp, xp, x);
- mpf_div_ui (t, xp, i);
+ mp_exp_t e;
- if (i % 2 == 0)
- mpf_sub (log, log, t);
- else
- mpf_add (log, log, t);
- }
-
- /* Add in the log (e^p) = p */
-
- if (p < 0)
- mpf_sub_ui (log, log, -p);
+ e = mpfr_get_z_exp (z, x);
+ if (e > 0)
+ mpz_mul_2exp (z, z, e);
else
- mpf_add_ui (log, log, p);
-
- mpf_clear (x);
- mpf_clear (xp);
- mpf_clear (t);
-
- mpf_set (*result, log);
- mpf_clear (log);
-}
-
-
-/* Calculate the common logarithm of arg. We use the natural
- logarithm of arg and of 10:
-
- log10(arg) = log(arg)/log(10) */
-
-void
-common_logarithm (mpf_t * arg, mpf_t * result)
-{
- mpf_t i10, log10;
-
- natural_logarithm (arg, result);
-
- mpf_init_set_ui (i10, 10);
- mpf_init (log10);
- natural_logarithm (&i10, &log10);
-
- mpf_div (*result, *result, log10);
- mpf_clear (i10);
- mpf_clear (log10);
+ mpz_tdiv_q_2exp (z, z, -e);
+ if (mpfr_sgn (x) < 0)
+ mpz_neg (z, z);
}
-/* Calculate exp(arg).
-
- We use a reduction of the form
-
- x = Nln2 + r
- Then we obtain exp(r) from the Maclaurin series.
- exp(x) is then recovered from the identity
-
- exp(x) = 2^N*exp(r). */
+/* Set the model number precision by the requested KIND. */
void
-exponential (mpf_t * arg, mpf_t * result)
+gfc_set_model_kind (int kind)
{
- mpf_t two, ln2, power, q, r, num, denom, term, x, xp;
- int i;
- long n;
- unsigned long p, mp;
-
-
- mpf_init_set (x, *arg);
-
- if (mpf_cmp_ui (x, 0) == 0)
- {
- mpf_set_ui (*result, 1);
- }
- else if (mpf_cmp_ui (x, 1) == 0)
- {
- mpf_set (*result, e);
- }
- else
- {
- mpf_init_set_ui (two, 2);
- mpf_init (ln2);
- mpf_init (q);
- mpf_init (r);
- mpf_init (power);
- mpf_init (term);
-
- natural_logarithm (&two, &ln2);
-
- mpf_div (q, x, ln2);
- mpf_floor (power, q);
- mpf_mul (q, power, ln2);
- mpf_sub (r, x, q);
-
- mpf_init_set_ui (xp, 1);
- mpf_init_set_ui (num, 1);
- mpf_init_set_ui (denom, 1);
-
- for (i = 1; i <= GFC_REAL_BITS + 10; i++)
+ switch (kind)
{
- mpf_mul (num, num, r);
- mpf_mul_ui (denom, denom, i);
- mpf_div (term, num, denom);
- mpf_add (xp, xp, term);
- }
-
- /* Reconstruction step */
- n = (long) mpf_get_d (power);
-
- if (n > 0)
- {
- p = (unsigned int) n;
- mpf_mul_2exp (*result, xp, p);
- }
- else
- {
- mp = (unsigned int) (-n);
- mpf_div_2exp (*result, xp, mp);
- }
-
- mpf_clear (two);
- mpf_clear (ln2);
- mpf_clear (q);
- mpf_clear (r);
- mpf_clear (power);
- mpf_clear (num);
- mpf_clear (denom);
- mpf_clear (term);
- mpf_clear (xp);
- }
-
- mpf_clear (x);
-}
-
-
-/* Calculate sin(arg).
-
- We use a reduction of the form
-
- x= N*2pi + r
-
- Then we obtain sin(r) from the Maclaurin series. */
-
-void
-sine (mpf_t * arg, mpf_t * result)
-{
- mpf_t factor, q, r, num, denom, term, x, xp;
- int i, sign;
-
- mpf_init_set (x, *arg);
-
- /* Special case (we do not treat multiples of pi due to roundoff issues) */
- if (mpf_cmp_ui (x, 0) == 0)
- {
- mpf_set_ui (*result, 0);
- }
- else
- {
- mpf_init (q);
- mpf_init (r);
- mpf_init (factor);
- mpf_init (term);
-
- mpf_div (q, x, two_pi);
- mpf_floor (factor, q);
- mpf_mul (q, factor, two_pi);
- mpf_sub (r, x, q);
-
- mpf_init_set_ui (xp, 0);
- mpf_init_set_ui (num, 1);
- mpf_init_set_ui (denom, 1);
-
- sign = -1;
- for (i = 1; i < GFC_REAL_BITS + 10; i++)
- {
- mpf_mul (num, num, r);
- mpf_mul_ui (denom, denom, i);
- if (i % 2 == 0)
- continue;
-
- sign = -sign;
- mpf_div (term, num, denom);
- if (sign > 0)
- mpf_add (xp, xp, term);
- else
- mpf_sub (xp, xp, term);
- }
-
- mpf_set (*result, xp);
-
- mpf_clear (q);
- mpf_clear (r);
- mpf_clear (factor);
- mpf_clear (num);
- mpf_clear (denom);
- mpf_clear (term);
- mpf_clear (xp);
- }
-
- mpf_clear (x);
-}
-
-
-/* Calculate cos(arg).
-
- Similar to sine. */
-
-void
-cosine (mpf_t * arg, mpf_t * result)
-{
- mpf_t factor, q, r, num, denom, term, x, xp;
- int i, sign;
-
- mpf_init_set (x, *arg);
-
- /* Special case (we do not treat multiples of pi due to roundoff issues) */
- if (mpf_cmp_ui (x, 0) == 0)
- {
- mpf_set_ui (*result, 1);
- }
- else
- {
- mpf_init (q);
- mpf_init (r);
- mpf_init (factor);
- mpf_init (term);
-
- mpf_div (q, x, two_pi);
- mpf_floor (factor, q);
- mpf_mul (q, factor, two_pi);
- mpf_sub (r, x, q);
-
- mpf_init_set_ui (xp, 1);
- mpf_init_set_ui (num, 1);
- mpf_init_set_ui (denom, 1);
-
- sign = 1;
- for (i = 1; i < GFC_REAL_BITS + 10; i++)
- {
- mpf_mul (num, num, r);
- mpf_mul_ui (denom, denom, i);
- if (i % 2 != 0)
- continue;
-
- sign = -sign;
- mpf_div (term, num, denom);
- if (sign > 0)
- mpf_add (xp, xp, term);
- else
- mpf_sub (xp, xp, term);
- }
- mpf_set (*result, xp);
-
- mpf_clear (q);
- mpf_clear (r);
- mpf_clear (factor);
- mpf_clear (num);
- mpf_clear (denom);
- mpf_clear (term);
- mpf_clear (xp);
+ case GFC_SP_KIND:
+ mpfr_set_default_prec (GFC_SP_PREC);
+ break;
+ case GFC_DP_KIND:
+ mpfr_set_default_prec (GFC_DP_PREC);
+ break;
+ case GFC_QP_KIND:
+ mpfr_set_default_prec (GFC_QP_PREC);
+ break;
+ default:
+ gfc_internal_error ("gfc_set_model_kind(): Bad model number");
}
-
- mpf_clear (x);
}
-/* Calculate atan(arg).
-
- Similar to sine but requires special handling for x near 1. */
+/* Set the model number precision from mpfr_t x. */
void
-arctangent (mpf_t * arg, mpf_t * result)
+gfc_set_model (mpfr_t x)
{
- mpf_t absval, convgu, convgl, num, term, x, xp;
- int i, sign;
-
- mpf_init_set (x, *arg);
-
- /* Special cases */
- if (mpf_cmp_ui (x, 0) == 0)
+ switch (mpfr_get_prec (x))
{
- mpf_set_ui (*result, 0);
- }
- else if (mpf_cmp_ui (x, 1) == 0)
- {
- mpf_init (num);
- mpf_div_ui (num, half_pi, 2);
- mpf_set (*result, num);
- mpf_clear (num);
- }
- else if (mpf_cmp_si (x, -1) == 0)
- {
- mpf_init (num);
- mpf_div_ui (num, half_pi, 2);
- mpf_neg (*result, num);
- mpf_clear (num);
- }
- else
- { /* General cases */
-
- mpf_init (absval);
- mpf_abs (absval, x);
-
- mpf_init_set_d (convgu, 1.5);
- mpf_init_set_d (convgl, 0.5);
- mpf_init_set_ui (num, 1);
- mpf_init (term);
-
- if (mpf_cmp (absval, convgl) < 0)
- {
- mpf_init_set_ui (xp, 0);
- sign = -1;
- for (i = 1; i < GFC_REAL_BITS + 10; i++)
- {
- mpf_mul (num, num, absval);
- if (i % 2 == 0)
- continue;
-
- sign = -sign;
- mpf_div_ui (term, num, i);
- if (sign > 0)
- mpf_add (xp, xp, term);
- else
- mpf_sub (xp, xp, term);
- }
- }
- else if (mpf_cmp (absval, convgu) >= 0)
- {
- mpf_init_set (xp, half_pi);
- sign = 1;
- for (i = 1; i < GFC_REAL_BITS + 10; i++)
- {
- mpf_div (num, num, absval);
- if (i % 2 == 0)
- continue;
-
- sign = -sign;
- mpf_div_ui (term, num, i);
- if (sign > 0)
- mpf_add (xp, xp, term);
- else
- mpf_sub (xp, xp, term);
- }
- }
- else
- {
- mpf_init_set_ui (xp, 0);
-
- mpf_sub_ui (num, absval, 1);
- mpf_add_ui (term, absval, 1);
- mpf_div (absval, num, term);
-
- mpf_set_ui (num, 1);
-
- sign = -1;
- for (i = 1; i < GFC_REAL_BITS + 10; i++)
- {
- mpf_mul (num, num, absval);
- if (i % 2 == 0)
- continue;
- sign = -sign;
- mpf_div_ui (term, num, i);
- if (sign > 0)
- mpf_add (xp, xp, term);
- else
- mpf_sub (xp, xp, term);
- }
-
- mpf_div_ui (term, half_pi, 2);
- mpf_add (xp, term, xp);
- }
-
- /* This makes sure to preserve the identity arctan(-x) = -arctan(x)
- and improves accuracy to boot. */
-
- if (mpf_cmp_ui (x, 0) > 0)
- mpf_set (*result, xp);
- else
- mpf_neg (*result, xp);
-
- mpf_clear (absval);
- mpf_clear (convgl);
- mpf_clear (convgu);
- mpf_clear (num);
- mpf_clear (term);
- mpf_clear (xp);
+ case GFC_SP_PREC:
+ mpfr_set_default_prec (GFC_SP_PREC);
+ break;
+ case GFC_DP_PREC:
+ mpfr_set_default_prec (GFC_DP_PREC);
+ break;
+ case GFC_QP_PREC:
+ mpfr_set_default_prec (GFC_QP_PREC);
+ break;
+ default:
+ gfc_internal_error ("gfc_set_model(): Bad model number");
}
- mpf_clear (x);
}
-
/* Calculate atan2 (y, x)
atan2(y, x) = atan(y/x) if x > 0,
*/
void
-arctangent2 (mpf_t * y, mpf_t * x, mpf_t * result)
+arctangent2 (mpfr_t y, mpfr_t x, mpfr_t result)
{
- mpf_t t;
+ int i;
+ mpfr_t t;
- mpf_init (t);
+ gfc_set_model (y);
+ mpfr_init (t);
- switch (mpf_sgn (*x))
+ i = mpfr_sgn(x);
+
+ if (i > 0)
{
- case 1:
- mpf_div (t, *y, *x);
- arctangent (&t, result);
- break;
- case -1:
- mpf_div (t, *y, *x);
- mpf_abs (t, t);
- arctangent (&t, &t);
- mpf_sub (*result, pi, t);
- if (mpf_sgn (*y) == -1)
- mpf_neg (*result, *result);
- break;
- case 0:
- if (mpf_sgn (*y) == 0)
- mpf_set_ui (*result, 0);
+ mpfr_div (t, y, x, GFC_RND_MODE);
+ mpfr_atan (result, t, GFC_RND_MODE);
+ }
+ else if (i < 0)
+ {
+ mpfr_const_pi (result, GFC_RND_MODE);
+ mpfr_div (t, y, x, GFC_RND_MODE);
+ mpfr_abs (t, t, GFC_RND_MODE);
+ mpfr_atan (t, t, GFC_RND_MODE);
+ mpfr_sub (result, result, t, GFC_RND_MODE);
+ if (mpfr_sgn (y) < 0)
+ mpfr_neg (result, result, GFC_RND_MODE);
+ }
+ else
+ {
+ if (mpfr_sgn (y) == 0)
+ mpfr_set_ui (result, 0, GFC_RND_MODE);
else
{
- mpf_set (*result, half_pi);
- if (mpf_sgn (*y) == -1)
- mpf_neg (*result, *result);
+ mpfr_const_pi (result, GFC_RND_MODE);
+ mpfr_div_ui (result, result, 2, GFC_RND_MODE);
+ if (mpfr_sgn (y) < 0)
+ mpfr_neg (result, result, GFC_RND_MODE);
}
- break;
}
- mpf_clear (t);
-}
-
-/* Calculate cosh(arg). */
-
-void
-hypercos (mpf_t * arg, mpf_t * result)
-{
- mpf_t neg, term1, term2, x, xp;
-
- mpf_init_set (x, *arg);
- mpf_init (neg);
- mpf_init (term1);
- mpf_init (term2);
- mpf_init (xp);
+ mpfr_clear (t);
- mpf_neg (neg, x);
-
- exponential (&x, &term1);
- exponential (&neg, &term2);
-
- mpf_add (xp, term1, term2);
- mpf_div_ui (*result, xp, 2);
-
- mpf_clear (neg);
- mpf_clear (term1);
- mpf_clear (term2);
- mpf_clear (x);
- mpf_clear (xp);
-}
-
-
-/* Calculate sinh(arg). */
-
-void
-hypersine (mpf_t * arg, mpf_t * result)
-{
- mpf_t neg, term1, term2, x, xp;
-
- mpf_init_set (x, *arg);
-
- mpf_init (neg);
- mpf_init (term1);
- mpf_init (term2);
- mpf_init (xp);
-
- mpf_neg (neg, x);
-
- exponential (&x, &term1);
- exponential (&neg, &term2);
-
- mpf_sub (xp, term1, term2);
- mpf_div_ui (*result, xp, 2);
-
- mpf_clear (neg);
- mpf_clear (term1);
- mpf_clear (term2);
- mpf_clear (x);
- mpf_clear (xp);
}
case ARITH_UNDERFLOW:
p = "Arithmetic underflow";
break;
+ case ARITH_NAN:
+ p = "Arithmetic NaN";
+ break;
case ARITH_DIV0:
p = "Division by zero";
break;
{
gfc_integer_info *int_info;
gfc_real_info *real_info;
- mpf_t a, b;
+ mpfr_t a, b, c;
mpz_t r;
- int i, n, limit;
-
- /* Set the default precision for GMP computations. */
- mpf_set_default_prec (GFC_REAL_BITS + 30);
-
- /* Calculate e, needed by the natural_logarithm() subroutine. */
- mpf_init (b);
- mpf_init_set_ui (e, 0);
- mpf_init_set_ui (a, 1);
-
- for (i = 1; i < 100; i++)
- {
- mpf_add (e, e, a);
- mpf_div_ui (a, a, i); /* 1/(i!) */
- }
-
- /* Calculate pi, 2pi, pi/2, and -pi/2, needed for trigonometric
- functions.
-
- We use the Bailey, Borwein and Plouffe formula:
-
- pi = \sum{n=0}^\infty (1/16)^n [4/(8n+1) - 2/(8n+4) - 1/(8n+5) - 1/(8n+6)]
-
- which gives about four bits per iteration. */
-
- mpf_init_set_ui (pi, 0);
-
- mpf_init (two_pi);
- mpf_init (half_pi);
-
- limit = (GFC_REAL_BITS / 4) + 10; /* (1/16)^n gives 4 bits per iteration */
-
- for (n = 0; n < limit; n++)
- {
- mpf_set_ui (b, 4);
- mpf_div_ui (b, b, 8 * n + 1); /* 4/(8n+1) */
-
- mpf_set_ui (a, 2);
- mpf_div_ui (a, a, 8 * n + 4); /* 2/(8n+4) */
- mpf_sub (b, b, a);
-
- mpf_set_ui (a, 1);
- mpf_div_ui (a, a, 8 * n + 5); /* 1/(8n+5) */
- mpf_sub (b, b, a);
-
- mpf_set_ui (a, 1);
- mpf_div_ui (a, a, 8 * n + 6); /* 1/(8n+6) */
- mpf_sub (b, b, a);
-
- mpf_set_ui (a, 16);
- mpf_pow_ui (a, a, n); /* 16^n */
-
- mpf_div (b, b, a);
+ int i;
- mpf_add (pi, pi, b);
- }
+ gfc_set_model_kind (GFC_QP_KIND);
- mpf_mul_ui (two_pi, pi, 2);
- mpf_div_ui (half_pi, pi, 2);
+ mpfr_init (a);
+ mpz_init (r);
/* Convert the minimum/maximum values for each kind into their
GNU MP representation. */
- mpz_init (r);
-
for (int_info = gfc_integer_kinds; int_info->kind != 0; int_info++)
{
/* Huge */
mpz_add_ui (int_info->max_int, int_info->max_int, 1);
/* Range */
- mpf_set_z (a, int_info->huge);
- common_logarithm (&a, &a);
- mpf_trunc (a, a);
- mpz_set_f (r, a);
+ mpfr_set_z (a, int_info->huge, GFC_RND_MODE);
+ mpfr_log10 (a, a, GFC_RND_MODE);
+ mpfr_trunc (a, a);
+ gfc_mpfr_to_mpz (r, a);
int_info->range = mpz_get_si (r);
}
- /* mpf_set_default_prec(GFC_REAL_BITS); */
+ mpfr_clear (a);
+
for (real_info = gfc_real_kinds; real_info->kind != 0; real_info++)
{
- /* Huge */
- mpf_set_ui (a, real_info->radix);
- mpf_set_ui (b, real_info->radix);
+ gfc_set_model_kind (real_info->kind);
- mpf_pow_ui (a, a, real_info->max_exponent);
- mpf_pow_ui (b, b, real_info->max_exponent - real_info->digits);
+ mpfr_init (a);
+ mpfr_init (b);
+ mpfr_init (c);
- mpf_init (real_info->huge);
- mpf_sub (real_info->huge, a, b);
+ /* huge(x) = (1 - b**(-p)) * b**(emax-1) * b */
+ /* a = 1 - b**(-p) */
+ mpfr_set_ui (a, 1, GFC_RND_MODE);
+ mpfr_set_ui (b, real_info->radix, GFC_RND_MODE);
+ mpfr_pow_si (b, b, -real_info->digits, GFC_RND_MODE);
+ mpfr_sub (a, a, b, GFC_RND_MODE);
- /* Tiny */
- mpf_set_ui (b, real_info->radix);
- mpf_pow_ui (b, b, 1 - real_info->min_exponent);
+ /* c = b**(emax-1) */
+ mpfr_set_ui (b, real_info->radix, GFC_RND_MODE);
+ mpfr_pow_ui (c, b, real_info->max_exponent - 1, GFC_RND_MODE);
- mpf_init (real_info->tiny);
- mpf_ui_div (real_info->tiny, 1, b);
+ /* a = a * c = (1 - b**(-p)) * b**(emax-1) */
+ mpfr_mul (a, a, c, GFC_RND_MODE);
- /* Epsilon */
- mpf_set_ui (b, real_info->radix);
- mpf_pow_ui (b, b, real_info->digits - 1);
+ /* a = (1 - b**(-p)) * b**(emax-1) * b */
+ mpfr_mul_ui (a, a, real_info->radix, GFC_RND_MODE);
- mpf_init (real_info->epsilon);
- mpf_ui_div (real_info->epsilon, 1, b);
+ mpfr_init (real_info->huge);
+ mpfr_set (real_info->huge, a, GFC_RND_MODE);
- /* Range */
- common_logarithm (&real_info->huge, &a);
- common_logarithm (&real_info->tiny, &b);
- mpf_neg (b, b);
+ /* tiny(x) = b**(emin-1) */
+ mpfr_set_ui (b, real_info->radix, GFC_RND_MODE);
+ mpfr_pow_si (b, b, real_info->min_exponent - 1, GFC_RND_MODE);
+
+ mpfr_init (real_info->tiny);
+ mpfr_set (real_info->tiny, b, GFC_RND_MODE);
+
+ /* epsilon(x) = b**(1-p) */
+ mpfr_set_ui (b, real_info->radix, GFC_RND_MODE);
+ mpfr_pow_si (b, b, 1 - real_info->digits, GFC_RND_MODE);
+
+ mpfr_init (real_info->epsilon);
+ mpfr_set (real_info->epsilon, b, GFC_RND_MODE);
- if (mpf_cmp (a, b) > 0)
- mpf_set (a, b); /* a = min(a, b) */
+ /* range(x) = int(min(log10(huge(x)), -log10(tiny)) */
+ mpfr_log10 (a, real_info->huge, GFC_RND_MODE);
+ mpfr_log10 (b, real_info->tiny, GFC_RND_MODE);
+ mpfr_neg (b, b, GFC_RND_MODE);
- mpf_trunc (a, a);
- mpz_set_f (r, a);
+ if (mpfr_cmp (a, b) > 0)
+ mpfr_set (a, b, GFC_RND_MODE); /* a = min(a, b) */
+
+ mpfr_trunc (a, a);
+ gfc_mpfr_to_mpz (r, a);
real_info->range = mpz_get_si (r);
- /* Precision */
- mpf_set_ui (a, real_info->radix);
- common_logarithm (&a, &a);
+ /* precision(x) = int((p - 1) * log10(b)) + k */
+ mpfr_set_ui (a, real_info->radix, GFC_RND_MODE);
+ mpfr_log10 (a, a, GFC_RND_MODE);
- mpf_mul_ui (a, a, real_info->digits - 1);
- mpf_trunc (a, a);
- mpz_set_f (r, a);
+ mpfr_mul_ui (a, a, real_info->digits - 1, GFC_RND_MODE);
+ mpfr_trunc (a, a);
+ gfc_mpfr_to_mpz (r, a);
real_info->precision = mpz_get_si (r);
/* If the radix is an integral power of 10, add one to the
for (i = 10; i <= real_info->radix; i *= 10)
if (i == real_info->radix)
real_info->precision++;
+
+ mpfr_clear (a);
+ mpfr_clear (b);
+ mpfr_clear (c);
}
mpz_clear (r);
- mpf_clear (a);
- mpf_clear (b);
}
gfc_integer_info *ip;
gfc_real_info *rp;
- mpf_clear (e);
-
- mpf_clear (pi);
- mpf_clear (half_pi);
- mpf_clear (two_pi);
-
for (ip = gfc_integer_kinds; ip->kind; ip++)
{
mpz_clear (ip->min_int);
for (rp = gfc_real_kinds; rp->kind; rp++)
{
- mpf_clear (rp->epsilon);
- mpf_clear (rp->huge);
- mpf_clear (rp->tiny);
+ mpfr_clear (rp->epsilon);
+ mpfr_clear (rp->huge);
+ mpfr_clear (rp->tiny);
}
}
ARITH_UNDERFLOW. */
static arith
-gfc_check_real_range (mpf_t p, int kind)
+gfc_check_real_range (mpfr_t p, int kind)
{
arith retval;
- mpf_t q;
+ mpfr_t q;
int i;
- mpf_init (q);
- mpf_abs (q, p);
-
i = validate_real (kind);
if (i == -1)
gfc_internal_error ("gfc_check_real_range(): Bad kind");
+ gfc_set_model (p);
+ mpfr_init (q);
+ mpfr_abs (q, p, GFC_RND_MODE);
+
retval = ARITH_OK;
- if (mpf_sgn (q) == 0)
+ if (mpfr_sgn (q) == 0)
goto done;
- if (mpf_cmp (q, gfc_real_kinds[i].huge) == 1)
+ if (mpfr_cmp (q, gfc_real_kinds[i].huge) > 0)
{
retval = ARITH_OVERFLOW;
goto done;
}
- if (mpf_cmp (q, gfc_real_kinds[i].tiny) == -1)
+ if (mpfr_cmp (q, gfc_real_kinds[i].tiny) < 0)
retval = ARITH_UNDERFLOW;
done:
- mpf_clear (q);
+ mpfr_clear (q);
return retval;
}
break;
case BT_REAL:
- mpf_init (result->value.real);
+ gfc_set_model_kind (kind);
+ mpfr_init (result->value.real);
break;
case BT_COMPLEX:
- mpf_init (result->value.complex.r);
- mpf_init (result->value.complex.i);
+ gfc_set_model_kind (kind);
+ mpfr_init (result->value.complex.r);
+ mpfr_init (result->value.complex.i);
break;
default:
/* Make sure a constant numeric expression is within the range for
- its type and kind. GMP is doing 130 bit arithmetic, so an UNDERFLOW
- is numerically zero for REAL(4) and REAL(8) types. Reset the value(s)
- to exactly 0 for UNDERFLOW. Note that there's also a gfc_check_range(),
+ its type and kind. Note that there's also a gfc_check_range(),
but that one deals with the intrinsic RANGE function. */
arith
case BT_REAL:
rc = gfc_check_real_range (e->value.real, e->ts.kind);
if (rc == ARITH_UNDERFLOW)
- mpf_set_ui (e->value.real, 0);
+ mpfr_set_ui (e->value.real, 0, GFC_RND_MODE);
break;
case BT_COMPLEX:
rc = gfc_check_real_range (e->value.complex.r, e->ts.kind);
if (rc == ARITH_UNDERFLOW)
- mpf_set_ui (e->value.complex.r, 0);
+ mpfr_set_ui (e->value.complex.r, 0, GFC_RND_MODE);
if (rc == ARITH_OK || rc == ARITH_UNDERFLOW)
{
rc = gfc_check_real_range (e->value.complex.i, e->ts.kind);
if (rc == ARITH_UNDERFLOW)
- mpf_set_ui (e->value.complex.i, 0);
+ mpfr_set_ui (e->value.complex.i, 0, GFC_RND_MODE);
}
break;
break;
case BT_REAL:
- mpf_neg (result->value.real, op1->value.real);
+ mpfr_neg (result->value.real, op1->value.real, GFC_RND_MODE);
break;
case BT_COMPLEX:
- mpf_neg (result->value.complex.r, op1->value.complex.r);
- mpf_neg (result->value.complex.i, op1->value.complex.i);
+ mpfr_neg (result->value.complex.r, op1->value.complex.r, GFC_RND_MODE);
+ mpfr_neg (result->value.complex.i, op1->value.complex.i, GFC_RND_MODE);
break;
default:
break;
case BT_REAL:
- mpf_add (result->value.real, op1->value.real, op2->value.real);
+ mpfr_add (result->value.real, op1->value.real, op2->value.real,
+ GFC_RND_MODE);
break;
case BT_COMPLEX:
- mpf_add (result->value.complex.r, op1->value.complex.r,
- op2->value.complex.r);
+ mpfr_add (result->value.complex.r, op1->value.complex.r,
+ op2->value.complex.r, GFC_RND_MODE);
- mpf_add (result->value.complex.i, op1->value.complex.i,
- op2->value.complex.i);
+ mpfr_add (result->value.complex.i, op1->value.complex.i,
+ op2->value.complex.i, GFC_RND_MODE);
break;
default:
break;
case BT_REAL:
- mpf_sub (result->value.real, op1->value.real, op2->value.real);
+ mpfr_sub (result->value.real, op1->value.real, op2->value.real,
+ GFC_RND_MODE);
break;
case BT_COMPLEX:
- mpf_sub (result->value.complex.r, op1->value.complex.r,
- op2->value.complex.r);
-
- mpf_sub (result->value.complex.i, op1->value.complex.i,
- op2->value.complex.i);
+ mpfr_sub (result->value.complex.r, op1->value.complex.r,
+ op2->value.complex.r, GFC_RND_MODE);
+ mpfr_sub (result->value.complex.i, op1->value.complex.i,
+ op2->value.complex.i, GFC_RND_MODE);
break;
default:
gfc_arith_times (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
{
gfc_expr *result;
- mpf_t x, y;
+ mpfr_t x, y;
arith rc;
result = gfc_constant_result (op1->ts.type, op1->ts.kind, &op1->where);
break;
case BT_REAL:
- mpf_mul (result->value.real, op1->value.real, op2->value.real);
+ mpfr_mul (result->value.real, op1->value.real, op2->value.real,
+ GFC_RND_MODE);
break;
case BT_COMPLEX:
- mpf_init (x);
- mpf_init (y);
- mpf_mul (x, op1->value.complex.r, op2->value.complex.r);
- mpf_mul (y, op1->value.complex.i, op2->value.complex.i);
- mpf_sub (result->value.complex.r, x, y);
+ /* FIXME: possible numericals problem. */
- mpf_mul (x, op1->value.complex.r, op2->value.complex.i);
- mpf_mul (y, op1->value.complex.i, op2->value.complex.r);
- mpf_add (result->value.complex.i, x, y);
+ gfc_set_model (op1->value.complex.r);
+ mpfr_init (x);
+ mpfr_init (y);
- mpf_clear (x);
- mpf_clear (y);
+ mpfr_mul (x, op1->value.complex.r, op2->value.complex.r, GFC_RND_MODE);
+ mpfr_mul (y, op1->value.complex.i, op2->value.complex.i, GFC_RND_MODE);
+ mpfr_sub (result->value.complex.r, x, y, GFC_RND_MODE);
+
+ mpfr_mul (x, op1->value.complex.r, op2->value.complex.i, GFC_RND_MODE);
+ mpfr_mul (y, op1->value.complex.i, op2->value.complex.r, GFC_RND_MODE);
+ mpfr_add (result->value.complex.i, x, y, GFC_RND_MODE);
+
+ mpfr_clear (x);
+ mpfr_clear (y);
break;
gfc_arith_divide (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
{
gfc_expr *result;
- mpf_t x, y, div;
+ mpfr_t x, y, div;
arith rc;
rc = ARITH_OK;
break;
case BT_REAL:
- if (mpf_sgn (op2->value.real) == 0)
+ /* FIXME: MPFR correctly generates NaN. This may not be needed. */
+ if (mpfr_sgn (op2->value.real) == 0)
{
rc = ARITH_DIV0;
break;
}
- mpf_div (result->value.real, op1->value.real, op2->value.real);
+ mpfr_div (result->value.real, op1->value.real, op2->value.real,
+ GFC_RND_MODE);
break;
case BT_COMPLEX:
- if (mpf_sgn (op2->value.complex.r) == 0
- && mpf_sgn (op2->value.complex.i) == 0)
+ /* FIXME: MPFR correctly generates NaN. This may not be needed. */
+ if (mpfr_sgn (op2->value.complex.r) == 0
+ && mpfr_sgn (op2->value.complex.i) == 0)
{
rc = ARITH_DIV0;
break;
}
- mpf_init (x);
- mpf_init (y);
- mpf_init (div);
+ gfc_set_model (op1->value.complex.r);
+ mpfr_init (x);
+ mpfr_init (y);
+ mpfr_init (div);
- mpf_mul (x, op2->value.complex.r, op2->value.complex.r);
- mpf_mul (y, op2->value.complex.i, op2->value.complex.i);
- mpf_add (div, x, y);
+ /* FIXME: possible numerical problems. */
+ mpfr_mul (x, op2->value.complex.r, op2->value.complex.r, GFC_RND_MODE);
+ mpfr_mul (y, op2->value.complex.i, op2->value.complex.i, GFC_RND_MODE);
+ mpfr_add (div, x, y, GFC_RND_MODE);
- mpf_mul (x, op1->value.complex.r, op2->value.complex.r);
- mpf_mul (y, op1->value.complex.i, op2->value.complex.i);
- mpf_add (result->value.complex.r, x, y);
- mpf_div (result->value.complex.r, result->value.complex.r, div);
+ mpfr_mul (x, op1->value.complex.r, op2->value.complex.r, GFC_RND_MODE);
+ mpfr_mul (y, op1->value.complex.i, op2->value.complex.i, GFC_RND_MODE);
+ mpfr_add (result->value.complex.r, x, y, GFC_RND_MODE);
+ mpfr_div (result->value.complex.r, result->value.complex.r, div,
+ GFC_RND_MODE);
- mpf_mul (x, op1->value.complex.i, op2->value.complex.r);
- mpf_mul (y, op1->value.complex.r, op2->value.complex.i);
- mpf_sub (result->value.complex.i, x, y);
- mpf_div (result->value.complex.i, result->value.complex.i, div);
+ mpfr_mul (x, op1->value.complex.i, op2->value.complex.r, GFC_RND_MODE);
+ mpfr_mul (y, op1->value.complex.r, op2->value.complex.i, GFC_RND_MODE);
+ mpfr_sub (result->value.complex.i, x, y, GFC_RND_MODE);
+ mpfr_div (result->value.complex.i, result->value.complex.i, div,
+ GFC_RND_MODE);
- mpf_clear (x);
- mpf_clear (y);
- mpf_clear (div);
+ mpfr_clear (x);
+ mpfr_clear (y);
+ mpfr_clear (div);
break;
static void
complex_reciprocal (gfc_expr * op)
{
- mpf_t mod, a, result_r, result_i;
-
- mpf_init (mod);
- mpf_init (a);
+ mpfr_t mod, a, re, im;
- mpf_mul (mod, op->value.complex.r, op->value.complex.r);
- mpf_mul (a, op->value.complex.i, op->value.complex.i);
- mpf_add (mod, mod, a);
+ gfc_set_model (op->value.complex.r);
+ mpfr_init (mod);
+ mpfr_init (a);
+ mpfr_init (re);
+ mpfr_init (im);
- mpf_init (result_r);
- mpf_div (result_r, op->value.complex.r, mod);
+ /* FIXME: another possible numerical problem. */
+ mpfr_mul (mod, op->value.complex.r, op->value.complex.r, GFC_RND_MODE);
+ mpfr_mul (a, op->value.complex.i, op->value.complex.i, GFC_RND_MODE);
+ mpfr_add (mod, mod, a, GFC_RND_MODE);
- mpf_init (result_i);
- mpf_neg (result_i, op->value.complex.i);
- mpf_div (result_i, result_i, mod);
+ mpfr_div (re, op->value.complex.r, mod, GFC_RND_MODE);
- mpf_set (op->value.complex.r, result_r);
- mpf_set (op->value.complex.i, result_i);
+ mpfr_neg (im, op->value.complex.i, GFC_RND_MODE);
+ mpfr_div (im, im, mod, GFC_RND_MODE);
- mpf_clear (result_r);
- mpf_clear (result_i);
+ mpfr_set (op->value.complex.r, re, GFC_RND_MODE);
+ mpfr_set (op->value.complex.i, im, GFC_RND_MODE);
- mpf_clear (mod);
- mpf_clear (a);
+ mpfr_clear (re);
+ mpfr_clear (im);
+ mpfr_clear (mod);
+ mpfr_clear (a);
}
static void
complex_pow_ui (gfc_expr * base, int power, gfc_expr * result)
{
- mpf_t temp_r, temp_i, a;
+ mpfr_t re, im, a;
- mpf_set_ui (result->value.complex.r, 1);
- mpf_set_ui (result->value.complex.i, 0);
+ gfc_set_model (base->value.complex.r);
+ mpfr_init (re);
+ mpfr_init (im);
+ mpfr_init (a);
- mpf_init (temp_r);
- mpf_init (temp_i);
- mpf_init (a);
+ mpfr_set_ui (result->value.complex.r, 1, GFC_RND_MODE);
+ mpfr_set_ui (result->value.complex.i, 0, GFC_RND_MODE);
for (; power > 0; power--)
{
- mpf_mul (temp_r, base->value.complex.r, result->value.complex.r);
- mpf_mul (a, base->value.complex.i, result->value.complex.i);
- mpf_sub (temp_r, temp_r, a);
+ mpfr_mul (re, base->value.complex.r, result->value.complex.r,
+ GFC_RND_MODE);
+ mpfr_mul (a, base->value.complex.i, result->value.complex.i,
+ GFC_RND_MODE);
+ mpfr_sub (re, re, a, GFC_RND_MODE);
- mpf_mul (temp_i, base->value.complex.r, result->value.complex.i);
- mpf_mul (a, base->value.complex.i, result->value.complex.r);
- mpf_add (temp_i, temp_i, a);
+ mpfr_mul (im, base->value.complex.r, result->value.complex.i,
+ GFC_RND_MODE);
+ mpfr_mul (a, base->value.complex.i, result->value.complex.r,
+ GFC_RND_MODE);
+ mpfr_add (im, im, a, GFC_RND_MODE);
- mpf_set (result->value.complex.r, temp_r);
- mpf_set (result->value.complex.i, temp_i);
+ mpfr_set (result->value.complex.r, re, GFC_RND_MODE);
+ mpfr_set (result->value.complex.i, im, GFC_RND_MODE);
}
- mpf_clear (temp_r);
- mpf_clear (temp_i);
- mpf_clear (a);
+ mpfr_clear (re);
+ mpfr_clear (im);
+ mpfr_clear (a);
}
int power, apower;
gfc_expr *result;
mpz_t unity_z;
- mpf_t unity_f;
+ mpfr_t unity_f;
arith rc;
rc = ARITH_OK;
rc = ARITH_0TO0;
else
mpz_set_ui (result->value.integer, 1);
-
break;
case BT_REAL:
- if (mpf_sgn (op1->value.real) == 0)
+ if (mpfr_sgn (op1->value.real) == 0)
rc = ARITH_0TO0;
else
- mpf_set_ui (result->value.real, 1);
-
+ mpfr_set_ui (result->value.real, 1, GFC_RND_MODE);
break;
case BT_COMPLEX:
- if (mpf_sgn (op1->value.complex.r) == 0
- && mpf_sgn (op1->value.complex.i) == 0)
+ if (mpfr_sgn (op1->value.complex.r) == 0
+ && mpfr_sgn (op1->value.complex.i) == 0)
rc = ARITH_0TO0;
else
{
- mpf_set_ui (result->value.complex.r, 1);
- mpf_set_ui (result->value.complex.i, 0);
+ mpfr_set_ui (result->value.complex.r, 1, GFC_RND_MODE);
+ mpfr_set_ui (result->value.complex.i, 0, GFC_RND_MODE);
}
break;
gfc_internal_error ("gfc_arith_power(): Bad base");
}
}
-
- if (power != 0)
+ else
{
apower = power;
if (power < 0)
break;
case BT_REAL:
- mpf_pow_ui (result->value.real, op1->value.real, apower);
+ mpfr_pow_ui (result->value.real, op1->value.real, apower,
+ GFC_RND_MODE);
if (power < 0)
{
- mpf_init_set_ui (unity_f, 1);
- mpf_div (result->value.real, unity_f, result->value.real);
- mpf_clear (unity_f);
+ gfc_set_model (op1->value.real);
+ mpfr_init (unity_f);
+ mpfr_set_ui (unity_f, 1, GFC_RND_MODE);
+ mpfr_div (result->value.real, unity_f, result->value.real,
+ GFC_RND_MODE);
+ mpfr_clear (unity_f);
}
-
break;
case BT_COMPLEX:
complex_pow_ui (op1, apower, result);
if (power < 0)
complex_reciprocal (result);
-
break;
default:
break;
case BT_REAL:
- rc = mpf_cmp (op1->value.real, op2->value.real);
+ rc = mpfr_cmp (op1->value.real, op2->value.real);
break;
case BT_CHARACTER:
compare_complex (gfc_expr * op1, gfc_expr * op2)
{
- return (mpf_cmp (op1->value.complex.r, op2->value.complex.r) == 0
- && mpf_cmp (op1->value.complex.i, op2->value.complex.i) == 0);
+ return (mpfr_cmp (op1->value.complex.r, op2->value.complex.r) == 0
+ && mpfr_cmp (op1->value.complex.i, op2->value.complex.i) == 0);
}
const char *t;
e = gfc_constant_result (BT_REAL, kind, where);
- /* a leading plus is allowed, but not by mpf_set_str */
+ /* A leading plus is allowed in Fortran, but not by mpfr_set_str */
if (buffer[0] == '+')
t = buffer + 1;
else
t = buffer;
- mpf_set_str (e->value.real, t, 10);
+ mpfr_set_str (e->value.real, t, 10, GFC_RND_MODE);
return e;
}
gfc_expr *e;
e = gfc_constant_result (BT_COMPLEX, kind, &real->where);
- mpf_set (e->value.complex.r, real->value.real);
- mpf_set (e->value.complex.i, imag->value.real);
+ mpfr_set (e->value.complex.r, real->value.real, GFC_RND_MODE);
+ mpfr_set (e->value.complex.i, imag->value.real, GFC_RND_MODE);
return e;
}
result = gfc_constant_result (BT_REAL, kind, &src->where);
- mpf_set_z (result->value.real, src->value.integer);
+ mpfr_set_z (result->value.real, src->value.integer, GFC_RND_MODE);
if ((rc = gfc_check_real_range (result->value.real, kind)) != ARITH_OK)
{
result = gfc_constant_result (BT_COMPLEX, kind, &src->where);
- mpf_set_z (result->value.complex.r, src->value.integer);
- mpf_set_ui (result->value.complex.i, 0);
+ mpfr_set_z (result->value.complex.r, src->value.integer, GFC_RND_MODE);
+ mpfr_set_ui (result->value.complex.i, 0, GFC_RND_MODE);
if ((rc = gfc_check_real_range (result->value.complex.r, kind)) != ARITH_OK)
{
result = gfc_constant_result (BT_INTEGER, kind, &src->where);
- mpz_set_f (result->value.integer, src->value.real);
+ gfc_mpfr_to_mpz (result->value.integer, src->value.real);
if ((rc = gfc_check_integer_range (result->value.integer, kind))
!= ARITH_OK)
result = gfc_constant_result (BT_REAL, kind, &src->where);
- mpf_set (result->value.real, src->value.real);
+ mpfr_set (result->value.real, src->value.real, GFC_RND_MODE);
rc = gfc_check_real_range (result->value.real, kind);
{
if (gfc_option.warn_underflow)
gfc_warning ("%s at %L", gfc_arith_error (rc), &src->where);
- mpf_set_ui(result->value.real, 0);
+ mpfr_set_ui(result->value.real, 0, GFC_RND_MODE);
}
else if (rc != ARITH_OK)
{
result = gfc_constant_result (BT_COMPLEX, kind, &src->where);
- mpf_set (result->value.complex.r, src->value.real);
- mpf_set_ui (result->value.complex.i, 0);
+ mpfr_set (result->value.complex.r, src->value.real, GFC_RND_MODE);
+ mpfr_set_ui (result->value.complex.i, 0, GFC_RND_MODE);
rc = gfc_check_real_range (result->value.complex.r, kind);
{
if (gfc_option.warn_underflow)
gfc_warning ("%s at %L", gfc_arith_error (rc), &src->where);
- mpf_set_ui(result->value.complex.r, 0);
+ mpfr_set_ui(result->value.complex.r, 0, GFC_RND_MODE);
}
else if (rc != ARITH_OK)
{
result = gfc_constant_result (BT_INTEGER, kind, &src->where);
- mpz_set_f (result->value.integer, src->value.complex.r);
+ gfc_mpfr_to_mpz(result->value.integer, src->value.complex.r);
if ((rc = gfc_check_integer_range (result->value.integer, kind))
!= ARITH_OK)
result = gfc_constant_result (BT_REAL, kind, &src->where);
- mpf_set (result->value.real, src->value.complex.r);
+ mpfr_set (result->value.real, src->value.complex.r, GFC_RND_MODE);
rc = gfc_check_real_range (result->value.real, kind);
{
if (gfc_option.warn_underflow)
gfc_warning ("%s at %L", gfc_arith_error (rc), &src->where);
- mpf_set_ui(result->value.real, 0);
+ mpfr_set_ui(result->value.real, 0, GFC_RND_MODE);
}
if (rc != ARITH_OK)
{
result = gfc_constant_result (BT_COMPLEX, kind, &src->where);
- mpf_set (result->value.complex.r, src->value.complex.r);
- mpf_set (result->value.complex.i, src->value.complex.i);
+ mpfr_set (result->value.complex.r, src->value.complex.r, GFC_RND_MODE);
+ mpfr_set (result->value.complex.i, src->value.complex.i, GFC_RND_MODE);
rc = gfc_check_real_range (result->value.complex.r, kind);
{
if (gfc_option.warn_underflow)
gfc_warning ("%s at %L", gfc_arith_error (rc), &src->where);
- mpf_set_ui(result->value.complex.r, 0);
+ mpfr_set_ui(result->value.complex.r, 0, GFC_RND_MODE);
}
else if (rc != ARITH_OK)
{
{
if (gfc_option.warn_underflow)
gfc_warning ("%s at %L", gfc_arith_error (rc), &src->where);
- mpf_set_ui(result->value.complex.i, 0);
+ mpfr_set_ui(result->value.complex.i, 0, GFC_RND_MODE);
}
else if (rc != ARITH_OK)
{
#include "gfortran.h"
-/* Constants calculated during initialization. */
-extern mpf_t pi, half_pi, two_pi, e;
-
-/* Calculate mathematically interesting functions. */
-void natural_logarithm (mpf_t *, mpf_t *);
-void common_logarithm (mpf_t *, mpf_t *);
-void exponential (mpf_t *, mpf_t *);
-void sine (mpf_t *, mpf_t *);
-void cosine (mpf_t *, mpf_t *);
-void arctangent (mpf_t *, mpf_t *);
-void arctangent2 (mpf_t *, mpf_t *, mpf_t *);
-void hypercos (mpf_t *, mpf_t *);
-void hypersine (mpf_t *, mpf_t *);
+/* MPFR does not have mpfr_atan2(), which needs to return the principle
+ value of atan2(). MPFR also does not have the conversion of a mpfr_t
+ to a mpz_t, so declare a function for this as well. */
+
+void arctangent2 (mpfr_t, mpfr_t, mpfr_t);
+void gfc_mpfr_to_mpz(mpz_t, mpfr_t);
+void gfc_set_model_kind (int);
+void gfc_set_model (mpfr_t);
/* Return a constant result of a given type and kind, with locus. */
gfc_expr *gfc_constant_result (bt, int, locus *);
break;
case BT_REAL:
- mpf_out_str (stdout, 10, 0, p->value.real);
+ mpfr_out_str (stdout, 10, 0, p->value.real, GFC_RND_MODE);
if (p->ts.kind != gfc_default_real_kind ())
gfc_status ("_%d", p->ts.kind);
break;
case BT_COMPLEX:
gfc_status ("(complex ");
- mpf_out_str (stdout, 10, 0, p->value.complex.r);
+ mpfr_out_str (stdout, 10, 0, p->value.complex.r, GFC_RND_MODE);
if (p->ts.kind != gfc_default_complex_kind ())
gfc_status ("_%d", p->ts.kind);
gfc_status (" ");
- mpf_out_str (stdout, 10, 0, p->value.complex.i);
+ mpfr_out_str (stdout, 10, 0, p->value.complex.i, GFC_RND_MODE);
if (p->ts.kind != gfc_default_complex_kind ())
gfc_status ("_%d", p->ts.kind);
break;
case BT_REAL:
- mpf_clear (e->value.real);
+ mpfr_clear (e->value.real);
break;
case BT_CHARACTER:
break;
case BT_COMPLEX:
- mpf_clear (e->value.complex.r);
- mpf_clear (e->value.complex.i);
+ mpfr_clear (e->value.complex.r);
+ mpfr_clear (e->value.complex.i);
break;
default:
break;
case BT_REAL:
- mpf_init_set (q->value.real, p->value.real);
+ gfc_set_model_kind (q->ts.kind);
+ mpfr_init (q->value.real);
+ mpfr_set (q->value.real, p->value.real, GFC_RND_MODE);
break;
case BT_COMPLEX:
- mpf_init_set (q->value.complex.r, p->value.complex.r);
- mpf_init_set (q->value.complex.i, p->value.complex.i);
+ gfc_set_model_kind (q->ts.kind);
+ mpfr_init (q->value.complex.r);
+ mpfr_init (q->value.complex.i);
+ mpfr_set (q->value.complex.r, p->value.complex.r, GFC_RND_MODE);
+ mpfr_set (q->value.complex.i, p->value.complex.i, GFC_RND_MODE);
break;
case BT_CHARACTER:
/* Major control parameters. */
#define GFC_MAX_SYMBOL_LEN 63
-#define GFC_REAL_BITS 100 /* Number of bits in g95's floating point numbers. */
#define GFC_MAX_LINE 132 /* Characters beyond this are not seen. */
#define GFC_MAX_DIMENSIONS 7 /* Maximum dimensions in an array. */
#define GFC_LETTERS 26 /* Number of letters in the alphabet. */
/* Arithmetic results. */
typedef enum
-{ ARITH_OK = 1, ARITH_OVERFLOW, ARITH_UNDERFLOW,
+{ ARITH_OK = 1, ARITH_OVERFLOW, ARITH_UNDERFLOW, ARITH_NAN,
ARITH_DIV0, ARITH_0TO0, ARITH_INCOMMENSURATE
}
arith;
EXPR_ARRAY An array constructor. */
#include <gmp.h>
+#include <mpfr.h>
+#define GFC_RND_MODE GMP_RNDN
typedef struct gfc_expr
{
union
{
- mpz_t integer;
- mpf_t real;
int logical;
+ mpz_t integer;
+
+ mpfr_t real;
struct
{
- mpf_t r, i;
+ mpfr_t r, i;
}
complex;
int kind, radix, digits, min_exponent, max_exponent;
int range, precision;
- mpf_t epsilon, huge, tiny;
+ mpfr_t epsilon, huge, tiny;
}
gfc_real_info;
/* simplify.c */
void gfc_simplify_init_1 (void);
-void gfc_simplify_done_1 (void);
/* match.c -- FIXME */
void gfc_free_iterator (gfc_iterator *, int);
gfc_check_rand, NULL, NULL,
i, BT_INTEGER, 4, 0);
+ /* Compatibility with HP FORTRAN 77/iX Reference. Note, rand() and
+ ran() use slightly different shoddy multiplicative congruential
+ PRNG. */
+ make_alias ("ran");
+
make_generic ("rand", GFC_ISYM_RAND);
add_sym_1 ("range", 0, 1, BT_INTEGER, di,
gfc_scanner_done_1 ();
gfc_intrinsic_done_1 ();
- gfc_simplify_done_1 ();
gfc_iresolve_done_1 ();
gfc_arith_done_1 ();
}
#include <time.h>
#include "gfortran.h"
+#include "arith.h"
#include "match.h"
#include "parse.h" /* FIXME */
tail->next = new;
tail = new;
- /* See what kind of interface we're dealing with. Asusume it is
+ /* See what kind of interface we're dealing with. Assume it is
not an operator. */
new->operator = INTRINSIC_NONE;
if (gfc_match_generic_spec (&type, name, &operator) == MATCH_ERROR)
static void
-mio_gmp_real (mpf_t * real)
+mio_gmp_real (mpfr_t * real)
{
mp_exp_t exponent;
char *p;
if (parse_atom () != ATOM_STRING)
bad_module ("Expected real string");
- mpf_init (*real);
- mpf_set_str (*real, atom_string, -16);
+ mpfr_init (*real);
+ mpfr_set_str (*real, atom_string, 16, GFC_RND_MODE);
gfc_free (atom_string);
}
else
{
- p = mpf_get_str (NULL, &exponent, 16, 0, *real);
+ p = mpfr_get_str (NULL, &exponent, 16, 0, *real, GFC_RND_MODE);
atom_string = gfc_getmem (strlen (p) + 20);
sprintf (atom_string, "0.%s@%ld", p, exponent);
break;
case BT_REAL:
+ gfc_set_model_kind (e->ts.kind);
mio_gmp_real (&e->value.real);
break;
case BT_COMPLEX:
+ gfc_set_model_kind (e->ts.kind);
mio_gmp_real (&e->value.complex.r);
mio_gmp_real (&e->value.complex.i);
break;
/* Primary expression subroutines
- Copyright (C) 2000, 2001, 2002 Free Software Foundation, Inc.
+ Copyright (C) 2000, 2001, 2002, 2004 Free Software Foundation, Inc.
Contributed by Andy Vaught
This file is part of GNU G95.
buffer = alloca (count + 1);
memset (buffer, '\0', count + 1);
- /* Hack for mpf_init_set_str(). */
+ /* Hack for mpfr_set_str(). */
p = buffer;
while (count > 0)
{
case ARITH_UNDERFLOW:
if (gfc_option.warn_underflow)
gfc_warning ("Real constant underflows its kind at %C");
- mpf_set_ui(e->value.real, 0);
+ mpfr_set_ui (e->value.real, 0, GFC_RND_MODE);
break;
default:
buffer = alloca (count + 1);
memset (buffer, '\0', count + 1);
- /* Hack for mpf_init_set_str(). */
+ /* Hack for mpfr_set_str(). */
p = buffer;
while (count > 0)
{
c = gfc_next_char ();
- if (c == 'd')
+ if (c == 'd' || c == 'q')
c = 'e';
*p++ = c;
count--;
#include "arith.h"
#include "intrinsic.h"
-static mpf_t mpf_zero, mpf_half, mpf_one;
-static mpz_t mpz_zero;
-
gfc_expr gfc_bad_expr;
gfc_simplify_abs (gfc_expr * e)
{
gfc_expr *result;
- mpf_t a, b;
+ mpfr_t a, b;
if (e->expr_type != EXPR_CONSTANT)
return NULL;
case BT_REAL:
result = gfc_constant_result (BT_REAL, e->ts.kind, &e->where);
- mpf_abs (result->value.real, e->value.real);
+ mpfr_abs (result->value.real, e->value.real, GFC_RND_MODE);
result = range_check (result, "ABS");
break;
case BT_COMPLEX:
result = gfc_constant_result (BT_REAL, e->ts.kind, &e->where);
- mpf_init (a);
- mpf_mul (a, e->value.complex.r, e->value.complex.r);
-
- mpf_init (b);
- mpf_mul (b, e->value.complex.i, e->value.complex.i);
-
- mpf_add (a, a, b);
- mpf_sqrt (result->value.real, a);
+ gfc_set_model_kind (e->ts.kind);
+ mpfr_init (a);
+ mpfr_init (b);
+ /* FIXME: Possible numerical problems. */
+ mpfr_mul (a, e->value.complex.r, e->value.complex.r, GFC_RND_MODE);
+ mpfr_mul (b, e->value.complex.i, e->value.complex.i, GFC_RND_MODE);
+ mpfr_add (a, a, b, GFC_RND_MODE);
+ mpfr_sqrt (result->value.real, a, GFC_RND_MODE);
- mpf_clear (a);
- mpf_clear (b);
+ mpfr_clear (a);
+ mpfr_clear (b);
result = range_check (result, "CABS");
break;
gfc_simplify_acos (gfc_expr * x)
{
gfc_expr *result;
- mpf_t negative, square, term;
if (x->expr_type != EXPR_CONSTANT)
return NULL;
- if (mpf_cmp_si (x->value.real, 1) > 0 || mpf_cmp_si (x->value.real, -1) < 0)
+ if (mpfr_cmp_si (x->value.real, 1) > 0 || mpfr_cmp_si (x->value.real, -1) < 0)
{
gfc_error ("Argument of ACOS at %L must be between -1 and 1",
&x->where);
result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
- if (mpf_cmp_si (x->value.real, 1) == 0)
- {
- mpf_set_ui (result->value.real, 0);
- return range_check (result, "ACOS");
- }
-
- if (mpf_cmp_si (x->value.real, -1) == 0)
- {
- mpf_set (result->value.real, pi);
- return range_check (result, "ACOS");
- }
-
- mpf_init (negative);
- mpf_init (square);
- mpf_init (term);
-
- mpf_pow_ui (square, x->value.real, 2);
- mpf_ui_sub (term, 1, square);
- mpf_sqrt (term, term);
- mpf_div (term, x->value.real, term);
- mpf_neg (term, term);
- arctangent (&term, &negative);
- mpf_add (result->value.real, half_pi, negative);
-
- mpf_clear (negative);
- mpf_clear (square);
- mpf_clear (term);
+ mpfr_acos (result->value.real, x->value.real, GFC_RND_MODE);
return range_check (result, "ACOS");
}
return NULL;
result = gfc_constant_result (BT_REAL, e->ts.kind, &e->where);
- mpf_set (result->value.real, e->value.complex.i);
+ mpfr_set (result->value.real, e->value.complex.i, GFC_RND_MODE);
return range_check (result, "AIMAG");
}
rtrunc = gfc_copy_expr (e);
- mpf_trunc (rtrunc->value.real, e->value.real);
+ mpfr_trunc (rtrunc->value.real, e->value.real);
result = gfc_real2real (rtrunc, kind);
gfc_free_expr (rtrunc);
rtrunc = gfc_copy_expr (e);
- mpf_trunc (rtrunc->value.real, e->value.real);
+ mpfr_trunc (rtrunc->value.real, e->value.real);
result = gfc_real2real (rtrunc, gfc_default_double_kind ());
gfc_free_expr (rtrunc);
{
gfc_expr *rtrunc, *result;
int kind, cmp;
+ mpfr_t half;
kind = get_kind (BT_REAL, k, "ANINT", e->ts.kind);
if (kind == -1)
rtrunc = gfc_copy_expr (e);
- cmp = mpf_cmp_ui (e->value.real, 0);
+ cmp = mpfr_cmp_ui (e->value.real, 0);
+
+ gfc_set_model_kind (kind);
+ mpfr_init (half);
+ mpfr_set_str (half, "0.5", 10, GFC_RND_MODE);
if (cmp > 0)
{
- mpf_add (rtrunc->value.real, e->value.real, mpf_half);
- mpf_trunc (result->value.real, rtrunc->value.real);
+ mpfr_add (rtrunc->value.real, e->value.real, half, GFC_RND_MODE);
+ mpfr_trunc (result->value.real, rtrunc->value.real);
}
else if (cmp < 0)
{
- mpf_sub (rtrunc->value.real, e->value.real, mpf_half);
- mpf_trunc (result->value.real, rtrunc->value.real);
+ mpfr_sub (rtrunc->value.real, e->value.real, half, GFC_RND_MODE);
+ mpfr_trunc (result->value.real, rtrunc->value.real);
}
else
- mpf_set_ui (result->value.real, 0);
+ mpfr_set_ui (result->value.real, 0, GFC_RND_MODE);
gfc_free_expr (rtrunc);
+ mpfr_clear (half);
return range_check (result, "ANINT");
}
{
gfc_expr *rtrunc, *result;
int cmp;
+ mpfr_t half;
if (e->expr_type != EXPR_CONSTANT)
return NULL;
rtrunc = gfc_copy_expr (e);
- cmp = mpf_cmp_ui (e->value.real, 0);
+ cmp = mpfr_cmp_ui (e->value.real, 0);
+
+ gfc_set_model_kind (gfc_default_double_kind ());
+ mpfr_init (half);
+ mpfr_set_str (half, "0.5", 10, GFC_RND_MODE);
if (cmp > 0)
{
- mpf_add (rtrunc->value.real, e->value.real, mpf_half);
- mpf_trunc (result->value.real, rtrunc->value.real);
+ mpfr_add (rtrunc->value.real, e->value.real, half, GFC_RND_MODE);
+ mpfr_trunc (result->value.real, rtrunc->value.real);
}
else if (cmp < 0)
{
- mpf_sub (rtrunc->value.real, e->value.real, mpf_half);
- mpf_trunc (result->value.real, rtrunc->value.real);
+ mpfr_sub (rtrunc->value.real, e->value.real, half, GFC_RND_MODE);
+ mpfr_trunc (result->value.real, rtrunc->value.real);
}
else
- mpf_set_ui (result->value.real, 0);
+ mpfr_set_ui (result->value.real, 0, GFC_RND_MODE);
gfc_free_expr (rtrunc);
+ mpfr_clear (half);
return range_check (result, "DNINT");
}
gfc_simplify_asin (gfc_expr * x)
{
gfc_expr *result;
- mpf_t negative, square, term;
if (x->expr_type != EXPR_CONSTANT)
return NULL;
- if (mpf_cmp_si (x->value.real, 1) > 0 || mpf_cmp_si (x->value.real, -1) < 0)
+ if (mpfr_cmp_si (x->value.real, 1) > 0 || mpfr_cmp_si (x->value.real, -1) < 0)
{
gfc_error ("Argument of ASIN at %L must be between -1 and 1",
&x->where);
result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
- if (mpf_cmp_si (x->value.real, 1) == 0)
- {
- mpf_set (result->value.real, half_pi);
- return range_check (result, "ASIN");
- }
-
- if (mpf_cmp_si (x->value.real, -1) == 0)
- {
- mpf_init (negative);
- mpf_neg (negative, half_pi);
- mpf_set (result->value.real, negative);
- mpf_clear (negative);
- return range_check (result, "ASIN");
- }
-
- mpf_init (square);
- mpf_init (term);
-
- mpf_pow_ui (square, x->value.real, 2);
- mpf_ui_sub (term, 1, square);
- mpf_sqrt (term, term);
- mpf_div (term, x->value.real, term);
- arctangent (&term, &result->value.real);
-
- mpf_clear (square);
- mpf_clear (term);
+ mpfr_asin(result->value.real, x->value.real, GFC_RND_MODE);
return range_check (result, "ASIN");
}
result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
- arctangent (&x->value.real, &result->value.real);
+ mpfr_atan(result->value.real, x->value.real, GFC_RND_MODE);
return range_check (result, "ATAN");
result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
-
- if (mpf_sgn (y->value.real) == 0 && mpf_sgn (x->value.real) == 0)
+ if (mpfr_sgn (y->value.real) == 0 && mpfr_sgn (x->value.real) == 0)
{
gfc_error
- ("If first argument of ATAN2 %L is zero, the second argument "
+ ("If first argument of ATAN2 %L is zero, then the second argument "
"must not be zero", &x->where);
gfc_free_expr (result);
return &gfc_bad_expr;
}
- arctangent2 (&y->value.real, &x->value.real, &result->value.real);
+ arctangent2 (y->value.real, x->value.real, result->value.real);
return range_check (result, "ATAN2");
ceil = gfc_copy_expr (e);
- mpf_ceil (ceil->value.real, e->value.real);
- mpz_set_f (result->value.integer, ceil->value.real);
+ mpfr_ceil (ceil->value.real, e->value.real);
+ gfc_mpfr_to_mpz(result->value.integer, ceil->value.real);
gfc_free_expr (ceil);
result = gfc_constant_result (BT_COMPLEX, kind, &x->where);
- mpf_set_ui (result->value.complex.i, 0);
+ mpfr_set_ui (result->value.complex.i, 0, GFC_RND_MODE);
switch (x->ts.type)
{
case BT_INTEGER:
- mpf_set_z (result->value.complex.r, x->value.integer);
+ mpfr_set_z (result->value.complex.r, x->value.integer, GFC_RND_MODE);
break;
case BT_REAL:
- mpf_set (result->value.complex.r, x->value.real);
+ mpfr_set (result->value.complex.r, x->value.real, GFC_RND_MODE);
break;
case BT_COMPLEX:
- mpf_set (result->value.complex.r, x->value.complex.r);
- mpf_set (result->value.complex.i, x->value.complex.i);
+ mpfr_set (result->value.complex.r, x->value.complex.r, GFC_RND_MODE);
+ mpfr_set (result->value.complex.i, x->value.complex.i, GFC_RND_MODE);
break;
default:
switch (y->ts.type)
{
case BT_INTEGER:
- mpf_set_z (result->value.complex.i, y->value.integer);
+ mpfr_set_z (result->value.complex.i, y->value.integer, GFC_RND_MODE);
break;
case BT_REAL:
- mpf_set (result->value.complex.i, y->value.real);
+ mpfr_set (result->value.complex.i, y->value.real, GFC_RND_MODE);
break;
default:
return NULL;
result = gfc_copy_expr (e);
- mpf_neg (result->value.complex.i, result->value.complex.i);
+ mpfr_neg (result->value.complex.i, result->value.complex.i, GFC_RND_MODE);
return range_check (result, "CONJG");
}
gfc_simplify_cos (gfc_expr * x)
{
gfc_expr *result;
- mpf_t xp, xq;
+ mpfr_t xp, xq;
if (x->expr_type != EXPR_CONSTANT)
return NULL;
switch (x->ts.type)
{
case BT_REAL:
- cosine (&x->value.real, &result->value.real);
+ mpfr_cos (result->value.real, x->value.real, GFC_RND_MODE);
break;
case BT_COMPLEX:
- mpf_init (xp);
- mpf_init (xq);
+ gfc_set_model_kind (x->ts.kind);
+ mpfr_init (xp);
+ mpfr_init (xq);
- cosine (&x->value.complex.r, &xp);
- hypercos (&x->value.complex.i, &xq);
- mpf_mul (result->value.complex.r, xp, xq);
+ mpfr_cos (xp, x->value.complex.r, GFC_RND_MODE);
+ mpfr_cosh (xq, x->value.complex.i, GFC_RND_MODE);
+ mpfr_mul(result->value.complex.r, xp, xq, GFC_RND_MODE);
- sine (&x->value.complex.r, &xp);
- hypersine (&x->value.complex.i, &xq);
- mpf_mul (xp, xp, xq);
- mpf_neg (result->value.complex.i, xp);
+ mpfr_sin (xp, x->value.complex.r, GFC_RND_MODE);
+ mpfr_sinh (xq, x->value.complex.i, GFC_RND_MODE);
+ mpfr_mul (xp, xp, xq, GFC_RND_MODE);
+ mpfr_neg (result->value.complex.i, xp, GFC_RND_MODE );
- mpf_clear (xp);
- mpf_clear (xq);
+ mpfr_clear (xp);
+ mpfr_clear (xq);
break;
default:
gfc_internal_error ("in gfc_simplify_cos(): Bad type");
result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
- hypercos (&x->value.real, &result->value.real);
+ mpfr_cosh (result->value.real, x->value.real, GFC_RND_MODE);
return range_check (result, "COSH");
}
if (mpz_cmp (x->value.integer, y->value.integer) > 0)
mpz_sub (result->value.integer, x->value.integer, y->value.integer);
else
- mpz_set (result->value.integer, mpz_zero);
+ mpz_set_ui (result->value.integer, 0);
break;
case BT_REAL:
- if (mpf_cmp (x->value.real, y->value.real) > 0)
- mpf_sub (result->value.real, x->value.real, y->value.real);
+ if (mpfr_cmp (x->value.real, y->value.real) > 0)
+ mpfr_sub (result->value.real, x->value.real, y->value.real, GFC_RND_MODE);
else
- mpf_set (result->value.real, mpf_zero);
+ mpfr_set_ui (result->value.real, 0, GFC_RND_MODE);
break;
gfc_expr *
gfc_simplify_dprod (gfc_expr * x, gfc_expr * y)
{
- gfc_expr *mult1, *mult2, *result;
+ gfc_expr *a1, *a2, *result;
if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
return NULL;
result =
gfc_constant_result (BT_REAL, gfc_default_double_kind (), &x->where);
- mult1 = gfc_real2real (x, gfc_default_double_kind ());
- mult2 = gfc_real2real (y, gfc_default_double_kind ());
+ a1 = gfc_real2real (x, gfc_default_double_kind ());
+ a2 = gfc_real2real (y, gfc_default_double_kind ());
- mpf_mul (result->value.real, mult1->value.real, mult2->value.real);
+ mpfr_mul (result->value.real, a1->value.real, a2->value.real, GFC_RND_MODE);
- gfc_free_expr (mult1);
- gfc_free_expr (mult2);
+ gfc_free_expr (a1);
+ gfc_free_expr (a2);
return range_check (result, "DPROD");
}
result = gfc_constant_result (BT_REAL, e->ts.kind, &e->where);
- mpf_set (result->value.real, gfc_real_kinds[i].epsilon);
+ mpfr_set (result->value.real, gfc_real_kinds[i].epsilon, GFC_RND_MODE);
return range_check (result, "EPSILON");
}
gfc_simplify_exp (gfc_expr * x)
{
gfc_expr *result;
- mpf_t xp, xq;
- double ln2, absval, rhuge;
+ mpfr_t xp, xq;
if (x->expr_type != EXPR_CONSTANT)
return NULL;
result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
- /* Exactitude doesn't matter here */
- ln2 = .6931472;
- rhuge = ln2 * mpz_get_d (gfc_integer_kinds[0].huge);
-
switch (x->ts.type)
{
case BT_REAL:
- absval = mpf_get_d (x->value.real);
- if (absval < 0)
- absval = -absval;
- if (absval > rhuge)
- {
- /* Underflow (set arg to zero) if x is negative and its
- magnitude is greater than the maximum C long int times
- ln2, because the exponential method in arith.c will fail
- for such values. */
- if (mpf_cmp_ui (x->value.real, 0) < 0)
- {
- if (pedantic == 1)
- gfc_warning_now
- ("Argument of EXP at %L is negative and too large, "
- "setting result to zero", &x->where);
- mpf_set_ui (result->value.real, 0);
- return range_check (result, "EXP");
- }
- /* Overflow if magnitude of x is greater than C long int
- huge times ln2. */
- else
- {
- gfc_error ("Argument of EXP at %L too large", &x->where);
- gfc_free_expr (result);
- return &gfc_bad_expr;
- }
- }
- exponential (&x->value.real, &result->value.real);
+ mpfr_exp(result->value.real, x->value.real, GFC_RND_MODE);
break;
case BT_COMPLEX:
- /* Using Euler's formula. */
- absval = mpf_get_d (x->value.complex.r);
- if (absval < 0)
- absval = -absval;
- if (absval > rhuge)
- {
- if (mpf_cmp_ui (x->value.complex.r, 0) < 0)
- {
- if (pedantic == 1)
- gfc_warning_now
- ("Real part of argument of EXP at %L is negative "
- "and too large, setting result to zero", &x->where);
-
- mpf_set_ui (result->value.complex.r, 0);
- mpf_set_ui (result->value.complex.i, 0);
- return range_check (result, "EXP");
- }
- else
- {
- gfc_error ("Real part of argument of EXP at %L too large",
- &x->where);
- gfc_free_expr (result);
- return &gfc_bad_expr;
- }
- }
- mpf_init (xp);
- mpf_init (xq);
- exponential (&x->value.complex.r, &xq);
- cosine (&x->value.complex.i, &xp);
- mpf_mul (result->value.complex.r, xq, xp);
- sine (&x->value.complex.i, &xp);
- mpf_mul (result->value.complex.i, xq, xp);
- mpf_clear (xp);
- mpf_clear (xq);
+ gfc_set_model_kind (x->ts.kind);
+ mpfr_init (xp);
+ mpfr_init (xq);
+ mpfr_exp (xq, x->value.complex.r, GFC_RND_MODE);
+ mpfr_cos (xp, x->value.complex.i, GFC_RND_MODE);
+ mpfr_mul (result->value.complex.r, xq, xp, GFC_RND_MODE);
+ mpfr_sin (xp, x->value.complex.i, GFC_RND_MODE);
+ mpfr_mul (result->value.complex.i, xq, xp, GFC_RND_MODE);
+ mpfr_clear (xp);
+ mpfr_clear (xq);
break;
default:
return range_check (result, "EXP");
}
-
+/* FIXME: MPFR should be able to do this better */
gfc_expr *
gfc_simplify_exponent (gfc_expr * x)
{
- mpf_t i2, absv, ln2, lnx;
+ mpfr_t i2, absv, ln2, lnx, zero;
gfc_expr *result;
if (x->expr_type != EXPR_CONSTANT)
result = gfc_constant_result (BT_INTEGER, gfc_default_integer_kind (),
&x->where);
- if (mpf_cmp (x->value.real, mpf_zero) == 0)
+ gfc_set_model (x->value.real);
+ mpfr_init (zero);
+ mpfr_set_ui (zero, 0, GFC_RND_MODE);
+
+ if (mpfr_cmp (x->value.real, zero) == 0)
{
mpz_set_ui (result->value.integer, 0);
+ mpfr_clear (zero);
return result;
}
- mpf_init_set_ui (i2, 2);
- mpf_init (absv);
- mpf_init (ln2);
- mpf_init (lnx);
+ mpfr_init (i2);
+ mpfr_init (absv);
+ mpfr_init (ln2);
+ mpfr_init (lnx);
+
+ mpfr_set_ui (i2, 2, GFC_RND_MODE);
- natural_logarithm (&i2, &ln2);
+ mpfr_log (ln2, i2, GFC_RND_MODE);
+ mpfr_abs (absv, x->value.real, GFC_RND_MODE);
+ mpfr_log (lnx, absv, GFC_RND_MODE);
- mpf_abs (absv, x->value.real);
- natural_logarithm (&absv, &lnx);
+ mpfr_div (lnx, lnx, ln2, GFC_RND_MODE);
+ mpfr_trunc (lnx, lnx);
+ mpfr_add_ui (lnx, lnx, 1, GFC_RND_MODE);
- mpf_div (lnx, lnx, ln2);
- mpf_trunc (lnx, lnx);
- mpf_add_ui (lnx, lnx, 1);
- mpz_set_f (result->value.integer, lnx);
+ gfc_mpfr_to_mpz (result->value.integer, lnx);
- mpf_clear (i2);
- mpf_clear (ln2);
- mpf_clear (lnx);
- mpf_clear (absv);
+ mpfr_clear (i2);
+ mpfr_clear (ln2);
+ mpfr_clear (lnx);
+ mpfr_clear (absv);
+ mpfr_clear (zero);
return range_check (result, "EXPONENT");
}
gfc_simplify_floor (gfc_expr * e, gfc_expr * k)
{
gfc_expr *result;
- mpf_t floor;
+ mpfr_t floor;
int kind;
kind = get_kind (BT_REAL, k, "FLOOR", gfc_default_real_kind ());
result = gfc_constant_result (BT_INTEGER, kind, &e->where);
- mpf_init (floor);
- mpf_floor (floor, e->value.real);
- mpz_set_f (result->value.integer, floor);
- mpf_clear (floor);
+ gfc_set_model_kind (kind);
+ mpfr_init (floor);
+ mpfr_floor (floor, e->value.real);
+
+ gfc_mpfr_to_mpz (result->value.integer, floor);
+
+ mpfr_clear (floor);
return range_check (result, "FLOOR");
}
gfc_simplify_fraction (gfc_expr * x)
{
gfc_expr *result;
- mpf_t i2, absv, ln2, lnx, pow2;
+ mpfr_t i2, absv, ln2, lnx, pow2, zero;
unsigned long exp2;
if (x->expr_type != EXPR_CONSTANT)
result = gfc_constant_result (BT_REAL, x->ts.kind, &x->where);
- if (mpf_cmp (x->value.real, mpf_zero) == 0)
+ gfc_set_model_kind (x->ts.kind);
+ mpfr_init (zero);
+ mpfr_set_ui (zero, 0, GFC_RND_MODE);
+
+ if (mpfr_cmp (x->value.real, zero) == 0)
{
- mpf_set (result->value.real, mpf_zero);
+ mpfr_set (result->value.real, zero, GFC_RND_MODE);
+ mpfr_clear (zero);
return result;
}
- mpf_init_set_ui (i2, 2);
- mpf_init (absv);
- mpf_init (ln2);
- mpf_init (lnx);
- mpf_init (pow2);
+ mpfr_init (i2);
+ mpfr_init (absv);
+ mpfr_init (ln2);
+ mpfr_init (lnx);
+ mpfr_init (pow2);
- natural_logarithm (&i2, &ln2);
+ mpfr_set_ui (i2, 2, GFC_RND_MODE);
- mpf_abs (absv, x->value.real);
- natural_logarithm (&absv, &lnx);
+ mpfr_log (ln2, i2, GFC_RND_MODE);
+ mpfr_abs (absv, x->value.real, GFC_RND_MODE);
+ mpfr_log (lnx, absv, GFC_RND_MODE);
- mpf_div (lnx, lnx, ln2);
- mpf_trunc (lnx, lnx);
- mpf_add_ui (lnx, lnx, 1);
+ mpfr_div (lnx, lnx, ln2, GFC_RND_MODE);
+ mpfr_trunc (lnx, lnx);
+ mpfr_add_ui (lnx, lnx, 1, GFC_RND_MODE);
- exp2 = (unsigned long) mpf_get_d (lnx);
- mpf_pow_ui (pow2, i2, exp2);
+ exp2 = (unsigned long) mpfr_get_d (lnx, GFC_RND_MODE);
+ mpfr_pow_ui (pow2, i2, exp2, GFC_RND_MODE);
- mpf_div (result->value.real, absv, pow2);
+ mpfr_div (result->value.real, absv, pow2, GFC_RND_MODE);
- mpf_clear (i2);
- mpf_clear (ln2);
- mpf_clear (absv);
- mpf_clear (lnx);
- mpf_clear (pow2);
+ mpfr_clear (i2);
+ mpfr_clear (ln2);
+ mpfr_clear (absv);
+ mpfr_clear (lnx);
+ mpfr_clear (pow2);
+ mpfr_clear (zero);
return range_check (result, "FRACTION");
}
break;
case BT_REAL:
- mpf_set (result->value.real, gfc_real_kinds[i].huge);
+ mpfr_set (result->value.real, gfc_real_kinds[i].huge, GFC_RND_MODE);
break;
bad_type:
case BT_REAL:
rtrunc = gfc_copy_expr (e);
- mpf_trunc (rtrunc->value.real, e->value.real);
- mpz_set_f (result->value.integer, rtrunc->value.real);
+ mpfr_trunc (rtrunc->value.real, e->value.real);
+ gfc_mpfr_to_mpz (result->value.integer, rtrunc->value.real);
gfc_free_expr (rtrunc);
break;
case BT_COMPLEX:
rpart = gfc_complex2real (e, kind);
rtrunc = gfc_copy_expr (rpart);
- mpf_trunc (rtrunc->value.real, rpart->value.real);
- mpz_set_f (result->value.integer, rtrunc->value.real);
+ mpfr_trunc (rtrunc->value.real, rpart->value.real);
+ gfc_mpfr_to_mpz (result->value.integer, rtrunc->value.real);
gfc_free_expr (rpart);
gfc_free_expr (rtrunc);
break;
rtrunc = gfc_copy_expr (e);
- mpf_trunc (rtrunc->value.real, e->value.real);
- mpz_set_f (result->value.integer, rtrunc->value.real);
+ mpfr_trunc (rtrunc->value.real, e->value.real);
+ gfc_mpfr_to_mpz (result->value.integer, rtrunc->value.real);
gfc_free_expr (rtrunc);
return range_check (result, "IFIX");
rtrunc = gfc_copy_expr (e);
- mpf_trunc (rtrunc->value.real, e->value.real);
- mpz_set_f (result->value.integer, rtrunc->value.real);
+ mpfr_trunc (rtrunc->value.real, e->value.real);
+ gfc_mpfr_to_mpz (result->value.integer, rtrunc->value.real);
gfc_free_expr (rtrunc);
return range_check (result, "IDINT");
gfc_simplify_log (gfc_expr * x)
{
gfc_expr *result;
- mpf_t xr, xi;
+ mpfr_t xr, xi, zero;
if (x->expr_type != EXPR_CONSTANT)
return NULL;
result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
+ gfc_set_model_kind (x->ts.kind);
+ mpfr_init (zero);
+ mpfr_set_ui (zero, 0, GFC_RND_MODE);
+
switch (x->ts.type)
{
case BT_REAL:
- if (mpf_cmp (x->value.real, mpf_zero) <= 0)
+ if (mpfr_cmp (x->value.real, zero) <= 0)
{
gfc_error
("Argument of LOG at %L cannot be less than or equal to zero",
&x->where);
gfc_free_expr (result);
+ mpfr_clear (zero);
return &gfc_bad_expr;
}
- natural_logarithm (&x->value.real, &result->value.real);
+ mpfr_log(result->value.real, x->value.real, GFC_RND_MODE);
+ mpfr_clear (zero);
break;
case BT_COMPLEX:
- if ((mpf_cmp (x->value.complex.r, mpf_zero) == 0)
- && (mpf_cmp (x->value.complex.i, mpf_zero) == 0))
+ if ((mpfr_cmp (x->value.complex.r, zero) == 0)
+ && (mpfr_cmp (x->value.complex.i, zero) == 0))
{
gfc_error ("Complex argument of LOG at %L cannot be zero",
&x->where);
gfc_free_expr (result);
+ mpfr_clear (zero);
return &gfc_bad_expr;
}
- mpf_init (xr);
- mpf_init (xi);
+ mpfr_init (xr);
+ mpfr_init (xi);
- arctangent2 (&x->value.complex.i, &x->value.complex.r,
- &result->value.complex.i);
+ arctangent2 (x->value.complex.i, x->value.complex.r,
+ result->value.complex.i);
- mpf_mul (xr, x->value.complex.r, x->value.complex.r);
- mpf_mul (xi, x->value.complex.i, x->value.complex.i);
- mpf_add (xr, xr, xi);
- mpf_sqrt (xr, xr);
- natural_logarithm (&xr, &result->value.complex.r);
+ mpfr_mul (xr, x->value.complex.r, x->value.complex.r, GFC_RND_MODE);
+ mpfr_mul (xi, x->value.complex.i, x->value.complex.i, GFC_RND_MODE);
+ mpfr_add (xr, xr, xi, GFC_RND_MODE);
+ mpfr_sqrt (xr, xr, GFC_RND_MODE);
+ mpfr_log (result->value.complex.r, xr, GFC_RND_MODE);
- mpf_clear (xr);
- mpf_clear (xi);
+ mpfr_clear (xr);
+ mpfr_clear (xi);
+ mpfr_clear (zero);
break;
gfc_simplify_log10 (gfc_expr * x)
{
gfc_expr *result;
+ mpfr_t zero;
if (x->expr_type != EXPR_CONSTANT)
return NULL;
- if (mpf_cmp (x->value.real, mpf_zero) <= 0)
+ gfc_set_model_kind (x->ts.kind);
+ mpfr_init (zero);
+ mpfr_set_ui (zero, 0, GFC_RND_MODE);
+
+ if (mpfr_cmp (x->value.real, zero) <= 0)
{
gfc_error
("Argument of LOG10 at %L cannot be less than or equal to zero",
&x->where);
+ mpfr_clear (zero);
return &gfc_bad_expr;
}
result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
- common_logarithm (&x->value.real, &result->value.real);
+ mpfr_log10 (result->value.real, x->value.real, GFC_RND_MODE);
+ mpfr_clear (zero);
return range_check (result, "LOG10");
}
break;
case BT_REAL:
- if (mpf_cmp (arg->expr->value.real, extremum->expr->value.real) *
+ if (mpfr_cmp (arg->expr->value.real, extremum->expr->value.real) *
sign > 0)
- mpf_set (extremum->expr->value.real, arg->expr->value.real);
+ mpfr_set (extremum->expr->value.real, arg->expr->value.real,
+ GFC_RND_MODE);
break;
gfc_simplify_mod (gfc_expr * a, gfc_expr * p)
{
gfc_expr *result;
- mpf_t quot, iquot, term;
+ mpfr_t quot, iquot, term;
if (a->expr_type != EXPR_CONSTANT || p->expr_type != EXPR_CONSTANT)
return NULL;
break;
case BT_REAL:
- if (mpf_cmp_ui (p->value.real, 0) == 0)
+ if (mpfr_cmp_ui (p->value.real, 0) == 0)
{
/* Result is processor-dependent. */
gfc_error ("Second argument of MOD at %L is zero", &p->where);
return &gfc_bad_expr;
}
- mpf_init (quot);
- mpf_init (iquot);
- mpf_init (term);
+ gfc_set_model_kind (a->ts.kind);
+ mpfr_init (quot);
+ mpfr_init (iquot);
+ mpfr_init (term);
- mpf_div (quot, a->value.real, p->value.real);
- mpf_trunc (iquot, quot);
- mpf_mul (term, iquot, p->value.real);
- mpf_sub (result->value.real, a->value.real, term);
+ mpfr_div (quot, a->value.real, p->value.real, GFC_RND_MODE);
+ mpfr_trunc (iquot, quot);
+ mpfr_mul (term, iquot, p->value.real, GFC_RND_MODE);
+ mpfr_sub (result->value.real, a->value.real, term, GFC_RND_MODE);
- mpf_clear (quot);
- mpf_clear (iquot);
- mpf_clear (term);
+ mpfr_clear (quot);
+ mpfr_clear (iquot);
+ mpfr_clear (term);
break;
default:
gfc_simplify_modulo (gfc_expr * a, gfc_expr * p)
{
gfc_expr *result;
- mpf_t quot, iquot, term;
+ mpfr_t quot, iquot, term;
if (a->expr_type != EXPR_CONSTANT || p->expr_type != EXPR_CONSTANT)
return NULL;
break;
case BT_REAL:
- if (mpf_cmp_ui (p->value.real, 0) == 0)
+ if (mpfr_cmp_ui (p->value.real, 0) == 0)
{
/* Result is processor-dependent. */
gfc_error ("Second argument of MODULO at %L is zero", &p->where);
return &gfc_bad_expr;
}
- mpf_init (quot);
- mpf_init (iquot);
- mpf_init (term);
+ gfc_set_model_kind (a->ts.kind);
+ mpfr_init (quot);
+ mpfr_init (iquot);
+ mpfr_init (term);
- mpf_div (quot, a->value.real, p->value.real);
- mpf_floor (iquot, quot);
- mpf_mul (term, iquot, p->value.real);
+ mpfr_div (quot, a->value.real, p->value.real, GFC_RND_MODE);
+ mpfr_floor (iquot, quot);
+ mpfr_mul (term, iquot, p->value.real, GFC_RND_MODE);
- mpf_clear (quot);
- mpf_clear (iquot);
- mpf_clear (term);
+ mpfr_clear (quot);
+ mpfr_clear (iquot);
+ mpfr_clear (term);
- mpf_sub (result->value.real, a->value.real, term);
+ mpfr_sub (result->value.real, a->value.real, term, GFC_RND_MODE);
break;
default:
result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
- val = mpf_get_d (x->value.real);
+ val = mpfr_get_d (x->value.real, GFC_RND_MODE);
p = gfc_real_kinds[k].digits;
eps = 1.;
/* TODO we should make sure that 'float' matches kind 4 */
match_float = gfc_real_kinds[k].kind == 4;
- if (mpf_cmp_ui (s->value.real, 0) > 0)
+ if (mpfr_cmp_ui (s->value.real, 0) > 0)
{
if (match_float)
{
rval = (float) val;
rval = rval + eps;
- mpf_set_d (result->value.real, rval);
+ mpfr_set_d (result->value.real, rval, GFC_RND_MODE);
}
else
{
val = val + eps;
- mpf_set_d (result->value.real, val);
+ mpfr_set_d (result->value.real, val, GFC_RND_MODE);
}
}
- else if (mpf_cmp_ui (s->value.real, 0) < 0)
+ else if (mpfr_cmp_ui (s->value.real, 0) < 0)
{
if (match_float)
{
rval = (float) val;
rval = rval - eps;
- mpf_set_d (result->value.real, rval);
+ mpfr_set_d (result->value.real, rval, GFC_RND_MODE);
}
else
{
val = val - eps;
- mpf_set_d (result->value.real, val);
+ mpfr_set_d (result->value.real, val, GFC_RND_MODE);
}
}
else
{
gfc_expr *rtrunc, *itrunc, *result;
int kind, cmp;
+ mpfr_t half;
kind = get_kind (BT_INTEGER, k, name, gfc_default_integer_kind ());
if (kind == -1)
rtrunc = gfc_copy_expr (e);
itrunc = gfc_copy_expr (e);
- cmp = mpf_cmp_ui (e->value.real, 0);
+ cmp = mpfr_cmp_ui (e->value.real, 0);
+
+ gfc_set_model (e->value.real);
+ mpfr_init (half);
+ mpfr_set_str (half, "0.5", 10, GFC_RND_MODE);
if (cmp > 0)
{
- mpf_add (rtrunc->value.real, e->value.real, mpf_half);
- mpf_trunc (itrunc->value.real, rtrunc->value.real);
+ mpfr_add (rtrunc->value.real, e->value.real, half, GFC_RND_MODE);
+ mpfr_trunc (itrunc->value.real, rtrunc->value.real);
}
else if (cmp < 0)
{
- mpf_sub (rtrunc->value.real, e->value.real, mpf_half);
- mpf_trunc (itrunc->value.real, rtrunc->value.real);
+ mpfr_sub (rtrunc->value.real, e->value.real, half, GFC_RND_MODE);
+ mpfr_trunc (itrunc->value.real, rtrunc->value.real);
}
else
- mpf_set_ui (itrunc->value.real, 0);
+ mpfr_set_ui (itrunc->value.real, 0, GFC_RND_MODE);
- mpz_set_f (result->value.integer, itrunc->value.real);
+ gfc_mpfr_to_mpz (result->value.integer, itrunc->value.real);
gfc_free_expr (itrunc);
gfc_free_expr (rtrunc);
+ mpfr_clear (half);
return range_check (result, name);
}
gfc_simplify_rrspacing (gfc_expr * x)
{
gfc_expr *result;
- mpf_t i2, absv, ln2, lnx, frac, pow2;
+ mpfr_t i2, absv, ln2, lnx, frac, pow2, zero;
unsigned long exp2;
int i, p;
p = gfc_real_kinds[i].digits;
- if (mpf_cmp (x->value.real, mpf_zero) == 0)
+ gfc_set_model_kind (x->ts.kind);
+ mpfr_init (zero);
+ mpfr_set_ui (zero, 0, GFC_RND_MODE);
+
+ if (mpfr_cmp (x->value.real, zero) == 0)
{
- mpf_ui_div (result->value.real, 1, gfc_real_kinds[i].tiny);
+ mpfr_ui_div (result->value.real, 1, gfc_real_kinds[i].tiny, GFC_RND_MODE);
+ mpfr_clear (zero);
return result;
}
- mpf_init_set_ui (i2, 2);
- mpf_init (ln2);
- mpf_init (absv);
- mpf_init (lnx);
- mpf_init (frac);
- mpf_init (pow2);
+ mpfr_init (i2);
+ mpfr_init (ln2);
+ mpfr_init (absv);
+ mpfr_init (lnx);
+ mpfr_init (frac);
+ mpfr_init (pow2);
- natural_logarithm (&i2, &ln2);
+ mpfr_set_ui (i2, 2, GFC_RND_MODE);
- mpf_abs (absv, x->value.real);
- natural_logarithm (&absv, &lnx);
+ mpfr_log (ln2, i2, GFC_RND_MODE);
+ mpfr_abs (absv, x->value.real, GFC_RND_MODE);
+ mpfr_log (lnx, absv, GFC_RND_MODE);
- mpf_div (lnx, lnx, ln2);
- mpf_trunc (lnx, lnx);
- mpf_add_ui (lnx, lnx, 1);
+ mpfr_div (lnx, lnx, ln2, GFC_RND_MODE);
+ mpfr_trunc (lnx, lnx);
+ mpfr_add_ui (lnx, lnx, 1, GFC_RND_MODE);
- exp2 = (unsigned long) mpf_get_d (lnx);
- mpf_pow_ui (pow2, i2, exp2);
- mpf_div (frac, absv, pow2);
+ exp2 = (unsigned long) mpfr_get_d (lnx, GFC_RND_MODE);
+ mpfr_pow_ui (pow2, i2, exp2, GFC_RND_MODE);
+ mpfr_div (frac, absv, pow2, GFC_RND_MODE);
exp2 = (unsigned long) p;
- mpf_mul_2exp (result->value.real, frac, exp2);
+ mpfr_mul_2exp (result->value.real, frac, exp2, GFC_RND_MODE);
- mpf_clear (i2);
- mpf_clear (ln2);
- mpf_clear (absv);
- mpf_clear (lnx);
- mpf_clear (frac);
- mpf_clear (pow2);
+ mpfr_clear (i2);
+ mpfr_clear (ln2);
+ mpfr_clear (absv);
+ mpfr_clear (lnx);
+ mpfr_clear (frac);
+ mpfr_clear (pow2);
+ mpfr_clear (zero);
return range_check (result, "RRSPACING");
}
gfc_simplify_scale (gfc_expr * x, gfc_expr * i)
{
int k, neg_flag, power, exp_range;
- mpf_t scale, radix;
+ mpfr_t scale, radix;
gfc_expr *result;
if (x->expr_type != EXPR_CONSTANT || i->expr_type != EXPR_CONSTANT)
result = gfc_constant_result (BT_REAL, x->ts.kind, &x->where);
- if (mpf_sgn (x->value.real) == 0)
+ if (mpfr_sgn (x->value.real) == 0)
{
- mpf_set_ui (result->value.real, 0);
+ mpfr_set_ui (result->value.real, 0, GFC_RND_MODE);
return result;
}
power = -power;
}
- mpf_init_set_ui (radix, gfc_real_kinds[k].radix);
- mpf_init (scale);
- mpf_pow_ui (scale, radix, power);
+ gfc_set_model_kind (x->ts.kind);
+ mpfr_init (scale);
+ mpfr_init (radix);
+ mpfr_set_ui (radix, gfc_real_kinds[k].radix, GFC_RND_MODE);
+ mpfr_pow_ui (scale, radix, power, GFC_RND_MODE);
if (neg_flag)
- mpf_div (result->value.real, x->value.real, scale);
+ mpfr_div (result->value.real, x->value.real, scale, GFC_RND_MODE);
else
- mpf_mul (result->value.real, x->value.real, scale);
+ mpfr_mul (result->value.real, x->value.real, scale, GFC_RND_MODE);
- mpf_clear (scale);
- mpf_clear (radix);
+ mpfr_clear (scale);
+ mpfr_clear (radix);
return range_check (result, "SCALE");
}
gfc_simplify_set_exponent (gfc_expr * x, gfc_expr * i)
{
gfc_expr *result;
- mpf_t i2, ln2, absv, lnx, pow2, frac;
+ mpfr_t i2, ln2, absv, lnx, pow2, frac, zero;
unsigned long exp2;
if (x->expr_type != EXPR_CONSTANT || i->expr_type != EXPR_CONSTANT)
result = gfc_constant_result (BT_REAL, x->ts.kind, &x->where);
- if (mpf_cmp (x->value.real, mpf_zero) == 0)
+ gfc_set_model_kind (x->ts.kind);
+ mpfr_init (zero);
+ mpfr_set_ui (zero, 0, GFC_RND_MODE);
+
+ if (mpfr_cmp (x->value.real, zero) == 0)
{
- mpf_set (result->value.real, mpf_zero);
+ mpfr_set (result->value.real, zero, GFC_RND_MODE);
+ mpfr_clear (zero);
return result;
}
- mpf_init_set_ui (i2, 2);
- mpf_init (ln2);
- mpf_init (absv);
- mpf_init (lnx);
- mpf_init (pow2);
- mpf_init (frac);
+ mpfr_init (i2);
+ mpfr_init (ln2);
+ mpfr_init (absv);
+ mpfr_init (lnx);
+ mpfr_init (pow2);
+ mpfr_init (frac);
- natural_logarithm (&i2, &ln2);
+ mpfr_set_ui (i2, 2, GFC_RND_MODE);
+ mpfr_log (ln2, i2, GFC_RND_MODE);
- mpf_abs (absv, x->value.real);
- natural_logarithm (&absv, &lnx);
+ mpfr_abs (absv, x->value.real, GFC_RND_MODE);
+ mpfr_log (lnx, absv, GFC_RND_MODE);
- mpf_div (lnx, lnx, ln2);
- mpf_trunc (lnx, lnx);
- mpf_add_ui (lnx, lnx, 1);
+ mpfr_div (lnx, lnx, ln2, GFC_RND_MODE);
+ mpfr_trunc (lnx, lnx);
+ mpfr_add_ui (lnx, lnx, 1, GFC_RND_MODE);
/* Old exponent value, and fraction. */
- exp2 = (unsigned long) mpf_get_d (lnx);
- mpf_pow_ui (pow2, i2, exp2);
+ exp2 = (unsigned long) mpfr_get_d (lnx, GFC_RND_MODE);
+ mpfr_pow_ui (pow2, i2, exp2, GFC_RND_MODE);
- mpf_div (frac, absv, pow2);
+ mpfr_div (frac, absv, pow2, GFC_RND_MODE);
/* New exponent. */
exp2 = (unsigned long) mpz_get_d (i->value.integer);
- mpf_mul_2exp (result->value.real, frac, exp2);
+ mpfr_mul_2exp (result->value.real, frac, exp2, GFC_RND_MODE);
- mpf_clear (i2);
- mpf_clear (ln2);
- mpf_clear (absv);
- mpf_clear (lnx);
- mpf_clear (pow2);
- mpf_clear (frac);
+ mpfr_clear (i2);
+ mpfr_clear (ln2);
+ mpfr_clear (absv);
+ mpfr_clear (lnx);
+ mpfr_clear (pow2);
+ mpfr_clear (frac);
+ mpfr_clear (zero);
return range_check (result, "SET_EXPONENT");
}
case BT_REAL:
/* TODO: Handle -0.0 and +0.0 correctly on machines that support
it. */
- mpf_abs (result->value.real, x->value.real);
- if (mpf_sgn (y->value.integer) < 0)
- mpf_neg (result->value.real, result->value.real);
+ mpfr_abs (result->value.real, x->value.real, GFC_RND_MODE);
+ if (mpfr_sgn (y->value.real) < 0)
+ mpfr_neg (result->value.real, result->value.real, GFC_RND_MODE);
break;
gfc_simplify_sin (gfc_expr * x)
{
gfc_expr *result;
- mpf_t xp, xq;
+ mpfr_t xp, xq;
if (x->expr_type != EXPR_CONSTANT)
return NULL;
switch (x->ts.type)
{
case BT_REAL:
- sine (&x->value.real, &result->value.real);
+ mpfr_sin (result->value.real, x->value.real, GFC_RND_MODE);
break;
case BT_COMPLEX:
- mpf_init (xp);
- mpf_init (xq);
+ gfc_set_model (x->value.real);
+ mpfr_init (xp);
+ mpfr_init (xq);
- sine (&x->value.complex.r, &xp);
- hypercos (&x->value.complex.i, &xq);
- mpf_mul (result->value.complex.r, xp, xq);
+ mpfr_sin (xp, x->value.complex.r, GFC_RND_MODE);
+ mpfr_cosh (xq, x->value.complex.i, GFC_RND_MODE);
+ mpfr_mul (result->value.complex.r, xp, xq, GFC_RND_MODE);
- cosine (&x->value.complex.r, &xp);
- hypersine (&x->value.complex.i, &xq);
- mpf_mul (result->value.complex.i, xp, xq);
+ mpfr_cos (xp, x->value.complex.r, GFC_RND_MODE);
+ mpfr_sinh (xq, x->value.complex.i, GFC_RND_MODE);
+ mpfr_mul (result->value.complex.i, xp, xq, GFC_RND_MODE);
- mpf_clear (xp);
- mpf_clear (xq);
+ mpfr_clear (xp);
+ mpfr_clear (xq);
break;
default:
result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
- hypersine (&x->value.real, &result->value.real);
+ mpfr_sinh(result->value.real, x->value.real, GFC_RND_MODE);
return range_check (result, "SINH");
}
gfc_simplify_spacing (gfc_expr * x)
{
gfc_expr *result;
- mpf_t i1, i2, ln2, absv, lnx;
+ mpfr_t i1, i2, ln2, absv, lnx, zero;
long diff;
unsigned long exp2;
int i, p;
result = gfc_constant_result (BT_REAL, x->ts.kind, &x->where);
- if (mpf_cmp (x->value.real, mpf_zero) == 0)
+ gfc_set_model_kind (x->ts.kind);
+ mpfr_init (zero);
+ mpfr_set_ui (zero, 0, GFC_RND_MODE);
+
+ if (mpfr_cmp (x->value.real, zero) == 0)
{
- mpf_set (result->value.real, gfc_real_kinds[i].tiny);
+ mpfr_set (result->value.real, gfc_real_kinds[i].tiny, GFC_RND_MODE);
+ mpfr_clear (zero);
return result;
}
- mpf_init_set_ui (i1, 1);
- mpf_init_set_ui (i2, 2);
- mpf_init (ln2);
- mpf_init (absv);
- mpf_init (lnx);
+ mpfr_init (i1);
+ mpfr_init (i2);
+ mpfr_init (ln2);
+ mpfr_init (absv);
+ mpfr_init (lnx);
- natural_logarithm (&i2, &ln2);
+ mpfr_set_ui (i1, 1, GFC_RND_MODE);
+ mpfr_set_ui (i2, 2, GFC_RND_MODE);
- mpf_abs (absv, x->value.real);
- natural_logarithm (&absv, &lnx);
+ mpfr_log (ln2, i2, GFC_RND_MODE);
+ mpfr_abs (absv, x->value.real, GFC_RND_MODE);
+ mpfr_log (lnx, absv, GFC_RND_MODE);
- mpf_div (lnx, lnx, ln2);
- mpf_trunc (lnx, lnx);
- mpf_add_ui (lnx, lnx, 1);
+ mpfr_div (lnx, lnx, ln2, GFC_RND_MODE);
+ mpfr_trunc (lnx, lnx);
+ mpfr_add_ui (lnx, lnx, 1, GFC_RND_MODE);
- diff = (long) mpf_get_d (lnx) - (long) p;
+ diff = (long) mpfr_get_d (lnx, GFC_RND_MODE) - (long) p;
if (diff >= 0)
{
exp2 = (unsigned) diff;
- mpf_mul_2exp (result->value.real, i1, exp2);
+ mpfr_mul_2exp (result->value.real, i1, exp2, GFC_RND_MODE);
}
else
{
diff = -diff;
exp2 = (unsigned) diff;
- mpf_div_2exp (result->value.real, i1, exp2);
+ mpfr_div_2exp (result->value.real, i1, exp2, GFC_RND_MODE);
}
- mpf_clear (i1);
- mpf_clear (i2);
- mpf_clear (ln2);
- mpf_clear (absv);
- mpf_clear (lnx);
+ mpfr_clear (i1);
+ mpfr_clear (i2);
+ mpfr_clear (ln2);
+ mpfr_clear (absv);
+ mpfr_clear (lnx);
+ mpfr_clear (zero);
- if (mpf_cmp (result->value.real, gfc_real_kinds[i].tiny) < 0)
- mpf_set (result->value.real, gfc_real_kinds[i].tiny);
+ if (mpfr_cmp (result->value.real, gfc_real_kinds[i].tiny) < 0)
+ mpfr_set (result->value.real, gfc_real_kinds[i].tiny, GFC_RND_MODE);
return range_check (result, "SPACING");
}
gfc_simplify_sqrt (gfc_expr * e)
{
gfc_expr *result;
- mpf_t ac, ad, s, t, w;
+ mpfr_t ac, ad, s, t, w;
if (e->expr_type != EXPR_CONSTANT)
return NULL;
switch (e->ts.type)
{
case BT_REAL:
- if (mpf_cmp_si (e->value.real, 0) < 0)
+ if (mpfr_cmp_si (e->value.real, 0) < 0)
goto negative_arg;
- mpf_sqrt (result->value.real, e->value.real);
+ mpfr_sqrt (result->value.real, e->value.real, GFC_RND_MODE);
break;
/* Formula taken from Numerical Recipes to avoid over- and
underflow. */
- mpf_init (ac);
- mpf_init (ad);
- mpf_init (s);
- mpf_init (t);
- mpf_init (w);
+ gfc_set_model (e->value.real);
+ mpfr_init (ac);
+ mpfr_init (ad);
+ mpfr_init (s);
+ mpfr_init (t);
+ mpfr_init (w);
- if (mpf_cmp_ui (e->value.complex.r, 0) == 0
- && mpf_cmp_ui (e->value.complex.i, 0) == 0)
+ if (mpfr_cmp_ui (e->value.complex.r, 0) == 0
+ && mpfr_cmp_ui (e->value.complex.i, 0) == 0)
{
- mpf_set_ui (result->value.complex.r, 0);
- mpf_set_ui (result->value.complex.i, 0);
+ mpfr_set_ui (result->value.complex.r, 0, GFC_RND_MODE);
+ mpfr_set_ui (result->value.complex.i, 0, GFC_RND_MODE);
break;
}
- mpf_abs (ac, e->value.complex.r);
- mpf_abs (ad, e->value.complex.i);
+ mpfr_abs (ac, e->value.complex.r, GFC_RND_MODE);
+ mpfr_abs (ad, e->value.complex.i, GFC_RND_MODE);
- if (mpf_cmp (ac, ad) >= 0)
+ if (mpfr_cmp (ac, ad) >= 0)
{
- mpf_div (t, e->value.complex.i, e->value.complex.r);
- mpf_mul (t, t, t);
- mpf_add_ui (t, t, 1);
- mpf_sqrt (t, t);
- mpf_add_ui (t, t, 1);
- mpf_div_ui (t, t, 2);
- mpf_sqrt (t, t);
- mpf_sqrt (s, ac);
- mpf_mul (w, s, t);
+ mpfr_div (t, e->value.complex.i, e->value.complex.r, GFC_RND_MODE);
+ mpfr_mul (t, t, t, GFC_RND_MODE);
+ mpfr_add_ui (t, t, 1, GFC_RND_MODE);
+ mpfr_sqrt (t, t, GFC_RND_MODE);
+ mpfr_add_ui (t, t, 1, GFC_RND_MODE);
+ mpfr_div_ui (t, t, 2, GFC_RND_MODE);
+ mpfr_sqrt (t, t, GFC_RND_MODE);
+ mpfr_sqrt (s, ac, GFC_RND_MODE);
+ mpfr_mul (w, s, t, GFC_RND_MODE);
}
else
{
- mpf_div (s, e->value.complex.r, e->value.complex.i);
- mpf_mul (t, s, s);
- mpf_add_ui (t, t, 1);
- mpf_sqrt (t, t);
- mpf_abs (s, s);
- mpf_add (t, t, s);
- mpf_div_ui (t, t, 2);
- mpf_sqrt (t, t);
- mpf_sqrt (s, ad);
- mpf_mul (w, s, t);
+ mpfr_div (s, e->value.complex.r, e->value.complex.i, GFC_RND_MODE);
+ mpfr_mul (t, s, s, GFC_RND_MODE);
+ mpfr_add_ui (t, t, 1, GFC_RND_MODE);
+ mpfr_sqrt (t, t, GFC_RND_MODE);
+ mpfr_abs (s, s, GFC_RND_MODE);
+ mpfr_add (t, t, s, GFC_RND_MODE);
+ mpfr_div_ui (t, t, 2, GFC_RND_MODE);
+ mpfr_sqrt (t, t, GFC_RND_MODE);
+ mpfr_sqrt (s, ad, GFC_RND_MODE);
+ mpfr_mul (w, s, t, GFC_RND_MODE);
}
- if (mpf_cmp_ui (w, 0) != 0 && mpf_cmp_ui (e->value.complex.r, 0) >= 0)
+ if (mpfr_cmp_ui (w, 0) != 0 && mpfr_cmp_ui (e->value.complex.r, 0) >= 0)
{
- mpf_mul_ui (t, w, 2);
- mpf_div (result->value.complex.i, e->value.complex.i, t);
- mpf_set (result->value.complex.r, w);
+ mpfr_mul_ui (t, w, 2, GFC_RND_MODE);
+ mpfr_div (result->value.complex.i, e->value.complex.i, t, GFC_RND_MODE);
+ mpfr_set (result->value.complex.r, w, GFC_RND_MODE);
}
- else if (mpf_cmp_ui (w, 0) != 0
- && mpf_cmp_ui (e->value.complex.r, 0) < 0
- && mpf_cmp_ui (e->value.complex.i, 0) >= 0)
+ else if (mpfr_cmp_ui (w, 0) != 0
+ && mpfr_cmp_ui (e->value.complex.r, 0) < 0
+ && mpfr_cmp_ui (e->value.complex.i, 0) >= 0)
{
- mpf_mul_ui (t, w, 2);
- mpf_div (result->value.complex.r, e->value.complex.i, t);
- mpf_set (result->value.complex.i, w);
+ mpfr_mul_ui (t, w, 2, GFC_RND_MODE);
+ mpfr_div (result->value.complex.r, e->value.complex.i, t, GFC_RND_MODE);
+ mpfr_set (result->value.complex.i, w, GFC_RND_MODE);
}
- else if (mpf_cmp_ui (w, 0) != 0
- && mpf_cmp_ui (e->value.complex.r, 0) < 0
- && mpf_cmp_ui (e->value.complex.i, 0) < 0)
+ else if (mpfr_cmp_ui (w, 0) != 0
+ && mpfr_cmp_ui (e->value.complex.r, 0) < 0
+ && mpfr_cmp_ui (e->value.complex.i, 0) < 0)
{
- mpf_mul_ui (t, w, 2);
- mpf_div (result->value.complex.r, ad, t);
- mpf_neg (w, w);
- mpf_set (result->value.complex.i, w);
+ mpfr_mul_ui (t, w, 2, GFC_RND_MODE);
+ mpfr_div (result->value.complex.r, ad, t, GFC_RND_MODE);
+ mpfr_neg (w, w, GFC_RND_MODE);
+ mpfr_set (result->value.complex.i, w, GFC_RND_MODE);
}
else
gfc_internal_error ("invalid complex argument of SQRT at %L",
&e->where);
- mpf_clear (s);
- mpf_clear (t);
- mpf_clear (ac);
- mpf_clear (ad);
- mpf_clear (w);
+ mpfr_clear (s);
+ mpfr_clear (t);
+ mpfr_clear (ac);
+ mpfr_clear (ad);
+ mpfr_clear (w);
break;
gfc_expr *
gfc_simplify_tan (gfc_expr * x)
{
- gfc_expr *result;
- mpf_t mpf_sin, mpf_cos, mag_cos;
int i;
+ gfc_expr *result;
if (x->expr_type != EXPR_CONSTANT)
return NULL;
result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
- mpf_init (mpf_sin);
- mpf_init (mpf_cos);
- mpf_init (mag_cos);
- sine (&x->value.real, &mpf_sin);
- cosine (&x->value.real, &mpf_cos);
- mpf_abs (mag_cos, mpf_cos);
- if (mpf_cmp_ui (mag_cos, 0) == 0)
- {
- gfc_error ("Tangent undefined at %L", &x->where);
- mpf_clear (mpf_sin);
- mpf_clear (mpf_cos);
- mpf_clear (mag_cos);
- gfc_free_expr (result);
- return &gfc_bad_expr;
- }
- else if (mpf_cmp (mag_cos, gfc_real_kinds[i].tiny) < 0)
- {
- gfc_error ("Tangent cannot be accurately evaluated at %L", &x->where);
- mpf_clear (mpf_sin);
- mpf_clear (mpf_cos);
- mpf_clear (mag_cos);
- gfc_free_expr (result);
- return &gfc_bad_expr;
- }
- else
- {
- mpf_div (result->value.real, mpf_sin, mpf_cos);
- mpf_clear (mpf_sin);
- mpf_clear (mpf_cos);
- mpf_clear (mag_cos);
- }
+ mpfr_tan (result->value.real, x->value.real, GFC_RND_MODE);
return range_check (result, "TAN");
}
gfc_simplify_tanh (gfc_expr * x)
{
gfc_expr *result;
- mpf_t xp, xq;
if (x->expr_type != EXPR_CONSTANT)
return NULL;
result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
- mpf_init (xp);
- mpf_init (xq);
-
- hypersine (&x->value.real, &xq);
- hypercos (&x->value.real, &xp);
-
- mpf_div (result->value.real, xq, xp);
-
- mpf_clear (xp);
- mpf_clear (xq);
+ mpfr_tanh (result->value.real, x->value.real, GFC_RND_MODE);
return range_check (result, "TANH");
gfc_internal_error ("gfc_simplify_error(): Bad kind");
result = gfc_constant_result (BT_REAL, e->ts.kind, &e->where);
- mpf_set (result->value.real, gfc_real_kinds[i].tiny);
+ mpfr_set (result->value.real, gfc_real_kinds[i].tiny, GFC_RND_MODE);
return result;
}
gfc_simplify_init_1 (void)
{
- mpf_init_set_str (mpf_zero, "0.0", 10);
- mpf_init_set_str (mpf_half, "0.5", 10);
- mpf_init_set_str (mpf_one, "1.0", 10);
- mpz_init_set_str (mpz_zero, "0", 10);
-
invert_table (ascii_table, xascii_table);
}
-
-
-void
-gfc_simplify_done_1 (void)
-{
-
- mpf_clear (mpf_zero);
- mpf_clear (mpf_half);
- mpf_clear (mpf_one);
- mpz_clear (mpz_zero);
-}
/* Converts a real constant into backend form. Uses an intermediate string
representation. */
tree
-gfc_conv_mpf_to_tree (mpf_t f, int kind)
+gfc_conv_mpfr_to_tree (mpfr_t f, int kind)
{
tree res;
tree type;
}
assert (gfc_real_kinds[n].kind);
- assert (gfc_real_kinds[n].radix == 2);
-
n = MAX (abs (gfc_real_kinds[n].min_exponent),
abs (gfc_real_kinds[n].max_exponent));
-#if 0
- edigits = 2 + (int) (log (n) / log (gfc_real_kinds[n].radix));
-#endif
+
edigits = 1;
while (n > 0)
{
edigits += 3;
}
+ if (kind == gfc_default_double_kind())
+ p = mpfr_get_str (NULL, &exp, 10, 17, f, GFC_RND_MODE);
+ else
+ p = mpfr_get_str (NULL, &exp, 10, 8, f, GFC_RND_MODE);
- p = mpf_get_str (NULL, &exp, 10, 0, f);
/* We also have one minus sign, "e", "." and a null terminator. */
q = (char *) gfc_getmem (strlen (p) + edigits + 4);
type = gfc_get_real_type (kind);
res = build_real (type, REAL_VALUE_ATOF (q, TYPE_MODE (type)));
+
gfc_free (q);
gfc_free (p);
return gfc_conv_mpz_to_tree (expr->value.integer, expr->ts.kind);
case BT_REAL:
- return gfc_conv_mpf_to_tree (expr->value.real, expr->ts.kind);
+ return gfc_conv_mpfr_to_tree (expr->value.real, expr->ts.kind);
case BT_LOGICAL:
return build_int_2 (expr->value.logical, 0);
case BT_COMPLEX:
{
- tree real = gfc_conv_mpf_to_tree (expr->value.complex.r,
+ tree real = gfc_conv_mpfr_to_tree (expr->value.complex.r,
expr->ts.kind);
- tree imag = gfc_conv_mpf_to_tree (expr->value.complex.i,
+ tree imag = gfc_conv_mpfr_to_tree (expr->value.complex.i,
expr->ts.kind);
return build_complex (NULL_TREE, real, imag);
tree gfc_conv_mpz_to_tree (mpz_t, int);
/* Returns a REAL_CST. */
-tree gfc_conv_mpf_to_tree (mpf_t, int);
+tree gfc_conv_mpfr_to_tree (mpfr_t, int);
/* Build a tree for a constant. Must be an EXPR_CONSTANT gfc_expr.
For CHARACTER literal constants, the caller still has to set the
#include "real.h"
#include "tree-gimple.h"
#include "flags.h"
-#include <gmp.h>
#include <assert.h>
#include "gfortran.h"
+#include "arith.h"
#include "intrinsic.h"
#include "trans.h"
#include "trans-const.h"
tree arg;
tree tmp;
tree cond;
- mpf_t huge;
+ mpfr_t huge;
int n;
int kind;
arg = gfc_evaluate_now (arg, &se->pre);
/* Test if the value is too large to handle sensibly. */
- mpf_init (huge);
+ gfc_set_model_kind (kind);
+ mpfr_init (huge);
n = gfc_validate_kind (BT_INTEGER, kind);
- mpf_set_z (huge, gfc_integer_kinds[n].huge);
- tmp = gfc_conv_mpf_to_tree (huge, kind);
+ mpfr_set_z (huge, gfc_integer_kinds[n].huge, GFC_RND_MODE);
+ tmp = gfc_conv_mpfr_to_tree (huge, kind);
cond = build (LT_EXPR, boolean_type_node, arg, tmp);
- mpf_neg (huge, huge);
- tmp = gfc_conv_mpf_to_tree (huge, kind);
+ mpfr_neg (huge, huge, GFC_RND_MODE);
+ tmp = gfc_conv_mpfr_to_tree (huge, kind);
tmp = build (GT_EXPR, boolean_type_node, arg, tmp);
cond = build (TRUTH_AND_EXPR, boolean_type_node, cond, tmp);
itype = gfc_get_int_type (kind);
tmp = build_fix_expr (&se->pre, arg, itype, op);
tmp = convert (type, tmp);
se->expr = build (COND_EXPR, type, cond, tmp, arg);
+ mpfr_clear (huge);
}
tree zero;
tree test;
tree test2;
- mpf_t huge;
+ mpfr_t huge;
int n;
arg = gfc_conv_intrinsic_function_args (se, expr);
tmp = build (RDIV_EXPR, type, arg, arg2);
/* Test if the value is too large to handle sensibly. */
- mpf_init (huge);
+ gfc_set_model_kind (expr->ts.kind);
+ mpfr_init (huge);
n = gfc_validate_kind (BT_INTEGER, expr->ts.kind);
- mpf_set_z (huge, gfc_integer_kinds[n].huge);
- test = gfc_conv_mpf_to_tree (huge, expr->ts.kind);
+ mpfr_set_z (huge, gfc_integer_kinds[n].huge, GFC_RND_MODE);
+ test = gfc_conv_mpfr_to_tree (huge, expr->ts.kind);
test2 = build (LT_EXPR, boolean_type_node, tmp, test);
- mpf_neg (huge, huge);
- test = gfc_conv_mpf_to_tree (huge, expr->ts.kind);
+ mpfr_neg (huge, huge, GFC_RND_MODE);
+ test = gfc_conv_mpfr_to_tree (huge, expr->ts.kind);
test = build (GT_EXPR, boolean_type_node, tmp, test);
test2 = build (TRUTH_AND_EXPR, boolean_type_node, test, test2);
tmp = build (COND_EXPR, type, test2, tmp, arg);
tmp = build (MULT_EXPR, type, tmp, arg2);
se->expr = build (MINUS_EXPR, type, arg, tmp);
+ mpfr_clear (huge);
break;
default:
switch (arrayexpr->ts.type)
{
case BT_REAL:
- tmp = gfc_conv_mpf_to_tree (gfc_real_kinds[n].huge, arrayexpr->ts.kind);
+ tmp = gfc_conv_mpfr_to_tree (gfc_real_kinds[n].huge, arrayexpr->ts.kind);
break;
case BT_INTEGER:
switch (expr->ts.type)
{
case BT_REAL:
- tmp = gfc_conv_mpf_to_tree (gfc_real_kinds[n].huge, expr->ts.kind);
+ tmp = gfc_conv_mpfr_to_tree (gfc_real_kinds[n].huge, expr->ts.kind);
break;
case BT_INTEGER:
OSDN Git Service
