/* See flt.h */
/* Last edited on 2003-01-02 06:08:49 by stolfi */

#include "flt.h"
#include <js.h>
#include <ioprotos.h>
#include <values.h>
#include <stdlib.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;
    affirm(x >= Zero, "flt_sqrt: argument is negative");
    *zp = sqrt(x);
    ROUND_DOWN;
    zlo = sqrt(x);
    ROUND_UP;
    zhi = sqrt(x);
    UPDATE_ERR;
  }

void flt_exp(Float x, FloatP zp, FloatP errp)
  { Float zhi, zlo;
    affirm(x >= Zero, "flt_sqrt: argument is negative");
    *zp = exp(x);
    ROUND_DOWN;
    zlo = exp(x);
    ROUND_UP;
    zhi = exp(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)
  {
    return FLT_RANDOM();
  }

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 = (rand() & 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
  }