/* See flt.h */

/****************************************************************************/
/* (C) Copyright 1993 Universidade Estadual de Campinas (UNICAMP)           */
/*                    Campinas, SP, Brazil                                  */
/*                                                                          */
/* This file can be freely distributed, modified, and used for any          */
/*   non-commercial purpose, provided that this copyright and authorship    */
/*   notice be included in any copy or derived version of this file.        */
/*                                                                          */
/* DISCLAIMER: This software is offered ``as is'', without any guarantee    */
/*   as to fitness for any particular purpose.  Neither the copyright       */
/*   holder nor the authors or their employers can be held responsible for  */
/*   any damages that may result from its use.                              */
/****************************************************************************/

#include "flt.h"
#include <js.h>
#include <ioprotos.h>
#include <values.h>
#include <floatingpoint.h>
#include <math.h>
#include <stdio.h>

char *ROUND_BUFP = NULL;

Float MaxFloat;
Float Infinity;

void flt_init(void)
  {
    Infinity = HUGE_VAL;
    MaxFloat = MAXFLOAT;
    ROUND_NEAR; 
    /* fprintf(stderr, "flt_init: fp_precision = %8x  fsr = %8x\n", fp_precision, flt_get_fsr()); */
  }

void flt_print (FILE *f, Float x)
  {
    fprintf(f, FLT_FMT_SPEC, x);
  }
  
/*
  There should be a simpler and faster way to do these: */

#define UPDATE_ERR                                                 \
    if (zlo != zhi)                                                \
      {                                                            \
	if ((zhi >= PlusInfinity) || (zlo <= MinusInfinity))       \
	  { *zp = PlusInfinity; *errp = PlusInfinity; return; }    \
	ROUND_UP;                                                  \
	{ Float d1 = zhi - (*zp);                                  \
	  Float d2 = (*zp) - zlo;                                  \
	  Float del = FMAX(d1, d2); /* Shouldn't overflow. */      \
	  *errp = *errp + del;  /* May overflow */                 \
	}                                                          \
      }
  
void flt_add(Float x, Float y, FloatP zp, FloatP errp)
  {
    Float zhi, zlo;
    *zp = x + y;         /* May overflow. */
    ROUND_DOWN;
    zlo = x + y;
    ROUND_UP;
    zhi = x + y;
    UPDATE_ERR
  }

void flt_sub(Float x, Float y, FloatP zp, FloatP errp)
  {
    Float zhi, zlo;
    *zp = x - y;         /* May overflow. */
    ROUND_DOWN;
    zlo = x - y;
    ROUND_UP;
    zhi = x - y;
    UPDATE_ERR
  }

void flt_mul(Float x, Float y, FloatP zp, FloatP errp)
  {
    Float zhi, zlo;
    *zp = x * y;
    ROUND_DOWN;
    zlo = x * y;
    ROUND_UP;
    zhi = x * y;
    UPDATE_ERR
  }

void flt_div(Float x, Float y, FloatP zp, FloatP errp)
  { Float zhi, zlo;
    *zp = x / y;         /* May overflow. */
    ROUND_DOWN;
    zlo = x / y;
    ROUND_UP;
    zhi = x / y;
    UPDATE_ERR
  }

void flt_inv(Float x, FloatP zp, FloatP errp)
  {
    Float zhi, zlo;
    *zp = One/x;
    ROUND_DOWN;
    zlo = One/x;
    ROUND_UP;
    zhi = One/x;
    UPDATE_ERR
  }

void flt_sqrt(Float x, FloatP zp, FloatP errp)
  {
    Float zhi, zlo;
    assert(x >= Zero, "flt_sqrt: argument is negative");
    *zp = sqrt(x);
    ROUND_DOWN;
    zlo = sqrt(x);
    ROUND_UP;
    zhi = sqrt(x);
    UPDATE_ERR;
  }

void flt_scale (
    Float x, Float alpha, 
    Float zeta,
    FloatP zp,
    FloatP errp
  )
  { double xa = ((double) x) * ((double) alpha); 
    /* With IEEE arithmetic, $xa$ should be exact        */
    /*   and without overflow, whatever the current rounding mode. */

    Float zhi, zlo;
    if (zeta == One)
      { *zp = (Float) xa;
	ROUND_DOWN;
	zlo = (Float) xa;
	ROUND_UP;
	zhi = (Float) xa;
      }
    else
      { *zp = ((Float) xa) / zeta;
        ROUND_DOWN;
        zlo = ((Float) xa) / zeta;
        ROUND_UP;
        zhi = ((Float) xa) / zeta;
      }
    UPDATE_ERR
  }

void flt_mix (
    Float x, Float alpha, 
    Float y, Float beta, 
    Float zeta,
    FloatP zp,
    FloatP errp
  )
  { double xa = ((double) x) * ((double) alpha); 
    double yb = ((double) y) * ((double) beta);  
    /* With IEEE arithmetic,  $xa$ and $yb$ should be exact        */
    /*   and without overflow, whatever the current rounding mode. */
    Float zhi, zlo;
    if (zeta == One)
      {  *zp = (Float)(xa + yb);
	ROUND_DOWN;
	zlo = (Float)(xa + yb);
	ROUND_UP;
	zhi = (Float)(xa + yb);
      }
    else
      { *zp = ((Float)(xa + yb)) / zeta;
	ROUND_DOWN;
	zlo = ((Float)(xa + yb)) / zeta;
	ROUND_UP;
	zhi = ((Float)(xa + yb)) / zeta;
      }
    UPDATE_ERR
  }
            
Float flt_random(void)
  {
    /* This code is not entirely correct for Float=float,
       because normalization may produce trailing zero bits 
       with probability slightly larger than 1/2.
       It is definitely wrong for Float=double. 
    */
    Float hi = ((random()&65535) + 0.0)/65536.0;
    Float lo = ((random()&65535) + 0.0)/65536.0;
    return(hi + lo/65536.0);
  }

Float flt_random_mag(int avg, int dev)
  {
    #define TwoToSixteen (65536.0)
    
    Float t = One;
    int exp;
    
    if (dev == 0)
      { exp = avg; }
    else
      { /* Compute exponent "exp": */
	int w = dev + dev + 1;
        unsigned mask = 1;
	/* find smallest all-ones mask that is no less than 2*dev + 1: */
	while ((mask & w) != w) mask = ((mask << 1) | 1);
	/* Generate random integer in range [0..2*dev]: */
	do { exp = (random() & mask); } while (exp >= w);
	/* Compute exponent: */
	exp = avg + (exp - dev);
      }
    
    /* Compute power: */
    { int j = 0;
      while (j + 16 <= exp) { t *= TwoToSixteen; j += 16; }
      while (j < exp) { t *= Two; j++; }
      while (j - 16 >= exp) { t /= TwoToSixteen; j -= 16; }
      while (j > exp) { t /= Two; j--; }
    }
    
    return(t);
    
    #undef TwoToSixteen
  }