/* This is -*- C -*- */
/* vim: set sw=2: */
/* $Id: gnan.c,v 1.9 2002/01/14 05:01:23 trow Exp $ */

/*
 * gnan.c
 *
 * Copyright (C) 2000 EMC Capital Management, Inc.
 *
 * Developed by Jon Trowbridge <trow@gnu.org>.
 *
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License as
 * published by the Free Software Foundation; either version 2 of the
 * License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
 * USA
 */

#include <config.h>
#include "gnan.h"

#include <string.h>
#include <signal.h>
#include <glib.h>
#include "guppi-debug.h"


union ieee754_big_endian_double {
  
  double d;

  struct {
    unsigned int signbit:1;
    unsigned int exp:11;
    unsigned int man0:20;
    unsigned int man1:32;
    
  } ieee;
};

union ieee754_little_endian_double {
  
  double d;

  struct {
    unsigned int man1:32;
    unsigned int man0:20;
    unsigned int exp:11;
    unsigned int signbit:1;
  } ieee;
};

/* little endian, but with big endian float word order */
union ieee754_mixed_endian_double {
  
  double d;

  struct {
    unsigned int man0:20;
    unsigned int exp:11;
    unsigned int signbit:1;
    unsigned int man1:32;
  } ieee;
};

#define ALPHA_TEST_VALUE 314159.314159
/* This is the bit-representation of ALPHA_TEST_VALUE under Alpha's
   non-IEEE FP. */
/* These values are WRONG */
static const guint alpha_test_bits[8] = {
  0x9b, 0xe5, 0xb2, 0x41, 0xbd, 0x2c, 0x13, 0x41
};
static const guint alpha_nan_bits_A[8] = {
  0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0xf8, 0x7f
};
static const guint alpha_nan_bits_B[8] = {
  0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0xf8, 0xff
};

const double G_NAN;

#ifndef NAN
static gboolean fake_isnan = FALSE;
#endif

#define FAKE_NAN 123.45678902468013579e-300

#define NUMBER 1.0
#define SIGNBIT 0
#define MAN0 0
#define MAN1 0
#define EXP 1023

gboolean
g_isnan (double x)
{
#ifndef NAN

  return fake_isnan ? (x == FAKE_NAN) : isnan (x);

#else

  return isnan (x);

#endif
}

gboolean
g_isinf (double x)
{
  /* FIXME: This needs to be made portable. */
  return isinf (x);
}

gboolean
g_finite (double x)
{
  /* FIXME: This needs to be made portable. */
  return finite (x);
}

void
gnan_init (void)
{
#ifdef NAN

  double d = NAN;

  /* If we have NAN, we just use it. */
  g_assert (isnan (NAN));
  memcpy ((double *)&G_NAN, &d, sizeof (double));

#else
  union ieee754_big_endian_double be;
  union ieee754_little_endian_double le;
  union ieee754_mixed_endian_double me;
  double d;

  /* Note: we declared NAN const so we have to be a bit careful about
     changing it.  Otherwise we outsmart the compiler and thus
     ourselves.  Not good.  */

  be.d = NUMBER;
  if (sizeof (double) == 8
      && be.ieee.signbit == SIGNBIT
      && be.ieee.man0 == MAN0
      && be.ieee.man1 == MAN1
      && be.ieee.exp == EXP) {
    
    be.ieee.signbit = 0;
    be.ieee.man0 = 1;
    be.ieee.man1 = 1;
    be.ieee.exp = 2047;

    memcpy ((double *)&G_NAN, &be.d, sizeof (double));
    if (isnan (G_NAN)) return;
  }

  le.d = NUMBER;
  if (sizeof (double) == 8
      && le.ieee.signbit == SIGNBIT
      && le.ieee.man0 == MAN0
      && le.ieee.man1 == MAN1
      && le.ieee.exp == EXP) {
    
    le.ieee.signbit = 0;
    le.ieee.man0 = 1;
    le.ieee.man1 = 1;
    le.ieee.exp = 2047;

    memcpy ((double *)&G_NAN, &le.d, sizeof (double));
    if (isnan (G_NAN)) return;
  }

  me.d = NUMBER;
  if (sizeof (double) == 8
      && me.ieee.signbit == SIGNBIT
      && me.ieee.man0 == MAN0
      && me.ieee.man1 == MAN1
      && me.ieee.exp == EXP) {
    
    me.ieee.signbit = 0;
    me.ieee.man0 = 1;
    me.ieee.man1 = 1;
    me.ieee.exp = 2047;

    memcpy ((double *)&G_NAN, &me.d, sizeof (double));
    if (isnan (G_NAN)) return;
  }

  /* Morten Welinder tells me that "Some people here argue that the
     specs guarantee a NaN to come out of this."  It is certainly
     worth a try... */
  {
    double a = -HUGE_VAL;
    double b = 0;
    *(double *)&G_NAN = a * b;
    if (isnan (G_NAN)) return;
  }
  
  /* This is abusive, but I need a quick and dirty way of
     handling Alpha FP right away... */
  d = ALPHA_TEST_VALUE;
  if (sizeof (double) == 8) {
    gint i;
    gboolean ok = TRUE;

    for (i=0; i<sizeof(double) && ok; ++i) {
      if (alpha_test_bits[i] != (guint)(((guchar *)&d)[i]))
	ok = FALSE;
    }

    if (ok) {
      memcpy ((gpointer)&G_NAN, alpha_nan_bits_A, 8);
      if (isnan (G_NAN)) return;
      memcpy ((gpointer)&G_NAN, alpha_nan_bits_B, 8);
      if (isnan (G_NAN)) return;
    }
  }

  /* This is crazy, but we are desparate at this point... */
  {
    struct sigaction fpe, old_fpe;

    fpe.sa_handler = SIG_IGN;

    sigaction (SIGFPE, &fpe, &old_fpe);
    d = sqrt (-1.0);
    sigaction (SIGFPE, &old_fpe, NULL);
    
    if (isnan (d)) {
      memcpy ((gpointer)&G_NAN, &d, sizeof (double));
      return;
    }
  }


  /* All else has failed, so we have to use the awful FAKE_NAN hack. */

  d = FAKE_NAN;
  fake_isnan = TRUE;
  memcpy ((gpointer)&G_NAN, &d, sizeof (double));

  g_assert (g_isnan (G_NAN));
  guppi_msg ("Using awful FAKE_NAN hack.");
  
#endif /* NAN */
}



/* $Id: gnan.c,v 1.9 2002/01/14 05:01:23 trow Exp $ */