[Insight-developers] vnl_math changes

Simon Warfield warfield at bwh.harvard.edu
Wed Aug 11 22:38:21 EDT 2004


I modified some of the vnl and vcl files in order to get the 
vnl_test_math to pass all of its test on the intel compiler.

There were a few  issues. One was to modify the code to be aware of the 
intel compiler, and conditionally compile the appropriate code.
Another was that the code in vnl_numeric_limits.cxx that defines e.g. 
NaN was using a bit pattern not correct for the Intel layout.
In particular, the Intel layout has an 80 bit long double in a 12 byte 
aligned structure, and so byte 9 needs to be set to 0x7f, not byte 11 as 
the code was doing.
Then, there are some issues with isnanl, isinfl etc. with the Intel 
compiler up to version 8.0.  Apparently these are resolved in version 
8.1 of the compiler, so the code calls the correct builtin functions for 
that or newer versions, and for 8.0 and earlier it makes a good try at 
working around the compiler problem (its pretty ugly in there).

On my system, with the changes in the attached files it now passes all 
of the tests in vnl_test_math.

-- 
Simon 

-------------- next part --------------
// This is core/vnl/vnl_numeric_limits.cxx
#ifdef VCL_NEEDS_PRAGMA_INTERFACE
#pragma implementation
#endif
//
// numeric_limits
// Author: Andrew W. Fitzgibbon, Oxford RRG
// Created: 28 Aug 96
// \deprecated in favour of vcl_limits.
// \verbatim
// Modifications
// IMS (Manchester) 21/10/2003: Deprecated - Can be deleted after VXL-1.1
// \endverbatim
//
//-----------------------------------------------------------------------------

#define vcl_deprecated_header_h_ // A hack to allow this vnl_numeric_limits.cxx to compile without complaint.
#include "vnl_numeric_limits.h"
#include <vxl_config.h> // for VXL_BIG_ENDIAN

// ----------------------------------------------------------------------
// Constants for int

#if !VCL_STATIC_CONST_INIT_INT_NO_DEFN
const bool vnl_numeric_limits<int>::is_specialized VCL_STATIC_CONST_INIT_INT_DEFN(true);
const int vnl_numeric_limits<int>::digits   VCL_STATIC_CONST_INIT_INT_DEFN(31);
const int vnl_numeric_limits<int>::digits10 VCL_STATIC_CONST_INIT_INT_DEFN(9);
const bool vnl_numeric_limits<int>::is_signed  VCL_STATIC_CONST_INIT_INT_DEFN(true);
const bool vnl_numeric_limits<int>::is_integer VCL_STATIC_CONST_INIT_INT_DEFN(true);
const bool vnl_numeric_limits<int>::is_exact   VCL_STATIC_CONST_INIT_INT_DEFN(true);
const int vnl_numeric_limits<int>::radix VCL_STATIC_CONST_INIT_INT_DEFN(2);
const int vnl_numeric_limits<int>::min_exponent   VCL_STATIC_CONST_INIT_INT_DEFN(-31);
const int vnl_numeric_limits<int>::min_exponent10 VCL_STATIC_CONST_INIT_INT_DEFN(-9);
const int vnl_numeric_limits<int>::max_exponent   VCL_STATIC_CONST_INIT_INT_DEFN(31);
const int vnl_numeric_limits<int>::max_exponent10 VCL_STATIC_CONST_INIT_INT_DEFN(9);
const bool vnl_numeric_limits<int>::has_infinity      VCL_STATIC_CONST_INIT_INT_DEFN(false);
const bool vnl_numeric_limits<int>::has_quiet_NaN     VCL_STATIC_CONST_INIT_INT_DEFN(false);
const bool vnl_numeric_limits<int>::has_signaling_NaN VCL_STATIC_CONST_INIT_INT_DEFN(false);
const bool vnl_numeric_limits<int>::has_denorm        VCL_STATIC_CONST_INIT_INT_DEFN(false);
const bool vnl_numeric_limits<int>::is_iec559  VCL_STATIC_CONST_INIT_INT_DEFN(false);
const bool vnl_numeric_limits<int>::is_bounded VCL_STATIC_CONST_INIT_INT_DEFN(true);
const bool vnl_numeric_limits<int>::is_modulo  VCL_STATIC_CONST_INIT_INT_DEFN(true);
const bool vnl_numeric_limits<int>::traps      VCL_STATIC_CONST_INIT_INT_DEFN(false);
const bool vnl_numeric_limits<int>::tinyness_before VCL_STATIC_CONST_INIT_INT_DEFN(false);
const vnl_float_round_style vnl_numeric_limits<int>::round_style VCL_STATIC_CONST_INIT_INT_DEFN(vnl_round_toward_zero);
#endif

// ----------------------------------------------------------------------
// Constants for long

#if !VCL_STATIC_CONST_INIT_INT_NO_DEFN
const bool vnl_numeric_limits<long>::is_specialized VCL_STATIC_CONST_INIT_INT_DEFN(true);
const int vnl_numeric_limits<long>::digits   VCL_STATIC_CONST_INIT_INT_DEFN(31);
const int vnl_numeric_limits<long>::digits10 VCL_STATIC_CONST_INIT_INT_DEFN(9);
const bool vnl_numeric_limits<long>::is_signed  VCL_STATIC_CONST_INIT_INT_DEFN(true);
const bool vnl_numeric_limits<long>::is_integer VCL_STATIC_CONST_INIT_INT_DEFN(true);
const bool vnl_numeric_limits<long>::is_exact   VCL_STATIC_CONST_INIT_INT_DEFN(true);
const int vnl_numeric_limits<long>::radix VCL_STATIC_CONST_INIT_INT_DEFN(2);
const int vnl_numeric_limits<long>::min_exponent   VCL_STATIC_CONST_INIT_INT_DEFN(-31);
const int vnl_numeric_limits<long>::min_exponent10 VCL_STATIC_CONST_INIT_INT_DEFN(-9);
const int vnl_numeric_limits<long>::max_exponent   VCL_STATIC_CONST_INIT_INT_DEFN(31);
const int vnl_numeric_limits<long>::max_exponent10 VCL_STATIC_CONST_INIT_INT_DEFN(9);
const bool vnl_numeric_limits<long>::has_infinity      VCL_STATIC_CONST_INIT_INT_DEFN(false);
const bool vnl_numeric_limits<long>::has_quiet_NaN     VCL_STATIC_CONST_INIT_INT_DEFN(false);
const bool vnl_numeric_limits<long>::has_signaling_NaN VCL_STATIC_CONST_INIT_INT_DEFN(false);
const bool vnl_numeric_limits<long>::has_denorm        VCL_STATIC_CONST_INIT_INT_DEFN(false);
const bool vnl_numeric_limits<long>::is_iec559  VCL_STATIC_CONST_INIT_INT_DEFN(false);
const bool vnl_numeric_limits<long>::is_bounded VCL_STATIC_CONST_INIT_INT_DEFN(true);
const bool vnl_numeric_limits<long>::is_modulo  VCL_STATIC_CONST_INIT_INT_DEFN(true);
const bool vnl_numeric_limits<long>::traps      VCL_STATIC_CONST_INIT_INT_DEFN(false);
const bool vnl_numeric_limits<long>::tinyness_before VCL_STATIC_CONST_INIT_INT_DEFN(false);
const vnl_float_round_style vnl_numeric_limits<long>::round_style VCL_STATIC_CONST_INIT_INT_DEFN(vnl_round_toward_zero);
#endif

// ----------------------------------------------------------------------
// Constants for unsigned long

#if !VCL_STATIC_CONST_INIT_INT_NO_DEFN
const bool vnl_numeric_limits<unsigned long>::is_specialized VCL_STATIC_CONST_INIT_INT_DEFN(true);
const int vnl_numeric_limits<unsigned long>::digits   VCL_STATIC_CONST_INIT_INT_DEFN(sizeof(unsigned long) * 8 );
const int vnl_numeric_limits<unsigned long>::digits10 VCL_STATIC_CONST_INIT_INT_DEFN( (digits * 301) / 1000 );
const bool vnl_numeric_limits<unsigned long>::is_signed  VCL_STATIC_CONST_INIT_INT_DEFN(false);
const bool vnl_numeric_limits<unsigned long>::is_integer VCL_STATIC_CONST_INIT_INT_DEFN(true);
const bool vnl_numeric_limits<unsigned long>::is_exact   VCL_STATIC_CONST_INIT_INT_DEFN(true);
const int vnl_numeric_limits<unsigned long>::radix VCL_STATIC_CONST_INIT_INT_DEFN(2);
const int vnl_numeric_limits<unsigned long>::min_exponent   VCL_STATIC_CONST_INIT_INT_DEFN(-31);
const int vnl_numeric_limits<unsigned long>::min_exponent10 VCL_STATIC_CONST_INIT_INT_DEFN(-9);
const int vnl_numeric_limits<unsigned long>::max_exponent   VCL_STATIC_CONST_INIT_INT_DEFN(31);
const int vnl_numeric_limits<unsigned long>::max_exponent10 VCL_STATIC_CONST_INIT_INT_DEFN(9);
const bool vnl_numeric_limits<unsigned long>::has_infinity      VCL_STATIC_CONST_INIT_INT_DEFN(false);
const bool vnl_numeric_limits<unsigned long>::has_quiet_NaN     VCL_STATIC_CONST_INIT_INT_DEFN(false);
const bool vnl_numeric_limits<unsigned long>::has_signaling_NaN VCL_STATIC_CONST_INIT_INT_DEFN(false);
const bool vnl_numeric_limits<unsigned long>::has_denorm        VCL_STATIC_CONST_INIT_INT_DEFN(false);
const bool vnl_numeric_limits<unsigned long>::is_iec559  VCL_STATIC_CONST_INIT_INT_DEFN(false);
const bool vnl_numeric_limits<unsigned long>::is_bounded VCL_STATIC_CONST_INIT_INT_DEFN(true);
const bool vnl_numeric_limits<unsigned long>::is_modulo  VCL_STATIC_CONST_INIT_INT_DEFN(true);
const bool vnl_numeric_limits<unsigned long>::traps      VCL_STATIC_CONST_INIT_INT_DEFN(false);
const bool vnl_numeric_limits<unsigned long>::tinyness_before VCL_STATIC_CONST_INIT_INT_DEFN(false);
const vnl_float_round_style vnl_numeric_limits<unsigned long>::round_style VCL_STATIC_CONST_INIT_INT_DEFN(vnl_round_toward_zero);
#endif

// ----------------------------------------------------------------------
// Constants for unsigned short

#if !VCL_STATIC_CONST_INIT_INT_NO_DEFN
const bool vnl_numeric_limits<unsigned short>::is_specialized VCL_STATIC_CONST_INIT_INT_DEFN(true);
const int vnl_numeric_limits<unsigned short>::digits   VCL_STATIC_CONST_INIT_INT_DEFN(sizeof(unsigned short) * 8 );
const int vnl_numeric_limits<unsigned short>::digits10 VCL_STATIC_CONST_INIT_INT_DEFN( (digits * 301) / 1000 );
const bool vnl_numeric_limits<unsigned short>::is_signed  VCL_STATIC_CONST_INIT_INT_DEFN(false);
const bool vnl_numeric_limits<unsigned short>::is_integer VCL_STATIC_CONST_INIT_INT_DEFN(true);
const bool vnl_numeric_limits<unsigned short>::is_exact   VCL_STATIC_CONST_INIT_INT_DEFN(true);
const int vnl_numeric_limits<unsigned short>::radix VCL_STATIC_CONST_INIT_INT_DEFN(2);
const int vnl_numeric_limits<unsigned short>::min_exponent   VCL_STATIC_CONST_INIT_INT_DEFN(-31);
const int vnl_numeric_limits<unsigned short>::min_exponent10 VCL_STATIC_CONST_INIT_INT_DEFN(-9);
const int vnl_numeric_limits<unsigned short>::max_exponent   VCL_STATIC_CONST_INIT_INT_DEFN(31);
const int vnl_numeric_limits<unsigned short>::max_exponent10 VCL_STATIC_CONST_INIT_INT_DEFN(9);
const bool vnl_numeric_limits<unsigned short>::has_infinity      VCL_STATIC_CONST_INIT_INT_DEFN(false);
const bool vnl_numeric_limits<unsigned short>::has_quiet_NaN     VCL_STATIC_CONST_INIT_INT_DEFN(false);
const bool vnl_numeric_limits<unsigned short>::has_signaling_NaN VCL_STATIC_CONST_INIT_INT_DEFN(false);
const bool vnl_numeric_limits<unsigned short>::has_denorm        VCL_STATIC_CONST_INIT_INT_DEFN(false);
const bool vnl_numeric_limits<unsigned short>::is_iec559  VCL_STATIC_CONST_INIT_INT_DEFN(false);
const bool vnl_numeric_limits<unsigned short>::is_bounded VCL_STATIC_CONST_INIT_INT_DEFN(true);
const bool vnl_numeric_limits<unsigned short>::is_modulo  VCL_STATIC_CONST_INIT_INT_DEFN(true);
const bool vnl_numeric_limits<unsigned short>::traps      VCL_STATIC_CONST_INIT_INT_DEFN(false);
const bool vnl_numeric_limits<unsigned short>::tinyness_before VCL_STATIC_CONST_INIT_INT_DEFN(false);
const vnl_float_round_style vnl_numeric_limits<unsigned short>::round_style VCL_STATIC_CONST_INIT_INT_DEFN(vnl_round_toward_zero);
#endif

// ----------------------------------------------------------------------
// Constants for short

#if !VCL_STATIC_CONST_INIT_INT_NO_DEFN
const bool vnl_numeric_limits<short>::is_specialized VCL_STATIC_CONST_INIT_INT_DEFN(true);
const int vnl_numeric_limits<short>::digits   VCL_STATIC_CONST_INIT_INT_DEFN(15);
const int vnl_numeric_limits<short>::digits10 VCL_STATIC_CONST_INIT_INT_DEFN(5);
const bool vnl_numeric_limits<short>::is_signed  VCL_STATIC_CONST_INIT_INT_DEFN(true);
const bool vnl_numeric_limits<short>::is_integer VCL_STATIC_CONST_INIT_INT_DEFN(true);
const bool vnl_numeric_limits<short>::is_exact   VCL_STATIC_CONST_INIT_INT_DEFN(true);
const int vnl_numeric_limits<short>::radix VCL_STATIC_CONST_INIT_INT_DEFN(2);
const int vnl_numeric_limits<short>::min_exponent   VCL_STATIC_CONST_INIT_INT_DEFN(-15);
const int vnl_numeric_limits<short>::min_exponent10 VCL_STATIC_CONST_INIT_INT_DEFN(-5);
const int vnl_numeric_limits<short>::max_exponent   VCL_STATIC_CONST_INIT_INT_DEFN(15);
const int vnl_numeric_limits<short>::max_exponent10 VCL_STATIC_CONST_INIT_INT_DEFN(5);
const bool vnl_numeric_limits<short>::has_infinity      VCL_STATIC_CONST_INIT_INT_DEFN(false);
const bool vnl_numeric_limits<short>::has_quiet_NaN     VCL_STATIC_CONST_INIT_INT_DEFN(false);
const bool vnl_numeric_limits<short>::has_signaling_NaN VCL_STATIC_CONST_INIT_INT_DEFN(false);
const bool vnl_numeric_limits<short>::has_denorm        VCL_STATIC_CONST_INIT_INT_DEFN(false);
const bool vnl_numeric_limits<short>::is_iec559  VCL_STATIC_CONST_INIT_INT_DEFN(false);
const bool vnl_numeric_limits<short>::is_bounded VCL_STATIC_CONST_INIT_INT_DEFN(true);
const bool vnl_numeric_limits<short>::is_modulo  VCL_STATIC_CONST_INIT_INT_DEFN(true);
const bool vnl_numeric_limits<short>::traps      VCL_STATIC_CONST_INIT_INT_DEFN(false);
const bool vnl_numeric_limits<short>::tinyness_before VCL_STATIC_CONST_INIT_INT_DEFN(false);
const vnl_float_round_style vnl_numeric_limits<short>::round_style VCL_STATIC_CONST_INIT_INT_DEFN(vnl_round_toward_zero);
#endif

// ----------------------------------------------------------------------
// Constants and functions for double

union vnl_numeric_limits_double_nan {
  double nan;
  unsigned char x[8];

  vnl_numeric_limits_double_nan() {
#if VXL_BIG_ENDIAN
    x[0] = 0x7f;
    x[1] = x[2] = x[3] = x[4] = x[5] = x[6] = x[7] = 0xff;
#else
    x[7] = 0x7f;
    x[0] = x[1] = x[2] = x[3] = x[4] = x[5] = x[6] = 0xff;
#endif
  }
};
static vnl_numeric_limits_double_nan dnan;

union vnl_numeric_limits_double_inf {
  double inf;
  unsigned char x[8];

  vnl_numeric_limits_double_inf() {
#ifdef __alpha__ // Alpha throws a floating exception when evaluating IEEE Inf
    x[7] = 0x7f; x[6] = 0xef;
    x[0] = x[1] = x[2] = x[3] = x[4] = x[5] = 0xff;
#elif VXL_BIG_ENDIAN
    x[0] = 0x7f; x[1] = 0xf0;
    x[2] = x[3] = x[4] = x[5] = x[6] = x[7] = 0x00;
#else
    x[7] = 0x7f; x[6] = 0xf0;
    x[0] = x[1] = x[2] = x[3] = x[4] = x[5] = 0x00;
#endif
  }
};
static vnl_numeric_limits_double_inf dinf;

double vnl_numeric_limits<double>::infinity()
{
  return dinf.inf;
}

double vnl_numeric_limits<double>::quiet_NaN()
{
  return dnan.nan;
}

double vnl_numeric_limits<double>::signaling_NaN()
{
  return quiet_NaN();
}

#if !VCL_STATIC_CONST_INIT_INT_NO_DEFN
const bool vnl_numeric_limits<double>::is_specialized VCL_STATIC_CONST_INIT_INT_DEFN(true);
const int  vnl_numeric_limits<double>::digits   VCL_STATIC_CONST_INIT_INT_DEFN(53);
const int  vnl_numeric_limits<double>::digits10 VCL_STATIC_CONST_INIT_INT_DEFN( 15);
const bool vnl_numeric_limits<double>::is_signed  VCL_STATIC_CONST_INIT_INT_DEFN(true);
const bool vnl_numeric_limits<double>::is_integer VCL_STATIC_CONST_INIT_INT_DEFN(false);
const bool vnl_numeric_limits<double>::is_exact   VCL_STATIC_CONST_INIT_INT_DEFN(false);
const int  vnl_numeric_limits<double>::radix VCL_STATIC_CONST_INIT_INT_DEFN(2);
const int  vnl_numeric_limits<double>::min_exponent   VCL_STATIC_CONST_INIT_INT_DEFN(-1021);
const int  vnl_numeric_limits<double>::min_exponent10 VCL_STATIC_CONST_INIT_INT_DEFN(-307);
const int  vnl_numeric_limits<double>::max_exponent   VCL_STATIC_CONST_INIT_INT_DEFN(1024);
const int  vnl_numeric_limits<double>::max_exponent10 VCL_STATIC_CONST_INIT_INT_DEFN(308);
const bool vnl_numeric_limits<double>::has_infinity      VCL_STATIC_CONST_INIT_INT_DEFN(true);
const bool vnl_numeric_limits<double>::has_quiet_NaN     VCL_STATIC_CONST_INIT_INT_DEFN(true);
const bool vnl_numeric_limits<double>::has_signaling_NaN VCL_STATIC_CONST_INIT_INT_DEFN(true);
const bool vnl_numeric_limits<double>::has_denorm        VCL_STATIC_CONST_INIT_INT_DEFN(true);
const bool vnl_numeric_limits<double>::is_iec559  VCL_STATIC_CONST_INIT_INT_DEFN(true);
const bool vnl_numeric_limits<double>::is_bounded VCL_STATIC_CONST_INIT_INT_DEFN(true);
const bool vnl_numeric_limits<double>::is_modulo  VCL_STATIC_CONST_INIT_INT_DEFN(false);
const bool vnl_numeric_limits<double>::traps      VCL_STATIC_CONST_INIT_INT_DEFN(true);
const bool vnl_numeric_limits<double>::tinyness_before VCL_STATIC_CONST_INIT_INT_DEFN(true);
const vnl_float_round_style vnl_numeric_limits<double>::round_style VCL_STATIC_CONST_INIT_INT_DEFN(vnl_round_to_nearest);
#endif


// ----------------------------------------------------------------------
// Constants and functions for long double

static const unsigned int szl = sizeof(long double);

union vnl_numeric_limits_long_double_nan {
  long double nan;
  unsigned char x[szl];

  vnl_numeric_limits_long_double_nan() {
    for (unsigned int i=0; i<szl; ++i) x[i] = 0xff;
#if VXL_BIG_ENDIAN
    x[0] = 0x7f;
#else
#ifdef __INTEL_COMPILER
  // 12 bytes for alignment, 10 bytes for the actual number.
  if (szl == 12) {
    x[szl-3] = 0x7f;
  }
#else
    x[szl-1] = 0x7f;
#endif
#endif
  }
};
static vnl_numeric_limits_long_double_nan ldnan;

union vnl_numeric_limits_long_double_inf {
  long double inf;
  unsigned char x[szl];

  vnl_numeric_limits_long_double_inf() {
    for (unsigned int i=0; i<szl; ++i) x[i] = 0x00;
#ifdef __alpha__ // Alpha throws a floating exception when evaluating IEEE Inf
    x[szl-1] = 0x7f; x[szl-2] = 0xef;
    for (unsigned int i=0; i<szl-2; ++i) x[i] = 0xff;
#elif VXL_BIG_ENDIAN
    x[0] = 0x7f; x[1] = 0xf0;
#else
# if defined(VCL_BORLAND)
    x[szl-1] = 0x7f; x[szl-2] = 0xff; x[szl-3] = 0x80;
# else
    x[szl-1] = 0x7f; x[szl-2] = 0xf0;
    if (szl == 12) // intel
      x[9]=x[11]=0x7f, x[8]=x[10]=0xff, x[7] = 0x80;
# endif
#endif
  }
};
static vnl_numeric_limits_long_double_inf ldinf;

long double vnl_numeric_limits<long double>::infinity()
{
  return ldinf.inf;
}

long double vnl_numeric_limits<long double>::quiet_NaN()
{
 // http://cch.loria.fr/documentation/IEEE754/numerical_comp_guide/ncg_math.doc.html#526
  return ldnan.nan;
}

long double vnl_numeric_limits<long double>::signaling_NaN()
{
  return quiet_NaN();
}

#if !VCL_STATIC_CONST_INIT_INT_NO_DEFN
const bool vnl_numeric_limits<long double>::is_specialized VCL_STATIC_CONST_INIT_INT_DEFN(true);
const int  vnl_numeric_limits<long double>::digits   VCL_STATIC_CONST_INIT_INT_DEFN((int)(85-10*szl+.75*szl*szl));
const int  vnl_numeric_limits<long double>::digits10 VCL_STATIC_CONST_INIT_INT_DEFN((int)(9-3.5*szl+.25*szl*szl-5));
// this is 15, 21, and 35 for sizes 8, 12, and 16.
const bool vnl_numeric_limits<long double>::is_signed  VCL_STATIC_CONST_INIT_INT_DEFN(true);
const bool vnl_numeric_limits<long double>::is_integer VCL_STATIC_CONST_INIT_INT_DEFN(false);
const bool vnl_numeric_limits<long double>::is_exact   VCL_STATIC_CONST_INIT_INT_DEFN(false);
const int  vnl_numeric_limits<long double>::radix VCL_STATIC_CONST_INIT_INT_DEFN(2);
const int  vnl_numeric_limits<long double>::min_exponent   VCL_STATIC_CONST_INIT_INT_DEFN(-1021);
const int  vnl_numeric_limits<long double>::min_exponent10 VCL_STATIC_CONST_INIT_INT_DEFN(-307);
const int  vnl_numeric_limits<long double>::max_exponent   VCL_STATIC_CONST_INIT_INT_DEFN(1024);
const int  vnl_numeric_limits<long double>::max_exponent10 VCL_STATIC_CONST_INIT_INT_DEFN(308);
const bool vnl_numeric_limits<long double>::has_infinity      VCL_STATIC_CONST_INIT_INT_DEFN(true);
const bool vnl_numeric_limits<long double>::has_quiet_NaN     VCL_STATIC_CONST_INIT_INT_DEFN(true);
const bool vnl_numeric_limits<long double>::has_signaling_NaN VCL_STATIC_CONST_INIT_INT_DEFN(true);
const bool vnl_numeric_limits<long double>::has_denorm        VCL_STATIC_CONST_INIT_INT_DEFN(true);
const bool vnl_numeric_limits<long double>::is_iec559  VCL_STATIC_CONST_INIT_INT_DEFN(true);
const bool vnl_numeric_limits<long double>::is_bounded VCL_STATIC_CONST_INIT_INT_DEFN(true);
const bool vnl_numeric_limits<long double>::is_modulo  VCL_STATIC_CONST_INIT_INT_DEFN(false);
const bool vnl_numeric_limits<long double>::traps      VCL_STATIC_CONST_INIT_INT_DEFN(true);
const bool vnl_numeric_limits<long double>::tinyness_before VCL_STATIC_CONST_INIT_INT_DEFN(true);
const vnl_float_round_style vnl_numeric_limits<long double>::round_style VCL_STATIC_CONST_INIT_INT_DEFN(vnl_round_to_nearest);
#endif

// ----------------------------------------------------------------------
// Constants and functions for float

union vnl_numeric_limits_float_nan {
  float nan;
  unsigned char x[4];

  vnl_numeric_limits_float_nan() {
#if VXL_BIG_ENDIAN
    x[0] = 0x7f; x[1] = x[2] = x[3] = 0xff;
#else
    x[3] = 0x7f; x[0] = x[1] = x[2] = 0xff;
#endif
  }
};
static vnl_numeric_limits_float_nan fnan;

union vnl_numeric_limits_float_inf {
  float inf;
  unsigned char x[4];

  vnl_numeric_limits_float_inf() {
#ifdef __alpha__ // Alpha throws a floating exception when evaluating IEEE Inf
    x[3] = 0x7f; x[2] = 0x7f; x[1] = x[0] = 0xff;
#elif VXL_BIG_ENDIAN
    x[0] = 0x7f; x[1] = 0x80; x[2] = x[3] = 0x00;
#else
    x[3] = 0x7f; x[2] = 0x80; x[1] = x[0] = 0x00;
#endif
  }
};
static vnl_numeric_limits_float_inf finf;

float vnl_numeric_limits<float>::infinity()
{
  return finf.inf;
}

float vnl_numeric_limits<float>::quiet_NaN()
{
  return fnan.nan;
}

float vnl_numeric_limits<float>::signaling_NaN()
{
  return quiet_NaN();
}

#if !VCL_STATIC_CONST_INIT_INT_NO_DEFN
const bool vnl_numeric_limits<float>::is_specialized VCL_STATIC_CONST_INIT_INT_DEFN(true);
const int  vnl_numeric_limits<float>::digits   VCL_STATIC_CONST_INIT_INT_DEFN(24);
const int  vnl_numeric_limits<float>::digits10 VCL_STATIC_CONST_INIT_INT_DEFN( 6);
const bool vnl_numeric_limits<float>::is_signed  VCL_STATIC_CONST_INIT_INT_DEFN(true);
const bool vnl_numeric_limits<float>::is_integer VCL_STATIC_CONST_INIT_INT_DEFN(false);
const bool vnl_numeric_limits<float>::is_exact   VCL_STATIC_CONST_INIT_INT_DEFN(false);
const int  vnl_numeric_limits<float>::radix VCL_STATIC_CONST_INIT_INT_DEFN(2);
const int  vnl_numeric_limits<float>::min_exponent   VCL_STATIC_CONST_INIT_INT_DEFN(-125);
const int  vnl_numeric_limits<float>::min_exponent10 VCL_STATIC_CONST_INIT_INT_DEFN(-37);
const int  vnl_numeric_limits<float>::max_exponent   VCL_STATIC_CONST_INIT_INT_DEFN(128);
const int  vnl_numeric_limits<float>::max_exponent10 VCL_STATIC_CONST_INIT_INT_DEFN(38);
const bool vnl_numeric_limits<float>::has_infinity      VCL_STATIC_CONST_INIT_INT_DEFN(true);
const bool vnl_numeric_limits<float>::has_quiet_NaN     VCL_STATIC_CONST_INIT_INT_DEFN(true);
const bool vnl_numeric_limits<float>::has_signaling_NaN VCL_STATIC_CONST_INIT_INT_DEFN(true);
const bool vnl_numeric_limits<float>::has_denorm        VCL_STATIC_CONST_INIT_INT_DEFN(false);
const bool vnl_numeric_limits<float>::is_iec559  VCL_STATIC_CONST_INIT_INT_DEFN(true);
const bool vnl_numeric_limits<float>::is_bounded VCL_STATIC_CONST_INIT_INT_DEFN(true);
const bool vnl_numeric_limits<float>::is_modulo  VCL_STATIC_CONST_INIT_INT_DEFN(false);
const bool vnl_numeric_limits<float>::traps      VCL_STATIC_CONST_INIT_INT_DEFN(true);
const bool vnl_numeric_limits<float>::tinyness_before VCL_STATIC_CONST_INIT_INT_DEFN(true);
const vnl_float_round_style vnl_numeric_limits<float>::round_style VCL_STATIC_CONST_INIT_INT_DEFN(vnl_round_to_nearest);
#endif
-------------- next part --------------
//-*- c++ -*-------------------------------------------------------------------
#ifndef vcl_compiler_h_
#define vcl_compiler_h_

//:
// \file
//
// It's much better to determine the compiler automatically here than to depend
// on command-line flags being set.

#if defined(__sgi) && !defined(__GNUC__)
# ifndef _COMPILER_VERSION
#  define VCL_SGI_CC_6
# else
#  if (_COMPILER_VERSION >= 700)
#   define VCL_SGI_CC_7
#  else
#   define VCL_SGI_CC_6
#  endif
#  if   (_COMPILER_VERSION >= 730)
#   define VCL_SGI_CC_730
#  elif (_COMPILER_VERSION >= 720)
#   define VCL_SGI_CC_720
#  endif
#  define VCL_SGI_CC
# endif
#endif

#if defined(__SUNPRO_CC)
# define VCL_SUNPRO_CC
# if (__SUNPRO_CC>=0x500)
#  define VCL_SUNPRO_CC_50
# else
#  undef VCL_SUNPRO_CC_50
# endif
# ifdef INSTANTIATE_TEMPLATES
#  define _RWSTD_COMPILE_INSTANTIATE
# endif
#endif

#if defined(__GNUC__)
# define VCL_GCC
# if (__GNUC__<=1)
#  error "forget it."
# elif (__GNUC__==2)
#  if (__GNUC_MINOR__>=100)
#   error "I need some help here."
#  elif (__GNUC_MINOR__>=95)
#   define VCL_GCC_295
#  elif (__GNUC_MINOR__>8)
#   define VCL_EGCS
#  elif (__GNUC_MINOR__>7)
#   define VCL_GCC_28
#  elif (__GNUC_MINOR__>6)
#   define VCL_GCC_27
#  endif
#  if (__GNUC_MINOR__>7)
#   define VCL_GCC_EGCS // so this is the union of EGCS, GCC_28 and GCC_295
#  endif
# elif (__GNUC__==3)
#  define VCL_GCC_30
#  if (__GNUC_MINOR__>0)
#   define VCL_GCC_31
#  endif
#  if (__GNUC_MINOR__>1)
#   define VCL_GCC_32
#  endif
# else
#  error "Dunno about this gcc"
# endif
#endif

#if defined(_WIN32) || defined(WIN32)
# define VCL_WIN32
# if defined(_MSC_VER)
#  define VCL_VC
#  if _MSC_VER >= 1300
#   define VCL_VC_DOTNET 1 // VC is at least version >= 7.0
#   if _MSC_VER >= 1310
#    define VCL_VC71 1     // .NET 2003 = Version 7.1
#   else
#    define VCL_VC70 1     // earlier .NET versions = Version 7.0
#   endif
#  elif _MSC_VER >= 1200   // last version before advent of .NET = Version 6.0
#   define VCL_VC60 1
#  else
#   define VCL_VC50 1
#  endif
# elif defined(__BORLANDC__)
#  define VCL_BORLAND
#  if __BORLANDC__ >= 0x0560
#   define VCL_BORLAND_56
#  elif __BORLANDC__ >= 0x0550
#   define VCL_BORLAND_55
#  endif
# endif
#endif

// win32 or vc++ ?
// awf hack alert:
#ifdef VCL_VC
#  ifdef VCL_VC60
#    pragma warning(disable:4786 4660 4661)
#    pragma warning(disable:4786 4660 4355 4390)
#  elif VCL_VC_DOTNET
// 4786: 'identifier' : identifier was truncated to 'number' characters in the debug information
// 4018: signed/unsigned mismatch
// 4146: unary minus operator applied to unsigned type, result still unsigned
// 4267: conversion related to size_t
// 4355: 'this' : used in base member initializer list
#    pragma warning(disable:4786 4018 4146 4267 4355)
#  endif
#endif

#if defined(__KCC) // KAI compiler
# define VCL_KAI
#endif

#if defined(__ICC) ||defined(__ECC) || defined(__INTEL_COMPILER) // Intel compiler?
# define VCL_ICC
#endif

#if defined(como4301) // Comeau C/C++ 4.3.0.1
# define VCL_COMO
#endif


#if defined(__MWERKS__)
// [sic]
# define VCL_METRO_WERKS
#endif

// include header files generated by configure.
#include <vcl_config_manual.h>
#include <vcl_config_compiler.h>
#include <vcl_config_headers.h>

// This *needs* to come after vcl_config_headers.h
#if defined(__GNUC__) && !defined(__INTEL_COMPILER)
# ifdef VCL_GCC_30
#  define GNU_LIBSTDCXX_V3 1
# elif !defined(GNU_LIBSTDCXX_V3) && defined(VCL_GCC_295) && VCL_CXX_HAS_HEADER_ISTREAM
// One difference between v2 and v3 is that the former has
// no <istream> header file whereas v3 has the lot.
#  define GNU_LIBSTDCXX_V3 1
# endif
#endif

// -------------------- default template parameters
#if VCL_CAN_DO_COMPLETE_DEFAULT_TYPE_PARAMETER
# define VCL_DFL_TYPE_PARAM_STLDECL(A, a) class A = a
#else
# define VCL_DFL_TYPE_PARAM_STLDECL(A, a) class A
#endif

#if VCL_CAN_DO_TEMPLATE_DEFAULT_TYPE_PARAMETER
# define VCL_DFL_TMPL_PARAM_STLDECL(A, a) class A = a
#else
# define VCL_DFL_TMPL_PARAM_STLDECL(A, a) class A
#endif

#define VCL_DFL_TMPL_ARG(classname) , classname

#if VCL_USE_NATIVE_STL
# define VCL_SUNPRO_ALLOCATOR_HACK(T) T VCL_SUNPRO_CLASS_SCOPE_HACK(std::allocator<T >)
#else
# define VCL_SUNPRO_ALLOCATOR_HACK(T) T // FIXME
#endif

//-------------------- template instantiation

// if the compiler doesn't understand "export", we just leave it out.
// gcc and SunPro 5.0 understand it, but they ignore it noisily.
#if !VCL_HAS_EXPORT || defined(VCL_EGCS) || defined(VCL_GCC_295) || defined(VCL_GCC_30) || defined(VCL_SUNPRO_CC_50)
# define export /* ignore */
#endif

#if VCL_NEEDS_INLINE_INSTANTIATION
# define VCL_INSTANTIATE_INLINE(symbol) template symbol
#else
# define VCL_INSTANTIATE_INLINE(symbol) /* */
#endif

// work-around to get template instantiation to work correctly with SunPro
// check flag to turn on inlining
#undef IUEi_STL_INLINE
#if defined(INLINE_EXPLICIT_FLAG) && defined(VCL_SUNPRO_CC) && defined(INSTANTIATE_TEMPLATES)
# define IUEi_STL_INLINE
#else
# define IUEi_STL_INLINE inline
#endif

//--------------------------------------------------------------------------------

#if VCL_FOR_SCOPE_HACK
# undef for
# define for if (false) { } else for
typedef int saw_VCL_FOR_SCOPE_HACK;
#endif

// fix to instantiate template functions
#define VCL_INSTANTIATE_NONINLINE(fn_decl) template fn_decl

// -------------------- handy macros

//: VCL_COMMA
//
// Handy for passing things with commas in them to CPP macros.  e.g.
// DO_MACRO(pair<A,B>) can be replaced by DO_MACRO(pair<A VCL_COMMA B>).
#define VCL_COMMA ,


//: VCL_VOID_RETURN
//
// VCL_VOID_RETURN is used as a return type where void is expected,
// as in return VCL_VOID_RETURN;
#define VCL_VOID_RETURN /*empty*/

#endif // vcl_compiler_h_
-------------- next part --------------
// This is core/vnl/vnl_math.cxx
#ifdef VCL_NEEDS_PRAGMA_INTERFACE
#pragma implementation
#endif
//:
// \file

#include "vnl_math.h"
#include <vxl_config.h>

#if defined(VCL_VC)
// I don't think we need this, because <ieeefp.h> is available -- fsm
# include <Float.h> // for 'isnan' and 'finite'
// # define isnan _isnan
# define finite _finite
# define finitef _finite
#ifndef finitel
# define finitel _finite
#endif
# define isnan _isnan
#elif VXL_IEEEFP_HAS_FINITE
# include <ieeefp.h>
# ifndef finitef
#  define finitef finite
# endif
# ifndef finitel
#  define finitel finite
# endif

#elif VXL_C_MATH_HAS_FINITE
# include <math.h>  // dont_vxl_filter: this is *not* supposed to be <cmath>
# if !defined(__alpha__) && !(VXL_C_MATH_HAS_FINITEF)// on Alpha, finitef() must be used for float args instead of finite()
#  define finitef finite
# endif
# if !(VXL_C_MATH_HAS_FINITEL)
#  define finitel finite
# endif

#elif defined(SYSV) && !defined(hppa)
// needed on platforms with finite() declared in strange places, e.g. on alpha
extern "C" int finite(double);
# ifdef __alpha__ // on Alpha, finitef() must be used for float args instead of finite()
extern "C" int finitef(float);
# else
#  define finitef finite
# endif

#elif defined(VCL_BORLAND) 
# include <math.h>
# include <float.h>
#else
# warning finite() is not declared on this platform
# define VNL_HAS_NO_FINITE
#endif

#ifdef VCL_SUNPRO_CC_50
# include <math.h> // dont_vxl_filter: no HUGE_VAL or isnan() in <cmath>
#endif

#if defined(__APPLE__)
# include <math.h>
# define isnan __isnan
#endif

//--------------------------------------------------------------------------------

#if !VCL_STATIC_CONST_INIT_FLOAT_NO_DEFN

// constants
const double vnl_math::e               VCL_STATIC_CONST_INIT_FLOAT_DEFN( 2.7182818284590452354  );
const double vnl_math::log2e           VCL_STATIC_CONST_INIT_FLOAT_DEFN( 1.4426950408889634074  );
const double vnl_math::log10e          VCL_STATIC_CONST_INIT_FLOAT_DEFN( 0.43429448190325182765 );
const double vnl_math::ln2             VCL_STATIC_CONST_INIT_FLOAT_DEFN( 0.69314718055994530942 );
const double vnl_math::ln10            VCL_STATIC_CONST_INIT_FLOAT_DEFN( 2.30258509299404568402 );
const double vnl_math::pi              VCL_STATIC_CONST_INIT_FLOAT_DEFN( 3.14159265358979323846 );
const double vnl_math::pi_over_2       VCL_STATIC_CONST_INIT_FLOAT_DEFN( 1.57079632679489661923 );
const double vnl_math::pi_over_4       VCL_STATIC_CONST_INIT_FLOAT_DEFN( 0.78539816339744830962 );
const double vnl_math::one_over_pi     VCL_STATIC_CONST_INIT_FLOAT_DEFN( 0.31830988618379067154 );
const double vnl_math::two_over_pi     VCL_STATIC_CONST_INIT_FLOAT_DEFN( 0.63661977236758134308 );
const double vnl_math::two_over_sqrtpi VCL_STATIC_CONST_INIT_FLOAT_DEFN( 1.12837916709551257390 );
const double vnl_math::sqrt2           VCL_STATIC_CONST_INIT_FLOAT_DEFN( 1.41421356237309504880 );
const double vnl_math::sqrt1_2         VCL_STATIC_CONST_INIT_FLOAT_DEFN( 0.70710678118654752440 );

// IEEE double machine precision
const double vnl_math::eps             VCL_STATIC_CONST_INIT_FLOAT_DEFN( 2.2204460492503131e-16 );
const double vnl_math::sqrteps         VCL_STATIC_CONST_INIT_FLOAT_DEFN( 1.490116119384766e-08  );

  //: IEEE single machine precision
const float vnl_math::float_eps        VCL_STATIC_CONST_INIT_FLOAT_DEFN( 1.192092896e-07f );
const float vnl_math::float_sqrteps    VCL_STATIC_CONST_INIT_FLOAT_DEFN( 3.4526698307e-4f );

#endif

//--------------------------------------------------------------------------------

#if defined(VCL_ICC)
#include <mathimf.h> // defines isnanf, isnan, and isnanl
//: Return true iff x is "Not a Number"
bool vnl_math_isnan(float x) { return isnanf(x); }
//: Return true iff x is "Not a Number"
bool vnl_math_isnan(double x) { return isnan(x); }
#if (__INTEL_COMPILER > 800)
// According to Intel, long double support is correct in icc 8.1.
//: Return true iff x is "Not a Number"
bool vnl_math_isnan(long double x) { return isnanl(x); }
#else
//: Return true iff x is "Not a Number"
bool vnl_math_isnan(long double x) { return isnan((double)x); }
#endif
#elif defined(VCL_BORLAND)
//: Return true iff x is "Not a Number"
bool vnl_math_isnan(float x) { return _isnan(x); }
//: Return true iff x is "Not a Number"
bool vnl_math_isnan(double x) { return _isnan(x); }
//: Return true iff x is "Not a Number"
bool vnl_math_isnan(long double x) { return _isnanl(x); }
#elif !defined(VNL_HAS_NO_FINITE) && !defined(VCL_SGI_CC_7) && !defined(__alpha__) && !defined(VCL_WIN32)
//: Return true iff x is "Not a Number"
bool vnl_math_isnan(float x) { return x != x; } // causes "floating exception" on alpha & sgi
//: Return true iff x is "Not a Number"
bool vnl_math_isnan(double x) { return x != x; }
//: Return true iff x is "Not a Number"
bool vnl_math_isnan(long double x) { return x != x; }
#else
// Auxiliary function to simplify notation
# ifndef DEBUG
static inline unsigned int bMp(void*x,unsigned int y,int p=0){return ((((unsigned int*)x)[p])&y);}
static inline bool bMe(void*x,unsigned int y,int p=0){return ((((unsigned int*)x)[p])&y)==y;}
# else
# include <vcl_iostream.h>
static inline unsigned int bMp(void* x, unsigned int y, int p=0) {
  unsigned char* v=(unsigned char*)x;
  vcl_cout<<int(v[4*p])<<' '<<int(v[4*p+1])<<' '<<int(v[4*p+2])<<' '<<int(v[4*p+3])<<" & ";
  v=(unsigned char*)(&y);
  vcl_cout<<int(v[0])<<' '<<int(v[1])<<' '<<int(v[2])<<' '<<int(v[3])<<" = ";
  unsigned int z = ((((unsigned int*)x)[p]) & y);
  v=(unsigned char*)(&z);
  vcl_cout<<int(v[0])<<' '<<int(v[1])<<' '<<int(v[2])<<' '<<int(v[3]);
  if (z == y) vcl_cout<<" ==";
  vcl_cout << '\n';
  return z;
}
static inline bool bMe(void* x, unsigned int y, int p=0) { return bMp(x,y,p) == y; }
# endif
# if VXL_BIG_ENDIAN
static const int sz_f = 0;
static const int sz_d = 0;
static const int sz_l = 0;
# else
static const int sz_f = sizeof(float)/sizeof(int) -1;
static const int sz_d = sizeof(double)/sizeof(int) -1;
static const int sz_l = sizeof(long double)/sizeof(int) -1;
# endif
// Assume IEEE floating point number representation
bool vnl_math_isnan( float x){return bMe(&x,0x7f800000L,sz_f)&&bMp(&x,0x007fffffL,sz_f);}
bool vnl_math_isnan(double x){return bMe(&x,0x7ff00000L,sz_d)&&bMp(&x,0x000fffffL,sz_d);}
bool vnl_math_isnan(long double x) {
  if (sizeof(long double) == 8) return bMe(&x,0x7ff00000L,sz_l) && bMp(&x,0x000fffffL,sz_l);
  else if (sizeof(long double) <= 12) return bMe(&x,0x4001ffffL,sz_l) && bMp(&x,0x40000000,sz_l-4);
  else return bMe(&x,0x7ff70000L,sz_l) && bMp(&x,0x0008ffffL,sz_l);
}
#endif

// fsm
// On linux noshared builds, with optimisation on, calling 'finite' within the
// scope of vnl_math causes vnl_math_isinf to be called. This blows the stack.
// Plausible theory : 'finite' is a preprocessor macro, defined in terms of a
// macro called 'isinf'.
#if defined(isinf)
# if defined(__GNUC__) || defined(VCL_METRO_WERKS) 
// I do not know if MW accepts #warning. Comment out the #undef if not.
// Ignore this for the INTEL_COMPILER icc
#if !defined(VCL_ICC)
#  warning macro isinf is defined
#  undef isinf
#endif
# else
// do not fail silently
#  error macro isinf is defined
# endif
#endif

#if defined(VCL_BORLAND)
//: Return true if x is neither NaN nor Inf.
bool vnl_math_isfinite(float x) { return _finite(x) != 0; }
//: Return true if x is neither NaN nor Inf.
bool vnl_math_isfinite(double x) { return _finite(x) != 0; }
//: Return true if x is neither NaN nor Inf.
bool vnl_math_isfinite(long double x) { return _finitel(x) != 0 && !_isnanl(x); }
#elif !defined(VNL_HAS_NO_FINITE) && !defined(VCL_ICC)
//: Return true if x is neither NaN nor Inf.
bool vnl_math_isfinite(float x) { return finitef(x) != 0; }
//: Return true if x is neither NaN nor Inf.
bool vnl_math_isfinite(double x) { return finite(x) != 0; }
//: Return true if x is neither NaN nor Inf.
bool vnl_math_isfinite(long double x) { return finitel(x) != 0; }
#elif defined(VCL_ICC)
#include <mathimf.h> // defines isfinite, isfinitef, isfinitel
//: Return true if x is neither NaN nor Inf.
bool vnl_math_isfinite(float x) { return isfinitef(x) != 0; }
//: Return true if x is neither NaN nor Inf.
bool vnl_math_isfinite(double x) { return isfinite(x) != 0; }
#if (__INTEL_COMPILER > 800)
// According to Intel, long double support is correct in icc 8.1.
//: Return true if x is neither NaN nor Inf.
bool vnl_math_isfinite(long double x) { return isfinitel(x) != 0; }
#else
//: Return true if x is neither NaN nor Inf.
bool vnl_math_isfinite(long double x) { return isfinite((double)x) != 0; }
#endif
#else
// Assume IEEE floating point number representation
bool vnl_math_isfinite(float x) { return !bMe(&x,0x7f800000L,sz_f) && bMp(&x,0x7fffffffL,sz_f) != 0x7f7fffffL; }
bool vnl_math_isfinite(double x) { return !bMe(&x,0x7ff00000L,sz_d); }
bool vnl_math_isfinite(long double x) {
  if (sizeof(long double) == 8) return !bMe(&x,0x7ff00000L,sz_l);
  else if (sizeof(long double) <= 12) return !bMe(&x,0xbfff7fffL,sz_l) && !bMe(&x,0x4001ffffL,sz_l);
  else return !bMe(&x,0x7ff70000L,sz_l);
}
#endif


#if defined(VCL_BORLAND)
//: Return true if x is inf
bool vnl_math_isinf(float x) { return !_finite(x) && !_isnan(x); }
//: Return true if x is inf
bool vnl_math_isinf(double x) { return !_finite(x) && !_isnan(x); }
//: Return true if x is inf
bool vnl_math_isinf(long double x) { return !_finitel(x) && !_isnanl(x); }
#elif !defined(VNL_HAS_NO_FINITE) && !defined(VCL_ICC)
//: Return true if x is inf
bool vnl_math_isinf(float x) { return !finitef(x) && !isnan(x); }
//: Return true if x is inf
bool vnl_math_isinf(double x) { return !finite(x) && !isnan(x); }
//: Return true if x is inf
bool vnl_math_isinf(long double x) { return !finitel(x) && !isnan(x); }
#elif defined(VCL_ICC)
#include <mathimf.h> // defines isfinite, isfinitef, isfinitel, isnan
// The Intel macros appear to find the right type (float, double, long double)
//: Return true if x is inf
bool vnl_math_isinf(float x) { return isinff(x) != 0; }
//: Return true if x is inf
bool vnl_math_isinf(double x) { return isinf(x) != 0; }
#if (__INTEL_COMPILER > 800)
// According to Intel, long double support is correct in icc 8.1.
//: Return true if x is inf
bool vnl_math_isinf(long double x) { return isinfl(x) != 0; }
#else
//: Return true if x is inf
bool vnl_math_isinf(long double x) { 
  if (vnl_math_isnan(x)) return false;
  if (x >= 0.0) {
    return ((x+1) == x);
  } else {
    return ((-x+1) == -x);
  }
}
#endif
#else
// Assume IEEE floating point number representation
bool vnl_math_isinf(float x) {return(bMe(&x,0x7f800000L,sz_f)&&!bMp(&x,0x007fffffL,sz_f))||bMp(&x,0x7fffffffL,sz_f)==0x7f7fffffL;}
bool vnl_math_isinf(double x) { return bMe(&x,0x7ff00000L,sz_d) && !bMp(&x,0x000fffffL,sz_d); }
bool vnl_math_isinf(long double x) {
  if (sizeof(long double) == 8) return bMe(&x,0x7ff00000L,sz_l) && !bMp(&x,0x000fffffL,sz_l);
  else if (sizeof(long double) <= 12) return (bMe(&x,0xbfff7fffL,sz_l)||bMe(&x,0x4001ffffL,sz_l))&&!bMp(&x,0x40000000,sz_l-4);
  else return bMe(&x,0x7ff70000L,sz_l) && !bMp(&x,0x0008ffffL,sz_l);
}
#endif

//----------------------------------------------------------------------

//: Type-accessible infinities for use in templates.
template <class T> T vnl_huge_val(T);
double vnl_huge_val(double) { return HUGE_VAL; }
float  vnl_huge_val(float)  { return (float)HUGE_VAL; }
#ifdef _INT_64BIT_
long int vnl_huge_val(long int) { return 0x7fffffffffffffffL; }
int    vnl_huge_val(int)    { return 0x7fffffffffffffffL; }
#else
int    vnl_huge_val(int)    { return 0x7fffffff; }
#endif
short  vnl_huge_val(short)  { return 0x7fff; }
char   vnl_huge_val(char)   { return 0x7f; }

//----------------------------------------------------------------------


More information about the Insight-developers mailing list