G95 Floating point arithmetics without GMP?

[Steven Bosscher <s dot bosscher at student dot tudelft dot nl>]

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));
}