Skip to content

R
riscv-gcc-1
Commit f8e566e5 by Steven G. Kargl Committed by Paul Brook
Options

arith.c: Add #define for model numbers. 

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.

From-SVN: r85652
parent 1b4ed0bc
Showing with 890 additions and 1334 deletions
gcc/fortran/ChangeLog
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>
......
gcc/fortran/arith.c
......@@ -32,8 +32,6 @@ Software Foundation, 59 Temple Place - Suite 330, Boston, MA
#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.
......@@ -69,9 +67,31 @@ gfc_logical_info gfc_logical_kinds[] = {
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)
};
......@@ -82,440 +102,67 @@ gfc_real_info gfc_real_kinds[] = {
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--;
}
mp_exp_t e;
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);
if (i % 2 == 0)
mpf_sub (log, log, t);
e = mpfr_get_z_exp (z, x);
if (e > 0)
mpz_mul_2exp (z, z, e);
else
mpf_add (log, log, t);
}
/* Add in the log (e^p) = p */
if (p < 0)
mpf_sub_ui (log, log, -p);
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
switch (kind)
{
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++)
{
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)
{
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)
switch (mpfr_get_prec (x))
{
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,
......@@ -525,97 +172,46 @@ 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);
mpf_neg (neg, x);
exponential (&x, &term1);
exponential (&neg, &term2);
mpf_add (xp, term1, term2);
mpf_div_ui (*result, xp, 2);
mpfr_clear (t);
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);
}
......@@ -638,6 +234,9 @@ gfc_arith_error (arith code)
case ARITH_UNDERFLOW:
p = "Arithmetic underflow";
break;
case ARITH_NAN:
p = "Arithmetic NaN";
break;
case ARITH_DIV0:
p = "Division by zero";
break;
......@@ -662,72 +261,17 @@ gfc_arith_init_1 (void)
{
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 */
......@@ -751,59 +295,76 @@ gfc_arith_init_1 (void)
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);
if (mpf_cmp (a, b) > 0)
mpf_set (a, b); /* a = min(a, b) */
mpfr_init (real_info->epsilon);
mpfr_set (real_info->epsilon, b, GFC_RND_MODE);
mpf_trunc (a, a);
mpz_set_f (r, a);
/* 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);
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
......@@ -811,11 +372,13 @@ gfc_arith_init_1 (void)
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);
}
......@@ -827,12 +390,6 @@ gfc_arith_done_1 (void)
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);
......@@ -842,9 +399,9 @@ gfc_arith_done_1 (void)
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);
}
}
......@@ -1022,34 +579,35 @@ gfc_check_integer_range (mpz_t p, int kind)
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;
}
......@@ -1081,12 +639,14 @@ gfc_constant_result (bt type, int kind, locus * where)
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:
......@@ -1173,9 +733,7 @@ gfc_arith_neqv (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
/* 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
......@@ -1192,18 +750,18 @@ gfc_range_check (gfc_expr * e)
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;
......@@ -1244,12 +802,12 @@ gfc_arith_uminus (gfc_expr * op1, gfc_expr ** resultp)
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:
......@@ -1289,15 +847,16 @@ gfc_arith_plus (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
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:
......@@ -1337,16 +896,16 @@ gfc_arith_minus (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
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:
......@@ -1375,7 +934,7 @@ static arith
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);
......@@ -1387,23 +946,28 @@ gfc_arith_times (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
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. */
gfc_set_model (op1->value.complex.r);
mpfr_init (x);
mpfr_init (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);
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);
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);
mpf_clear (x);
mpf_clear (y);
mpfr_clear (x);
mpfr_clear (y);
break;
......@@ -1433,7 +997,7 @@ static arith
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;
......@@ -1454,44 +1018,51 @@ gfc_arith_divide (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
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;
......@@ -1523,30 +1094,31 @@ gfc_arith_divide (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
static void
complex_reciprocal (gfc_expr * op)
{
mpf_t mod, a, result_r, result_i;
mpfr_t mod, a, re, im;
mpf_init (mod);
mpf_init (a);
gfc_set_model (op->value.complex.r);
mpfr_init (mod);
mpfr_init (a);
mpfr_init (re);
mpfr_init (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);
/* 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_r);
mpf_div (result_r, op->value.complex.r, mod);
mpfr_div (re, op->value.complex.r, mod, GFC_RND_MODE);
mpf_init (result_i);
mpf_neg (result_i, op->value.complex.i);
mpf_div (result_i, result_i, mod);
mpfr_neg (im, op->value.complex.i, GFC_RND_MODE);
mpfr_div (im, im, mod, GFC_RND_MODE);
mpf_set (op->value.complex.r, result_r);
mpf_set (op->value.complex.i, result_i);
mpfr_set (op->value.complex.r, re, GFC_RND_MODE);
mpfr_set (op->value.complex.i, im, GFC_RND_MODE);
mpf_clear (result_r);
mpf_clear (result_i);
mpf_clear (mod);
mpf_clear (a);
mpfr_clear (re);
mpfr_clear (im);
mpfr_clear (mod);
mpfr_clear (a);
}
......@@ -1555,32 +1127,37 @@ complex_reciprocal (gfc_expr * op)
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);
}
......@@ -1592,7 +1169,7 @@ gfc_arith_power (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
int power, apower;
gfc_expr *result;
mpz_t unity_z;
mpf_t unity_f;
mpfr_t unity_f;
arith rc;
rc = ARITH_OK;
......@@ -1611,25 +1188,23 @@ gfc_arith_power (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
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;
......@@ -1638,8 +1213,7 @@ gfc_arith_power (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
gfc_internal_error ("gfc_arith_power(): Bad base");
}
}
if (power != 0)
else
{
apower = power;
if (power < 0)
......@@ -1661,22 +1235,24 @@ gfc_arith_power (gfc_expr * op1, gfc_expr * op2, gfc_expr ** resultp)
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:
......@@ -1748,7 +1324,7 @@ gfc_compare_expr (gfc_expr * op1, gfc_expr * op2)
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:
......@@ -1775,8 +1351,8 @@ static int
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);
}
......@@ -2544,12 +2120,12 @@ gfc_convert_real (const char *buffer, int kind, locus * where)
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;
}
......@@ -2564,8 +2140,8 @@ gfc_convert_complex (gfc_expr * real, gfc_expr * imag, int kind)
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;
}
......@@ -2621,7 +2197,7 @@ gfc_int2real (gfc_expr * src, int kind)
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)
{
......@@ -2644,8 +2220,8 @@ gfc_int2complex (gfc_expr * src, int kind)
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)
{
......@@ -2668,7 +2244,7 @@ gfc_real2int (gfc_expr * src, int kind)
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)
......@@ -2692,7 +2268,7 @@ gfc_real2real (gfc_expr * src, int kind)
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);
......@@ -2700,7 +2276,7 @@ gfc_real2real (gfc_expr * src, int 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)
{
......@@ -2723,8 +2299,8 @@ gfc_real2complex (gfc_expr * src, int kind)
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);
......@@ -2732,7 +2308,7 @@ gfc_real2complex (gfc_expr * src, int 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)
{
......@@ -2755,7 +2331,7 @@ gfc_complex2int (gfc_expr * src, int kind)
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)
......@@ -2779,7 +2355,7 @@ gfc_complex2real (gfc_expr * src, int kind)
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);
......@@ -2787,7 +2363,7 @@ gfc_complex2real (gfc_expr * src, int 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)
{
......@@ -2810,8 +2386,8 @@ gfc_complex2complex (gfc_expr * src, int kind)
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);
......@@ -2819,7 +2395,7 @@ gfc_complex2complex (gfc_expr * src, int 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)
{
......@@ -2834,7 +2410,7 @@ gfc_complex2complex (gfc_expr * src, int kind)
{
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)
{
......
gcc/fortran/arith.h
......@@ -24,19 +24,14 @@ Software Foundation, 59 Temple Place - Suite 330, Boston, MA
#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 *);
......
gcc/fortran/dump-parse-tree.c
......@@ -363,7 +363,7 @@ gfc_show_expr (gfc_expr * p)
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;
......@@ -388,13 +388,13 @@ gfc_show_expr (gfc_expr * p)
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);
......
gcc/fortran/expr.c
......@@ -154,7 +154,7 @@ free_expr0 (gfc_expr * e)
break;
case BT_REAL:
mpf_clear (e->value.real);
mpfr_clear (e->value.real);
break;
case BT_CHARACTER:
......@@ -162,8 +162,8 @@ free_expr0 (gfc_expr * e)
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:
......@@ -365,12 +365,17 @@ gfc_copy_expr (gfc_expr * p)
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:
......
gcc/fortran/gfortran.h
......@@ -59,7 +59,6 @@ char *alloca ();
/* 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. */
......@@ -184,7 +183,7 @@ extern mstring intrinsic_operators[];
/* 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;
......@@ -930,6 +929,8 @@ gfc_intrinsic_sym;
EXPR_ARRAY An array constructor. */
#include <gmp.h>
#include <mpfr.h>
#define GFC_RND_MODE GMP_RNDN
typedef struct gfc_expr
{
......@@ -953,13 +954,14 @@ 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;
......@@ -1023,7 +1025,7 @@ typedef struct
int kind, radix, digits, min_exponent, max_exponent;
int range, precision;
mpf_t epsilon, huge, tiny;
mpfr_t epsilon, huge, tiny;
}
gfc_real_info;
......@@ -1555,7 +1557,6 @@ match gfc_intrinsic_sub_interface (gfc_code *, int);
/* simplify.c */
void gfc_simplify_init_1 (void);
void gfc_simplify_done_1 (void);
/* match.c -- FIXME */
void gfc_free_iterator (gfc_iterator *, int);
......
gcc/fortran/intrinsic.c
......@@ -1492,6 +1492,11 @@ add_functions (void)
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,
......
gcc/fortran/misc.c
......@@ -309,7 +309,6 @@ gfc_done_1 (void)
gfc_scanner_done_1 ();
gfc_intrinsic_done_1 ();
gfc_simplify_done_1 ();
gfc_iresolve_done_1 ();
gfc_arith_done_1 ();
}
......
gcc/fortran/module.c
......@@ -71,6 +71,7 @@ Software Foundation, 59 Temple Place - Suite 330, Boston, MA
#include <time.h>
#include "gfortran.h"
#include "arith.h"
#include "match.h"
#include "parse.h" /* FIXME */
......@@ -519,7 +520,7 @@ gfc_match_use (void)
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)
......@@ -2245,7 +2246,7 @@ mio_gmp_integer (mpz_t * integer)
static void
mio_gmp_real (mpf_t * real)
mio_gmp_real (mpfr_t * real)
{
mp_exp_t exponent;
char *p;
......@@ -2255,14 +2256,14 @@ mio_gmp_real (mpf_t * real)
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);
......@@ -2507,10 +2508,12 @@ mio_expr (gfc_expr ** ep)
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;
......
gcc/fortran/primary.c
/* 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.
......@@ -436,7 +436,7 @@ done:
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)
{
......@@ -497,7 +497,7 @@ done:
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:
......@@ -1076,12 +1076,12 @@ done:
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--;
......
gcc/fortran/simplify.c
......@@ -30,9 +30,6 @@ Software Foundation, 59 Temple Place - Suite 330, Boston, MA
#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;
......@@ -148,7 +145,7 @@ gfc_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;
......@@ -166,7 +163,7 @@ gfc_simplify_abs (gfc_expr * e)
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;
......@@ -174,17 +171,17 @@ gfc_simplify_abs (gfc_expr * e)
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;
......@@ -231,12 +228,11 @@ gfc_expr *
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);
......@@ -245,33 +241,7 @@ gfc_simplify_acos (gfc_expr * x)
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");
}
......@@ -370,7 +340,7 @@ gfc_simplify_aimag (gfc_expr * e)
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");
}
......@@ -391,7 +361,7 @@ gfc_simplify_aint (gfc_expr * e, gfc_expr * k)
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);
......@@ -410,7 +380,7 @@ gfc_simplify_dint (gfc_expr * e)
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);
......@@ -425,6 +395,7 @@ gfc_simplify_anint (gfc_expr * e, gfc_expr * k)
{
gfc_expr *rtrunc, *result;
int kind, cmp;
mpfr_t half;
kind = get_kind (BT_REAL, k, "ANINT", e->ts.kind);
if (kind == -1)
......@@ -437,22 +408,27 @@ gfc_simplify_anint (gfc_expr * e, gfc_expr * k)
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");
}
......@@ -463,6 +439,7 @@ gfc_simplify_dnint (gfc_expr * e)
{
gfc_expr *rtrunc, *result;
int cmp;
mpfr_t half;
if (e->expr_type != EXPR_CONSTANT)
return NULL;
......@@ -472,22 +449,27 @@ gfc_simplify_dnint (gfc_expr * e)
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");
}
......@@ -497,12 +479,11 @@ gfc_expr *
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);
......@@ -511,32 +492,7 @@ gfc_simplify_asin (gfc_expr * x)
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");
}
......@@ -552,7 +508,7 @@ gfc_simplify_atan (gfc_expr * x)
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");
......@@ -569,17 +525,16 @@ gfc_simplify_atan2 (gfc_expr * y, gfc_expr * x)
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");
......@@ -635,8 +590,8 @@ gfc_simplify_ceiling (gfc_expr * e, gfc_expr * k)
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);
......@@ -684,21 +639,21 @@ simplify_cmplx (const char *name, gfc_expr * x, gfc_expr * y, int kind)
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:
......@@ -710,11 +665,11 @@ simplify_cmplx (const char *name, gfc_expr * x, gfc_expr * y, int kind)
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:
......@@ -752,7 +707,7 @@ gfc_simplify_conjg (gfc_expr * e)
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");
}
......@@ -762,7 +717,7 @@ gfc_expr *
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;
......@@ -772,23 +727,24 @@ gfc_simplify_cos (gfc_expr * x)
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");
......@@ -809,7 +765,7 @@ gfc_simplify_cosh (gfc_expr * x)
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");
}
......@@ -902,15 +858,15 @@ gfc_simplify_dim (gfc_expr * x, gfc_expr * y)
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;
......@@ -925,7 +881,7 @@ gfc_simplify_dim (gfc_expr * x, gfc_expr * y)
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;
......@@ -933,13 +889,13 @@ gfc_simplify_dprod (gfc_expr * x, gfc_expr * y)
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");
}
......@@ -957,7 +913,7 @@ gfc_simplify_epsilon (gfc_expr * e)
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");
}
......@@ -967,86 +923,30 @@ gfc_expr *
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:
......@@ -1056,11 +956,11 @@ gfc_simplify_exp (gfc_expr * x)
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)
......@@ -1069,31 +969,39 @@ gfc_simplify_exponent (gfc_expr * x)
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);
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);
mpz_set_f (result->value.integer, lnx);
mpfr_div (lnx, lnx, ln2, GFC_RND_MODE);
mpfr_trunc (lnx, lnx);
mpfr_add_ui (lnx, lnx, 1, GFC_RND_MODE);
mpf_clear (i2);
mpf_clear (ln2);
mpf_clear (lnx);
mpf_clear (absv);
gfc_mpfr_to_mpz (result->value.integer, lnx);
mpfr_clear (i2);
mpfr_clear (ln2);
mpfr_clear (lnx);
mpfr_clear (absv);
mpfr_clear (zero);
return range_check (result, "EXPONENT");
}
......@@ -1116,7 +1024,7 @@ gfc_expr *
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 ());
......@@ -1128,10 +1036,13 @@ gfc_simplify_floor (gfc_expr * e, gfc_expr * k)
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");
}
......@@ -1141,7 +1052,7 @@ gfc_expr *
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)
......@@ -1149,37 +1060,44 @@ gfc_simplify_fraction (gfc_expr * x)
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");
}
......@@ -1204,7 +1122,7 @@ gfc_simplify_huge (gfc_expr * e)
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:
......@@ -1605,16 +1523,16 @@ gfc_simplify_int (gfc_expr * e, gfc_expr * k)
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;
......@@ -1642,8 +1560,8 @@ gfc_simplify_ifix (gfc_expr * e)
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");
......@@ -1663,8 +1581,8 @@ gfc_simplify_idint (gfc_expr * e)
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");
......@@ -2000,52 +1918,60 @@ gfc_expr *
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;
......@@ -2061,21 +1987,28 @@ gfc_expr *
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");
}
......@@ -2142,9 +2075,10 @@ simplify_min_max (gfc_expr * expr, int sign)
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;
......@@ -2235,7 +2169,7 @@ gfc_expr *
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;
......@@ -2256,7 +2190,7 @@ gfc_simplify_mod (gfc_expr * a, gfc_expr * p)
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);
......@@ -2264,18 +2198,19 @@ gfc_simplify_mod (gfc_expr * a, gfc_expr * p)
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:
......@@ -2290,7 +2225,7 @@ gfc_expr *
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;
......@@ -2313,7 +2248,7 @@ gfc_simplify_modulo (gfc_expr * a, gfc_expr * p)
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);
......@@ -2321,19 +2256,20 @@ gfc_simplify_modulo (gfc_expr * a, gfc_expr * p)
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:
......@@ -2376,7 +2312,7 @@ gfc_simplify_nearest (gfc_expr * x, gfc_expr * s)
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.;
......@@ -2387,32 +2323,32 @@ gfc_simplify_nearest (gfc_expr * x, gfc_expr * s)
/* 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
......@@ -2432,6 +2368,7 @@ simplify_nint (const char *name, gfc_expr * e, gfc_expr * k)
{
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)
......@@ -2445,25 +2382,30 @@ simplify_nint (const char *name, gfc_expr * e, gfc_expr * k)
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);
}
......@@ -2937,7 +2879,7 @@ gfc_expr *
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;
......@@ -2952,41 +2894,48 @@ gfc_simplify_rrspacing (gfc_expr * x)
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");
}
......@@ -2996,7 +2945,7 @@ gfc_expr *
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)
......@@ -3004,9 +2953,9 @@ gfc_simplify_scale (gfc_expr * x, gfc_expr * i)
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;
}
......@@ -3035,17 +2984,19 @@ gfc_simplify_scale (gfc_expr * x, gfc_expr * i)
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");
}
......@@ -3195,7 +3146,7 @@ gfc_expr *
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)
......@@ -3203,44 +3154,51 @@ gfc_simplify_set_exponent (gfc_expr * x, gfc_expr * i)
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");
}
......@@ -3352,9 +3310,9 @@ gfc_simplify_sign (gfc_expr * x, gfc_expr * y)
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;
......@@ -3370,7 +3328,7 @@ gfc_expr *
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;
......@@ -3380,23 +3338,24 @@ gfc_simplify_sin (gfc_expr * x)
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:
......@@ -3417,7 +3376,7 @@ gfc_simplify_sinh (gfc_expr * x)
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");
}
......@@ -3443,7 +3402,7 @@ gfc_expr *
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;
......@@ -3459,48 +3418,56 @@ gfc_simplify_spacing (gfc_expr * x)
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");
}
......@@ -3510,7 +3477,7 @@ gfc_expr *
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;
......@@ -3520,9 +3487,9 @@ gfc_simplify_sqrt (gfc_expr * e)
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;
......@@ -3530,82 +3497,83 @@ gfc_simplify_sqrt (gfc_expr * e)
/* 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;
......@@ -3625,9 +3593,8 @@ negative_arg:
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;
......@@ -3638,37 +3605,7 @@ gfc_simplify_tan (gfc_expr * x)
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");
}
......@@ -3678,23 +3615,13 @@ gfc_expr *
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");
......@@ -3712,7 +3639,7 @@ gfc_simplify_tiny (gfc_expr * e)
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;
}
......@@ -3988,21 +3915,5 @@ void
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);
}
gcc/fortran/trans-const.c
......@@ -234,7 +234,7 @@ gfc_conv_mpz_to_tree (mpz_t i, int kind)
/* 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;
......@@ -251,13 +251,9 @@ gfc_conv_mpf_to_tree (mpf_t f, int kind)
}
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)
{
......@@ -265,8 +261,11 @@ gfc_conv_mpf_to_tree (mpf_t f, int kind)
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);
......@@ -294,6 +293,7 @@ gfc_conv_mpf_to_tree (mpf_t f, int kind)
type = gfc_get_real_type (kind);
res = build_real (type, REAL_VALUE_ATOF (q, TYPE_MODE (type)));
gfc_free (q);
gfc_free (p);
......@@ -321,16 +321,16 @@ gfc_conv_constant_to_tree (gfc_expr * expr)
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);
......
gcc/fortran/trans-const.h
......@@ -23,7 +23,7 @@ Software Foundation, 59 Temple Place - Suite 330, Boston, MA
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
......
gcc/fortran/trans-intrinsic.c
......@@ -33,9 +33,9 @@ Software Foundation, 59 Temple Place - Suite 330, Boston, MA
#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"
......@@ -308,7 +308,7 @@ gfc_conv_intrinsic_aint (gfc_se * se, gfc_expr * expr, int op)
tree arg;
tree tmp;
tree cond;
mpf_t huge;
mpfr_t huge;
int n;
int kind;
......@@ -363,14 +363,15 @@ gfc_conv_intrinsic_aint (gfc_se * se, gfc_expr * expr, int op)
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);
......@@ -378,6 +379,7 @@ gfc_conv_intrinsic_aint (gfc_se * se, gfc_expr * expr, int op)
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);
}
......@@ -777,7 +779,7 @@ gfc_conv_intrinsic_mod (gfc_se * se, gfc_expr * expr, int modulo)
tree zero;
tree test;
tree test2;
mpf_t huge;
mpfr_t huge;
int n;
arg = gfc_conv_intrinsic_function_args (se, expr);
......@@ -799,14 +801,15 @@ gfc_conv_intrinsic_mod (gfc_se * se, gfc_expr * expr, int modulo)
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);
......@@ -816,6 +819,7 @@ gfc_conv_intrinsic_mod (gfc_se * se, gfc_expr * expr, int modulo)
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:
......@@ -1423,7 +1427,7 @@ gfc_conv_intrinsic_minmaxloc (gfc_se * se, gfc_expr * expr, int op)
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:
......@@ -1564,7 +1568,7 @@ gfc_conv_intrinsic_minmaxval (gfc_se * se, gfc_expr * expr, int op)
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:
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment