#ifndef _CEPHES_EMATH_H
#define _CEPHES_EMATH_H

/* This file is extracted from S L Moshier's  ioldoubl.c,
 * modified for use in MinGW 
 *
 * Extended precision arithmetic functions for long double I/O.
 * This program has been placed in the public domain.
 */

   

/*
 * Revision history:
 *
 *  5 Jan 84	PDP-11 assembly language version
 *  6 Dec 86	C language version
 * 30 Aug 88	100 digit version, improved rounding
 * 15 May 92    80-bit long double support
 *
 * Author:  S. L. Moshier.
 *
 * 6 Oct 02	Modified for MinGW by inlining utility routines,
 * 		removing global variables, and splitting out strtold
 *		from _IO_ldtoa and _IO_ldtostr.
 *  
 *		Danny Smith <dannysmith@users.sourceforge.net>
 * 
 */

  

/*							ieee.c
 *
 *    Extended precision IEEE binary floating point arithmetic routines
 *
 * Numbers are stored in C language as arrays of 16-bit unsigned
 * short integers.  The arguments of the routines are pointers to
 * the arrays.
 *
 *
 * External e type data structure, simulates Intel 8087 chip
 * temporary real format but possibly with a larger significand:
 *
 *	NE-1 significand words	(least significant word first,
 *				 most significant bit is normally set)
 *	exponent		(value = EXONE for 1.0,
 *				top bit is the sign)
 *
 *
 * Internal data structure of a number (a "word" is 16 bits):
 *
 * ei[0]	sign word	(0 for positive, 0xffff for negative)
 * ei[1]	biased __exponent	(value = EXONE for the number 1.0)
 * ei[2]	high guard word	(always zero after normalization)
 * ei[3]
 * to ei[NI-2]	significand	(NI-4 significand words,
 *				 most significant word first,
 *				 most significant bit is set)
 * ei[NI-1]	low guard word	(0x8000 bit is rounding place)
 *
 *
 *
 *		Routines for external format numbers
 *
 *	__asctoe64( string, &d )	ASCII string to long double
 *	__asctoeg( string, e, prec )	ASCII string to specified precision
 *	__e64toe( &d, e )		IEEE long double precision to e type
 *	__eadd( a, b, c )		c = b + a
 *	__eclear(e)			e = 0
 *	__ecmp (a, b)			Returns 1 if a > b, 0 if a == b,
 *					-1 if a < b, -2 if either a or b is a NaN.
 *	__ediv( a, b, c )		c = b / a
 *	__efloor( a, b )		truncate to integer, toward -infinity
 *	__efrexp( a, exp, s )		extract exponent and significand
 *	__eifrac( e, &l, frac )   	e to long integer and e type fraction
 *	__euifrac( e, &l, frac )  	e to unsigned long integer and e type fraction
 *	__einfin( e )			set e to infinity, leaving its sign alone
 *	__eldexp( a, n, b )		multiply by 2**n
 *	__emov( a, b )			b = a
 *	__emul( a, b, c )		c = b * a
 *	__eneg(e)			e = -e
 *	__eround( a, b )		b = nearest integer value to a
 *	__esub( a, b, c )		c = b - a
 *	__e24toasc( &f, str, n )	single to ASCII string, n digits after decimal
 *	__e53toasc( &d, str, n )	double to ASCII string, n digits after decimal
 *	__e64toasc( &d, str, n )	long double to ASCII string
 *	__etoasc( e, str, n )		e to ASCII string, n digits after decimal
 *	__etoe24( e, &f )		convert e type to IEEE single precision
 *	__etoe53( e, &d )		convert e type to IEEE double precision
 *	__etoe64( e, &d )		convert e type to IEEE long double precision
 *	__eisneg( e )             	1 if sign bit of e != 0, else 0
 *	__eisinf( e )             	1 if e has maximum exponent (non-IEEE)
 *					or is infinite (IEEE)
 *	__eisnan( e )             	1 if e is a NaN
 *	__esqrt( a, b )			b = square root of a
 *
 *
 *		Routines for internal format numbers
 *
 *	__eaddm( ai, bi )		add significands, bi = bi + ai
 *	__ecleaz(ei)		ei = 0
 *	__ecleazs(ei)		set ei = 0 but leave its sign alone
 *	__ecmpm( ai, bi )		compare significands, return 1, 0, or -1
 *	__edivm( ai, bi )		divide  significands, bi = bi / ai
 *	__emdnorm(ai,l,s,exp)	normalize and round off
 *	__emovi( a, ai )		convert external a to internal ai
 *	__emovo( ai, a )		convert internal ai to external a
 *	__emovz( ai, bi )		bi = ai, low guard word of bi = 0
 *	__emulm( ai, bi )		multiply significands, bi = bi * ai
 *	__enormlz(ei)		left-justify the significand
 *	__eshdn1( ai )		shift significand and guards down 1 bit
 *	__eshdn8( ai )		shift down 8 bits
 *	__eshdn6( ai )		shift down 16 bits
 *	__eshift( ai, n )		shift ai n bits up (or down if n < 0)
 *	__eshup1( ai )		shift significand and guards up 1 bit
 *	__eshup8( ai )		shift up 8 bits
 *	__eshup6( ai )		shift up 16 bits
 *	__esubm( ai, bi )		subtract significands, bi = bi - ai
 *
 *
 * The result is always normalized and rounded to NI-4 word precision
 * after each arithmetic operation.
 *
 * Exception flags are NOT fully supported.
 *
 * Define INFINITY in mconf.h for support of infinity; otherwise a
 * saturation arithmetic is implemented.
 *
 * Define NANS for support of Not-a-Number items; otherwise the
 * arithmetic will never produce a NaN output, and might be confused
 * by a NaN input.
 * If NaN's are supported, the output of ecmp(a,b) is -2 if
 * either a or b is a NaN. This means asking if(ecmp(a,b) < 0)
 * may not be legitimate. Use if(ecmp(a,b) == -1) for less-than
 * if in doubt.
 * Signaling NaN's are NOT supported; they are treated the same
 * as quiet NaN's.
 *
 * Denormals are always supported here where appropriate (e.g., not
 * for conversion to DEC numbers).
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <errno.h>
#include <math.h>
#include <locale.h>
#include <ctype.h>

#define alloca __builtin_alloca

/* Don't build non-ANSI _IO_ldtoa.  It is not thread safe. */ 
#ifndef USE_LDTOA
#define USE_LDTOA 0
#endif


 /* Number of 16 bit words in external x type format */
#define NE 6

 /* Number of 16 bit words in internal format */
#define NI (NE+3)

 /* Array offset to exponent */
#define E 1

 /* Array offset to high guard word */
#define M 2

 /* Number of bits of precision */
#define NBITS ((NI-4)*16)

 /* Maximum number of decimal digits in ASCII conversion
  * = NBITS*log10(2)
  */
#define NDEC (NBITS*8/27)

 /* The exponent of 1.0 */
#define EXONE (0x3fff)


#define  mtherr(x,y) 


extern long double strtold (const char * __restrict__ s, char ** __restrict__ se);
extern int __asctoe64(const char * __restrict__ ss,
	       short unsigned int * __restrict__ y);
extern void __emul(const short unsigned int *  a,
		 const short unsigned int *  b,
		 short unsigned int * c);
extern int __ecmp(const short unsigned int * __restrict__ a,
		const short unsigned int *  __restrict__ b);
extern int __enormlz(short unsigned int *x);
extern int __eshift(short unsigned int *x, int sc);
extern void __eaddm(const short unsigned int  *  __restrict__  x,
	          short unsigned int *  __restrict__  y);
extern void __esubm(const short unsigned int * __restrict__  x,
		  short unsigned int *  __restrict__ y);
extern void __emdnorm(short unsigned int *s, int lost, int subflg,
		      int exp, int rcntrl, const int rndprc);
extern void __toe64(short unsigned int *  __restrict__  a,
		  short unsigned int *  __restrict__  b);
extern int __edivm(short unsigned int *  __restrict__  den,
		 short unsigned int * __restrict__  num);
extern int __emulm(const short unsigned int *  __restrict__ a,
		 short unsigned int *  __restrict__ b);
extern void __emovi(const short unsigned int * __restrict__ a,
		  short unsigned int * __restrict__ b);
extern void __emovo(const short unsigned int * __restrict__ a,
		  short unsigned int * __restrict__ b);

#if USE_LDTOA

extern char * _IO_ldtoa(long double, int, int, int *, int *, char **);
extern void _IO_ldtostr(long double *x, char *string, int ndigs,
			int flags, char fmt);

extern void __eiremain(short unsigned int * __restrict__ den,
		       short unsigned int *__restrict__ num,
		       short unsigned int *__restrict__ equot);
extern void __efloor(short unsigned int *x, short unsigned int *y);
extern void __eadd1(const short unsigned int * __restrict__ a,
		  const short unsigned int * __restrict__ b,
		  short unsigned int * __restrict__ c,
		  int subflg);
extern void __esub(const short unsigned int *a, const short unsigned int *b,
		   short unsigned int *c);
extern void __ediv(const short unsigned int *a, const short unsigned int *b,
		   short unsigned int *c);
extern void __e64toe(short unsigned int *pe, short unsigned int *y);


#endif

static  __inline__ int __eisneg(const short unsigned int *x);
static  __inline__ int __eisinf(const short unsigned int *x);
static __inline__ int __eisnan(const short unsigned int *x);
static __inline__ int __eiszero(const short unsigned int *a);
static __inline__ void __emovz(register const short unsigned int * __restrict__ a,
			       register short unsigned int * __restrict__ b);
static __inline__ void __eclear(register short unsigned int *x);
static __inline__ void __ecleaz(register short unsigned int *xi);
static __inline__ void __ecleazs(register short unsigned int *xi);
static  __inline__ int __eiisinf(const short unsigned int *x);
static __inline__ int __eiisnan(const short unsigned int *x);
static __inline__ int __eiiszero(const short unsigned int *x);
static __inline__ void __enan_64(short unsigned int *nan);
static __inline__ void __enan_NBITS (short unsigned int *nan);
static __inline__ void __enan_NI16 (short unsigned int *nan);
static __inline__ void __einfin(register short unsigned int *x);
static __inline__ void __eneg(short unsigned int *x);
static __inline__ void __eshup1(register short unsigned int *x);
static __inline__ void __eshup8(register short unsigned int *x);
static __inline__ void __eshup6(register short unsigned int *x);
static __inline__ void __eshdn1(register short unsigned int *x);
static __inline__ void __eshdn8(register short unsigned int *x);
static __inline__ void __eshdn6(register short unsigned int *x);



/* Intel IEEE, low order words come first:
 */
#define IBMPC 1

/* Define 1 for ANSI C atan2() function
 * See atan.c and clog.c.
 */
#define ANSIC 1

/*define VOLATILE volatile*/
#define VOLATILE 

/* For 12-byte long doubles on an i386, pad a 16-bit short 0
 * to the end of real constants initialized by integer arrays.
 *
 * #define XPD 0,
 *
 * Otherwise, the type is 10 bytes long and XPD should be
 * defined blank.
 *
 * #define XPD
 */
#define XPD 0,
/* #define XPD */
#define NANS

/* NaN's require infinity support. */
#ifdef NANS
#ifndef INFINITY
#define INFINITY
#endif
#endif

/* This handles 64-bit long ints. */
#define LONGBITS (8 * sizeof(long))


#define NTEN 12
#define MAXP 4096

/*
; Clear out entire external format number.
;
; unsigned short x[];
; eclear( x );
*/

static __inline__ void __eclear(register short unsigned int *x)
{
  memset(x, 0, NE * sizeof(unsigned short));
}


/* Move external format number from a to b.
 *
 * emov( a, b );
 */

static __inline__ void __emov(register const short unsigned int * __restrict__ a,
			    register short unsigned int * __restrict__ b)
{
  memcpy(b, a, NE * sizeof(unsigned short));
}


/*
;	Negate external format number
;
;	unsigned short x[NE];
;	eneg( x );
*/

static __inline__ void __eneg(short unsigned int *x)
{

#ifdef NANS
if( __eisnan(x) )
	return;
#endif
x[NE-1] ^= 0x8000; /* Toggle the sign bit */
}


/* Return 1 if external format number is negative,
 * else return zero.
 */
static __inline__ int __eisneg(const short unsigned int *x)
{

#ifdef NANS
if( __eisnan(x) )
	return( 0 );
#endif
if( x[NE-1] & 0x8000 )
	return( 1 );
else
	return( 0 );
}


/* Return 1 if external format number has maximum possible exponent,
 * else return zero.
 */
static __inline__ int __eisinf(const short unsigned int *x)
{

if( (x[NE-1] & 0x7fff) == 0x7fff )
	{
#ifdef NANS
	if( __eisnan(x) )
		return( 0 );
#endif
	return( 1 );
	}
else
	return( 0 );
}

/* Check if e-type number is not a number.
 */
static __inline__ int __eisnan(const short unsigned int *x)
{
#ifdef NANS
int i;
/* NaN has maximum __exponent */
if( (x[NE-1] & 0x7fff) == 0x7fff )
/* ... and non-zero significand field. */
    for( i=0; i<NE-1; i++ )
	{
	if( *x++ != 0 )
		return (1);
	}
#endif
return (0);
}

/*
; Fill __entire number, including __exponent and significand, with
; largest possible number.  These programs implement a saturation
; value that is an ordinary, legal number.  A special value
; "infinity" may also be implemented; this would require tests
; for that value and implementation of special rules for arithmetic
; operations involving inifinity.
*/

static __inline__ void __einfin(register short unsigned int *x)
{
register int i;

#ifdef INFINITY
for( i=0; i<NE-1; i++ )
	*x++ = 0;
*x |= 32767;
#else
for( i=0; i<NE-1; i++ )
	*x++ = 0xffff;
*x |= 32766;
*(x-5) = 0;
#endif
}

/* Clear out internal format number.
 */

static __inline__ void __ecleaz(register short unsigned int *xi)
{
  memset(xi, 0, NI * sizeof(unsigned short));
}

/* same, but don't touch the sign. */

static __inline__ void __ecleazs(register short unsigned int *xi)
{
  ++xi;
  memset(xi, 0, (NI-1) * sizeof(unsigned short));
}



/* Move internal format number from a to b.
 */
static __inline__ void __emovz(register const short unsigned int * __restrict__ a,
			     register short unsigned int * __restrict__ b)
{
  memcpy(b, a, (NI-1) * sizeof(unsigned short));
  b[NI-1]=0;
}

/* Return nonzero if internal format number is a NaN.
 */

static __inline__ int __eiisnan (const short unsigned int *x)
{
int i;

if( (x[E] & 0x7fff) == 0x7fff )
	{
	for( i=M+1; i<NI; i++ )
		{
		if( x[i] != 0 )
			return(1);
		}
	}
return(0);
}

/* Return nonzero if external format number is zero. */

static __inline__ int
__eiszero(const short unsigned int * a)
{
if (*((long double*) a) == 0)
	return (1);
return (0);
}

/* Return nonzero if internal format number is zero. */

static __inline__ int
__eiiszero(const short unsigned int * ai)
{
  int i;
  /* skip the sign word */
  for( i=1; i<NI-1; i++ )
    {
      if( ai[i] != 0 )
        return (0);
    }
  return (1);
}


/* Return nonzero if internal format number is infinite. */

static __inline__ int 
__eiisinf (const unsigned short *x)
{

#ifdef NANS
  if (__eiisnan (x))
    return (0);
#endif
  if ((x[E] & 0x7fff) == 0x7fff)
    return (1);
  return (0);
}

/*
;	Compare significands of numbers in internal format.
;	Guard words are included in the comparison.
;
;	unsigned short a[NI], b[NI];
;	cmpm( a, b );
;
;	for the significands:
;	returns	+1 if a > b
;		 0 if a == b
;		-1 if a < b
*/
static __inline__ int __ecmpm(register const short unsigned int * __restrict__ a,
			    register const short unsigned int *  __restrict__ b)
{
int i;

a += M; /* skip up to significand area */
b += M;
for( i=M; i<NI; i++ )
	{
	if( *a++ != *b++ )
		goto difrnt;
	}
return(0);

difrnt:
if( *(--a) > *(--b) )
	return(1);
else
	return(-1);
}


/*
;	Shift significand down by 1 bit
*/

static __inline__ void __eshdn1(register short unsigned int *x)
{
register unsigned short bits;
int i;

x += M;	/* point to significand area */

bits = 0;
for( i=M; i<NI; i++ )
	{
	if( *x & 1 )
		bits |= 1;
	*x >>= 1;
	if( bits & 2 )
		*x |= 0x8000;
	bits <<= 1;
	++x;
	}	
}

/*
;	Shift significand up by 1 bit
*/

static __inline__ void __eshup1(register short unsigned int *x)
{
register unsigned short bits;
int i;

x += NI-1;
bits = 0;

for( i=M; i<NI; i++ )
	{
	if( *x & 0x8000 )
		bits |= 1;
	*x <<= 1;
	if( bits & 2 )
		*x |= 1;
	bits <<= 1;
	--x;
	}
}



/*
;	Shift significand down by 8 bits
*/

static __inline__ void __eshdn8(register short unsigned int *x)
{
register unsigned short newbyt, oldbyt;
int i;

x += M;
oldbyt = 0;
for( i=M; i<NI; i++ )
	{
	newbyt = *x << 8;
	*x >>= 8;
	*x |= oldbyt;
	oldbyt = newbyt;
	++x;
	}
}

/*
;	Shift significand up by 8 bits
*/

static __inline__ void __eshup8(register short unsigned int *x)
{
int i;
register unsigned short newbyt, oldbyt;

x += NI-1;
oldbyt = 0;

for( i=M; i<NI; i++ )
	{
	newbyt = *x >> 8;
	*x <<= 8;
	*x |= oldbyt;
	oldbyt = newbyt;
	--x;
	}
}

/*
;	Shift significand up by 16 bits
*/

static __inline__ void __eshup6(register short unsigned int *x)
{
int i;
register unsigned short *p;

p = x + M;
x += M + 1;

for( i=M; i<NI-1; i++ )
	*p++ = *x++;

*p = 0;
}

/*
;	Shift significand down by 16 bits
*/

static __inline__ void __eshdn6(register short unsigned int *x)
{
int i;
register unsigned short *p;

x += NI-1;
p = x + 1;

for( i=M; i<NI-1; i++ )
	*(--p) = *(--x);

*(--p) = 0;
}

/*
;	Add significands
;	x + y replaces y
*/

static __inline__ void __enan_64(unsigned short* nan)
{

  int i;
  for( i=0; i<3; i++ )
    *nan++ = 0;
  *nan++ = 0xc000;
  *nan++ = 0x7fff;
  *nan = 0;
  return;
}

static __inline__ void __enan_NBITS(unsigned short* nan)
{
  int i; 
  for( i=0; i<NE-2; i++ )
    *nan++ = 0;
  *nan++ = 0xc000;
  *nan = 0x7fff;
  return;
}

static __inline__ void __enan_NI16(unsigned short* nan)
{
  int i; 
  *nan++ = 0;
  *nan++ = 0x7fff;
  *nan++ = 0;
  *nan++ = 0xc000;
  for( i=4; i<NI; i++ )
    *nan++ = 0;
  return;
}

#endif /* _CEPHES_EMATH_H */