Re: REAL(16) ...on x86/x86-64

[Steve Kargl <sgk at troutmask dot apl dot washington dot edu>]

On Mon, Aug 25, 2008 at 07:09:12PM -0700, Jerry DeLisle wrote:
> Steve Kargl wrote:
> >On Mon, Aug 25, 2008 at 03:15:42PM -0700, Steve Kargl wrote:
> >>On Mon, Aug 25, 2008 at 11:50:16PM +0200, Dominique Dhumieres wrote:
> >>>>REAL(16) needs to be done in software -- on x86, x86-64 -- as it is not
> >>>>supported in hardware; if you want to use more than REAL(8) on x86,
> >>>>x86-64 you can use REAL(10).  Or you use a system such as PowerPC which
> >>>>supports REAL(16) in silicon.
> >>>As far as I can tell, the ifc implementation is a "real" one where the 
> >>>full 128 bits are used to code the real number.  So far i did not have 
> >>>the time to play with it.
> >>>
> >>>Now concerning gfortran and since gcc requires gmp and mpfr, how 
> >>>difficult (efficient) would it be to use these libraries to implement 
> >>>REAL(xxx)?
> >>>
> >>Getting +, -, *, and / working is almost trivial.  The hard part is 
> >>getting REAL(16) working in all the fun corners of the standard 
> >>as well as getting IO working.  Here's a start where one may need
> >>to distinguish between hardware FP and software emulation.
> >>
> >>
> >
> >With the patch I sent earlier with an off-by-1 fix in exponents, I get
> >
> >     Kind: 4 8 10 16
> >Precision: 6 15 18 33
> >   Digits: 24 53 64 113
> >    Radix: 2 2 2 2
> >Min. exp.: -125 -1021 -16381 -16381
> >Max. exp.: 128 1024 16384 16384
> >
> >on x86_64-*-freebsd.
> >
> This looks like a nice start.  We already have support for real(16) I/O 
> working so that will be a minor piece.

I think you'ld need to adapt the IO for software emulated FP. 
Everything would remain an mpfr_t type, so you (this is a
generic you) need to transfer mpfr_t during IO.

> There exists float-128 routines within the gcc libraries now.  We need to 
> figure out how to adapt/configure and use these.

I had the impression some would like to see real.c die.

> Steve, will you carry this patch forward to get +,-,*, and / working?
 
I don't have any plans to carry this any further.  It was just a pointer
to others where to start.  I'll note that if the REAL(16) stuff is 
completely implemented within framework of mpfr, then one can implement
-freal-bits=255,15.  This then can be used to define a real type
with 255 bits of precision and an exponent range of 2**15-3:2**15-1.
In trans-type, one does

  if (option.real_bits) then
  else
     what I already posted.
  end if

> In the meantime, is there someone knowledgeable of the flt-128 library that 
> can guide us in this area.  Maybe between Steve and I we can get it working 
> with some mentoring on the configuration stuff. (target gcc 4.5)

Well, my idea (based on Dominique's comment) is to keep everything
in a mpfr_t type.  For example, in the following code,

   real(16) x, y
   x = 1._16
   y = 2._16
   x = x + y
   end

The binary operator + ultimately gets translated to 

   mpfr_t tmp;
   mpfr_init(tmp, 113);
   mpfr_add(tmp, x, y, RND);
   mpfr_set(x, tmp, RND);
   mpfr_clear(tmp)

BTW, constant folding already works.

   program probe
   real(16) d, e
   d = 1.e0_16 + 3.14e1_16
   e = 2. / 3._16
   d =  4._16 * 3.3_16
   e = 42. - 8.4_16
   d = - 3._16
   d = sin(1.2345_16)
   a = sin(1.2345)
   end program probe

troutmask:sgk[204] gfc4x -o z -fdump-tree-original a.f90
troutmask:sgk[205] tail a.f90.003t.original 
  d = 3.23999999999999999999999999999999987674048355921690540442e+1;
  e = 6.666666666666666666666666666666666345678342602127357824e-1;
  d = 1.31999999999999999999999999999999993837024177960845270221e+1;
  e = 3.36000000000000000000000000000000012325951644078309459558e+1;
  d = -3.0e+0;
  d = 9.43983323944511166155703552452566891228124607439110465785e-1;
  a = 9.439833164215087890625e-1;

Now, for something like a very basic binary operator, one needs
to convert a = b + c into similar mpfr code.  This can be done in
trans-expr.c(gfc_conv_expr_op) by short circuiting the logic. 

  if (expr->value.op.op1->ts.kind == 16 || expr->value.op.op2->ts.kind == 16)
    {
      gfc_error_now ("Tsk! Tsk!  Need to convert op into mpfr call");
    }

-- 
Steve