import newlib-2000-02-17 snapshot
This commit is contained in:
Christopher Faylor
97
newlib/libm/common/Makefile.am
Normal file
97
newlib/libm/common/Makefile.am
Normal file
@ -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
|
||||
341
newlib/libm/common/Makefile.in
Normal file
341
newlib/libm/common/Makefile.in
Normal file
@ -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:
|
||||
12
newlib/libm/common/common.tex
Normal file
12
newlib/libm/common/common.tex
Normal file
@ -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
|
||||
261
newlib/libm/common/fdlibm.h
Normal file
261
newlib/libm/common/fdlibm.h
Normal file
@ -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)
|
||||
123
newlib/libm/common/s_cbrt.c
Normal file
123
newlib/libm/common/s_cbrt.c
Normal file
@ -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 */
|
||||
82
newlib/libm/common/s_copysign.c
Normal file
82
newlib/libm/common/s_copysign.c
Normal file
@ -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 */
|
||||
271
newlib/libm/common/s_expm1.c
Normal file
271
newlib/libm/common/s_expm1.c
Normal file
@ -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 */
|
||||
35
newlib/libm/common/s_finite.c
Normal file
35
newlib/libm/common/s_finite.c
Normal file
@ -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 */
|
||||
92
newlib/libm/common/s_ilogb.c
Normal file
92
newlib/libm/common/s_ilogb.c
Normal file
@ -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 */
|
||||
48
newlib/libm/common/s_infinity.c
Normal file
48
newlib/libm/common/s_infinity.c
Normal file
@ -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 */
|
||||
35
newlib/libm/common/s_lib_ver.c
Normal file
35
newlib/libm/common/s_lib_ver.c
Normal file
@ -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
|
||||
217
newlib/libm/common/s_log1p.c
Normal file
217
newlib/libm/common/s_log1p.c
Normal file
@ -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 */
|
||||
42
newlib/libm/common/s_logb.c
Normal file
42
newlib/libm/common/s_logb.c
Normal file
@ -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 */
|
||||
123
newlib/libm/common/s_matherr.c
Normal file
123
newlib/libm/common/s_matherr.c
Normal file
@ -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;
|
||||
}
|
||||
131
newlib/libm/common/s_modf.c
Normal file
131
newlib/libm/common/s_modf.c
Normal file
@ -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 */
|
||||
47
newlib/libm/common/s_nan.c
Normal file
47
newlib/libm/common/s_nan.c
Normal file
@ -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 */
|
||||
121
newlib/libm/common/s_nextafter.c
Normal file
121
newlib/libm/common/s_nextafter.c
Normal file
@ -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 */
|
||||
90
newlib/libm/common/s_rint.c
Normal file
90
newlib/libm/common/s_rint.c
Normal file
@ -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 */
|
||||
|
||||
|
||||
|
||||
103
newlib/libm/common/s_scalbn.c
Normal file
103
newlib/libm/common/s_scalbn.c
Normal file
@ -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 */
|
||||
93
newlib/libm/common/sf_cbrt.c
Normal file
93
newlib/libm/common/sf_cbrt.c
Normal file
@ -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) */
|
||||
50
newlib/libm/common/sf_copysign.c
Normal file
50
newlib/libm/common/sf_copysign.c
Normal file
@ -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) */
|
||||
146
newlib/libm/common/sf_expm1.c
Normal file
146
newlib/libm/common/sf_expm1.c
Normal file
@ -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) */
|
||||
47
newlib/libm/common/sf_finite.c
Normal file
47
newlib/libm/common/sf_finite.c
Normal file
@ -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) */
|
||||
53
newlib/libm/common/sf_ilogb.c
Normal file
53
newlib/libm/common/sf_ilogb.c
Normal file
@ -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) */
|
||||
23
newlib/libm/common/sf_infinity.c
Normal file
23
newlib/libm/common/sf_infinity.c
Normal file
@ -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) */
|
||||
121
newlib/libm/common/sf_log1p.c
Normal file
121
newlib/libm/common/sf_log1p.c
Normal file
@ -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) */
|
||||
48
newlib/libm/common/sf_logb.c
Normal file
48
newlib/libm/common/sf_logb.c
Normal file
@ -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) */
|
||||
73
newlib/libm/common/sf_modf.c
Normal file
73
newlib/libm/common/sf_modf.c
Normal file
@ -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) */
|
||||
23
newlib/libm/common/sf_nan.c
Normal file
23
newlib/libm/common/sf_nan.c
Normal file
@ -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) */
|
||||
79
newlib/libm/common/sf_nextafter.c
Normal file
79
newlib/libm/common/sf_nextafter.c
Normal file
@ -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) */
|
||||
81
newlib/libm/common/sf_rint.c
Normal file
81
newlib/libm/common/sf_rint.c
Normal file
@ -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) */
|
||||
79
newlib/libm/common/sf_scalbn.c
Normal file
79
newlib/libm/common/sf_scalbn.c
Normal file
@ -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) */
|
||||
Reference in New Issue
Block a user