ldtoa.c

Newlib 嵌入式 C库 标准实现代码

 /* Extended precision arithmetic functions for long double I/O.
  * This program has been placed in the public domain.
  */

#include <_ansi.h>
#include <reent.h>
#include <string.h>
#include <stdlib.h>
#include "mprec.h"

/* These are the externally visible entries. */
/* linux name:  long double _IO_strtold (char *, char **); */
long double _strtold (char *, char **);
char * _ldtoa_r (struct _reent *, long double, int, int, int *, int *, char **);
int    _ldcheck (long double *);
#if 0
void _IO_ldtostr(long double *, char *, int, int, char);
#endif

 /* Number of 16 bit words in external x type format */
 #define NE 10

 /* 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)

 /* Maximum exponent digits - base 10 */
 #define MAX_EXP_DIGITS 5

/* Control structure for long double conversion including rounding precision values.
 * rndprc can be set to 80 (if NE=6), 64, 56, 53, or 24 bits.
 */
typedef struct
{
  int rlast;
  int rndprc;
  int rw;
  int re;
  int outexpon;
  unsigned short rmsk;
  unsigned short rmbit;
  unsigned short rebit;
  unsigned short rbit[NI];
  unsigned short equot[NI];
} LDPARMS;

static void esub(short unsigned int *a, short unsigned int *b, short unsigned int *c, LDPARMS *ldp);
static void emul(short unsigned int *a, short unsigned int *b, short unsigned int *c, LDPARMS *ldp);
static void ediv(short unsigned int *a, short unsigned int *b, short unsigned int *c, LDPARMS *ldp);
static int ecmp(short unsigned int *a, short unsigned int *b);
static int enormlz(short unsigned int *x);
static int eshift(short unsigned int *x, int sc);
static void eshup1(register short unsigned int *x);
static void eshup8(register short unsigned int *x);
static void eshup6(register short unsigned int *x);
static void eshdn1(register short unsigned int *x);
static void eshdn8(register short unsigned int *x);
static void eshdn6(register short unsigned int *x);
static void eneg(short unsigned int *x);
static void emov(register short unsigned int *a, register short unsigned int *b);
static void eclear(register short unsigned int *x);
static void einfin(register short unsigned int *x, register LDPARMS *ldp);
static void efloor(short unsigned int *x, short unsigned int *y, LDPARMS *ldp);
static void etoasc(short unsigned int *x, char *string, int ndigs, int outformat, LDPARMS *ldp);

union uconv
{
  unsigned short pe;
  long double d;
};

#if LDBL_MANT_DIG == 24
static void e24toe(short unsigned int *pe, short unsigned int *y, LDPARMS *ldp);
#elif LDBL_MANT_DIG == 53
static void e53toe(short unsigned int *pe, short unsigned int *y, LDPARMS *ldp);
#elif LDBL_MANT_DIG == 64
static void e64toe(short unsigned int *pe, short unsigned int *y, LDPARMS *ldp);
#else
static void e113toe(short unsigned int *pe, short unsigned int *y, LDPARMS *ldp);
#endif

/*							econst.c	*/
/*  e type constants used by high precision check routines */

#if NE == 10
/* 0.0 */
static unsigned short ezero[NE] =
 {0x0000, 0x0000, 0x0000, 0x0000,
  0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,};

/* 1.0E0 */
static unsigned short eone[NE] =
 {0x0000, 0x0000, 0x0000, 0x0000,
  0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3fff,};

#else

/* 0.0 */
static unsigned short ezero[NE] = {
0, 0000000,0000000,0000000,0000000,0000000,};
/* 1.0E0 */
static unsigned short eone[NE] = {
0, 0000000,0000000,0000000,0100000,0x3fff,};

#endif

/* Debugging routine for displaying errors */
#ifdef DEBUG
/* Notice: the order of appearance of the following
 * messages is bound to the error codes defined
 * in mconf.h.
 */
static char *ermsg[7] = {
"unknown",      /* error code 0 */
"domain",       /* error code 1 */
"singularity",  /* et seq.      */
"overflow",
"underflow",
"total loss of precision",
"partial loss of precision"
};
#define mtherr(name, code) printf( "\n%s %s error\n", name, ermsg[code] );
#else
#define mtherr(name, code)
#endif

/*							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
 *
 *	asctoe( string, e )	ASCII string to extended double e type
 *	asctoe64( string, &d )	ASCII string to long double
 *	asctoe53( string, &d )	ASCII string to double
 *	asctoe24( string, &f )	ASCII string to single
 *	asctoeg( string, e, prec, ldp ) ASCII string to specified precision
 *	e24toe( &f, e, ldp )	IEEE single precision to e type
 *	e53toe( &d, e, ldp )	IEEE double precision to e type
 *	e64toe( &d, e, ldp )	IEEE long double precision to e type
 *	e113toe( &d, e, ldp )	IEEE long double precision to e type
 *	eabs(e)			absolute value
 *	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, ldp )	c = b / a
 *	efloor( a, b, ldp )	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, ldp )	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, ldp )	c = b * a
 *	eneg(e)			e = -e
 *	eround( a, b )		b = nearest integer value to a
 *	esub( a, b, c, ldp )	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,fmt,ldp)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
 *	ltoe( &l, e )		long (32 bit) integer to e type
 *	ultoe( &l, e )		unsigned long (32 bit) integer to e type
 *      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, ldp )	divide  significands, bi = bi / ai
 *	emdnorm(ai,l,s,exp,ldp) normalize and round off
 *	emovi( a, ai )		convert external a to internal ai
 *	emovo( ai, a, ldp )	convert internal ai to external a
 *	emovz( ai, bi )		bi = ai, low guard word of bi = 0
 *	emulm( ai, bi, ldp )	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 USE_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).
 */

/*
 * 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
 * 22 Nov 00    Revised to fit into newlib by Jeff Johnston <[e-mail]>
 *
 * Author:  S. L. Moshier.
 *
 * Copyright (c) 1984,2000 S.L. Moshier
 *
 * Permission to use, copy, modify, and distribute this software for any
 * purpose without fee is hereby granted, provided that this entire notice
 * is included in all copies of any software which is or includes a copy
 * or modification of this software and in all copies of the supporting
 * documentation for such software.
 *
 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
 * WARRANTY.  IN PARTICULAR,  THE AUTHOR MAKES NO REPRESENTATION
 * OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY OF THIS
 * SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
 *
 */

#include <stdio.h>
/* #include "\usr\include\stdio.h" */
/*#include "ehead.h"*/
/*#include "mconf.h"*/
/*							mconf.h
 *
 *	Common include file for math routines
 *
 *
 *
 * SYNOPSIS:
 *
 * #include "mconf.h"
 *
 *
 *
 * DESCRIPTION:
 *
 * This file contains definitions for error codes that are
 * passed to the common error handling routine mtherr()
 * (which see).
 *
 * The file also includes a conditional assembly definition
 * for the type of computer arithmetic (IEEE, DEC, Motorola
 * IEEE, or UNKnown).
 *
 * For Digital Equipment PDP-11 and VAX computers, certain
 * IBM systems, and others that use numbers with a 56-bit
 * significand, the symbol DEC should be defined.  In this
 * mode, most floating point constants are given as arrays
 * of octal integers to eliminate decimal to binary conversion
 * errors that might be introduced by the compiler.
 *
 * For computers, such as IBM PC, that follow the IEEE 
 * Standard for Binary Floating Point Arithmetic (ANSI/IEEE
 * Std 754-1985), the symbol IBMPC should be defined.  These
 * numbers have 53-bit significands.  In this mode, constants
 * are provided as arrays of hexadecimal 16 bit integers.
 *
 * To accommodate other types of computer arithmetic, all
 * constants are also provided in a normal decimal radix
 * which one can hope are correctly converted to a suitable
 * format by the available C language compiler.  To invoke
 * this mode, the symbol UNK is defined.
 *
 * An important difference among these modes is a predefined
 * set of machine arithmetic constants for each.  The numbers
 * MACHEP (the machine roundoff error), MAXNUM (largest number
 * represented), and several other parameters are preset by
 * the configuration symbol.  Check the file const.c to
 * ensure that these values are correct for your computer.
 *
 * For ANSI C compatibility, define ANSIC equal to 1.  Currently
 * this affects only the atan2() function and others that use it.
 */

/* Constant definitions for math error conditions
 */

#define DOMAIN		1	/* argument domain error */
#define SING		2	/* argument singularity */
#define OVERFLOW	3	/* overflow range error */
#define UNDERFLOW	4	/* underflow range error */
#define TLOSS		5	/* total loss of precision */
#define PLOSS		6	/* partial loss of precision */

#define EDOM		33
#define ERANGE		34

typedef struct
	{
	double r;
	double i;
	}cmplx;

/* Type of computer arithmetic */

#ifndef DEC
#ifdef __IEEE_LITTLE_ENDIAN
#define IBMPC 1
#else  /* !__IEEE_LITTLE_ENDIAN */
#define MIEEE 1
#endif /* !__IEEE_LITTLE_ENDIAN */
#endif /* !DEC */

/* Define 1 for ANSI C atan2() function
 * See atan.c and clog.c.
 */
#define ANSIC 1

/*define VOLATILE volatile*/
#define VOLATILE 

#define NANS
#define USE_INFINITY

/* NaN's require infinity support. */
#ifdef NANS
#ifndef INFINITY
#define USE_INFINITY
#endif
#endif

/* This handles 64-bit long ints. */
#define LONGBITS (8 * sizeof(long))


static void eaddm(short unsigned int *x, short unsigned int *y);
static void esubm(short unsigned int *x, short unsigned int *y);
static void emdnorm(short unsigned int *s, int lost, int subflg, long int exp, int rcntrl, LDPARMS *ldp);
static int  asctoeg(char *ss, short unsigned int *y, int oprec, LDPARMS *ldp);
static void enan(short unsigned int *nan, int size);
#if LDBL_MANT_DIG == 24
static void toe24(short unsigned int *x, short unsigned int *y);
#elif LDBL_MANT_DIG == 53
static void toe53(short unsigned int *x, short unsigned int *y);
#elif LDBL_MANT_DIG == 64
static void toe64(short unsigned int *a, short unsigned int *b);
#else
static void toe113(short unsigned int *a, short unsigned int *b);
#endif
static void eiremain(short unsigned int *den, short unsigned int *num, LDPARMS *ldp);
static int ecmpm(register short unsigned int *a, register short unsigned int *b);
static int edivm(short unsigned int *den, short unsigned int *num, LDPARMS *ldp);
static int emulm(short unsigned int *a, short unsigned int *b, LDPARMS *ldp);
static int eisneg(short unsigned int *x);
static int eisinf(short unsigned int *x);
static void emovi(short unsigned int *a, short unsigned int *b);
static void emovo(short unsigned int *a, short unsigned int *b, LDPARMS *ldp);
static void emovz(register short unsigned int *a, register short unsigned int *b);
static void ecleaz(register short unsigned int *xi);
static void eadd1(short unsigned int *a, short unsigned int *b, short unsigned int *c, int subflg, LDPARMS *ldp);
static int eisnan(short unsigned int *x);
static int eiisnan(short unsigned int *x);

#ifdef DEC
static void etodec(), todec(), dectoe();
#endif

/*
; Clear out entire external format number.
;
; unsigned short x[];
; eclear( x );
*/

static void eclear(register short unsigned int *x)
{
register int i;

for( i=0; i<NE; i++ )
	*x++ = 0;
}



/* Move external format number from a to b.
 *
 * emov( a, b );
 */

static void emov(register short unsigned int *a, register short unsigned int *b)
{
register int i;

for( i=0; i<NE; i++ )
	*b++ = *a++;
}


/*
;	Negate external format number
;
;	unsigned short x[NE];
;	eneg( x );
*/

static 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 int eisneg(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 int eisinf(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 int eisnan(short unsigned int *x)
{

#ifdef NANS
int i;
/* NaN has maximum exponent */
if( (x[NE-1] & 0x7fff) != 0x7fff )
	return (0);
/* ... 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 void einfin(register short unsigned int *x, register LDPARMS *ldp)
{
register int i;

#ifdef USE_INFINITY
for( i=0; i<NE-1; i++ )
	*x++ = 0;
*x |= 32767;
ldp = ldp;
#else
for( i=0; i<NE-1; i++ )
	*x++ = 0xffff;
*x |= 32766;
if( ldp->rndprc < NBITS )
	{
	if (ldp->rndprc == 113)
		{
		*(x - 9) = 0;
		*(x - 8) = 0;
		}
	if( ldp->rndprc == 64 )
		{
		*(x-5) = 0;
		}
	if( ldp->rndprc == 53 )
		{
		*(x-4) = 0xf800;
		}
	else
		{
		*(x-4) = 0;
		*(x-3) = 0;
		*(x-2) = 0xff00;
		}
	}
#endif
}

/* Move in external format number,
 * converting it to internal format.
 */
static void emovi(short unsigned int *a, short unsigned int *b)
{
register unsigned short *p, *q;
int i;

q = b;
p = a + (NE-1);	/* point to last word of external number */
/* get the sign bit */
if( *p & 0x8000 )
	*q++ = 0xffff;
else
	*q++ = 0;
/* get the exponent */
*q = *p--;
*q++ &= 0x7fff;	/* delete the sign bit */
#ifdef USE_INFINITY
if( (*(q-1) & 0x7fff) == 0x7fff )
	{
#ifdef NANS
	if( eisnan(a) )
		{
		*q++ = 0;
		for( i=3; i<NI; i++ )
			*q++ = *p--;
		return;
		}
#endif
	for( i=2; i<NI; i++ )
		*q++ = 0;
	return;
	}
#endif
/* clear high guard word */
*q++ = 0;
/* move in the significand */
for( i=0; i<NE-1; i++ )
	*q++ = *p--;
/* clear low guard word */
*q = 0;
}


/* Move internal format number out,
 * converting it to external format.
 */
static void emovo(short unsigned int *a, short unsigned int *b, LDPARMS *ldp)
{
register unsigned short *p, *q;
unsigned short i;

p = a;
q = b + (NE-1); /* point to output exponent */
/* combine sign and exponent */
i = *p++;
if( i )
	*q-- = *p++ | 0x8000;
else
	*q-- = *p++;
#ifdef USE_INFINITY
if( *(p-1) == 0x7fff )
	{
#ifdef NANS
	if( eiisnan(a) )
		{
		enan( b, NBITS );
		return;
		}
#endif
	einfin(b, ldp);
	return;
	}
#endif
/* skip over guard word */
++p;
/* move the significand */
for( i=0; i<NE-1; i++ )
	*q-- = *p++;
}


/* Clear out internal format number.
 */

static void ecleaz(register short unsigned int *xi)
{
register int i;

for( i=0; i<NI; i++ )
	*xi++ = 0;
}

/* same, but don't touch the sign. */

static void ecleazs(register short unsigned int *xi)
{
register int i;

++xi;
for(i=0; i<NI-1; i++)
	*xi++ = 0;
}