This is the mail archive of the
mailing list for the GCC project.
G95 Floating point arithmetics without GMP?
- From: Steven Bosscher <s dot bosscher at student dot tudelft dot nl>
- To: gcc at gcc dot gnu dot org
- Date: 08 May 2002 20:46:52 +0200
- Subject: G95 Floating point arithmetics without GMP?
I'm trying to figure out if/how we can do G95's arithmetics without GMP.
Without GMP, it would be much earier to rewrite g95 to use GCC's 'tree'
structure in the parser, instead of a dozen++ different "native"
structures. So, I've started playing with GCC REAL_ARITHMETICS
I've rewritten the first to functions of arith.c to use GCC middle-end
stuff. But REAL_ARITHMETICS doesn't provide all the arithmetics
operations we need.
For example, g95 can simplify exponentiation expressions (A**B) but
since C doesn't have that operator, REAL_ARITHMETICS doesn't allow me to
use eldexp() in real.c (there's no POWER_EXPR tree code).
eldexp(x,pwr2,y) returns Y = X * 2^PWR2, which is exactly what I need.
Is it safe to declare this function as extern and just use it, or should
I add an extra macro to real.h (say, REAL_VALUE_2EXP, or something like
Another problem I ran into, is that there's no function in real.c to
calculate log(x). G95 uses McLaurin series for log, and also for sine,
cosine, and arctangent. I've tried to do this with REAL_ARITHMETICS, but
I'm not sure if this is correct. I'm not really familiar with floating
So I'd like your comments: can I do this? Better ideas?
---- excerpt from my arith.c ----
/* Ugly hack to make it work: x = x+0. Strictly speaking, this is in
fact not true for all REAL_VALUE_TYPEs, but we do it anyway. */
#define REAL_COPY(from, to) \
REAL_VALUE_ARITH (to, PLUS_EXPR, &(from), &(dconst0))
/* Compute the natural logarithm of ARG, which must be a constant
expression of a REAL_TYPE.
We use the series:
ln(x) = (x-1) - (x-1)^/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. */
natural_logarithm (tree arg)
REAL_VALUE_TYPE x, t, xp, log, tmp;
HOST_WIDE_INT i, p;
type = TREE_TYPE (arg);
x = TREE_REAL_CST (fold (arg));
/* If the argument is not a number, forget about trying to
get a sane log value for it. */
if (REAL_VALUE_ISNAN (x))
/* Get the argument into the range 0.5 (1/2) to 1.5 (3/2) by
successive multiplications or divisions by e. */
p = 0;
REAL_ARITHMETIC (t, RDIV_EXPR, dconst1, dconst2); /* t = 0.5 */
while (REAL_VALUES_LESS (x, t)) /* while (x < .5) x=x*e */
REAL_ARITHMETIC (x, MULT_EXPR, type, x, e);
REAL_VALUE_FROM_INT (tmp, 3, 0, TYPE_MODE (type));
REAL_ARITHMETIC (t, RDIV_EXPR, tmp, dconst2); /* t = 1.5 */
while (!REAL_VALUES_LESS (x, t)) /* while not (x < 1.5) x=x/e */
REAL_ARITHMETIC (x, RDIV_EXPR, x, e);
/* Compute the series. */
REAL_COPY (dconst0 log);
REAL_COPY (dconst1, xp);
REAL_ARITHMETIC (x, MINUS_EXPR, x, xp);
for (i = 1; i < TYPE_PRECISION (type); i++) /* Is this sane? */
REAL_VALUE_FROM_INT (tmp, i, 0, TYPE_MODE (type));
REAL_ARITHMETIC (xp, MULT_EXPR, xp, x); /* xp = x^i */
REAL_ARITHMETIC (t, RDIV_EXPR, xp, tmp); /* t = (x^i)/i */
REAL_ARITHMETIC (log, (i % 2 == 0) ?
MINUS_EXPR : PLUS_EXPR, log, t); /* sum */
/* Sanity check */
if REAL_VALUE_ISNAN (log)
abort(); /* Series did not converge??? */
/* Build the return tree and add the log (e^p) = p. */
if (p >= 0)
REAL_VALUE_FROM_INT (tmp, p, 0, TYPE_MODE (type));
REAL_VALUE_FROM_INT (tmp, p, -1, TYPE_MODE (type));
return fold (build (PLUS_EXPR, type, build_real(type,log), tmp));
/* Compute the common logarithm of ARG, which must be a REAL_TYPE. */
common_logarithm (tree arg)
tree type, log10;
type = TREE_TYPE (arg);
return fold (build (RDIV_EXPR, type, natural_logarithm (arg), log10));