FW: lo_ieee_init: unrecognized floating point format! on ARM platform

[John W. Eaton]

On 22-Sep-2005, Simon Pickering wrote:

| >From looking at the code, it looks like there's a set of test routines which
| try to work out the fp format for whatever platform the code is
| generated/run on(?) and 4 sets of test data against which this is compared
| to determine which fp type it matches. I'll have to have a look through this
| code which is used to determine the fp type; the numbers hard-coded in the
| oct_mach_info::flt_fmt_* types are presumably pre-calculated somehow (on a
| known platform from the fp type test code I guess

Yes, but the constants in that file can be determined by knowing the
floating point format.  The four numbers used for comparison are

   xmin    - the smallest non-vanishing normalized floating-point
             power of the radix, i.e.,  xmin = FLOAT(ibeta)**minexp
             (d1mach(1))

   xmax    - the largest finite floating-point number.  In
             particular  xmax = (1.0-epsneg)*FLOAT(ibeta)**maxexp
             Note - on some machines  xmax  will be only the
             second, or perhaps third, largest number, being
             too small by 1 or 2 units in the last digit of
             the significand.
             (d1mach(2))

   epsneg  - A small positive floating-point number such that
             1.0-epsneg .NE. 1.0. In particular, if ibeta = 2
             or  IRND = 0, epsneg = FLOAT(ibeta)**negeps.
             Otherwise,  epsneg = (ibeta**negeps)/2.  Because
             negeps is bounded below by -(it+3), epsneg may not
             be the smallest number that can alter 1.0 by
             subtraction.
             (d1mach(3))

   eps     - the smallest positive floating-point number such
             that  1.0+eps .NE. 1.0. In particular, if either
             ibeta = 2  or  IRND = 0, eps = FLOAT(ibeta)**machep.
             Otherwise,  eps = (FLOAT(ibeta)**machep)/2
             (d1mach(4))

See d1mach.f and machar.c files in the libcruft/misc source
directory.

| Here are the results in decimal with 'format long':
| 
| In decimal:
| Eps =     2.22044604925031e-16
| Realmin = 2.22507385850720e-308
| Realmax = 1.79769313486232e+308

These and the bit patterns for them seem to be correct, so I'm
guessing that the epsneg calculation is wrong somehow.  Maybe this is
a compiler problem, or a bug in machar.c?

Looking a little closer, I see that these numbers are defined in
Octave using the C macros DBL_MAX, DBL_MIN, and DBL_EPSILON (in
src/data.cc).  But the comparison that is failing is based on computed
values.  So what happens if you try comiling the attached file with

  mkoctfile d1mach

and then at the Octave prompt,

  format bit
  [xmin, xmax, epsneg, eps] = d1mach?

On your system (apparently IEEE big endian), I think the results
should be

  xmin = 0000000000000000000000000000000000000000000100000000000000000000
  xmax = 1111111111111111111111111111111101111111111011111111111111111111
  epsneg = 0000000000000000000000000000000000111100101000000000000000000000
  eps = 0000000000000000000000000000000000111100101100000000000000000000

If they are not (and I'd guess not, since the comparison is failing
when Octave starts), then you'll need to find out why the computation
of these constants is failing in machar.c.  My guess is that compiling
that routine with -ffloat-store will avoid the problem, but that is
just a guess.

jwe

#include <octave/oct.h>
#include <octave/f77-fcn.h>
#include <octave/oct-map.h>

extern "C"
{
  F77_RET_T
  F77_FUNC (machar, MACHAR) (double& xmin, double& xmax, double& epsneg,
                             double& eps, double& log10_ibeta);
}

DEFUN_DLD (d1mach, args, , "d1mach")
{
  octave_value_list retval;

  double xmin, xmax, epsneg, eps, log10_ibeta;

  F77_FUNC (machar, MACHAR) (xmin, xmax, epsneg, eps, log10_ibeta);
  
  retval(4) = log10_ibeta;
  retval(3) = eps;
  retval(2) = epsneg;
  retval(1) = xmax;
  retval(0) = xmin;

  return retval;
}