/* Some extra floating-point tools, used by ia.c and aa.c  */

/****************************************************************************/
/* (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.                              */
/****************************************************************************/

#ifndef FLT_H
#define FLT_H

#include <stdio.h>
#include <sys/ieeefp.h>

/* The floating-point type used by the operations below.: */
typedef
  float Float;

typedef
  float *FloatP;

/* Default $printf$ parameters for type Float: */

#define FLT_FMT_DIGITS 8
#define FLT_FMT_WIDTH (FLT_FMT_DIGITS + 6)
#define FLT_FMT_DECIMALS (FLT_FMT_DIGITS - 1)
#define FLT_FMT_SPEC "%14.7e"

/* Common constants for type Float: */

#define One     (1.0)
#define Zero    (0.0)
#define Half    (0.5)
#define Quarter (0.25)
#define Two     (2.0)
#define Three   (3.0)
#define Four    (4.0)
#define Eight   (8.0)
#define Sixteen (16.0)

extern Float MaxFloat;
extern Float Infinity;

#define MinusInfinity (-Infinity)
#define PlusInfinity Infinity

/* Basic operations for type Float */
#define FABS(x)   ((x) > Zero ? (x) : -(x))
#define FMAX(x,y) ((x) > (y) ? (x) : (y))
#define FMIN(x,y) ((x) > (y) ? (y) : (x))

/* Procedures: */

void flt_init (void);
  /* Initializes the constants (Infinity, MaxFloat, MinusInfinity, etc.) */
  /* This routine MUST be called at least once before any other routine below. */

void flt_print (FILE *f, Float x);
  /* Prints $x$ on $f$, using "%e" format with full precision. */

/* Error-monitoring arithmetic routines */
/* The routines below perform *zp = x OP y, where OP is the appropriate */
/*   arithmetic operation, using the current rounding mode. */
/* They also add to *errp a quantity d that is an upper bound to */
/*   the error made in the computation of *zp. */
/* They will return *errp = PlusInfinity in case of overflow. */
/* WARNING: they may change the current rounding mode. */

void flt_add(Float x, Float y, FloatP zp, FloatP errp);
  /* Computes $*zp = x + y$.  */

void flt_sub(Float x, Float y, FloatP zp, FloatP errp);
  /* Computes $*zp = x - y$.  */

void flt_mul(Float x, Float y, FloatP zp, FloatP errp);
  /* Computes $*zp = x * y$. */

void flt_div(Float x, Float y, FloatP zp, FloatP errp);
  /* Computes $*zp = x / y$. */

void flt_inv(Float x, FloatP zp, FloatP errp);
  /* Computes $*zp = 1/x$. */

void flt_sqrt(Float x, FloatP zp, FloatP errp);
  /* Computes $*zp = sqrt(x)$. */

void flt_scale (
    Float x, Float alpha, 
    Float zeta,
    FloatP zp,
    FloatP errp
  );
  /* Computes $*zp = \alpha x / zeta$. */

void flt_mix (
    Float x, Float alpha, 
    Float y, Float beta, 
    Float zeta,
    FloatP zp,
    FloatP errp
  );
  /* Computes $*zp = ( \alpha x + \beta y ) / \zeta$. */

Float flt_random(void);
  /* A random Float value in [0 __ 1).
     The client must have called $srandom(<seed>)$. */

Float flt_random_mag(int avg, int dev);
  /* Two raised to a random exponent, uniform between 
     $avg-dev$ and $avg+dev$.  (The result may be infinity). 
     The client must have called $srandom(<seed>)$. */

/* Setting the IEEE rounding mode bits on a SPARC: */

#define ROUND_DOWN  _flt_round_down()
#define ROUND_UP    _flt_round_up()
#define ROUND_NEAR  _flt_round_near()
#define ROUND_ZERO  _flt_round_zero()

void _flt_round_down (void);
void _flt_round_up   (void);
void _flt_round_near (void);
void _flt_round_zero (void);
  /* These routines set the IEEE FP rounding direction. */
  /* They should be a lot faster than "ieee_flags". */

int flt_get_fsr(void);
  /* Returns the FP status register. (For debugging) */

#endif