Re: Anyone working on 64-bit octave

[Thomas D. Dean]

On 02/17/12 12:25, Jordi Gutiérrez Hermoso wrote:
> On 17 February 2012 15:12, Matyas Sustik<address@hidden>  wrote:
> > > I am looking at building BLAS and LAPACK with
> > >
> > >     gfortran -fdefault-integer-8 and -default-real-8
> > >
> > > This makes double precision 128 bits.
> >
> > My search turned up the following that may help:
> >
> > http://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=2&t=2739
> >
> > If you got it working, let me (us) know...
>
> Note that even if you do get this working, essentially all of the
> Octave source assumes double precision. There are no quadruple
> precision types in Octave.

More investigation shows this to be a really hard problem.  The state of 
computing with household PC's is not equal to the task.  Maybe in a few 
years when compilers catch up with hardware.

sizeof(long double) returns 16(bytes) == 128-bits.  But, the gcc libs, 
etc., change this into 80-bits precision.

The __float128 type will cause gcc46 to issue 128-bit instructions. 
But, to do i/o, special function calls must be made.  And, libquadmath 
linked.

I fear this would lead to non-standard variable definitions in the C 
code and lead to a non-supportable situation when compilers catch up.

On FreeBSD and most likely any system that uses gcc46 or earlier:

As it turns out, gcc46 does not support the full precision of long 
double.  Or, is it just because the port lang/gcc46 uses the system 
libc, which is gcc 4.2.1?

#include <quadmath.h>
#include <stdio.h>
int main() {
  __float128 f128;
  long double ld;
  char s[512];

  f128 = (__float128)1.0;
  f128 += ((__float128)2.0) * FLT128_EPSILON;
  ld = (long double)f128;

  quadmath_snprintf(s, sizeof s, "%.70Qe", f128);
  printf("%s\n%.70Le\n",s,ld);
  if (ld != f128) printf("not equal\n");

  return 0;
}

> env | grep LD
LD_LIBRARY_PATH=/usr/local/lib:/usr/local/lib/gcc46:/lib
> which gcc46
/usr/local/bin/gcc46
> locate libquadmath.so
/usr/local/lib/gcc46/libquadmath.so
/usr/local/lib/gcc46/libquadmath.so.0

> gcc46 -march=corei7 -msse2 -Wall quadtest.c -o quadtest -lquadmath
> ./quadtest
1.0000000000000000000000000000000003851859888774471706111955885169854637e+00
1.0000000000000000000000000000000000000000000000000000000000000000000000e+00
not equal

Oops!

The conversion from __float128 to long double uses 80-bit floating point 
registers.

In reality, long double does not provide 16 bytes precision as reported 
by sizeof(), but, is 80-bits, like the original 386 FPU's.

gcc docs, in discussing 128-bit floats:

Notice that neither of these options enable any extra precision over the 
x87 standard of 80 bits for a `long double'.

Tom Dean