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. */
}
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;
/* Reset emin and emax. */
mpfr_set_emin (emin);
mpfr_set_emax (emax);
-#endif
/* Copy sign if needed. */
if (mpfr_sgn (p) < 0)
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");
}
int i;
gfc_expr *result;
-#if defined(GFC_MPFR_TOO_OLD)
- mpfr_t tmp;
-#endif
-
if (x->expr_type != EXPR_CONSTANT)
return NULL;
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");
}
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);
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;
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");
}
}
-#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)
{
return range_check (result, "RRSPACING");
}
-#endif
+
gfc_expr *
gfc_simplify_scale (gfc_expr * x, gfc_expr * i)
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)
{
return range_check (result, "SPACING");
}
-#endif
+
gfc_expr *
gfc_simplify_sqrt (gfc_expr * e)