This is the mail archive of the
`gcc@gcc.gnu.org`
mailing list for the GCC project.

Index Nav: | [Date Index] [Subject Index] [Author Index] [Thread Index] | |
---|---|---|

Message Nav: | [Date Prev] [Date Next] | [Thread Prev] [Thread Next] |

Other format: | [Raw text] |

*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?

Hi all, 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 (mainline). 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 that)? 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 point arithmetics. So I'd like your comments: can I do this? Better ideas? Greetz Steven ---- 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. */ tree natural_logarithm (tree arg) { tree type; 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)) return arg; /* 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); p--; } 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); p++; } /* 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 */ } #ifdef G95_DEBUG /* Sanity check */ if REAL_VALUE_ISNAN (log) abort(); /* Series did not converge??? */ #endif /* Build the return tree and add the log (e^p) = p. */ if (p >= 0) REAL_VALUE_FROM_INT (tmp, p, 0, TYPE_MODE (type)); else 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. */ tree common_logarithm (tree arg) { tree type, log10; type = TREE_TYPE (arg); log10 = natural_logarithm (build_real_from_int_cst(type, build_int_2(10,0))); return fold (build (RDIV_EXPR, type, natural_logarithm (arg), log10)); }

**Follow-Ups**:**Re: G95 Floating point arithmetics without GMP?***From:*Zack Weinberg

Index Nav: | [Date Index] [Subject Index] [Author Index] [Thread Index] | |
---|---|---|

Message Nav: | [Date Prev] [Date Next] | [Thread Prev] [Thread Next] |