public inbox for gcc@gcc.gnu.org
 help / color / mirror / Atom feed
* G95 Floating point arithmetics without GMP?
@ 2002-05-08 13:05 Steven Bosscher
  2002-05-08 13:11 ` Zack Weinberg
  0 siblings, 1 reply; 8+ messages in thread
From: Steven Bosscher @ 2002-05-08 13:05 UTC (permalink / raw)
  To: gcc

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




^ permalink raw reply	[flat|nested] 8+ messages in thread
end of thread, other threads:[~2002-05-09 11:20 UTC | newest]

Thread overview: 8+ messages (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2002-05-08 13:05 G95 Floating point arithmetics without GMP? Steven Bosscher
2002-05-08 13:11 ` Zack Weinberg
2002-05-08 13:18   ` Steven Bosscher
2002-05-08 13:38     ` Richard Henderson
2002-05-08 14:19     ` law
2002-05-08 16:15       ` Joseph S. Myers
2002-05-08 18:10         ` law
2002-05-09  5:35           ` Toon Moene

This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).