diff options
Diffstat (limited to 'newlib/libm/common')
32 files changed, 3187 insertions, 0 deletions
diff --git a/newlib/libm/common/Makefile.am b/newlib/libm/common/Makefile.am new file mode 100644 index 000000000..57b3ff64d --- /dev/null +++ b/newlib/libm/common/Makefile.am @@ -0,0 +1,97 @@ +## Process this file with automake to generate Makefile.in + +AUTOMAKE_OPTIONS = cygnus + +INCLUDES = $(NEWLIB_CFLAGS) $(CROSS_CFLAGS) $(TARGET_CFLAGS) + +noinst_LIBRARIES = lib.a + +src = s_finite.c s_copysign.c s_modf.c s_scalbn.c \ + s_cbrt.c s_expm1.c s_ilogb.c \ + s_infinity.c s_log1p.c s_nan.c s_nextafter.c \ + s_rint.c s_logb.c s_matherr.c s_lib_ver.c + +fsrc = sf_finite.c sf_copysign.c sf_modf.c sf_scalbn.c \ + sf_cbrt.c sf_expm1.c sf_ilogb.c \ + sf_infinity.c sf_log1p.c sf_nan.c sf_nextafter.c \ + sf_rint.c sf_logb.c + +lib_a_SOURCES = $(src) $(fsrc) + +chobj = scbrt.def scopysign.def sexpm1.def silogb.def \ + sinfinity.def slog1p.def smatherr.def smodf.def \ + snan.def snextafter.def sscalbn.def + +SUFFIXES = .def + +CHEW = ../../doc/makedoc -f $(srcdir)/../../doc/doc.str + +.c.def: + $(CHEW) < $< > $*.def 2> $*.ref + touch stmp-def + +TARGETDOC = ../tmp.texi + +doc: $(chobj) + cat $(srcdir)/common.tex >> $(TARGETDOC) + +CLEANFILES = $(chobj) *.ref + +# Texinfo does not appear to support underscores in file names, so we +# name the .def files without underscores. + +smodf.def: s_modf.c + $(CHEW) < $(srcdir)/s_modf.c >$@ 2>/dev/null + touch stmp-def + +scopysign.def: s_copysign.c + $(CHEW) < $(srcdir)/s_copysign.c >$@ 2>/dev/null + touch stmp-def + +sscalbn.def: s_scalbn.c + $(CHEW) < $(srcdir)/s_scalbn.c >$@ 2>/dev/null + touch stmp-def + +scbrt.def: s_cbrt.c + $(CHEW) < $(srcdir)/s_cbrt.c >$@ 2>/dev/null + touch stmp-def + +serf.def: s_erf.c + $(CHEW) < $(srcdir)/s_serf.c >$@ 2>/dev/null + touch stmp-def + +sexpn.def: s_expm.c + $(CHEW) < $(srcdir)/s_expn.c >$@ 2>/dev/null + touch stmp-def + +sexpm1.def: s_expm1.c + $(CHEW) < $(srcdir)/s_expm1.c >$@ 2>/dev/null + touch stmp-def + +silogb.def: s_ilogb.c + $(CHEW) < $(srcdir)/s_ilogb.c >$@ 2>/dev/null + touch stmp-def + +sinfinity.def: s_infinity.c + $(CHEW) < $(srcdir)/s_infinity.c >$@ 2>/dev/null + touch stmp-def + +slog1p.def: s_log1p.c + $(CHEW) < $(srcdir)/s_log1p.c >$@ 2>/dev/null + touch stmp-def + +smatherr.def: s_matherr.c + $(CHEW) < $(srcdir)/s_matherr.c >$@ 2>/dev/null + touch stmp-def + +snan.def: s_nan.c + $(CHEW) < $(srcdir)/s_nan.c >$@ 2>/dev/null + touch stmp-def + +snextafter.def: s_nextafter.c + $(CHEW) < $(srcdir)/s_nextafter.c >$@ 2>/dev/null + touch stmp-def + +# A partial dependency list. + +$(lib_a_OBJECTS): $(srcdir)/../../libc/include/math.h fdlibm.h diff --git a/newlib/libm/common/Makefile.in b/newlib/libm/common/Makefile.in new file mode 100644 index 000000000..3ba380d8f --- /dev/null +++ b/newlib/libm/common/Makefile.in @@ -0,0 +1,341 @@ +# Makefile.in generated automatically by automake 1.3b from Makefile.am + +# Copyright (C) 1994, 1995, 1996, 1997, 1998 Free Software Foundation, Inc. +# This Makefile.in is free software; the Free Software Foundation +# gives unlimited permission to copy and/or distribute it, +# with or without modifications, as long as this notice is preserved. + +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY, to the extent permitted by law; without +# even the implied warranty of MERCHANTABILITY or FITNESS FOR A +# PARTICULAR PURPOSE. + + +SHELL = @SHELL@ + +srcdir = @srcdir@ +top_srcdir = @top_srcdir@ +VPATH = @srcdir@ +prefix = @prefix@ +exec_prefix = @exec_prefix@ + +bindir = @bindir@ +sbindir = @sbindir@ +libexecdir = @libexecdir@ +datadir = @datadir@ +sysconfdir = @sysconfdir@ +sharedstatedir = @sharedstatedir@ +localstatedir = @localstatedir@ +libdir = @libdir@ +infodir = @infodir@ +mandir = @mandir@ +includedir = @includedir@ +oldincludedir = /usr/include + +DESTDIR = + +pkgdatadir = $(datadir)/@PACKAGE@ +pkglibdir = $(libdir)/@PACKAGE@ +pkgincludedir = $(includedir)/@PACKAGE@ + +top_builddir = .. + +ACLOCAL = @ACLOCAL@ +AUTOCONF = @AUTOCONF@ +AUTOMAKE = @AUTOMAKE@ +AUTOHEADER = @AUTOHEADER@ + +INSTALL = @INSTALL@ +INSTALL_PROGRAM = @INSTALL_PROGRAM@ +INSTALL_DATA = @INSTALL_DATA@ +INSTALL_SCRIPT = @INSTALL_SCRIPT@ +transform = @program_transform_name@ + +NORMAL_INSTALL = : +PRE_INSTALL = : +POST_INSTALL = : +NORMAL_UNINSTALL = : +PRE_UNINSTALL = : +POST_UNINSTALL = : +host_alias = @host_alias@ +host_triplet = @host@ +AR = @AR@ +AS = @AS@ +CC = @CC@ +CPP = @CPP@ +EXEEXT = @EXEEXT@ +MAINT = @MAINT@ +MAKEINFO = @MAKEINFO@ +NEWLIB_CFLAGS = @NEWLIB_CFLAGS@ +PACKAGE = @PACKAGE@ +RANLIB = @RANLIB@ +VERSION = @VERSION@ +mach_add_objs = @mach_add_objs@ +machine_dir = @machine_dir@ +newlib_basedir = @newlib_basedir@ +sys_dir = @sys_dir@ + +AUTOMAKE_OPTIONS = cygnus + +INCLUDES = $(NEWLIB_CFLAGS) $(CROSS_CFLAGS) $(TARGET_CFLAGS) + +noinst_LIBRARIES = lib.a + +src = s_finite.c s_copysign.c s_modf.c s_scalbn.c \ + s_cbrt.c s_expm1.c s_ilogb.c \ + s_infinity.c s_log1p.c s_nan.c s_nextafter.c \ + s_rint.c s_logb.c s_matherr.c s_lib_ver.c + +fsrc = sf_finite.c sf_copysign.c sf_modf.c sf_scalbn.c \ + sf_cbrt.c sf_expm1.c sf_ilogb.c \ + sf_infinity.c sf_log1p.c sf_nan.c sf_nextafter.c \ + sf_rint.c sf_logb.c + +lib_a_SOURCES = $(src) $(fsrc) + +chobj = scbrt.def scopysign.def sexpm1.def silogb.def \ + sinfinity.def slog1p.def smatherr.def smodf.def \ + snan.def snextafter.def sscalbn.def + +SUFFIXES = .def + +CHEW = ../../doc/makedoc -f $(srcdir)/../../doc/doc.str + +TARGETDOC = ../tmp.texi + +CLEANFILES = $(chobj) *.ref +mkinstalldirs = $(SHELL) $(top_srcdir)/../../mkinstalldirs +CONFIG_CLEAN_FILES = +LIBRARIES = $(noinst_LIBRARIES) + + +DEFS = @DEFS@ -I. -I$(srcdir) +CPPFLAGS = @CPPFLAGS@ +LDFLAGS = @LDFLAGS@ +LIBS = @LIBS@ +lib_a_LIBADD = +lib_a_OBJECTS = s_finite.o s_copysign.o s_modf.o s_scalbn.o s_cbrt.o \ +s_expm1.o s_ilogb.o s_infinity.o s_log1p.o s_nan.o s_nextafter.o \ +s_rint.o s_logb.o s_matherr.o s_lib_ver.o sf_finite.o sf_copysign.o \ +sf_modf.o sf_scalbn.o sf_cbrt.o sf_expm1.o sf_ilogb.o sf_infinity.o \ +sf_log1p.o sf_nan.o sf_nextafter.o sf_rint.o sf_logb.o +CFLAGS = @CFLAGS@ +COMPILE = $(CC) $(DEFS) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) +LINK = $(CC) $(AM_CFLAGS) $(CFLAGS) $(LDFLAGS) -o $@ +DIST_COMMON = Makefile.am Makefile.in + + +DISTFILES = $(DIST_COMMON) $(SOURCES) $(HEADERS) $(TEXINFOS) $(EXTRA_DIST) + +TAR = tar +GZIP = --best +SOURCES = $(lib_a_SOURCES) +OBJECTS = $(lib_a_OBJECTS) + +all: Makefile $(LIBRARIES) + +.SUFFIXES: +.SUFFIXES: .S .c .def .o .s +$(srcdir)/Makefile.in: @MAINT@ Makefile.am $(top_srcdir)/configure.in $(ACLOCAL_M4) + cd $(top_srcdir) && $(AUTOMAKE) --cygnus common/Makefile + +Makefile: $(srcdir)/Makefile.in $(top_builddir)/config.status + cd $(top_builddir) \ + && CONFIG_FILES=$(subdir)/$@ CONFIG_HEADERS= $(SHELL) ./config.status + + +mostlyclean-noinstLIBRARIES: + +clean-noinstLIBRARIES: + -test -z "$(noinst_LIBRARIES)" || rm -f $(noinst_LIBRARIES) + +distclean-noinstLIBRARIES: + +maintainer-clean-noinstLIBRARIES: + +.c.o: + $(COMPILE) -c $< + +.s.o: + $(COMPILE) -c $< + +.S.o: + $(COMPILE) -c $< + +mostlyclean-compile: + -rm -f *.o core *.core + +clean-compile: + +distclean-compile: + -rm -f *.tab.c + +maintainer-clean-compile: + +lib.a: $(lib_a_OBJECTS) $(lib_a_DEPENDENCIES) + -rm -f lib.a + $(AR) cru lib.a $(lib_a_OBJECTS) $(lib_a_LIBADD) + $(RANLIB) lib.a + +tags: TAGS + +ID: $(HEADERS) $(SOURCES) $(LISP) + here=`pwd` && cd $(srcdir) \ + && mkid -f$$here/ID $(SOURCES) $(HEADERS) $(LISP) + +TAGS: $(HEADERS) $(SOURCES) $(TAGS_DEPENDENCIES) $(LISP) + tags=; \ + here=`pwd`; \ + list='$(SOURCES) $(HEADERS)'; \ + unique=`for i in $$list; do echo $$i; done | \ + awk ' { files[$$0] = 1; } \ + END { for (i in files) print i; }'`; \ + test -z "$(ETAGS_ARGS)$$unique$(LISP)$$tags" \ + || (cd $(srcdir) && etags $(ETAGS_ARGS) $$tags $$unique $(LISP) -o $$here/TAGS) + +mostlyclean-tags: + +clean-tags: + +distclean-tags: + -rm -f TAGS ID + +maintainer-clean-tags: + +distdir = $(top_builddir)/$(PACKAGE)-$(VERSION)/$(subdir) + +subdir = common + +distdir: $(DISTFILES) + @for file in $(DISTFILES); do \ + if test -f $$file; then d=.; else d=$(srcdir); fi; \ + test -f $(distdir)/$$file \ + || ln $$d/$$file $(distdir)/$$file 2> /dev/null \ + || cp -p $$d/$$file $(distdir)/$$file; \ + done +info: +dvi: +check: +installcheck: +install-info: +install-exec: + @$(NORMAL_INSTALL) + +install-data: + @$(NORMAL_INSTALL) + +install: install-exec install-data all + @: + +uninstall: + +install-strip: + $(MAKE) $(AM_MAKEFLAGS) INSTALL_PROGRAM='$(INSTALL_PROGRAM) -s' INSTALL_SCRIPT='$(INSTALL_PROGRAM)' install +installdirs: + + +mostlyclean-generic: + +clean-generic: + -test -z "$(CLEANFILES)" || rm -f $(CLEANFILES) + +distclean-generic: + -rm -f Makefile $(CONFIG_CLEAN_FILES) + -rm -f config.cache config.log stamp-h stamp-h[0-9]* + +maintainer-clean-generic: +mostlyclean: mostlyclean-noinstLIBRARIES mostlyclean-compile \ + mostlyclean-tags mostlyclean-generic + +clean: clean-noinstLIBRARIES clean-compile clean-tags clean-generic \ + mostlyclean + +distclean: distclean-noinstLIBRARIES distclean-compile distclean-tags \ + distclean-generic clean + -rm -f config.status + +maintainer-clean: maintainer-clean-noinstLIBRARIES \ + maintainer-clean-compile maintainer-clean-tags \ + maintainer-clean-generic distclean + @echo "This command is intended for maintainers to use;" + @echo "it deletes files that may require special tools to rebuild." + +.PHONY: mostlyclean-noinstLIBRARIES distclean-noinstLIBRARIES \ +clean-noinstLIBRARIES maintainer-clean-noinstLIBRARIES \ +mostlyclean-compile distclean-compile clean-compile \ +maintainer-clean-compile tags mostlyclean-tags distclean-tags \ +clean-tags maintainer-clean-tags distdir info dvi installcheck \ +install-info install-exec install-data install uninstall all \ +installdirs mostlyclean-generic distclean-generic clean-generic \ +maintainer-clean-generic clean mostlyclean distclean maintainer-clean + + +.c.def: + $(CHEW) < $< > $*.def 2> $*.ref + touch stmp-def + +doc: $(chobj) + cat $(srcdir)/common.tex >> $(TARGETDOC) + +# Texinfo does not appear to support underscores in file names, so we +# name the .def files without underscores. + +smodf.def: s_modf.c + $(CHEW) < $(srcdir)/s_modf.c >$@ 2>/dev/null + touch stmp-def + +scopysign.def: s_copysign.c + $(CHEW) < $(srcdir)/s_copysign.c >$@ 2>/dev/null + touch stmp-def + +sscalbn.def: s_scalbn.c + $(CHEW) < $(srcdir)/s_scalbn.c >$@ 2>/dev/null + touch stmp-def + +scbrt.def: s_cbrt.c + $(CHEW) < $(srcdir)/s_cbrt.c >$@ 2>/dev/null + touch stmp-def + +serf.def: s_erf.c + $(CHEW) < $(srcdir)/s_serf.c >$@ 2>/dev/null + touch stmp-def + +sexpn.def: s_expm.c + $(CHEW) < $(srcdir)/s_expn.c >$@ 2>/dev/null + touch stmp-def + +sexpm1.def: s_expm1.c + $(CHEW) < $(srcdir)/s_expm1.c >$@ 2>/dev/null + touch stmp-def + +silogb.def: s_ilogb.c + $(CHEW) < $(srcdir)/s_ilogb.c >$@ 2>/dev/null + touch stmp-def + +sinfinity.def: s_infinity.c + $(CHEW) < $(srcdir)/s_infinity.c >$@ 2>/dev/null + touch stmp-def + +slog1p.def: s_log1p.c + $(CHEW) < $(srcdir)/s_log1p.c >$@ 2>/dev/null + touch stmp-def + +smatherr.def: s_matherr.c + $(CHEW) < $(srcdir)/s_matherr.c >$@ 2>/dev/null + touch stmp-def + +snan.def: s_nan.c + $(CHEW) < $(srcdir)/s_nan.c >$@ 2>/dev/null + touch stmp-def + +snextafter.def: s_nextafter.c + $(CHEW) < $(srcdir)/s_nextafter.c >$@ 2>/dev/null + touch stmp-def + +# A partial dependency list. + +$(lib_a_OBJECTS): $(srcdir)/../../libc/include/math.h fdlibm.h + +# Tell versions [3.59,3.63) of GNU make to not export all variables. +# Otherwise a system limit (for SysV at least) may be exceeded. +.NOEXPORT: diff --git a/newlib/libm/common/common.tex b/newlib/libm/common/common.tex new file mode 100644 index 000000000..0e8a9863f --- /dev/null +++ b/newlib/libm/common/common.tex @@ -0,0 +1,12 @@ +@page +@include common/scbrt.def +@include common/scopysign.def +@include common/sexpm1.def +@include common/silogb.def +@include common/sinfinity.def +@include common/slog1p.def +@include common/smatherr.def +@include common/smodf.def +@include common/snan.def +@include common/snextafter.def +@include common/sscalbn.def diff --git a/newlib/libm/common/fdlibm.h b/newlib/libm/common/fdlibm.h new file mode 100644 index 000000000..752d24246 --- /dev/null +++ b/newlib/libm/common/fdlibm.h @@ -0,0 +1,261 @@ + +/* @(#)fdlibm.h 5.1 93/09/24 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* CYGNUS LOCAL: Include files. */ +#include <math.h> +#include <machine/ieeefp.h> + +/* CYGNUS LOCAL: Default to XOPEN_MODE. */ +#define _XOPEN_MODE + +#ifdef __STDC__ +#define __P(p) p +#else +#define __P(p) () +#endif + +#define HUGE ((float)3.40282346638528860e+38) + +/* + * set X_TLOSS = pi*2**52, which is possibly defined in <values.h> + * (one may replace the following line by "#include <values.h>") + */ + +#define X_TLOSS 1.41484755040568800000e+16 + +/* Functions that are not documented, and are not in <math.h>. */ + +extern double logb __P((double)); +#ifdef _SCALB_INT +extern double scalb __P((double, int)); +#else +extern double scalb __P((double, double)); +#endif +extern double significand __P((double)); + +/* ieee style elementary functions */ +extern double __ieee754_sqrt __P((double)); +extern double __ieee754_acos __P((double)); +extern double __ieee754_acosh __P((double)); +extern double __ieee754_log __P((double)); +extern double __ieee754_atanh __P((double)); +extern double __ieee754_asin __P((double)); +extern double __ieee754_atan2 __P((double,double)); +extern double __ieee754_exp __P((double)); +extern double __ieee754_cosh __P((double)); +extern double __ieee754_fmod __P((double,double)); +extern double __ieee754_pow __P((double,double)); +extern double __ieee754_lgamma_r __P((double,int *)); +extern double __ieee754_gamma_r __P((double,int *)); +extern double __ieee754_log10 __P((double)); +extern double __ieee754_sinh __P((double)); +extern double __ieee754_hypot __P((double,double)); +extern double __ieee754_j0 __P((double)); +extern double __ieee754_j1 __P((double)); +extern double __ieee754_y0 __P((double)); +extern double __ieee754_y1 __P((double)); +extern double __ieee754_jn __P((int,double)); +extern double __ieee754_yn __P((int,double)); +extern double __ieee754_remainder __P((double,double)); +extern __int32_t __ieee754_rem_pio2 __P((double,double*)); +#ifdef _SCALB_INT +extern double __ieee754_scalb __P((double,int)); +#else +extern double __ieee754_scalb __P((double,double)); +#endif + +/* fdlibm kernel function */ +extern double __kernel_standard __P((double,double,int)); +extern double __kernel_sin __P((double,double,int)); +extern double __kernel_cos __P((double,double)); +extern double __kernel_tan __P((double,double,int)); +extern int __kernel_rem_pio2 __P((double*,double*,int,int,int,const __int32_t*)); + +/* Undocumented float functions. */ +extern float logbf __P((float)); +#ifdef _SCALB_INT +extern float scalbf __P((float, int)); +#else +extern float scalbf __P((float, float)); +#endif +extern float significandf __P((float)); + +/* ieee style elementary float functions */ +extern float __ieee754_sqrtf __P((float)); +extern float __ieee754_acosf __P((float)); +extern float __ieee754_acoshf __P((float)); +extern float __ieee754_logf __P((float)); +extern float __ieee754_atanhf __P((float)); +extern float __ieee754_asinf __P((float)); +extern float __ieee754_atan2f __P((float,float)); +extern float __ieee754_expf __P((float)); +extern float __ieee754_coshf __P((float)); +extern float __ieee754_fmodf __P((float,float)); +extern float __ieee754_powf __P((float,float)); +extern float __ieee754_lgammaf_r __P((float,int *)); +extern float __ieee754_gammaf_r __P((float,int *)); +extern float __ieee754_log10f __P((float)); +extern float __ieee754_sinhf __P((float)); +extern float __ieee754_hypotf __P((float,float)); +extern float __ieee754_j0f __P((float)); +extern float __ieee754_j1f __P((float)); +extern float __ieee754_y0f __P((float)); +extern float __ieee754_y1f __P((float)); +extern float __ieee754_jnf __P((int,float)); +extern float __ieee754_ynf __P((int,float)); +extern float __ieee754_remainderf __P((float,float)); +extern __int32_t __ieee754_rem_pio2f __P((float,float*)); +#ifdef _SCALB_INT +extern float __ieee754_scalbf __P((float,int)); +#else +extern float __ieee754_scalbf __P((float,float)); +#endif + +/* float versions of fdlibm kernel functions */ +extern float __kernel_sinf __P((float,float,int)); +extern float __kernel_cosf __P((float,float)); +extern float __kernel_tanf __P((float,float,int)); +extern int __kernel_rem_pio2f __P((float*,float*,int,int,int,const __int32_t*)); + +/* The original code used statements like + n0 = ((*(int*)&one)>>29)^1; * index of high word * + ix0 = *(n0+(int*)&x); * high word of x * + ix1 = *((1-n0)+(int*)&x); * low word of x * + to dig two 32 bit words out of the 64 bit IEEE floating point + value. That is non-ANSI, and, moreover, the gcc instruction + scheduler gets it wrong. We instead use the following macros. + Unlike the original code, we determine the endianness at compile + time, not at run time; I don't see much benefit to selecting + endianness at run time. */ + +#ifndef __IEEE_BIG_ENDIAN +#ifndef __IEEE_LITTLE_ENDIAN + #error Must define endianness +#endif +#endif + +/* A union which permits us to convert between a double and two 32 bit + ints. */ + +#ifdef __IEEE_BIG_ENDIAN + +typedef union +{ + double value; + struct + { + __uint32_t msw; + __uint32_t lsw; + } parts; +} ieee_double_shape_type; + +#endif + +#ifdef __IEEE_LITTLE_ENDIAN + +typedef union +{ + double value; + struct + { + __uint32_t lsw; + __uint32_t msw; + } parts; +} ieee_double_shape_type; + +#endif + +/* Get two 32 bit ints from a double. */ + +#define EXTRACT_WORDS(ix0,ix1,d) \ +do { \ + ieee_double_shape_type ew_u; \ + ew_u.value = (d); \ + (ix0) = ew_u.parts.msw; \ + (ix1) = ew_u.parts.lsw; \ +} while (0) + +/* Get the more significant 32 bit int from a double. */ + +#define GET_HIGH_WORD(i,d) \ +do { \ + ieee_double_shape_type gh_u; \ + gh_u.value = (d); \ + (i) = gh_u.parts.msw; \ +} while (0) + +/* Get the less significant 32 bit int from a double. */ + +#define GET_LOW_WORD(i,d) \ +do { \ + ieee_double_shape_type gl_u; \ + gl_u.value = (d); \ + (i) = gl_u.parts.lsw; \ +} while (0) + +/* Set a double from two 32 bit ints. */ + +#define INSERT_WORDS(d,ix0,ix1) \ +do { \ + ieee_double_shape_type iw_u; \ + iw_u.parts.msw = (ix0); \ + iw_u.parts.lsw = (ix1); \ + (d) = iw_u.value; \ +} while (0) + +/* Set the more significant 32 bits of a double from an int. */ + +#define SET_HIGH_WORD(d,v) \ +do { \ + ieee_double_shape_type sh_u; \ + sh_u.value = (d); \ + sh_u.parts.msw = (v); \ + (d) = sh_u.value; \ +} while (0) + +/* Set the less significant 32 bits of a double from an int. */ + +#define SET_LOW_WORD(d,v) \ +do { \ + ieee_double_shape_type sl_u; \ + sl_u.value = (d); \ + sl_u.parts.lsw = (v); \ + (d) = sl_u.value; \ +} while (0) + +/* A union which permits us to convert between a float and a 32 bit + int. */ + +typedef union +{ + float value; + __uint32_t word; +} ieee_float_shape_type; + +/* Get a 32 bit int from a float. */ + +#define GET_FLOAT_WORD(i,d) \ +do { \ + ieee_float_shape_type gf_u; \ + gf_u.value = (d); \ + (i) = gf_u.word; \ +} while (0) + +/* Set a float from a 32 bit int. */ + +#define SET_FLOAT_WORD(d,i) \ +do { \ + ieee_float_shape_type sf_u; \ + sf_u.word = (i); \ + (d) = sf_u.value; \ +} while (0) diff --git a/newlib/libm/common/s_cbrt.c b/newlib/libm/common/s_cbrt.c new file mode 100644 index 000000000..95185d0fa --- /dev/null +++ b/newlib/libm/common/s_cbrt.c @@ -0,0 +1,123 @@ + +/* @(#)s_cbrt.c 5.1 93/09/24 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + * + */ + +/* +FUNCTION + <<cbrt>>, <<cbrtf>>---cube root + +INDEX + cbrt +INDEX + cbrtf + +ANSI_SYNOPSIS + #include <math.h> + double cbrt(double <[x]>); + float cbrtf(float <[x]>); + +TRAD_SYNOPSIS + #include <math.h> + double cbrt(<[x]>); + float cbrtf(<[x]>); + +DESCRIPTION + <<cbrt>> computes the cube root of the argument. + +RETURNS + The cube root is returned. + +PORTABILITY + <<cbrt>> is in System V release 4. <<cbrtf>> is an extension. +*/ + +#include "fdlibm.h" + +#ifndef _DOUBLE_IS_32BITS + +/* cbrt(x) + * Return cube root of x + */ +#ifdef __STDC__ +static const __uint32_t +#else +static __uint32_t +#endif + B1 = 715094163, /* B1 = (682-0.03306235651)*2**20 */ + B2 = 696219795; /* B2 = (664-0.03306235651)*2**20 */ + +#ifdef __STDC__ +static const double +#else +static double +#endif +C = 5.42857142857142815906e-01, /* 19/35 = 0x3FE15F15, 0xF15F15F1 */ +D = -7.05306122448979611050e-01, /* -864/1225 = 0xBFE691DE, 0x2532C834 */ +E = 1.41428571428571436819e+00, /* 99/70 = 0x3FF6A0EA, 0x0EA0EA0F */ +F = 1.60714285714285720630e+00, /* 45/28 = 0x3FF9B6DB, 0x6DB6DB6E */ +G = 3.57142857142857150787e-01; /* 5/14 = 0x3FD6DB6D, 0xB6DB6DB7 */ + +#ifdef __STDC__ + double cbrt(double x) +#else + double cbrt(x) + double x; +#endif +{ + __int32_t hx; + double r,s,t=0.0,w; + __uint32_t sign; + __uint32_t high,low; + + GET_HIGH_WORD(hx,x); + sign=hx&0x80000000; /* sign= sign(x) */ + hx ^=sign; + if(hx>=0x7ff00000) return(x+x); /* cbrt(NaN,INF) is itself */ + GET_LOW_WORD(low,x); + if((hx|low)==0) + return(x); /* cbrt(0) is itself */ + + SET_HIGH_WORD(x,hx); /* x <- |x| */ + /* rough cbrt to 5 bits */ + if(hx<0x00100000) /* subnormal number */ + {SET_HIGH_WORD(t,0x43500000); /* set t= 2**54 */ + t*=x; GET_HIGH_WORD(high,t); SET_HIGH_WORD(t,high/3+B2); + } + else + SET_HIGH_WORD(t,hx/3+B1); + + + /* new cbrt to 23 bits, may be implemented in single precision */ + r=t*t/x; + s=C+r*t; + t*=G+F/(s+E+D/s); + + /* chopped to 20 bits and make it larger than cbrt(x) */ + GET_HIGH_WORD(high,t); + INSERT_WORDS(t,high+0x00000001,0); + + + /* one step newton iteration to 53 bits with error less than 0.667 ulps */ + s=t*t; /* t*t is exact */ + r=x/s; + w=t+t; + r=(r-t)/(w+r); /* r-s is exact */ + t=t+t*r; + + /* retore the sign bit */ + GET_HIGH_WORD(high,t); + SET_HIGH_WORD(t,high|sign); + return(t); +} + +#endif /* _DOUBLE_IS_32BITS */ diff --git a/newlib/libm/common/s_copysign.c b/newlib/libm/common/s_copysign.c new file mode 100644 index 000000000..bfc546db5 --- /dev/null +++ b/newlib/libm/common/s_copysign.c @@ -0,0 +1,82 @@ + +/* @(#)s_copysign.c 5.1 93/09/24 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* +FUNCTION +<<copysign>>, <<copysignf>>---sign of <[y]>, magnitude of <[x]> + +INDEX + copysign +INDEX + copysignf + +ANSI_SYNOPSIS + #include <math.h> + double copysign (double <[x]>, double <[y]>); + float copysignf (float <[x]>, float <[y]>); + +TRAD_SYNOPSIS + #include <math.h> + double copysign (<[x]>, <[y]>) + double <[x]>; + double <[y]>; + + float copysignf (<[x]>, <[y]>) + float <[x]>; + float <[y]>; + +DESCRIPTION +<<copysign>> constructs a number with the magnitude (absolute value) +of its first argument, <[x]>, and the sign of its second argument, +<[y]>. + +<<copysignf>> does the same thing; the two functions differ only in +the type of their arguments and result. + +RETURNS +<<copysign>> returns a <<double>> with the magnitude of +<[x]> and the sign of <[y]>. +<<copysignf>> returns a <<float>> with the magnitude of +<[x]> and the sign of <[y]>. + +PORTABILITY +<<copysign>> is not required by either ANSI C or the System V Interface +Definition (Issue 2). + +*/ + +/* + * copysign(double x, double y) + * copysign(x,y) returns a value with the magnitude of x and + * with the sign bit of y. + */ + +#include "fdlibm.h" + +#ifndef _DOUBLE_IS_32BITS + +#ifdef __STDC__ + double copysign(double x, double y) +#else + double copysign(x,y) + double x,y; +#endif +{ + __uint32_t hx,hy; + GET_HIGH_WORD(hx,x); + GET_HIGH_WORD(hy,y); + SET_HIGH_WORD(x,(hx&0x7fffffff)|(hy&0x80000000)); + return x; +} + +#endif /* _DOUBLE_IS_32BITS */ diff --git a/newlib/libm/common/s_expm1.c b/newlib/libm/common/s_expm1.c new file mode 100644 index 000000000..bc0e2f23e --- /dev/null +++ b/newlib/libm/common/s_expm1.c @@ -0,0 +1,271 @@ + +/* @(#)s_expm1.c 5.1 93/09/24 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* +FUNCTION + <<expm1>>, <<expm1f>>---exponential minus 1 +INDEX + expm1 +INDEX + expm1f + +ANSI_SYNOPSIS + #include <math.h> + double expm1(double <[x]>); + float expm1f(float <[x]>); + +TRAD_SYNOPSIS + #include <math.h> + double expm1(<[x]>); + double <[x]>; + + float expm1f(<[x]>); + float <[x]>; + +DESCRIPTION + <<expm1>> and <<expm1f>> calculate the exponential of <[x]> + and subtract 1, that is, + @ifinfo + e raised to the power <[x]> minus 1 (where e + @end ifinfo + @tex + $e^x - 1$ (where $e$ + @end tex + is the base of the natural system of logarithms, approximately + 2.71828). The result is accurate even for small values of + <[x]>, where using <<exp(<[x]>)-1>> would lose many + significant digits. + +RETURNS + e raised to the power <[x]>, minus 1. + +PORTABILITY + Neither <<expm1>> nor <<expm1f>> is required by ANSI C or by + the System V Interface Definition (Issue 2). +*/ + +/* expm1(x) + * Returns exp(x)-1, the exponential of x minus 1. + * + * Method + * 1. Argument reduction: + * Given x, find r and integer k such that + * + * x = k*ln2 + r, |r| <= 0.5*ln2 ~ 0.34658 + * + * Here a correction term c will be computed to compensate + * the error in r when rounded to a floating-point number. + * + * 2. Approximating expm1(r) by a special rational function on + * the interval [0,0.34658]: + * Since + * r*(exp(r)+1)/(exp(r)-1) = 2+ r^2/6 - r^4/360 + ... + * we define R1(r*r) by + * r*(exp(r)+1)/(exp(r)-1) = 2+ r^2/6 * R1(r*r) + * That is, + * R1(r**2) = 6/r *((exp(r)+1)/(exp(r)-1) - 2/r) + * = 6/r * ( 1 + 2.0*(1/(exp(r)-1) - 1/r)) + * = 1 - r^2/60 + r^4/2520 - r^6/100800 + ... + * We use a special Reme algorithm on [0,0.347] to generate + * a polynomial of degree 5 in r*r to approximate R1. The + * maximum error of this polynomial approximation is bounded + * by 2**-61. In other words, + * R1(z) ~ 1.0 + Q1*z + Q2*z**2 + Q3*z**3 + Q4*z**4 + Q5*z**5 + * where Q1 = -1.6666666666666567384E-2, + * Q2 = 3.9682539681370365873E-4, + * Q3 = -9.9206344733435987357E-6, + * Q4 = 2.5051361420808517002E-7, + * Q5 = -6.2843505682382617102E-9; + * (where z=r*r, and the values of Q1 to Q5 are listed below) + * with error bounded by + * | 5 | -61 + * | 1.0+Q1*z+...+Q5*z - R1(z) | <= 2 + * | | + * + * expm1(r) = exp(r)-1 is then computed by the following + * specific way which minimize the accumulation rounding error: + * 2 3 + * r r [ 3 - (R1 + R1*r/2) ] + * expm1(r) = r + --- + --- * [--------------------] + * 2 2 [ 6 - r*(3 - R1*r/2) ] + * + * To compensate the error in the argument reduction, we use + * expm1(r+c) = expm1(r) + c + expm1(r)*c + * ~ expm1(r) + c + r*c + * Thus c+r*c will be added in as the correction terms for + * expm1(r+c). Now rearrange the term to avoid optimization + * screw up: + * ( 2 2 ) + * ({ ( r [ R1 - (3 - R1*r/2) ] ) } r ) + * expm1(r+c)~r - ({r*(--- * [--------------------]-c)-c} - --- ) + * ({ ( 2 [ 6 - r*(3 - R1*r/2) ] ) } 2 ) + * ( ) + * + * = r - E + * 3. Scale back to obtain expm1(x): + * From step 1, we have + * expm1(x) = either 2^k*[expm1(r)+1] - 1 + * = or 2^k*[expm1(r) + (1-2^-k)] + * 4. Implementation notes: + * (A). To save one multiplication, we scale the coefficient Qi + * to Qi*2^i, and replace z by (x^2)/2. + * (B). To achieve maximum accuracy, we compute expm1(x) by + * (i) if x < -56*ln2, return -1.0, (raise inexact if x!=inf) + * (ii) if k=0, return r-E + * (iii) if k=-1, return 0.5*(r-E)-0.5 + * (iv) if k=1 if r < -0.25, return 2*((r+0.5)- E) + * else return 1.0+2.0*(r-E); + * (v) if (k<-2||k>56) return 2^k(1-(E-r)) - 1 (or exp(x)-1) + * (vi) if k <= 20, return 2^k((1-2^-k)-(E-r)), else + * (vii) return 2^k(1-((E+2^-k)-r)) + * + * Special cases: + * expm1(INF) is INF, expm1(NaN) is NaN; + * expm1(-INF) is -1, and + * for finite argument, only expm1(0)=0 is exact. + * + * Accuracy: + * according to an error analysis, the error is always less than + * 1 ulp (unit in the last place). + * + * Misc. info. + * For IEEE double + * if x > 7.09782712893383973096e+02 then expm1(x) overflow + * + * Constants: + * The hexadecimal values are the intended ones for the following + * constants. The decimal values may be used, provided that the + * compiler will convert from decimal to binary accurately enough + * to produce the hexadecimal values shown. + */ + +#include "fdlibm.h" + +#ifndef _DOUBLE_IS_32BITS + +#ifdef __STDC__ +static const double +#else +static double +#endif +one = 1.0, +huge = 1.0e+300, +tiny = 1.0e-300, +o_threshold = 7.09782712893383973096e+02,/* 0x40862E42, 0xFEFA39EF */ +ln2_hi = 6.93147180369123816490e-01,/* 0x3fe62e42, 0xfee00000 */ +ln2_lo = 1.90821492927058770002e-10,/* 0x3dea39ef, 0x35793c76 */ +invln2 = 1.44269504088896338700e+00,/* 0x3ff71547, 0x652b82fe */ + /* scaled coefficients related to expm1 */ +Q1 = -3.33333333333331316428e-02, /* BFA11111 111110F4 */ +Q2 = 1.58730158725481460165e-03, /* 3F5A01A0 19FE5585 */ +Q3 = -7.93650757867487942473e-05, /* BF14CE19 9EAADBB7 */ +Q4 = 4.00821782732936239552e-06, /* 3ED0CFCA 86E65239 */ +Q5 = -2.01099218183624371326e-07; /* BE8AFDB7 6E09C32D */ + +#ifdef __STDC__ + double expm1(double x) +#else + double expm1(x) + double x; +#endif +{ + double y,hi,lo,c,t,e,hxs,hfx,r1; + __int32_t k,xsb; + __uint32_t hx; + + GET_HIGH_WORD(hx,x); + xsb = hx&0x80000000; /* sign bit of x */ + if(xsb==0) y=x; else y= -x; /* y = |x| */ + hx &= 0x7fffffff; /* high word of |x| */ + + /* filter out huge and non-finite argument */ + if(hx >= 0x4043687A) { /* if |x|>=56*ln2 */ + if(hx >= 0x40862E42) { /* if |x|>=709.78... */ + if(hx>=0x7ff00000) { + __uint32_t low; + GET_LOW_WORD(low,x); + if(((hx&0xfffff)|low)!=0) + return x+x; /* NaN */ + else return (xsb==0)? x:-1.0;/* exp(+-inf)={inf,-1} */ + } + if(x > o_threshold) return huge*huge; /* overflow */ + } + if(xsb!=0) { /* x < -56*ln2, return -1.0 with inexact */ + if(x+tiny<0.0) /* raise inexact */ + return tiny-one; /* return -1 */ + } + } + + /* argument reduction */ + if(hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */ + if(hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */ + if(xsb==0) + {hi = x - ln2_hi; lo = ln2_lo; k = 1;} + else + {hi = x + ln2_hi; lo = -ln2_lo; k = -1;} + } else { + k = invln2*x+((xsb==0)?0.5:-0.5); + t = k; + hi = x - t*ln2_hi; /* t*ln2_hi is exact here */ + lo = t*ln2_lo; + } + x = hi - lo; + c = (hi-x)-lo; + } + else if(hx < 0x3c900000) { /* when |x|<2**-54, return x */ + t = huge+x; /* return x with inexact flags when x!=0 */ + return x - (t-(huge+x)); + } + else k = 0; + + /* x is now in primary range */ + hfx = 0.5*x; + hxs = x*hfx; + r1 = one+hxs*(Q1+hxs*(Q2+hxs*(Q3+hxs*(Q4+hxs*Q5)))); + t = 3.0-r1*hfx; + e = hxs*((r1-t)/(6.0 - x*t)); + if(k==0) return x - (x*e-hxs); /* c is 0 */ + else { + e = (x*(e-c)-c); + e -= hxs; + if(k== -1) return 0.5*(x-e)-0.5; + if(k==1) + if(x < -0.25) return -2.0*(e-(x+0.5)); + else return one+2.0*(x-e); + if (k <= -2 || k>56) { /* suffice to return exp(x)-1 */ + __uint32_t high; + y = one-(e-x); + GET_HIGH_WORD(high,y); + SET_HIGH_WORD(y,high+(k<<20)); /* add k to y's exponent */ + return y-one; + } + t = one; + if(k<20) { + __uint32_t high; + SET_HIGH_WORD(t,0x3ff00000 - (0x200000>>k)); /* t=1-2^-k */ + y = t-(e-x); + GET_HIGH_WORD(high,y); + SET_HIGH_WORD(y,high+(k<<20)); /* add k to y's exponent */ + } else { + __uint32_t high; + SET_HIGH_WORD(t,((0x3ff-k)<<20)); /* 2^-k */ + y = x-(e+t); + y += one; + GET_HIGH_WORD(high,y); + SET_HIGH_WORD(y,high+(k<<20)); /* add k to y's exponent */ + } + } + return y; +} + +#endif /* _DOUBLE_IS_32BITS */ diff --git a/newlib/libm/common/s_finite.c b/newlib/libm/common/s_finite.c new file mode 100644 index 000000000..17f4e849b --- /dev/null +++ b/newlib/libm/common/s_finite.c @@ -0,0 +1,35 @@ + +/* @(#)s_finite.c 5.1 93/09/24 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* + * finite(x) returns 1 is x is finite, else 0; + * no branching! + */ + +#include "fdlibm.h" + +#ifndef _DOUBLE_IS_32BITS + +#ifdef __STDC__ + int finite(double x) +#else + int finite(x) + double x; +#endif +{ + __int32_t hx; + GET_HIGH_WORD(hx,x); + return (int)((__uint32_t)((hx&0x7fffffff)-0x7ff00000)>>31); +} + +#endif /* _DOUBLE_IS_32BITS */ diff --git a/newlib/libm/common/s_ilogb.c b/newlib/libm/common/s_ilogb.c new file mode 100644 index 000000000..4e3e69f12 --- /dev/null +++ b/newlib/libm/common/s_ilogb.c @@ -0,0 +1,92 @@ + +/* @(#)s_ilogb.c 5.1 93/09/24 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* +FUNCTION + <<ilogb>>, <<ilogbf>>---get exponent of floating point number +INDEX + ilogb +INDEX + ilogbf + +ANSI_SYNOPSIS + #include <math.h> + int ilogb(double <[val]>); + int ilogbf(float <[val]>); + +TRAD_SYNOPSIS + #include <math.h> + int ilogb(<[val]>) + double <[val]>; + + int ilogbf(<[val]>) + float <[val]>; + + +DESCRIPTION + + All non zero, normal numbers can be described as <[m]> * + 2**<[p]>. <<ilogb>> and <<ilogbf>> examine the argument + <[val]>, and return <[p]>. The functions <<frexp>> and + <<frexpf>> are similar to <<ilogb>> and <<ilogbf>>, but also + return <[m]>. + +RETURNS + +<<ilogb>> and <<ilogbf>> return the power of two used to form the +floating point argument. If <[val]> is <<0>>, they return <<- +INT_MAX>> (<<INT_MAX>> is defined in limits.h). If <[val]> is +infinite, or NaN, they return <<INT_MAX>>. + +PORTABILITY + Neither <<ilogb>> nor <<ilogbf>> is required by ANSI C or by + the System V Interface Definition (Issue 2). */ + +/* ilogb(double x) + * return the binary exponent of non-zero x + * ilogb(0) = 0x80000001 + * ilogb(inf/NaN) = 0x7fffffff (no signal is raised) + */ + +#include "fdlibm.h" +#include <limits.h> + +#ifndef _DOUBLE_IS_32BITS + +#ifdef __STDC__ + int ilogb(double x) +#else + int ilogb(x) + double x; +#endif +{ + __int32_t hx,lx,ix; + + EXTRACT_WORDS(hx,lx,x); + hx &= 0x7fffffff; + if(hx<0x00100000) { + if((hx|lx)==0) + return - INT_MAX; /* ilogb(0) = 0x80000001 */ + else /* subnormal x */ + if(hx==0) { + for (ix = -1043; lx>0; lx<<=1) ix -=1; + } else { + for (ix = -1022,hx<<=11; hx>0; hx<<=1) ix -=1; + } + return ix; + } + else if (hx<0x7ff00000) return (hx>>20)-1023; + else return INT_MAX; +} + +#endif /* _DOUBLE_IS_32BITS */ diff --git a/newlib/libm/common/s_infinity.c b/newlib/libm/common/s_infinity.c new file mode 100644 index 000000000..6508216bd --- /dev/null +++ b/newlib/libm/common/s_infinity.c @@ -0,0 +1,48 @@ +/* + * infinity () returns the representation of infinity. + * Added by Cygnus Support. + */ + +/* +FUNCTION + <<infinity>>, <<infinityf>>---representation of infinity + +INDEX + infinity +INDEX + infinityf + +ANSI_SYNOPSIS + #include <math.h> + double infinity(void); + float infinityf(void); + +TRAD_SYNOPSIS + #include <math.h> + double infinity(); + float infinityf(); + + +DESCRIPTION + <<infinity>> and <<infinityf>> return the special number IEEE + infinity in double and single precision arithmetic + respectivly. + +QUICKREF + infinity - pure + +*/ + +#include "fdlibm.h" + +#ifndef _DOUBLE_IS_32BITS + + double infinity() +{ + double x; + + INSERT_WORDS(x,0x7ff00000,0); + return x; +} + +#endif /* _DOUBLE_IS_32BITS */ diff --git a/newlib/libm/common/s_lib_ver.c b/newlib/libm/common/s_lib_ver.c new file mode 100644 index 000000000..d42f62404 --- /dev/null +++ b/newlib/libm/common/s_lib_ver.c @@ -0,0 +1,35 @@ + +/* @(#)s_lib_ver.c 5.1 93/09/24 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* + * MACRO for standards + */ + +#include "fdlibm.h" + +/* + * define and initialize _LIB_VERSION + */ +#ifdef _POSIX_MODE +_CONST _LIB_VERSION_TYPE _LIB_VERSION = _POSIX_; +#else +#ifdef _XOPEN_MODE +_CONST _LIB_VERSION_TYPE _LIB_VERSION = _XOPEN_; +#else +#ifdef _SVID3_MODE +_CONST _LIB_VERSION_TYPE _LIB_VERSION = _SVID_; +#else /* default _IEEE_MODE */ +_CONST _LIB_VERSION_TYPE _LIB_VERSION = _IEEE_; +#endif +#endif +#endif diff --git a/newlib/libm/common/s_log1p.c b/newlib/libm/common/s_log1p.c new file mode 100644 index 000000000..3c3d49733 --- /dev/null +++ b/newlib/libm/common/s_log1p.c @@ -0,0 +1,217 @@ + +/* @(#)s_log1p.c 5.1 93/09/24 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* +FUNCTION +<<log1p>>, <<log1pf>>---log of <<1 + <[x]>>> + +INDEX + log1p +INDEX + log1pf + +ANSI_SYNOPSIS + #include <math.h> + double log1p(double <[x]>); + float log1pf(float <[x]>); + +TRAD_SYNOPSIS + #include <math.h> + double log1p(<[x]>) + double <[x]>; + + float log1pf(<[x]>) + float <[x]>; + +DESCRIPTION +<<log1p>> calculates +@tex +$ln(1+x)$, +@end tex +the natural logarithm of <<1+<[x]>>>. You can use <<log1p>> rather +than `<<log(1+<[x]>)>>' for greater precision when <[x]> is very +small. + +<<log1pf>> calculates the same thing, but accepts and returns +<<float>> values rather than <<double>>. + +RETURNS +<<log1p>> returns a <<double>>, the natural log of <<1+<[x]>>>. +<<log1pf>> returns a <<float>>, the natural log of <<1+<[x]>>>. + +PORTABILITY +Neither <<log1p>> nor <<log1pf>> is required by ANSI C or by the System V +Interface Definition (Issue 2). + +*/ + +/* double log1p(double x) + * + * Method : + * 1. Argument Reduction: find k and f such that + * 1+x = 2^k * (1+f), + * where sqrt(2)/2 < 1+f < sqrt(2) . + * + * Note. If k=0, then f=x is exact. However, if k!=0, then f + * may not be representable exactly. In that case, a correction + * term is need. Let u=1+x rounded. Let c = (1+x)-u, then + * log(1+x) - log(u) ~ c/u. Thus, we proceed to compute log(u), + * and add back the correction term c/u. + * (Note: when x > 2**53, one can simply return log(x)) + * + * 2. Approximation of log1p(f). + * Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s) + * = 2s + 2/3 s**3 + 2/5 s**5 + ....., + * = 2s + s*R + * We use a special Reme algorithm on [0,0.1716] to generate + * a polynomial of degree 14 to approximate R The maximum error + * of this polynomial approximation is bounded by 2**-58.45. In + * other words, + * 2 4 6 8 10 12 14 + * R(z) ~ Lp1*s +Lp2*s +Lp3*s +Lp4*s +Lp5*s +Lp6*s +Lp7*s + * (the values of Lp1 to Lp7 are listed in the program) + * and + * | 2 14 | -58.45 + * | Lp1*s +...+Lp7*s - R(z) | <= 2 + * | | + * Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2. + * In order to guarantee error in log below 1ulp, we compute log + * by + * log1p(f) = f - (hfsq - s*(hfsq+R)). + * + * 3. Finally, log1p(x) = k*ln2 + log1p(f). + * = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo))) + * Here ln2 is split into two floating point number: + * ln2_hi + ln2_lo, + * where n*ln2_hi is always exact for |n| < 2000. + * + * Special cases: + * log1p(x) is NaN with signal if x < -1 (including -INF) ; + * log1p(+INF) is +INF; log1p(-1) is -INF with signal; + * log1p(NaN) is that NaN with no signal. + * + * Accuracy: + * according to an error analysis, the error is always less than + * 1 ulp (unit in the last place). + * + * Constants: + * The hexadecimal values are the intended ones for the following + * constants. The decimal values may be used, provided that the + * compiler will convert from decimal to binary accurately enough + * to produce the hexadecimal values shown. + * + * Note: Assuming log() return accurate answer, the following + * algorithm can be used to compute log1p(x) to within a few ULP: + * + * u = 1+x; + * if(u==1.0) return x ; else + * return log(u)*(x/(u-1.0)); + * + * See HP-15C Advanced Functions Handbook, p.193. + */ + +#include "fdlibm.h" + +#ifndef _DOUBLE_IS_32BITS + +#ifdef __STDC__ +static const double +#else +static double +#endif +ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */ +ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */ +two54 = 1.80143985094819840000e+16, /* 43500000 00000000 */ +Lp1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */ +Lp2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */ +Lp3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */ +Lp4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */ +Lp5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */ +Lp6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */ +Lp7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */ + +#ifdef __STDC__ +static const double zero = 0.0; +#else +static double zero = 0.0; +#endif + +#ifdef __STDC__ + double log1p(double x) +#else + double log1p(x) + double x; +#endif +{ + double hfsq,f,c,s,z,R,u; + __int32_t k,hx,hu,ax; + + GET_HIGH_WORD(hx,x); + ax = hx&0x7fffffff; + + k = 1; + if (hx < 0x3FDA827A) { /* x < 0.41422 */ + if(ax>=0x3ff00000) { /* x <= -1.0 */ + if(x==-1.0) return -two54/zero; /* log1p(-1)=+inf */ + else return (x-x)/(x-x); /* log1p(x<-1)=NaN */ + } + if(ax<0x3e200000) { /* |x| < 2**-29 */ + if(two54+x>zero /* raise inexact */ + &&ax<0x3c900000) /* |x| < 2**-54 */ + return x; + else + return x - x*x*0.5; + } + if(hx>0||hx<=((__int32_t)0xbfd2bec3)) { + k=0;f=x;hu=1;} /* -0.2929<x<0.41422 */ + } + if (hx >= 0x7ff00000) return x+x; + if(k!=0) { + if(hx<0x43400000) { + u = 1.0+x; + GET_HIGH_WORD(hu,u); + k = (hu>>20)-1023; + c = (k>0)? 1.0-(u-x):x-(u-1.0);/* correction term */ + c /= u; + } else { + u = x; + GET_HIGH_WORD(hu,u); + k = (hu>>20)-1023; + c = 0; + } + hu &= 0x000fffff; + if(hu<0x6a09e) { + SET_HIGH_WORD(u,hu|0x3ff00000); /* normalize u */ + } else { + k += 1; + SET_HIGH_WORD(u,hu|0x3fe00000); /* normalize u/2 */ + hu = (0x00100000-hu)>>2; + } + f = u-1.0; + } + hfsq=0.5*f*f; + if(hu==0) { /* |f| < 2**-20 */ + if(f==zero) if(k==0) return zero; + else {c += k*ln2_lo; return k*ln2_hi+c;} + R = hfsq*(1.0-0.66666666666666666*f); + if(k==0) return f-R; else + return k*ln2_hi-((R-(k*ln2_lo+c))-f); + } + s = f/(2.0+f); + z = s*s; + R = z*(Lp1+z*(Lp2+z*(Lp3+z*(Lp4+z*(Lp5+z*(Lp6+z*Lp7)))))); + if(k==0) return f-(hfsq-s*(hfsq+R)); else + return k*ln2_hi-((hfsq-(s*(hfsq+R)+(k*ln2_lo+c)))-f); +} + +#endif /* _DOUBLE_IS_32BITS */ diff --git a/newlib/libm/common/s_logb.c b/newlib/libm/common/s_logb.c new file mode 100644 index 000000000..dbddd1994 --- /dev/null +++ b/newlib/libm/common/s_logb.c @@ -0,0 +1,42 @@ + +/* @(#)s_logb.c 5.1 93/09/24 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* + * double logb(x) + * IEEE 754 logb. Included to pass IEEE test suite. Not recommend. + * Use ilogb instead. + */ + +#include "fdlibm.h" + +#ifndef _DOUBLE_IS_32BITS + +#ifdef __STDC__ + double logb(double x) +#else + double logb(x) + double x; +#endif +{ + __int32_t lx,ix; + EXTRACT_WORDS(ix,lx,x); + ix &= 0x7fffffff; /* high |x| */ + if((ix|lx)==0) return -1.0/fabs(x); + if(ix>=0x7ff00000) return x*x; + if((ix>>=20)==0) /* IEEE 754 logb */ + return -1022.0; + else + return (double) (ix-1023); +} + +#endif /* _DOUBLE_IS_32BITS */ diff --git a/newlib/libm/common/s_matherr.c b/newlib/libm/common/s_matherr.c new file mode 100644 index 000000000..58e242834 --- /dev/null +++ b/newlib/libm/common/s_matherr.c @@ -0,0 +1,123 @@ + +/* @(#)s_matherr.c 5.1 93/09/24 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* + +FUNCTION + <<matherr>>---modifiable math error handler + +INDEX + matherr + +ANSI_SYNOPSIS + #include <math.h> + int matherr(struct exception *<[e]>); + +TRAD_SYNOPSIS + #include <math.h> + int matherr(*<[e]>) + struct exception *<[e]>; + +DESCRIPTION +<<matherr>> is called whenever a math library function generates an error. +You can replace <<matherr>> by your own subroutine to customize +error treatment. The customized <<matherr>> must return 0 if +it fails to resolve the error, and non-zero if the error is resolved. + +When <<matherr>> returns a nonzero value, no error message is printed +and the value of <<errno>> is not modified. You can accomplish either +or both of these things in your own <<matherr>> using the information +passed in the structure <<*<[e]>>>. + +This is the <<exception>> structure (defined in `<<math.h>>'): +. struct exception { +. int type; +. char *name; +. double arg1, arg2, retval; +. int err; +. }; + +The members of the exception structure have the following meanings: +o+ +o type +The type of mathematical error that occured; macros encoding error +types are also defined in `<<math.h>>'. + +o name +a pointer to a null-terminated string holding the +name of the math library function where the error occurred. + +o arg1, arg2 +The arguments which caused the error. + +o retval +The error return value (what the calling function will return). + +o err +If set to be non-zero, this is the new value assigned to <<errno>>. +o- + +The error types defined in `<<math.h>>' represent possible mathematical +errors as follows: + +o+ +o DOMAIN +An argument was not in the domain of the function; e.g. <<log(-1.0)>>. + +o SING +The requested calculation would result in a singularity; e.g. <<pow(0.0,-2.0)>> + +o OVERFLOW +A calculation would produce a result too large to represent; e.g. +<<exp(1000.0)>>. + +o UNDERFLOW +A calculation would produce a result too small to represent; e.g. +<<exp(-1000.0)>>. + +o TLOSS +Total loss of precision. The result would have no significant digits; +e.g. <<sin(10e70)>>. + +o PLOSS +Partial loss of precision. +o- + + +RETURNS +The library definition for <<matherr>> returns <<0>> in all cases. + +You can change the calling function's result from a customized <<matherr>> +by modifying <<e->retval>>, which propagates backs to the caller. + +If <<matherr>> returns <<0>> (indicating that it was not able to resolve +the error) the caller sets <<errno>> to an appropriate value, and prints +an error message. + +PORTABILITY +<<matherr>> is not ANSI C. +*/ + +#include "fdlibm.h" + +#ifdef __STDC__ + int matherr(struct exception *x) +#else + int matherr(x) + struct exception *x; +#endif +{ + int n=0; + if(x->arg1!=x->arg1) return 0; + return n; +} diff --git a/newlib/libm/common/s_modf.c b/newlib/libm/common/s_modf.c new file mode 100644 index 000000000..01151397d --- /dev/null +++ b/newlib/libm/common/s_modf.c @@ -0,0 +1,131 @@ + +/* @(#)s_modf.c 5.1 93/09/24 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* +FUNCTION + <<modf>>, <<modff>>---split fractional and integer parts + +INDEX + modf +INDEX + modff + +ANSI_SYNOPSIS + #include <math.h> + double modf(double <[val]>, double *<[ipart]>); + float modff(float <[val]>, float *<[ipart]>); + +TRAD_SYNOPSIS + #include <math.h> + double modf(<[val]>, <[ipart]>) + double <[val]>; + double *<[ipart]>; + + float modff(<[val]>, <[ipart]>) + float <[val]>; + float *<[ipart]>; + +DESCRIPTION + <<modf>> splits the double <[val]> apart into an integer part + and a fractional part, returning the fractional part and + storing the integer part in <<*<[ipart]>>>. No rounding + whatsoever is done; the sum of the integer and fractional + parts is guaranteed to be exactly equal to <[val]>. That + is, if . <[realpart]> = modf(<[val]>, &<[intpart]>); then + `<<<[realpart]>+<[intpart]>>>' is the same as <[val]>. + <<modff>> is identical, save that it takes and returns + <<float>> rather than <<double>> values. + +RETURNS + The fractional part is returned. Each result has the same + sign as the supplied argument <[val]>. + +PORTABILITY + <<modf>> is ANSI C. <<modff>> is an extension. + +QUICKREF + modf ansi pure + modff - pure + +*/ + +/* + * modf(double x, double *iptr) + * return fraction part of x, and return x's integral part in *iptr. + * Method: + * Bit twiddling. + * + * Exception: + * No exception. + */ + +#include "fdlibm.h" + +#ifndef _DOUBLE_IS_32BITS + +#ifdef __STDC__ +static const double one = 1.0; +#else +static double one = 1.0; +#endif + +#ifdef __STDC__ + double modf(double x, double *iptr) +#else + double modf(x, iptr) + double x,*iptr; +#endif +{ + __int32_t i0,i1,j0; + __uint32_t i; + EXTRACT_WORDS(i0,i1,x); + j0 = ((i0>>20)&0x7ff)-0x3ff; /* exponent of x */ + if(j0<20) { /* integer part in high x */ + if(j0<0) { /* |x|<1 */ + INSERT_WORDS(*iptr,i0&0x80000000,0); /* *iptr = +-0 */ + return x; + } else { + i = (0x000fffff)>>j0; + if(((i0&i)|i1)==0) { /* x is integral */ + __uint32_t high; + *iptr = x; + GET_HIGH_WORD(high,x); + INSERT_WORDS(x,high&0x80000000,0); /* return +-0 */ + return x; + } else { + INSERT_WORDS(*iptr,i0&(~i),0); + return x - *iptr; + } + } + } else if (j0>51) { /* no fraction part */ + __uint32_t high; + *iptr = x*one; + GET_HIGH_WORD(high,x); + INSERT_WORDS(x,high&0x80000000,0); /* return +-0 */ + return x; + } else { /* fraction part in low x */ + i = ((__uint32_t)(0xffffffff))>>(j0-20); + if((i1&i)==0) { /* x is integral */ + __uint32_t high; + *iptr = x; + GET_HIGH_WORD(high,x); + INSERT_WORDS(x,high&0x80000000,0); /* return +-0 */ + return x; + } else { + INSERT_WORDS(*iptr,i0,i1&(~i)); + return x - *iptr; + } + } +} + +#endif /* _DOUBLE_IS_32BITS */ diff --git a/newlib/libm/common/s_nan.c b/newlib/libm/common/s_nan.c new file mode 100644 index 000000000..f06242647 --- /dev/null +++ b/newlib/libm/common/s_nan.c @@ -0,0 +1,47 @@ +/* + * nan () returns a nan. + * Added by Cygnus Support. + */ + +/* +FUNCTION + <<nan>>, <<nanf>>---representation of infinity + +INDEX + nan +INDEX + nanf + +ANSI_SYNOPSIS + #include <math.h> + double nan(void); + float nanf(void); + +TRAD_SYNOPSIS + #include <math.h> + double nan(); + float nanf(); + + +DESCRIPTION + <<nan>> and <<nanf>> return an IEEE NaN (Not a Number) in + double and single precision arithmetic respectivly. + +QUICKREF + nan - pure + +*/ + +#include "fdlibm.h" + +#ifndef _DOUBLE_IS_32BITS + + double nan() +{ + double x; + + INSERT_WORDS(x,0x7ff80000,0); + return x; +} + +#endif /* _DOUBLE_IS_32BITS */ diff --git a/newlib/libm/common/s_nextafter.c b/newlib/libm/common/s_nextafter.c new file mode 100644 index 000000000..82eb8e3ed --- /dev/null +++ b/newlib/libm/common/s_nextafter.c @@ -0,0 +1,121 @@ + +/* @(#)s_nextafter.c 5.1 93/09/24 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* +FUNCTION + <<nextafter>>, <<nextafterf>>---get next number + +INDEX + nextafter +INDEX + nextafterf + +ANSI_SYNOPSIS + #include <math.h> + double nextafter(double <[val]>, double <[dir]>); + float nextafterf(float <[val]>, float <[dir]>); + +TRAD_SYNOPSIS + #include <math.h> + + double nextafter(<[val]>, <[dir]>) + double <[val]>; + double <[exp]>; + + float nextafter(<[val]>, <[dir]>) + float <[val]>; + float <[dir]>; + + +DESCRIPTION +<<nextafter>> returns the double) precision floating point number +closest to <[val]> in the direction toward <[dir]>. <<nextafterf>> +performs the same operation in single precision. For example, +<<nextafter(0.0,1.0)>> returns the smallest positive number which is +representable in double precision. + +RETURNS +Returns the next closest number to <[val]> in the direction toward +<[dir]>. + +PORTABILITY + Neither <<nextafter>> nor <<nextafterf>> is required by ANSI C + or by the System V Interface Definition (Issue 2). +*/ + +/* IEEE functions + * nextafter(x,y) + * return the next machine floating-point number of x in the + * direction toward y. + * Special cases: + */ + +#include "fdlibm.h" + +#ifndef _DOUBLE_IS_32BITS + +#ifdef __STDC__ + double nextafter(double x, double y) +#else + double nextafter(x,y) + double x,y; +#endif +{ + __int32_t hx,hy,ix,iy; + __uint32_t lx,ly; + + EXTRACT_WORDS(hx,lx,x); + EXTRACT_WORDS(hy,ly,y); + ix = hx&0x7fffffff; /* |x| */ + iy = hy&0x7fffffff; /* |y| */ + + if(((ix>=0x7ff00000)&&((ix-0x7ff00000)|lx)!=0) || /* x is nan */ + ((iy>=0x7ff00000)&&((iy-0x7ff00000)|ly)!=0)) /* y is nan */ + return x+y; + if(x==y) return x; /* x=y, return x */ + if((ix|lx)==0) { /* x == 0 */ + INSERT_WORDS(x,hy&0x80000000,1); /* return +-minsubnormal */ + y = x*x; + if(y==x) return y; else return x; /* raise underflow flag */ + } + if(hx>=0) { /* x > 0 */ + if(hx>hy||((hx==hy)&&(lx>ly))) { /* x > y, x -= ulp */ + if(lx==0) hx -= 1; + lx -= 1; + } else { /* x < y, x += ulp */ + lx += 1; + if(lx==0) hx += 1; + } + } else { /* x < 0 */ + if(hy>=0||hx>hy||((hx==hy)&&(lx>ly))){/* x < y, x -= ulp */ + if(lx==0) hx -= 1; + lx -= 1; + } else { /* x > y, x += ulp */ + lx += 1; + if(lx==0) hx += 1; + } + } + hy = hx&0x7ff00000; + if(hy>=0x7ff00000) return x+x; /* overflow */ + if(hy<0x00100000) { /* underflow */ + y = x*x; + if(y!=x) { /* raise underflow flag */ + INSERT_WORDS(y,hx,lx); + return y; + } + } + INSERT_WORDS(x,hx,lx); + return x; +} + +#endif /* _DOUBLE_IS_32BITS */ diff --git a/newlib/libm/common/s_rint.c b/newlib/libm/common/s_rint.c new file mode 100644 index 000000000..ba33e4786 --- /dev/null +++ b/newlib/libm/common/s_rint.c @@ -0,0 +1,90 @@ + +/* @(#)s_rint.c 5.1 93/09/24 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* + * rint(x) + * Return x rounded to integral value according to the prevailing + * rounding mode. + * Method: + * Using floating addition. + * Exception: + * Inexact flag raised if x not equal to rint(x). + */ + +#include "fdlibm.h" + +#ifndef _DOUBLE_IS_32BITS + +#ifdef __STDC__ +static const double +#else +static double +#endif +TWO52[2]={ + 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */ + -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */ +}; + +#ifdef __STDC__ + double rint(double x) +#else + double rint(x) + double x; +#endif +{ + __int32_t i0,j0,sx; + __uint32_t i,i1; + double t; + volatile double w; + EXTRACT_WORDS(i0,i1,x); + sx = (i0>>31)&1; + j0 = ((i0>>20)&0x7ff)-0x3ff; + if(j0<20) { + if(j0<0) { + if(((i0&0x7fffffff)|i1)==0) return x; + i1 |= (i0&0x0fffff); + i0 &= 0xfffe0000; + i0 |= ((i1|-i1)>>12)&0x80000; + SET_HIGH_WORD(x,i0); + w = TWO52[sx]+x; + t = w-TWO52[sx]; + GET_HIGH_WORD(i0,t); + SET_HIGH_WORD(t,(i0&0x7fffffff)|(sx<<31)); + return t; + } else { + i = (0x000fffff)>>j0; + if(((i0&i)|i1)==0) return x; /* x is integral */ + i>>=1; + if(((i0&i)|i1)!=0) { + if(j0==19) i1 = 0x40000000; else + i0 = (i0&(~i))|((0x20000)>>j0); + } + } + } else if (j0>51) { + if(j0==0x400) return x+x; /* inf or NaN */ + else return x; /* x is integral */ + } else { + i = ((__uint32_t)(0xffffffff))>>(j0-20); + if((i1&i)==0) return x; /* x is integral */ + i>>=1; + if((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-20)); + } + INSERT_WORDS(x,i0,i1); + w = TWO52[sx]+x; + return w-TWO52[sx]; +} + +#endif /* _DOUBLE_IS_32BITS */ + + + diff --git a/newlib/libm/common/s_scalbn.c b/newlib/libm/common/s_scalbn.c new file mode 100644 index 000000000..245888fc8 --- /dev/null +++ b/newlib/libm/common/s_scalbn.c @@ -0,0 +1,103 @@ + +/* @(#)s_scalbn.c 5.1 93/09/24 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* +FUNCTION +<<scalbn>>, <<scalbnf>>---scale by integer +INDEX + scalbn +INDEX + scalbnf + +ANSI_SYNOPSIS + #include <math.h> + double scalbn(double <[x]>, int <[y]>); + float scalbnf(float <[x]>, int <[y]>); + +TRAD_SYNOPSIS + #include <math.h> + double scalbn(<[x]>,<[y]>) + double <[x]>; + int <[y]>; + float scalbnf(<[x]>,<[y]>) + float <[x]>; + int <[y]>; + +DESCRIPTION +<<scalbn>> and <<scalbnf>> scale <[x]> by <[n]>, returning <[x]> times +2 to the power <[n]>. The result is computed by manipulating the +exponent, rather than by actually performing an exponentiation or +multiplication. + +RETURNS +<[x]> times 2 to the power <[n]>. + +PORTABILITY +Neither <<scalbn>> nor <<scalbnf>> is required by ANSI C or by the System V +Interface Definition (Issue 2). + +*/ + +/* + * scalbn (double x, int n) + * scalbn(x,n) returns x* 2**n computed by exponent + * manipulation rather than by actually performing an + * exponentiation or a multiplication. + */ + +#include "fdlibm.h" + +#ifndef _DOUBLE_IS_32BITS + +#ifdef __STDC__ +static const double +#else +static double +#endif +two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */ +twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */ +huge = 1.0e+300, +tiny = 1.0e-300; + +#ifdef __STDC__ + double scalbn (double x, int n) +#else + double scalbn (x,n) + double x; int n; +#endif +{ + __int32_t k,hx,lx; + EXTRACT_WORDS(hx,lx,x); + k = (hx&0x7ff00000)>>20; /* extract exponent */ + if (k==0) { /* 0 or subnormal x */ + if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */ + x *= two54; + GET_HIGH_WORD(hx,x); + k = ((hx&0x7ff00000)>>20) - 54; + if (n< -50000) return tiny*x; /*underflow*/ + } + if (k==0x7ff) return x+x; /* NaN or Inf */ + k = k+n; + if (k > 0x7fe) return huge*copysign(huge,x); /* overflow */ + if (k > 0) /* normal result */ + {SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); return x;} + if (k <= -54) + if (n > 50000) /* in case integer overflow in n+k */ + return huge*copysign(huge,x); /*overflow*/ + else return tiny*copysign(tiny,x); /*underflow*/ + k += 54; /* subnormal result */ + SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); + return x*twom54; +} + +#endif /* _DOUBLE_IS_32BITS */ diff --git a/newlib/libm/common/sf_cbrt.c b/newlib/libm/common/sf_cbrt.c new file mode 100644 index 000000000..c053d9548 --- /dev/null +++ b/newlib/libm/common/sf_cbrt.c @@ -0,0 +1,93 @@ +/* sf_cbrt.c -- float version of s_cbrt.c. + * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. + */ + +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + * + */ + +#include "fdlibm.h" + +/* cbrtf(x) + * Return cube root of x + */ +#ifdef __STDC__ +static const __uint32_t +#else +static __uint32_t +#endif + B1 = 709958130, /* B1 = (84+2/3-0.03306235651)*2**23 */ + B2 = 642849266; /* B2 = (76+2/3-0.03306235651)*2**23 */ + +#ifdef __STDC__ +static const float +#else +static float +#endif +C = 5.4285717010e-01, /* 19/35 = 0x3f0af8b0 */ +D = -7.0530611277e-01, /* -864/1225 = 0xbf348ef1 */ +E = 1.4142856598e+00, /* 99/70 = 0x3fb50750 */ +F = 1.6071428061e+00, /* 45/28 = 0x3fcdb6db */ +G = 3.5714286566e-01; /* 5/14 = 0x3eb6db6e */ + +#ifdef __STDC__ + float cbrtf(float x) +#else + float cbrtf(x) + float x; +#endif +{ + __int32_t hx; + float r,s,t; + __uint32_t sign; + __uint32_t high; + + GET_FLOAT_WORD(hx,x); + sign=hx&0x80000000; /* sign= sign(x) */ + hx ^=sign; + if(hx>=0x7f800000) return(x+x); /* cbrt(NaN,INF) is itself */ + if(hx==0) + return(x); /* cbrt(0) is itself */ + + SET_FLOAT_WORD(x,hx); /* x <- |x| */ + /* rough cbrt to 5 bits */ + if(hx<0x00800000) /* subnormal number */ + {SET_FLOAT_WORD(t,0x4b800000); /* set t= 2**24 */ + t*=x; GET_FLOAT_WORD(high,t); SET_FLOAT_WORD(t,high/3+B2); + } + else + SET_FLOAT_WORD(t,hx/3+B1); + + + /* new cbrt to 23 bits */ + r=t*t/x; + s=C+r*t; + t*=G+F/(s+E+D/s); + + /* retore the sign bit */ + GET_FLOAT_WORD(high,t); + SET_FLOAT_WORD(t,high|sign); + return(t); +} + +#ifdef _DOUBLE_IS_32BITS + +#ifdef __STDC__ + double cbrt(double x) +#else + double cbrt(x) + double x; +#endif +{ + return (double) cbrtf((float) x); +} + +#endif /* defined(_DOUBLE_IS_32BITS) */ diff --git a/newlib/libm/common/sf_copysign.c b/newlib/libm/common/sf_copysign.c new file mode 100644 index 000000000..f547c82ed --- /dev/null +++ b/newlib/libm/common/sf_copysign.c @@ -0,0 +1,50 @@ +/* sf_copysign.c -- float version of s_copysign.c. + * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. + */ + +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* + * copysignf(float x, float y) + * copysignf(x,y) returns a value with the magnitude of x and + * with the sign bit of y. + */ + +#include "fdlibm.h" + +#ifdef __STDC__ + float copysignf(float x, float y) +#else + float copysignf(x,y) + float x,y; +#endif +{ + __uint32_t ix,iy; + GET_FLOAT_WORD(ix,x); + GET_FLOAT_WORD(iy,y); + SET_FLOAT_WORD(x,(ix&0x7fffffff)|(iy&0x80000000)); + return x; +} + +#ifdef _DOUBLE_IS_32BITS + +#ifdef __STDC__ + double copysign(double x, double y) +#else + double copysign(x,y) + double x,y; +#endif +{ + return (double) copysignf((float) x, (float) y); +} + +#endif /* defined(_DOUBLE_IS_32BITS) */ diff --git a/newlib/libm/common/sf_expm1.c b/newlib/libm/common/sf_expm1.c new file mode 100644 index 000000000..1df6f80d2 --- /dev/null +++ b/newlib/libm/common/sf_expm1.c @@ -0,0 +1,146 @@ +/* sf_expm1.c -- float version of s_expm1.c. + * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. + */ + +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +#include "fdlibm.h" + +#ifdef __v810__ +#define const +#endif + +#ifdef __STDC__ +static const float +#else +static float +#endif +one = 1.0, +huge = 1.0e+30, +tiny = 1.0e-30, +o_threshold = 8.8721679688e+01,/* 0x42b17180 */ +ln2_hi = 6.9313812256e-01,/* 0x3f317180 */ +ln2_lo = 9.0580006145e-06,/* 0x3717f7d1 */ +invln2 = 1.4426950216e+00,/* 0x3fb8aa3b */ + /* scaled coefficients related to expm1 */ +Q1 = -3.3333335072e-02, /* 0xbd088889 */ +Q2 = 1.5873016091e-03, /* 0x3ad00d01 */ +Q3 = -7.9365076090e-05, /* 0xb8a670cd */ +Q4 = 4.0082177293e-06, /* 0x36867e54 */ +Q5 = -2.0109921195e-07; /* 0xb457edbb */ + +#ifdef __STDC__ + float expm1f(float x) +#else + float expm1f(x) + float x; +#endif +{ + float y,hi,lo,c,t,e,hxs,hfx,r1; + __int32_t k,xsb; + __uint32_t hx; + + GET_FLOAT_WORD(hx,x); + xsb = hx&0x80000000; /* sign bit of x */ + if(xsb==0) y=x; else y= -x; /* y = |x| */ + hx &= 0x7fffffff; /* high word of |x| */ + + /* filter out huge and non-finite argument */ + if(hx >= 0x4195b844) { /* if |x|>=27*ln2 */ + if(hx >= 0x42b17218) { /* if |x|>=88.721... */ + if(hx>0x7f800000) + return x+x; /* NaN */ + if(hx==0x7f800000) + return (xsb==0)? x:-1.0;/* exp(+-inf)={inf,-1} */ + if(x > o_threshold) return huge*huge; /* overflow */ + } + if(xsb!=0) { /* x < -27*ln2, return -1.0 with inexact */ + if(x+tiny<(float)0.0) /* raise inexact */ + return tiny-one; /* return -1 */ + } + } + + /* argument reduction */ + if(hx > 0x3eb17218) { /* if |x| > 0.5 ln2 */ + if(hx < 0x3F851592) { /* and |x| < 1.5 ln2 */ + if(xsb==0) + {hi = x - ln2_hi; lo = ln2_lo; k = 1;} + else + {hi = x + ln2_hi; lo = -ln2_lo; k = -1;} + } else { + k = invln2*x+((xsb==0)?(float)0.5:(float)-0.5); + t = k; + hi = x - t*ln2_hi; /* t*ln2_hi is exact here */ + lo = t*ln2_lo; + } + x = hi - lo; + c = (hi-x)-lo; + } + else if(hx < 0x33000000) { /* when |x|<2**-25, return x */ + t = huge+x; /* return x with inexact flags when x!=0 */ + return x - (t-(huge+x)); + } + else k = 0; + + /* x is now in primary range */ + hfx = (float)0.5*x; + hxs = x*hfx; + r1 = one+hxs*(Q1+hxs*(Q2+hxs*(Q3+hxs*(Q4+hxs*Q5)))); + t = (float)3.0-r1*hfx; + e = hxs*((r1-t)/((float)6.0 - x*t)); + if(k==0) return x - (x*e-hxs); /* c is 0 */ + else { + e = (x*(e-c)-c); + e -= hxs; + if(k== -1) return (float)0.5*(x-e)-(float)0.5; + if(k==1) + if(x < (float)-0.25) return -(float)2.0*(e-(x+(float)0.5)); + else return one+(float)2.0*(x-e); + if (k <= -2 || k>56) { /* suffice to return exp(x)-1 */ + __int32_t i; + y = one-(e-x); + GET_FLOAT_WORD(i,y); + SET_FLOAT_WORD(y,i+(k<<23)); /* add k to y's exponent */ + return y-one; + } + t = one; + if(k<23) { + __int32_t i; + SET_FLOAT_WORD(t,0x3f800000 - (0x1000000>>k)); /* t=1-2^-k */ + y = t-(e-x); + GET_FLOAT_WORD(i,y); + SET_FLOAT_WORD(y,i+(k<<23)); /* add k to y's exponent */ + } else { + __int32_t i; + SET_FLOAT_WORD(t,((0x7f-k)<<23)); /* 2^-k */ + y = x-(e+t); + y += one; + GET_FLOAT_WORD(i,y); + SET_FLOAT_WORD(y,i+(k<<23)); /* add k to y's exponent */ + } + } + return y; +} + +#ifdef _DOUBLE_IS_32BITS + +#ifdef __STDC__ + double expm1(double x) +#else + double expm1(x) + double x; +#endif +{ + return (double) expm1f((float) x); +} + +#endif /* defined(_DOUBLE_IS_32BITS) */ diff --git a/newlib/libm/common/sf_finite.c b/newlib/libm/common/sf_finite.c new file mode 100644 index 000000000..4c48f400f --- /dev/null +++ b/newlib/libm/common/sf_finite.c @@ -0,0 +1,47 @@ +/* sf_finite.c -- float version of s_finite.c. + * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. + */ + +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* + * finitef(x) returns 1 is x is finite, else 0; + * no branching! + */ + +#include "fdlibm.h" + +#ifdef __STDC__ + int finitef(float x) +#else + int finitef(x) + float x; +#endif +{ + __int32_t ix; + GET_FLOAT_WORD(ix,x); + return (int)((__uint32_t)((ix&0x7fffffff)-0x7f800000)>>31); +} + +#ifdef _DOUBLE_IS_32BITS + +#ifdef __STDC__ + int finite(double x) +#else + int finite(x) + double x; +#endif +{ + return finitef((float) x); +} + +#endif /* defined(_DOUBLE_IS_32BITS) */ diff --git a/newlib/libm/common/sf_ilogb.c b/newlib/libm/common/sf_ilogb.c new file mode 100644 index 000000000..ee65594b1 --- /dev/null +++ b/newlib/libm/common/sf_ilogb.c @@ -0,0 +1,53 @@ +/* sf_ilogb.c -- float version of s_ilogb.c. + * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. + */ + +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +#include "fdlibm.h" +#include <limits.h> + +#ifdef __STDC__ + int ilogbf(float x) +#else + int ilogbf(x) + float x; +#endif +{ + __int32_t hx,ix; + + GET_FLOAT_WORD(hx,x); + hx &= 0x7fffffff; + if(hx<0x00800000) { + if(hx==0) + return - INT_MAX; /* ilogb(0) = 0x80000001 */ + else /* subnormal x */ + for (ix = -126,hx<<=8; hx>0; hx<<=1) ix -=1; + return ix; + } + else if (hx<0x7f800000) return (hx>>23)-127; + else return INT_MAX; +} + +#ifdef _DOUBLE_IS_32BITS + +#ifdef __STDC__ + int ilogb(double x) +#else + int ilogb(x) + double x; +#endif +{ + return ilogbf((float) x); +} + +#endif /* defined(_DOUBLE_IS_32BITS) */ diff --git a/newlib/libm/common/sf_infinity.c b/newlib/libm/common/sf_infinity.c new file mode 100644 index 000000000..8722596c9 --- /dev/null +++ b/newlib/libm/common/sf_infinity.c @@ -0,0 +1,23 @@ +/* + * infinityf () returns the representation of infinity. + * Added by Cygnus Support. + */ + +#include "fdlibm.h" + + float infinityf() +{ + float x; + + SET_FLOAT_WORD(x,0x7f800000); + return x; +} + +#ifdef _DOUBLE_IS_32BITS + + double infinity() +{ + return (double) infinityf(); +} + +#endif /* defined(_DOUBLE_IS_32BITS) */ diff --git a/newlib/libm/common/sf_log1p.c b/newlib/libm/common/sf_log1p.c new file mode 100644 index 000000000..31c426fe7 --- /dev/null +++ b/newlib/libm/common/sf_log1p.c @@ -0,0 +1,121 @@ +/* sf_log1p.c -- float version of s_log1p.c. + * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. + */ + +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +#include "fdlibm.h" + +#ifdef __STDC__ +static const float +#else +static float +#endif +ln2_hi = 6.9313812256e-01, /* 0x3f317180 */ +ln2_lo = 9.0580006145e-06, /* 0x3717f7d1 */ +two25 = 3.355443200e+07, /* 0x4c000000 */ +Lp1 = 6.6666668653e-01, /* 3F2AAAAB */ +Lp2 = 4.0000000596e-01, /* 3ECCCCCD */ +Lp3 = 2.8571429849e-01, /* 3E924925 */ +Lp4 = 2.2222198546e-01, /* 3E638E29 */ +Lp5 = 1.8183572590e-01, /* 3E3A3325 */ +Lp6 = 1.5313838422e-01, /* 3E1CD04F */ +Lp7 = 1.4798198640e-01; /* 3E178897 */ + +#ifdef __STDC__ +static const float zero = 0.0; +#else +static float zero = 0.0; +#endif + +#ifdef __STDC__ + float log1pf(float x) +#else + float log1pf(x) + float x; +#endif +{ + float hfsq,f,c,s,z,R,u; + __int32_t k,hx,hu,ax; + + GET_FLOAT_WORD(hx,x); + ax = hx&0x7fffffff; + + k = 1; + if (hx < 0x3ed413d7) { /* x < 0.41422 */ + if(ax>=0x3f800000) { /* x <= -1.0 */ + if(x==(float)-1.0) return -two25/zero; /* log1p(-1)=+inf */ + else return (x-x)/(x-x); /* log1p(x<-1)=NaN */ + } + if(ax<0x31000000) { /* |x| < 2**-29 */ + if(two25+x>zero /* raise inexact */ + &&ax<0x24800000) /* |x| < 2**-54 */ + return x; + else + return x - x*x*(float)0.5; + } + if(hx>0||hx<=((__int32_t)0xbe95f61f)) { + k=0;f=x;hu=1;} /* -0.2929<x<0.41422 */ + } + if (hx >= 0x7f800000) return x+x; + if(k!=0) { + if(hx<0x5a000000) { + u = (float)1.0+x; + GET_FLOAT_WORD(hu,u); + k = (hu>>23)-127; + /* correction term */ + c = (k>0)? (float)1.0-(u-x):x-(u-(float)1.0); + c /= u; + } else { + u = x; + GET_FLOAT_WORD(hu,u); + k = (hu>>23)-127; + c = 0; + } + hu &= 0x007fffff; + if(hu<0x3504f7) { + SET_FLOAT_WORD(u,hu|0x3f800000);/* normalize u */ + } else { + k += 1; + SET_FLOAT_WORD(u,hu|0x3f000000); /* normalize u/2 */ + hu = (0x00800000-hu)>>2; + } + f = u-(float)1.0; + } + hfsq=(float)0.5*f*f; + if(hu==0) { /* |f| < 2**-20 */ + if(f==zero) if(k==0) return zero; + else {c += k*ln2_lo; return k*ln2_hi+c;} + R = hfsq*((float)1.0-(float)0.66666666666666666*f); + if(k==0) return f-R; else + return k*ln2_hi-((R-(k*ln2_lo+c))-f); + } + s = f/((float)2.0+f); + z = s*s; + R = z*(Lp1+z*(Lp2+z*(Lp3+z*(Lp4+z*(Lp5+z*(Lp6+z*Lp7)))))); + if(k==0) return f-(hfsq-s*(hfsq+R)); else + return k*ln2_hi-((hfsq-(s*(hfsq+R)+(k*ln2_lo+c)))-f); +} + +#ifdef _DOUBLE_IS_32BITS + +#ifdef __STDC__ + double log1p(double x) +#else + double log1p(x) + double x; +#endif +{ + return (double) log1pf((float) x); +} + +#endif /* defined(_DOUBLE_IS_32BITS) */ diff --git a/newlib/libm/common/sf_logb.c b/newlib/libm/common/sf_logb.c new file mode 100644 index 000000000..6e67637fe --- /dev/null +++ b/newlib/libm/common/sf_logb.c @@ -0,0 +1,48 @@ +/* sf_logb.c -- float version of s_logb.c. + * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. + */ + +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +#include "fdlibm.h" + +#ifdef __STDC__ + float logbf(float x) +#else + float logbf(x) + float x; +#endif +{ + __int32_t ix; + GET_FLOAT_WORD(ix,x); + ix &= 0x7fffffff; /* high |x| */ + if(ix==0) return (float)-1.0/fabsf(x); + if(ix>=0x7f800000) return x*x; + if((ix>>=23)==0) /* IEEE 754 logb */ + return -126.0; + else + return (float) (ix-127); +} + +#ifdef _DOUBLE_IS_32BITS + +#ifdef __STDC__ + double logb(double x) +#else + double logb(x) + double x; +#endif +{ + return (double) logbf((float) x); +} + +#endif /* defined(_DOUBLE_IS_32BITS) */ diff --git a/newlib/libm/common/sf_modf.c b/newlib/libm/common/sf_modf.c new file mode 100644 index 000000000..6c64e3fa0 --- /dev/null +++ b/newlib/libm/common/sf_modf.c @@ -0,0 +1,73 @@ +/* sf_modf.c -- float version of s_modf.c. + * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. + */ + +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +#include "fdlibm.h" + +#ifdef __STDC__ +static const float one = 1.0; +#else +static float one = 1.0; +#endif + +#ifdef __STDC__ + float modff(float x, float *iptr) +#else + float modff(x, iptr) + float x,*iptr; +#endif +{ + __int32_t i0,j0; + __uint32_t i; + GET_FLOAT_WORD(i0,x); + j0 = ((i0>>23)&0xff)-0x7f; /* exponent of x */ + if(j0<23) { /* integer part in x */ + if(j0<0) { /* |x|<1 */ + SET_FLOAT_WORD(*iptr,i0&0x80000000); /* *iptr = +-0 */ + return x; + } else { + i = (0x007fffff)>>j0; + if((i0&i)==0) { /* x is integral */ + __uint32_t ix; + *iptr = x; + GET_FLOAT_WORD(ix,x); + SET_FLOAT_WORD(x,ix&0x80000000); /* return +-0 */ + return x; + } else { + SET_FLOAT_WORD(*iptr,i0&(~i)); + return x - *iptr; + } + } + } else { /* no fraction part */ + __uint32_t ix; + *iptr = x*one; + GET_FLOAT_WORD(ix,x); + SET_FLOAT_WORD(x,ix&0x80000000); /* return +-0 */ + return x; + } +} + +#ifdef _DOUBLE_IS_32BITS + +#ifdef __STDC__ + double modf(double x, double *iptr) +#else + double modf(x, iptr) + double x,*iptr; +#endif +{ + return (double) modff((float) x, (float *) iptr); +} + +#endif /* defined(_DOUBLE_IS_32BITS) */ diff --git a/newlib/libm/common/sf_nan.c b/newlib/libm/common/sf_nan.c new file mode 100644 index 000000000..cb3e1cd03 --- /dev/null +++ b/newlib/libm/common/sf_nan.c @@ -0,0 +1,23 @@ +/* + * nanf () returns a nan. + * Added by Cygnus Support. + */ + +#include "fdlibm.h" + + float nanf() +{ + float x; + + SET_FLOAT_WORD(x,0x7fc00000); + return x; +} + +#ifdef _DOUBLE_IS_32BITS + + double nan() +{ + return (double) nanf(); +} + +#endif /* defined(_DOUBLE_IS_32BITS) */ diff --git a/newlib/libm/common/sf_nextafter.c b/newlib/libm/common/sf_nextafter.c new file mode 100644 index 000000000..cd938d339 --- /dev/null +++ b/newlib/libm/common/sf_nextafter.c @@ -0,0 +1,79 @@ +/* sf_nextafter.c -- float version of s_nextafter.c. + * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. + */ + +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +#include "fdlibm.h" + +#ifdef __STDC__ + float nextafterf(float x, float y) +#else + float nextafterf(x,y) + float x,y; +#endif +{ + __int32_t hx,hy,ix,iy; + + GET_FLOAT_WORD(hx,x); + GET_FLOAT_WORD(hy,y); + ix = hx&0x7fffffff; /* |x| */ + iy = hy&0x7fffffff; /* |y| */ + + if((ix>0x7f800000) || /* x is nan */ + (iy>0x7f800000)) /* y is nan */ + return x+y; + if(x==y) return x; /* x=y, return x */ + if(ix==0) { /* x == 0 */ + SET_FLOAT_WORD(x,(hy&0x80000000)|1);/* return +-minsubnormal */ + y = x*x; + if(y==x) return y; else return x; /* raise underflow flag */ + } + if(hx>=0) { /* x > 0 */ + if(hx>hy) { /* x > y, x -= ulp */ + hx -= 1; + } else { /* x < y, x += ulp */ + hx += 1; + } + } else { /* x < 0 */ + if(hy>=0||hx>hy){ /* x < y, x -= ulp */ + hx -= 1; + } else { /* x > y, x += ulp */ + hx += 1; + } + } + hy = hx&0x7f800000; + if(hy>=0x7f800000) return x+x; /* overflow */ + if(hy<0x00800000) { /* underflow */ + y = x*x; + if(y!=x) { /* raise underflow flag */ + SET_FLOAT_WORD(y,hx); + return y; + } + } + SET_FLOAT_WORD(x,hx); + return x; +} + +#ifdef _DOUBLE_IS_32BITS + +#ifdef __STDC__ + double nextafter(double x, double y) +#else + double nextafter(x,y) + double x,y; +#endif +{ + return (double) nextafterf((float) x, (float) x); +} + +#endif /* defined(_DOUBLE_IS_32BITS) */ diff --git a/newlib/libm/common/sf_rint.c b/newlib/libm/common/sf_rint.c new file mode 100644 index 000000000..d38080a5d --- /dev/null +++ b/newlib/libm/common/sf_rint.c @@ -0,0 +1,81 @@ +/* sf_rint.c -- float version of s_rint.c. + * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. + */ + +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +#include "fdlibm.h" + +#ifdef __STDC__ +static const float +#else +static float +#endif +TWO23[2]={ + 8.3886080000e+06, /* 0x4b000000 */ + -8.3886080000e+06, /* 0xcb000000 */ +}; + +#ifdef __STDC__ + float rintf(float x) +#else + float rintf(x) + float x; +#endif +{ + __int32_t i0,j0,sx; + __uint32_t i,i1; + float t; + volatile float w; + GET_FLOAT_WORD(i0,x); + sx = (i0>>31)&1; + j0 = ((i0>>23)&0xff)-0x7f; + if(j0<23) { + if(j0<0) { + if((i0&0x7fffffff)==0) return x; + i1 = (i0&0x07fffff); + i0 &= 0xfff00000; + i0 |= ((i1|-i1)>>9)&0x400000; + SET_FLOAT_WORD(x,i0); + w = TWO23[sx]+x; + t = w-TWO23[sx]; + GET_FLOAT_WORD(i0,t); + SET_FLOAT_WORD(t,(i0&0x7fffffff)|(sx<<31)); + return t; + } else { + i = (0x007fffff)>>j0; + if((i0&i)==0) return x; /* x is integral */ + i>>=1; + if((i0&i)!=0) i0 = (i0&(~i))|((0x100000)>>j0); + } + } else { + if(j0==0x80) return x+x; /* inf or NaN */ + else return x; /* x is integral */ + } + SET_FLOAT_WORD(x,i0); + w = TWO23[sx]+x; + return w-TWO23[sx]; +} + +#ifdef _DOUBLE_IS_32BITS + +#ifdef __STDC__ + double rint(double x) +#else + double rint(x) + double x; +#endif +{ + return (double) rintf((float) x); +} + +#endif /* defined(_DOUBLE_IS_32BITS) */ diff --git a/newlib/libm/common/sf_scalbn.c b/newlib/libm/common/sf_scalbn.c new file mode 100644 index 000000000..ee65f40ad --- /dev/null +++ b/newlib/libm/common/sf_scalbn.c @@ -0,0 +1,79 @@ +/* sf_scalbn.c -- float version of s_scalbn.c. + * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. + */ + +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +#include "fdlibm.h" +#include <limits.h> + +#if INT_MAX > 50000 +#define OVERFLOW_INT 50000 +#else +#define OVERFLOW_INT 30000 +#endif + +#ifdef __STDC__ +static const float +#else +static float +#endif +two25 = 3.355443200e+07, /* 0x4c000000 */ +twom25 = 2.9802322388e-08, /* 0x33000000 */ +huge = 1.0e+30, +tiny = 1.0e-30; + +#ifdef __STDC__ + float scalbnf (float x, int n) +#else + float scalbnf (x,n) + float x; int n; +#endif +{ + __int32_t k,ix; + GET_FLOAT_WORD(ix,x); + k = (ix&0x7f800000)>>23; /* extract exponent */ + if (k==0) { /* 0 or subnormal x */ + if ((ix&0x7fffffff)==0) return x; /* +-0 */ + x *= two25; + GET_FLOAT_WORD(ix,x); + k = ((ix&0x7f800000)>>23) - 25; + if (n< -50000) return tiny*x; /*underflow*/ + } + if (k==0xff) return x+x; /* NaN or Inf */ + k = k+n; + if (k > 0xfe) return huge*copysignf(huge,x); /* overflow */ + if (k > 0) /* normal result */ + {SET_FLOAT_WORD(x,(ix&0x807fffff)|(k<<23)); return x;} + if (k <= -25) + if (n > OVERFLOW_INT) /* in case integer overflow in n+k */ + return huge*copysignf(huge,x); /*overflow*/ + else return tiny*copysignf(tiny,x); /*underflow*/ + k += 25; /* subnormal result */ + SET_FLOAT_WORD(x,(ix&0x807fffff)|(k<<23)); + return x*twom25; +} + +#ifdef _DOUBLE_IS_32BITS + +#ifdef __STDC__ + double scalbn(double x, int n) +#else + double scalbn(x,n) + double x; + int n; +#endif +{ + return (double) scalbnf((float) x, n); +} + +#endif /* defined(_DOUBLE_IS_32BITS) */ |