G95 Floating point arithmetics without GMP?

G95 Floating point arithmetics without GMP?

[Steven Bosscher]

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

Re: G95 Floating point arithmetics without GMP?

[Zack Weinberg]

On Wed, May 08, 2002 at 08:46:52PM +0200, Steven Bosscher wrote:
> 
> 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).

You should add one.

In general, do not hesitate to add tree codes to tree.def when they
can reasonably express a language-independent concept.  POWER_EXPR is
one such.

> 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.

Hmm, it would be better to add such routines to real.c and have them
operate directly on the internal floating-point representation.  You'd
need LOG_EXPR, SIN_EXPR, COS_EXPR, ARCTAN_EXPR, but that's okay.

I can't help much with details of numerical analysis, sorry. 

zw

Re: G95 Floating point arithmetics without GMP?

[Steven Bosscher]

Op wo 08-05-2002, om 21:15 schreef Zack Weinberg:
> On Wed, May 08, 2002 at 08:46:52PM +0200, Steven Bosscher wrote:
> > 
> > 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).
> 
> You should add one.
> 
> In general, do not hesitate to add tree codes to tree.def when they
> can reasonably express a language-independent concept.  POWER_EXPR is
> one such.

I defined it in f95-tree.def because I thought it shouldn't be in
tree.def.

I think that right now GCC can generate code for all tree codes defined
in tree.def (right?). If I add POWER_EXPR, I expect we'd need to add
support for it in the back-end as well? For G95, I wanted to replace it
with a libcall.

Greetz
Steven

Re: G95 Floating point arithmetics without GMP?

[Richard Henderson]

On Wed, May 08, 2002 at 10:10:11PM +0200, Steven Bosscher wrote:
> If I add POWER_EXPR, I expect we'd need to add support for it in
> the back-end as well? For G95, I wanted to replace it with a libcall.

Well, you'd add support for it in the tree->rtl converter.
You'd do this by initializing an "optab" with the functions
to be called.

And really, this is exactly what you want, since you'd like
POWER_EXPR to survive through the early tree optimizations.


r~

Re: G95 Floating point arithmetics without GMP?

[law]

In message <1020888612.759.77.camel@steven>, Steven Bosscher writes:
 > Op wo 08-05-2002, om 21:15 schreef Zack Weinberg:
 > > On Wed, May 08, 2002 at 08:46:52PM +0200, Steven Bosscher wrote:
 > > > 
 > > > 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).
 > > 
 > > You should add one.
 > > 
 > > In general, do not hesitate to add tree codes to tree.def when they
 > > can reasonably express a language-independent concept.  POWER_EXPR is
 > > one such.
 > 
 > I defined it in f95-tree.def because I thought it shouldn't be in
 > tree.def.
It's a generic enough concept that it probably should migrate into
the toplevel tree.def.

 > I think that right now GCC can generate code for all tree codes defined
 > in tree.def (right?).
Not always.  It's dependent on the processor; if we don't have an instruction
for the operation and we can't open-code it using other primitives, then 
we call out to a library.

jeff

Re: G95 Floating point arithmetics without GMP?

[Joseph S. Myers]

On Wed, 8 May 2002 law@redhat.com wrote:

>  > > In general, do not hesitate to add tree codes to tree.def when they
>  > > can reasonably express a language-independent concept.  POWER_EXPR is
>  > > one such.
>  > 
>  > I defined it in f95-tree.def because I thought it shouldn't be in
>  > tree.def.
> It's a generic enough concept that it probably should migrate into
> the toplevel tree.def.

There used to be something called EXPON_EXPR as a placeholder in tree.def.

revision 1.41
date: 2001/07/15 02:16:34;  author: rth;  state: Exp;  lines: +0 -4
        * tree.def (EXPON_EXPR) remove. Never supported anyway.

        * f/com.c (ffecom_overlap_): Remove references to EXPON_EXPR.
        (ffecom_tree_canonize_ref_): Likewise.

        * java/check-init.c (check_init): Remove references to EXPON_EXPR.

-- 
Joseph S. Myers
jsm28@cam.ac.uk

Re: G95 Floating point arithmetics without GMP?

[law]

In message <Pine.LNX.4.33.0205082336590.12048-100000@kern.srcf.societies.cam.ac
.uk>, "Joseph S. Myers" writes:
 > On Wed, 8 May 2002 law@redhat.com wrote:
 > 
 > >  > > In general, do not hesitate to add tree codes to tree.def when they
 > >  > > can reasonably express a language-independent concept.  POWER_EXPR is
 > >  > > one such.
 > >  > 
 > >  > I defined it in f95-tree.def because I thought it shouldn't be in
 > >  > tree.def.
 > > It's a generic enough concept that it probably should migrate into
 > > the toplevel tree.def.
 > 
 > There used to be something called EXPON_EXPR as a placeholder in tree.def.
And your point is?

If this is something the g95 folks will need/use, then it really
belongs in tree.def.  It was probably removed because nobody actually
used it.

jeff

Re: G95 Floating point arithmetics without GMP?

[Toon Moene]

law@redhat.com wrote:

> If this is something the g95 folks will need/use, then it really
> belongs in tree.def.  It was probably removed because nobody actually
> used it.

It was probably removed based on the analysis I provided on March, 12,
2001.

	http://gcc.gnu.org/ml/gcc/2001-03/msg00412.html

If this analysis is wrong, the change might have to be reverted.

Please advise,

-- 
Toon Moene - mailto:toon@moene.indiv.nluug.nl - phoneto: +31 346 214290
Saturnushof 14, 3738 XG  Maartensdijk, The Netherlands
Maintainer, GNU Fortran 77: http://gcc.gnu.org/onlinedocs/g77_news.html
Join GNU Fortran 95: http://g95.sourceforge.net/ (under construction)