[PATCH,fortran] Removekludges to work with older versions of MPFR
Steve Kargl
sgk@troutmask.apl.washington.edu
Mon Oct 23 06:39:00 GMT 2006
This patch has been bootstrapped and regression tested
on i386-*-freebsd.
It removes all of the kludges that allowed gfortran to
mostly work with older versions of MPFR. The patch is
self explanatory, but I'll note that it removes about
270 lines of code.
2006-10-22 Steven G. Kargl <kargl@gcc.gnu.org>
* gfortran.h: Remove GFC_MPFR_TOO_OLD.
* arith.c (arctangent2): Remove function
(gfc_check_real_range): Remove subnormal kludge.
* arith.h: Remove arctangent2 prototype.
* simplify.c: (gfc_simplify_atan2): Remove use of arctangent2.
(gfc_simplify_exponent, gfc_simplify_log, gfc_simplify_nearest,
gfc_simplify_rrspacing, gfc_simplify_spacing): Remove mpfr kludges.
--
Steve
-------------- next part --------------
Index: gfortran.h
===================================================================
--- gfortran.h (revision 117959)
+++ gfortran.h (working copy)
@@ -1224,10 +1224,6 @@ gfc_intrinsic_sym;
#include <gmp.h>
#include <mpfr.h>
#define GFC_RND_MODE GMP_RNDN
-#undef GFC_MPFR_TOO_OLD
-#if MPFR_VERSION_MAJOR < 2 || (MPFR_VERSION_MAJOR == 2 && MPFR_VERSION_MINOR < 2)
-#define GFC_MPFR_TOO_OLD 1
-#endif
typedef struct gfc_expr
{
Index: arith.c
===================================================================
--- arith.c (revision 117959)
+++ arith.c (working copy)
@@ -75,57 +75,6 @@ gfc_set_model (mpfr_t x)
mpfr_set_default_prec (mpfr_get_prec (x));
}
-#if defined(GFC_MPFR_TOO_OLD)
-/* Calculate atan2 (y, x)
-
-atan2(y, x) = atan(y/x) if x > 0,
- sign(y)*(pi - atan(|y/x|)) if x < 0,
- 0 if x = 0 && y == 0,
- sign(y)*pi/2 if x = 0 && y != 0.
-*/
-
-void
-arctangent2 (mpfr_t y, mpfr_t x, mpfr_t result)
-{
- int i;
- mpfr_t t;
-
- gfc_set_model (y);
- mpfr_init (t);
-
- i = mpfr_sgn (x);
-
- if (i > 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
- {
- 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);
- }
- }
-
- mpfr_clear (t);
-}
-#endif
/* Given an arithmetic error code, return a pointer to a string that
explains the error. */
@@ -412,31 +361,6 @@ gfc_check_real_range (mpfr_t p, int kind
}
else if (mpfr_cmp (q, gfc_real_kinds[i].tiny) < 0)
{
-#if defined(GFC_MPFR_TOO_OLD)
- /* MPFR operates on a number with a given precision and enormous
- exponential range. To represent subnormal numbers, the exponent is
- allowed to become smaller than emin, but always retains the full
- precision. This code resets unused bits to 0 to alleviate
- rounding problems. Note, a future version of MPFR will have a
- mpfr_subnormalize() function, which handles this truncation in a
- more efficient and robust way. */
-
- int j, k;
- char *bin, *s;
- mp_exp_t e;
-
- bin = mpfr_get_str (NULL, &e, gfc_real_kinds[i].radix, 0, q, GMP_RNDN);
- k = gfc_real_kinds[i].digits - (gfc_real_kinds[i].min_exponent - e);
- for (j = k; j < gfc_real_kinds[i].digits; j++)
- bin[j] = '0';
- /* Need space for '0.', bin, 'E', and e */
- s = (char *) gfc_getmem (strlen(bin) + 10);
- sprintf (s, "0.%sE%d", bin, (int) e);
- mpfr_set_str (q, s, gfc_real_kinds[i].radix, GMP_RNDN);
-
- gfc_free (s);
- gfc_free (bin);
-#else
mp_exp_t emin, emax;
int en;
@@ -453,7 +377,6 @@ gfc_check_real_range (mpfr_t p, int kind
/* Reset emin and emax. */
mpfr_set_emin (emin);
mpfr_set_emax (emax);
-#endif
/* Copy sign if needed. */
if (mpfr_sgn (p) < 0)
Index: arith.h
===================================================================
--- arith.h (revision 117959)
+++ arith.h (working copy)
@@ -25,11 +25,9 @@ Software Foundation, 51 Franklin Street,
#include "gfortran.h"
-/* 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. */
+/* 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);
Index: simplify.c
===================================================================
--- simplify.c (revision 117959)
+++ simplify.c (working copy)
@@ -607,11 +607,7 @@ gfc_simplify_atan2 (gfc_expr * y, gfc_ex
return &gfc_bad_expr;
}
-#if defined(GFC_MPFR_TOO_OLD)
- arctangent2 (y->value.real, x->value.real, result->value.real);
-#else
mpfr_atan2 (result->value.real, y->value.real, x->value.real, GFC_RND_MODE);
-#endif
return range_check (result, "ATAN2");
}
@@ -1060,10 +1056,6 @@ gfc_simplify_exponent (gfc_expr * x)
int i;
gfc_expr *result;
-#if defined(GFC_MPFR_TOO_OLD)
- mpfr_t tmp;
-#endif
-
if (x->expr_type != EXPR_CONSTANT)
return NULL;
@@ -1078,27 +1070,8 @@ gfc_simplify_exponent (gfc_expr * x)
return result;
}
-#if defined(GFC_MPFR_TOO_OLD)
- /* PR fortran/28276 suffers from a buggy MPFR, and this block of code
- does not function correctly. */
- mpfr_init (tmp);
-
- mpfr_abs (tmp, x->value.real, GFC_RND_MODE);
- mpfr_log2 (tmp, tmp, GFC_RND_MODE);
-
- gfc_mpfr_to_mpz (result->value.integer, tmp);
-
- /* The model number for tiny(x) is b**(emin - 1) where b is the base and emin
- is the smallest exponent value. So, we need to add 1 if x is tiny(x). */
- i = gfc_validate_kind (x->ts.type, x->ts.kind, false);
- if (mpfr_cmp (x->value.real, gfc_real_kinds[i].tiny) == 0)
- mpz_add_ui (result->value.integer,result->value.integer, 1);
-
- mpfr_clear (tmp);
-#else
i = (int) mpfr_get_exp (x->value.real);
mpz_set_si (result->value.integer, i);
-#endif
return range_check (result, "EXPONENT");
}
@@ -2160,13 +2133,8 @@ gfc_simplify_log (gfc_expr * x)
mpfr_init (xr);
mpfr_init (xi);
-#if defined(GFC_MPFR_TOO_OLD)
- arctangent2 (x->value.complex.i, x->value.complex.r, result->value.complex.i);
-#else
mpfr_atan2 (result->value.complex.i, x->value.complex.i, x->value.complex.r,
GFC_RND_MODE);
-#endif
-
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);
@@ -2494,9 +2462,6 @@ gfc_simplify_nearest (gfc_expr * x, gfc_
gfc_expr *result;
mpfr_t tmp;
int sgn;
-#if defined(GFC_MPFR_TOO_OLD)
- int direction;
-#endif
if (x->expr_type != EXPR_CONSTANT || s->expr_type != EXPR_CONSTANT)
return NULL;
@@ -2510,60 +2475,11 @@ gfc_simplify_nearest (gfc_expr * x, gfc_
gfc_set_model_kind (x->ts.kind);
result = gfc_copy_expr (x);
-#if defined(GFC_MPFR_TOO_OLD)
-
- direction = mpfr_sgn (s->value.real);
- sgn = mpfr_sgn (x->value.real);
-
- if (sgn == 0)
- {
- int k = gfc_validate_kind (BT_REAL, x->ts.kind, 0);
-
- if (direction > 0)
- mpfr_add (result->value.real,
- x->value.real, gfc_real_kinds[k].subnormal, GFC_RND_MODE);
- else
- mpfr_sub (result->value.real,
- x->value.real, gfc_real_kinds[k].subnormal, GFC_RND_MODE);
- }
- else
- {
- if (sgn < 0)
- {
- direction = -direction;
- mpfr_neg (result->value.real, result->value.real, GFC_RND_MODE);
- }
-
- if (direction > 0)
- mpfr_add_one_ulp (result->value.real, GFC_RND_MODE);
- else
- {
- /* In this case the exponent can shrink, which makes us skip
- over one number because we subtract one ulp with the
- larger exponent. Thus we need to compensate for this. */
- mpfr_init_set (tmp, result->value.real, GFC_RND_MODE);
-
- mpfr_sub_one_ulp (result->value.real, GFC_RND_MODE);
- mpfr_add_one_ulp (result->value.real, GFC_RND_MODE);
-
- /* If we're back to where we started, the spacing is one
- ulp, and we get the correct result by subtracting. */
- if (mpfr_cmp (tmp, result->value.real) == 0)
- mpfr_sub_one_ulp (result->value.real, GFC_RND_MODE);
-
- mpfr_clear (tmp);
- }
-
- if (sgn < 0)
- mpfr_neg (result->value.real, result->value.real, GFC_RND_MODE);
- }
-#else
sgn = mpfr_sgn (s->value.real);
mpfr_init (tmp);
mpfr_set_inf (tmp, sgn);
mpfr_nexttoward (result->value.real, tmp);
mpfr_clear(tmp);
-#endif
return range_check (result, "NEAREST");
}
@@ -3112,57 +3028,6 @@ bad_reshape:
}
-#if defined(GFC_MPFR_TOO_OLD)
-gfc_expr *
-gfc_simplify_rrspacing (gfc_expr * x)
-{
- gfc_expr *result;
- mpfr_t absv, log2, exp, frac, pow2;
- int i, p;
-
- if (x->expr_type != EXPR_CONSTANT)
- return NULL;
-
- i = gfc_validate_kind (x->ts.type, x->ts.kind, false);
-
- result = gfc_constant_result (BT_REAL, x->ts.kind, &x->where);
-
- p = gfc_real_kinds[i].digits;
-
- gfc_set_model_kind (x->ts.kind);
-
- if (mpfr_sgn (x->value.real) == 0)
- {
- mpfr_set_ui (result->value.real, 0, GFC_RND_MODE);
- return result;
- }
-
- mpfr_init (log2);
- mpfr_init (absv);
- mpfr_init (frac);
- mpfr_init (pow2);
- mpfr_init (exp);
-
- mpfr_abs (absv, x->value.real, GFC_RND_MODE);
- mpfr_log2 (log2, absv, GFC_RND_MODE);
-
- mpfr_trunc (log2, log2);
- mpfr_add_ui (exp, log2, 1, GFC_RND_MODE);
-
- mpfr_ui_pow (pow2, 2, exp, GFC_RND_MODE);
- mpfr_div (frac, absv, pow2, GFC_RND_MODE);
-
- mpfr_mul_2exp (result->value.real, frac, (unsigned long)p, GFC_RND_MODE);
-
- mpfr_clear (log2);
- mpfr_clear (absv);
- mpfr_clear (frac);
- mpfr_clear (pow2);
- mpfr_clear (exp);
-
- return range_check (result, "RRSPACING");
-}
-#else
gfc_expr *
gfc_simplify_rrspacing (gfc_expr * x)
{
@@ -3195,7 +3060,7 @@ gfc_simplify_rrspacing (gfc_expr * x)
return range_check (result, "RRSPACING");
}
-#endif
+
gfc_expr *
gfc_simplify_scale (gfc_expr * x, gfc_expr * i)
@@ -3639,58 +3504,7 @@ gfc_simplify_sngl (gfc_expr * a)
return range_check (result, "SNGL");
}
-#if defined(GFC_MPFR_TOO_OLD)
-gfc_expr *
-gfc_simplify_spacing (gfc_expr * x)
-{
- gfc_expr *result;
- mpfr_t absv, log2;
- long diff;
- int i, p;
-
- if (x->expr_type != EXPR_CONSTANT)
- return NULL;
-
- i = gfc_validate_kind (x->ts.type, x->ts.kind, false);
- p = gfc_real_kinds[i].digits;
-
- result = gfc_constant_result (BT_REAL, x->ts.kind, &x->where);
-
- gfc_set_model_kind (x->ts.kind);
-
- /* Special case x = 0 and -0. */
- mpfr_init (absv);
- mpfr_abs (absv, x->value.real, GFC_RND_MODE);
- if (mpfr_sgn (absv) == 0)
- {
- mpfr_set (result->value.real, gfc_real_kinds[i].tiny, GFC_RND_MODE);
- return result;
- }
-
- mpfr_init (log2);
- mpfr_log2 (log2, absv, GFC_RND_MODE);
- mpfr_trunc (log2, log2);
-
- mpfr_add_ui (log2, log2, 1, GFC_RND_MODE);
-
- /* FIXME: We should be using mpfr_get_si here, but this function is
- not available with the version of mpfr distributed with gmp (as of
- 2004-09-17). Replace once mpfr has been imported into the gcc cvs
- tree. */
- diff = (long)mpfr_get_d (log2, GFC_RND_MODE) - (long)p;
- mpfr_set_ui (result->value.real, 1, GFC_RND_MODE);
- mpfr_mul_2si (result->value.real, result->value.real, diff, GFC_RND_MODE);
-
- mpfr_clear (log2);
- mpfr_clear (absv);
-
- 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");
-}
-#else
gfc_expr *
gfc_simplify_spacing (gfc_expr * x)
{
@@ -3727,7 +3541,7 @@ gfc_simplify_spacing (gfc_expr * x)
return range_check (result, "SPACING");
}
-#endif
+
gfc_expr *
gfc_simplify_sqrt (gfc_expr * e)
More information about the Fortran
mailing list