diff options
Diffstat (limited to 'gsl-1.9/sum')
| -rw-r--r-- | gsl-1.9/sum/ChangeLog | 84 | ||||
| -rw-r--r-- | gsl-1.9/sum/Makefile.am | 17 | ||||
| -rw-r--r-- | gsl-1.9/sum/Makefile.in | 544 | ||||
| -rw-r--r-- | gsl-1.9/sum/gsl_sum.h | 162 | ||||
| -rw-r--r-- | gsl-1.9/sum/levin_u.c | 255 | ||||
| -rw-r--r-- | gsl-1.9/sum/levin_utrunc.c | 202 | ||||
| -rw-r--r-- | gsl-1.9/sum/test.c | 234 | ||||
| -rw-r--r-- | gsl-1.9/sum/work_u.c | 94 | ||||
| -rw-r--r-- | gsl-1.9/sum/work_utrunc.c | 67 |
9 files changed, 1659 insertions, 0 deletions
diff --git a/gsl-1.9/sum/ChangeLog b/gsl-1.9/sum/ChangeLog new file mode 100644 index 0000000..400e925 --- /dev/null +++ b/gsl-1.9/sum/ChangeLog @@ -0,0 +1,84 @@ +Thu Dec 21 21:46:54 2000 Brian Gough <bjg@network-theory.co.uk> + + * changed err to abserr, to make it clear that it is an absolute + error in line with other modules. + +Thu Nov 2 20:08:14 2000 Brian Gough <bjg@network-theory.co.uk> + + * added support for a workspace so that the user does not have to + allocate memory + + * made the names of the functions consistent as either levin_u, + for the full u transform with error estimate, or levin_utrunc for + the u transform with only a truncation error estimate. + +Mon Apr 24 21:15:27 2000 Brian Gough <bjg@network-theory.co.uk> + + * gsl_sum.h: added #include <stdlib.h> for size_t + +Mon Nov 1 12:50:17 1999 Brian Gough <bjg@network-theory.co.uk> + + * test.c (main): added tests using series for the Euler constant + and eta(1/2) + + * resolved problems with spurious failures by replacing the + directly computed truncation error by an estimate which varies + more smoothly. I have used the average of the previous two values, + which seems to give a reliable estimate of the truncation + error. The direct evaluation of the truncation error sometimes + fluctuated wildly, due to cancellation effects. + +Thu Oct 28 12:05:47 1999 Brian Gough <bjg@network-theory.co.uk> + + * test.c: cleaned up tests, now find that everything works in + double-precision but not extended-precision where there are two + failures depending on the optimization level + + * levin_uerr.c (gsl_sum_levin_u_accel_minmax): changed loop + maximum from <=n to <n to avoid access of uninitialized memory, + need further correctness checks on algorithm though. + +Thu Oct 7 12:03:36 1999 Brian Gough <bjg@network-theory.co.uk> + + * levin_u.c levin_uerr.c: changed DBL_MAX to GSL_DBL_MAX since we + don't rely on DBL_MAX + +Sat Feb 6 20:35:26 1999 Brian Gough <bjg@netsci.freeserve.co.uk> + + * test.c: adjusted the precision check to allow for "infinite + accuracy" which occurs when two results agree to machine precision + +Thu Nov 19 13:10:19 1998 Brian Gough <bjg@vvv.lanl.gov> + + * added an n_used parameter to all routines which gives the number + of terms actually used + +Tue Nov 17 12:31:03 1998 Brian Gough <bjg@vvv.lanl.gov> + + * test.c: added #include <config.h> + + * renamed test_sum.c to test.c + + * renamed all the functions so that _with_derivs is now the + default and _trunc is the case of no error estimates from the + derivatives + + * test_sum.c (main): cleaned up tests + +Mon Nov 9 22:05:45 1998 Brian Gough <bjg@vvv.lanl.gov> + + * levin_u.c (gsl_sum_levin_u_accel_minmax): got rid of noise + variables since they aren't used in the non-derivative case + +1998-11-06 <bjg@ancho.lanl.gov> + + * test_sum.c: replace variable N by macro to avoid variable length + array warning + +Tue Oct 27 18:06:16 1998 Brian Gough <bjg@vvv.lanl.gov> + + * levin_u.c: added in noise but it does not seem to be giving the + right answer for the error estimate. The actual value for the + accelerated sum is correct though. Check toms/602 for the original + algorithm. + diff --git a/gsl-1.9/sum/Makefile.am b/gsl-1.9/sum/Makefile.am new file mode 100644 index 0000000..df348aa --- /dev/null +++ b/gsl-1.9/sum/Makefile.am @@ -0,0 +1,17 @@ +noinst_LTLIBRARIES = libgslsum.la + +pkginclude_HEADERS = gsl_sum.h + +INCLUDES= -I$(top_builddir) + +libgslsum_la_SOURCES = levin_u.c levin_utrunc.c work_u.c work_utrunc.c + +TESTS = $(check_PROGRAMS) + +check_PROGRAMS = test + +test_LDADD = libgslsum.la ../ieee-utils/libgslieeeutils.la ../err/libgslerr.la ../test/libgsltest.la ../sys/libgslsys.la ../utils/libutils.la + +test_SOURCES = test.c + + diff --git a/gsl-1.9/sum/Makefile.in b/gsl-1.9/sum/Makefile.in new file mode 100644 index 0000000..198d27b --- /dev/null +++ b/gsl-1.9/sum/Makefile.in @@ -0,0 +1,544 @@ +# Makefile.in generated by automake 1.9.6 from Makefile.am. +# @configure_input@ + +# Copyright (C) 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2001, 2002, +# 2003, 2004, 2005 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. + +@SET_MAKE@ + + +srcdir = @srcdir@ +top_srcdir = @top_srcdir@ +VPATH = @srcdir@ +pkgdatadir = $(datadir)/@PACKAGE@ +pkglibdir = $(libdir)/@PACKAGE@ +pkgincludedir = $(includedir)/@PACKAGE@ +top_builddir = .. +am__cd = CDPATH="$${ZSH_VERSION+.}$(PATH_SEPARATOR)" && cd +INSTALL = @INSTALL@ +install_sh_DATA = $(install_sh) -c -m 644 +install_sh_PROGRAM = $(install_sh) -c +install_sh_SCRIPT = $(install_sh) -c +INSTALL_HEADER = $(INSTALL_DATA) +transform = $(program_transform_name) +NORMAL_INSTALL = : +PRE_INSTALL = : +POST_INSTALL = : +NORMAL_UNINSTALL = : +PRE_UNINSTALL = : +POST_UNINSTALL = : +build_triplet = @build@ +host_triplet = @host@ +check_PROGRAMS = test$(EXEEXT) +subdir = sum +DIST_COMMON = $(pkginclude_HEADERS) $(srcdir)/Makefile.am \ + $(srcdir)/Makefile.in ChangeLog +ACLOCAL_M4 = $(top_srcdir)/aclocal.m4 +am__aclocal_m4_deps = $(top_srcdir)/configure.ac +am__configure_deps = $(am__aclocal_m4_deps) $(CONFIGURE_DEPENDENCIES) \ + $(ACLOCAL_M4) +mkinstalldirs = $(SHELL) $(top_srcdir)/mkinstalldirs +CONFIG_HEADER = $(top_builddir)/config.h +CONFIG_CLEAN_FILES = +LTLIBRARIES = $(noinst_LTLIBRARIES) +libgslsum_la_LIBADD = +am_libgslsum_la_OBJECTS = levin_u.lo levin_utrunc.lo work_u.lo \ + work_utrunc.lo +libgslsum_la_OBJECTS = $(am_libgslsum_la_OBJECTS) +am_test_OBJECTS = test.$(OBJEXT) +test_OBJECTS = $(am_test_OBJECTS) +test_DEPENDENCIES = libgslsum.la ../ieee-utils/libgslieeeutils.la \ + ../err/libgslerr.la ../test/libgsltest.la ../sys/libgslsys.la \ + ../utils/libutils.la +DEFAULT_INCLUDES = -I. -I$(srcdir) -I$(top_builddir) +depcomp = +am__depfiles_maybe = +COMPILE = $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) \ + $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) +LTCOMPILE = $(LIBTOOL) --tag=CC --mode=compile $(CC) $(DEFS) \ + $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) \ + $(AM_CFLAGS) $(CFLAGS) +CCLD = $(CC) +LINK = $(LIBTOOL) --tag=CC --mode=link $(CCLD) $(AM_CFLAGS) $(CFLAGS) \ + $(AM_LDFLAGS) $(LDFLAGS) -o $@ +SOURCES = $(libgslsum_la_SOURCES) $(test_SOURCES) +DIST_SOURCES = $(libgslsum_la_SOURCES) $(test_SOURCES) +am__vpath_adj_setup = srcdirstrip=`echo "$(srcdir)" | sed 's|.|.|g'`; +am__vpath_adj = case $$p in \ + $(srcdir)/*) f=`echo "$$p" | sed "s|^$$srcdirstrip/||"`;; \ + *) f=$$p;; \ + esac; +am__strip_dir = `echo $$p | sed -e 's|^.*/||'`; +am__installdirs = "$(DESTDIR)$(pkgincludedir)" +pkgincludeHEADERS_INSTALL = $(INSTALL_HEADER) +HEADERS = $(pkginclude_HEADERS) +ETAGS = etags +CTAGS = ctags +DISTFILES = $(DIST_COMMON) $(DIST_SOURCES) $(TEXINFOS) $(EXTRA_DIST) +ACLOCAL = @ACLOCAL@ +AMTAR = @AMTAR@ +AR = @AR@ +AUTOCONF = @AUTOCONF@ +AUTOHEADER = @AUTOHEADER@ +AUTOMAKE = @AUTOMAKE@ +AWK = @AWK@ +CC = @CC@ +CFLAGS = @CFLAGS@ +CPP = @CPP@ +CPPFLAGS = @CPPFLAGS@ +CYGPATH_W = @CYGPATH_W@ +DEFS = @DEFS@ +ECHO = @ECHO@ +ECHO_C = @ECHO_C@ +ECHO_N = @ECHO_N@ +ECHO_T = @ECHO_T@ +EGREP = @EGREP@ +EXEEXT = @EXEEXT@ +GSL_CFLAGS = @GSL_CFLAGS@ +GSL_LIBS = @GSL_LIBS@ +GSL_LT_CBLAS_VERSION = @GSL_LT_CBLAS_VERSION@ +GSL_LT_VERSION = @GSL_LT_VERSION@ +HAVE_AIX_IEEE_INTERFACE = @HAVE_AIX_IEEE_INTERFACE@ +HAVE_DARWIN86_IEEE_INTERFACE = @HAVE_DARWIN86_IEEE_INTERFACE@ +HAVE_DARWIN_IEEE_INTERFACE = @HAVE_DARWIN_IEEE_INTERFACE@ +HAVE_EXTENDED_PRECISION_REGISTERS = @HAVE_EXTENDED_PRECISION_REGISTERS@ +HAVE_FREEBSD_IEEE_INTERFACE = @HAVE_FREEBSD_IEEE_INTERFACE@ +HAVE_GNUM68K_IEEE_INTERFACE = @HAVE_GNUM68K_IEEE_INTERFACE@ +HAVE_GNUPPC_IEEE_INTERFACE = @HAVE_GNUPPC_IEEE_INTERFACE@ +HAVE_GNUSPARC_IEEE_INTERFACE = @HAVE_GNUSPARC_IEEE_INTERFACE@ +HAVE_GNUX86_IEEE_INTERFACE = @HAVE_GNUX86_IEEE_INTERFACE@ +HAVE_HPUX11_IEEE_INTERFACE = @HAVE_HPUX11_IEEE_INTERFACE@ +HAVE_HPUX_IEEE_INTERFACE = @HAVE_HPUX_IEEE_INTERFACE@ +HAVE_IEEE_COMPARISONS = @HAVE_IEEE_COMPARISONS@ +HAVE_IEEE_DENORMALS = @HAVE_IEEE_DENORMALS@ +HAVE_INLINE = @HAVE_INLINE@ +HAVE_IRIX_IEEE_INTERFACE = @HAVE_IRIX_IEEE_INTERFACE@ +HAVE_NETBSD_IEEE_INTERFACE = @HAVE_NETBSD_IEEE_INTERFACE@ +HAVE_OPENBSD_IEEE_INTERFACE = @HAVE_OPENBSD_IEEE_INTERFACE@ +HAVE_OS2EMX_IEEE_INTERFACE = @HAVE_OS2EMX_IEEE_INTERFACE@ +HAVE_PRINTF_LONGDOUBLE = @HAVE_PRINTF_LONGDOUBLE@ +HAVE_SOLARIS_IEEE_INTERFACE = @HAVE_SOLARIS_IEEE_INTERFACE@ +HAVE_SUNOS4_IEEE_INTERFACE = @HAVE_SUNOS4_IEEE_INTERFACE@ +HAVE_TRU64_IEEE_INTERFACE = @HAVE_TRU64_IEEE_INTERFACE@ +INSTALL_DATA = @INSTALL_DATA@ +INSTALL_PROGRAM = @INSTALL_PROGRAM@ +INSTALL_SCRIPT = @INSTALL_SCRIPT@ +INSTALL_STRIP_PROGRAM = @INSTALL_STRIP_PROGRAM@ +LDFLAGS = @LDFLAGS@ +LIBOBJS = @LIBOBJS@ +LIBS = @LIBS@ +LIBTOOL = @LIBTOOL@ +LN_S = @LN_S@ +LTLIBOBJS = @LTLIBOBJS@ +MAINT = @MAINT@ +MAINTAINER_MODE_FALSE = @MAINTAINER_MODE_FALSE@ +MAINTAINER_MODE_TRUE = @MAINTAINER_MODE_TRUE@ +MAKEINFO = @MAKEINFO@ +OBJEXT = @OBJEXT@ +PACKAGE = @PACKAGE@ +PACKAGE_BUGREPORT = @PACKAGE_BUGREPORT@ +PACKAGE_NAME = @PACKAGE_NAME@ +PACKAGE_STRING = @PACKAGE_STRING@ +PACKAGE_TARNAME = @PACKAGE_TARNAME@ +PACKAGE_VERSION = @PACKAGE_VERSION@ +PATH_SEPARATOR = @PATH_SEPARATOR@ +RANLIB = @RANLIB@ +RELEASED = @RELEASED@ +SET_MAKE = @SET_MAKE@ +SHELL = @SHELL@ +STRIP = @STRIP@ +VERSION = @VERSION@ +ac_ct_AR = @ac_ct_AR@ +ac_ct_CC = @ac_ct_CC@ +ac_ct_RANLIB = @ac_ct_RANLIB@ +ac_ct_STRIP = @ac_ct_STRIP@ +am__leading_dot = @am__leading_dot@ +am__tar = @am__tar@ +am__untar = @am__untar@ +bindir = @bindir@ +build = @build@ +build_alias = @build_alias@ +build_cpu = @build_cpu@ +build_os = @build_os@ +build_vendor = @build_vendor@ +datadir = @datadir@ +exec_prefix = @exec_prefix@ +host = @host@ +host_alias = @host_alias@ +host_cpu = @host_cpu@ +host_os = @host_os@ +host_vendor = @host_vendor@ +includedir = @includedir@ +infodir = @infodir@ +install_sh = @install_sh@ +libdir = @libdir@ +libexecdir = @libexecdir@ +localstatedir = @localstatedir@ +mandir = @mandir@ +mkdir_p = @mkdir_p@ +oldincludedir = @oldincludedir@ +prefix = @prefix@ +program_transform_name = @program_transform_name@ +sbindir = @sbindir@ +sharedstatedir = @sharedstatedir@ +sysconfdir = @sysconfdir@ +target_alias = @target_alias@ +noinst_LTLIBRARIES = libgslsum.la +pkginclude_HEADERS = gsl_sum.h +INCLUDES = -I$(top_builddir) +libgslsum_la_SOURCES = levin_u.c levin_utrunc.c work_u.c work_utrunc.c +TESTS = $(check_PROGRAMS) +test_LDADD = libgslsum.la ../ieee-utils/libgslieeeutils.la ../err/libgslerr.la ../test/libgsltest.la ../sys/libgslsys.la ../utils/libutils.la +test_SOURCES = test.c +all: all-am + +.SUFFIXES: +.SUFFIXES: .c .lo .o .obj +$(srcdir)/Makefile.in: @MAINTAINER_MODE_TRUE@ $(srcdir)/Makefile.am $(am__configure_deps) + @for dep in $?; do \ + case '$(am__configure_deps)' in \ + *$$dep*) \ + cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh \ + && exit 0; \ + exit 1;; \ + esac; \ + done; \ + echo ' cd $(top_srcdir) && $(AUTOMAKE) --gnu --ignore-deps sum/Makefile'; \ + cd $(top_srcdir) && \ + $(AUTOMAKE) --gnu --ignore-deps sum/Makefile +.PRECIOUS: Makefile +Makefile: $(srcdir)/Makefile.in $(top_builddir)/config.status + @case '$?' in \ + *config.status*) \ + cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh;; \ + *) \ + echo ' cd $(top_builddir) && $(SHELL) ./config.status $(subdir)/$@ $(am__depfiles_maybe)'; \ + cd $(top_builddir) && $(SHELL) ./config.status $(subdir)/$@ $(am__depfiles_maybe);; \ + esac; + +$(top_builddir)/config.status: $(top_srcdir)/configure $(CONFIG_STATUS_DEPENDENCIES) + cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh + +$(top_srcdir)/configure: @MAINTAINER_MODE_TRUE@ $(am__configure_deps) + cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh +$(ACLOCAL_M4): @MAINTAINER_MODE_TRUE@ $(am__aclocal_m4_deps) + cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh + +clean-noinstLTLIBRARIES: + -test -z "$(noinst_LTLIBRARIES)" || rm -f $(noinst_LTLIBRARIES) + @list='$(noinst_LTLIBRARIES)'; for p in $$list; do \ + dir="`echo $$p | sed -e 's|/[^/]*$$||'`"; \ + test "$$dir" != "$$p" || dir=.; \ + echo "rm -f \"$${dir}/so_locations\""; \ + rm -f "$${dir}/so_locations"; \ + done +libgslsum.la: $(libgslsum_la_OBJECTS) $(libgslsum_la_DEPENDENCIES) + $(LINK) $(libgslsum_la_LDFLAGS) $(libgslsum_la_OBJECTS) $(libgslsum_la_LIBADD) $(LIBS) + +clean-checkPROGRAMS: + @list='$(check_PROGRAMS)'; for p in $$list; do \ + f=`echo $$p|sed 's/$(EXEEXT)$$//'`; \ + echo " rm -f $$p $$f"; \ + rm -f $$p $$f ; \ + done +test$(EXEEXT): $(test_OBJECTS) $(test_DEPENDENCIES) + @rm -f test$(EXEEXT) + $(LINK) $(test_LDFLAGS) $(test_OBJECTS) $(test_LDADD) $(LIBS) + +mostlyclean-compile: + -rm -f *.$(OBJEXT) + +distclean-compile: + -rm -f *.tab.c + +.c.o: + $(COMPILE) -c $< + +.c.obj: + $(COMPILE) -c `$(CYGPATH_W) '$<'` + +.c.lo: + $(LTCOMPILE) -c -o $@ $< + +mostlyclean-libtool: + -rm -f *.lo + +clean-libtool: + -rm -rf .libs _libs + +distclean-libtool: + -rm -f libtool +uninstall-info-am: +install-pkgincludeHEADERS: $(pkginclude_HEADERS) + @$(NORMAL_INSTALL) + test -z "$(pkgincludedir)" || $(mkdir_p) "$(DESTDIR)$(pkgincludedir)" + @list='$(pkginclude_HEADERS)'; for p in $$list; do \ + if test -f "$$p"; then d=; else d="$(srcdir)/"; fi; \ + f=$(am__strip_dir) \ + echo " $(pkgincludeHEADERS_INSTALL) '$$d$$p' '$(DESTDIR)$(pkgincludedir)/$$f'"; \ + $(pkgincludeHEADERS_INSTALL) "$$d$$p" "$(DESTDIR)$(pkgincludedir)/$$f"; \ + done + +uninstall-pkgincludeHEADERS: + @$(NORMAL_UNINSTALL) + @list='$(pkginclude_HEADERS)'; for p in $$list; do \ + f=$(am__strip_dir) \ + echo " rm -f '$(DESTDIR)$(pkgincludedir)/$$f'"; \ + rm -f "$(DESTDIR)$(pkgincludedir)/$$f"; \ + done + +ID: $(HEADERS) $(SOURCES) $(LISP) $(TAGS_FILES) + list='$(SOURCES) $(HEADERS) $(LISP) $(TAGS_FILES)'; \ + unique=`for i in $$list; do \ + if test -f "$$i"; then echo $$i; else echo $(srcdir)/$$i; fi; \ + done | \ + $(AWK) ' { files[$$0] = 1; } \ + END { for (i in files) print i; }'`; \ + mkid -fID $$unique +tags: TAGS + +TAGS: $(HEADERS) $(SOURCES) $(TAGS_DEPENDENCIES) \ + $(TAGS_FILES) $(LISP) + tags=; \ + here=`pwd`; \ + list='$(SOURCES) $(HEADERS) $(LISP) $(TAGS_FILES)'; \ + unique=`for i in $$list; do \ + if test -f "$$i"; then echo $$i; else echo $(srcdir)/$$i; fi; \ + done | \ + $(AWK) ' { files[$$0] = 1; } \ + END { for (i in files) print i; }'`; \ + if test -z "$(ETAGS_ARGS)$$tags$$unique"; then :; else \ + test -n "$$unique" || unique=$$empty_fix; \ + $(ETAGS) $(ETAGSFLAGS) $(AM_ETAGSFLAGS) $(ETAGS_ARGS) \ + $$tags $$unique; \ + fi +ctags: CTAGS +CTAGS: $(HEADERS) $(SOURCES) $(TAGS_DEPENDENCIES) \ + $(TAGS_FILES) $(LISP) + tags=; \ + here=`pwd`; \ + list='$(SOURCES) $(HEADERS) $(LISP) $(TAGS_FILES)'; \ + unique=`for i in $$list; do \ + if test -f "$$i"; then echo $$i; else echo $(srcdir)/$$i; fi; \ + done | \ + $(AWK) ' { files[$$0] = 1; } \ + END { for (i in files) print i; }'`; \ + test -z "$(CTAGS_ARGS)$$tags$$unique" \ + || $(CTAGS) $(CTAGSFLAGS) $(AM_CTAGSFLAGS) $(CTAGS_ARGS) \ + $$tags $$unique + +GTAGS: + here=`$(am__cd) $(top_builddir) && pwd` \ + && cd $(top_srcdir) \ + && gtags -i $(GTAGS_ARGS) $$here + +distclean-tags: + -rm -f TAGS ID GTAGS GRTAGS GSYMS GPATH tags + +check-TESTS: $(TESTS) + @failed=0; all=0; xfail=0; xpass=0; skip=0; \ + srcdir=$(srcdir); export srcdir; \ + list='$(TESTS)'; \ + if test -n "$$list"; then \ + for tst in $$list; do \ + if test -f ./$$tst; then dir=./; \ + elif test -f $$tst; then dir=; \ + else dir="$(srcdir)/"; fi; \ + if $(TESTS_ENVIRONMENT) $${dir}$$tst; then \ + all=`expr $$all + 1`; \ + case " $(XFAIL_TESTS) " in \ + *" $$tst "*) \ + xpass=`expr $$xpass + 1`; \ + failed=`expr $$failed + 1`; \ + echo "XPASS: $$tst"; \ + ;; \ + *) \ + echo "PASS: $$tst"; \ + ;; \ + esac; \ + elif test $$? -ne 77; then \ + all=`expr $$all + 1`; \ + case " $(XFAIL_TESTS) " in \ + *" $$tst "*) \ + xfail=`expr $$xfail + 1`; \ + echo "XFAIL: $$tst"; \ + ;; \ + *) \ + failed=`expr $$failed + 1`; \ + echo "FAIL: $$tst"; \ + ;; \ + esac; \ + else \ + skip=`expr $$skip + 1`; \ + echo "SKIP: $$tst"; \ + fi; \ + done; \ + if test "$$failed" -eq 0; then \ + if test "$$xfail" -eq 0; then \ + banner="All $$all tests passed"; \ + else \ + banner="All $$all tests behaved as expected ($$xfail expected failures)"; \ + fi; \ + else \ + if test "$$xpass" -eq 0; then \ + banner="$$failed of $$all tests failed"; \ + else \ + banner="$$failed of $$all tests did not behave as expected ($$xpass unexpected passes)"; \ + fi; \ + fi; \ + dashes="$$banner"; \ + skipped=""; \ + if test "$$skip" -ne 0; then \ + skipped="($$skip tests were not run)"; \ + test `echo "$$skipped" | wc -c` -le `echo "$$banner" | wc -c` || \ + dashes="$$skipped"; \ + fi; \ + report=""; \ + if test "$$failed" -ne 0 && test -n "$(PACKAGE_BUGREPORT)"; then \ + report="Please report to $(PACKAGE_BUGREPORT)"; \ + test `echo "$$report" | wc -c` -le `echo "$$banner" | wc -c` || \ + dashes="$$report"; \ + fi; \ + dashes=`echo "$$dashes" | sed s/./=/g`; \ + echo "$$dashes"; \ + echo "$$banner"; \ + test -z "$$skipped" || echo "$$skipped"; \ + test -z "$$report" || echo "$$report"; \ + echo "$$dashes"; \ + test "$$failed" -eq 0; \ + else :; fi + +distdir: $(DISTFILES) + @srcdirstrip=`echo "$(srcdir)" | sed 's|.|.|g'`; \ + topsrcdirstrip=`echo "$(top_srcdir)" | sed 's|.|.|g'`; \ + list='$(DISTFILES)'; for file in $$list; do \ + case $$file in \ + $(srcdir)/*) file=`echo "$$file" | sed "s|^$$srcdirstrip/||"`;; \ + $(top_srcdir)/*) file=`echo "$$file" | sed "s|^$$topsrcdirstrip/|$(top_builddir)/|"`;; \ + esac; \ + if test -f $$file || test -d $$file; then d=.; else d=$(srcdir); fi; \ + dir=`echo "$$file" | sed -e 's,/[^/]*$$,,'`; \ + if test "$$dir" != "$$file" && test "$$dir" != "."; then \ + dir="/$$dir"; \ + $(mkdir_p) "$(distdir)$$dir"; \ + else \ + dir=''; \ + fi; \ + if test -d $$d/$$file; then \ + if test -d $(srcdir)/$$file && test $$d != $(srcdir); then \ + cp -pR $(srcdir)/$$file $(distdir)$$dir || exit 1; \ + fi; \ + cp -pR $$d/$$file $(distdir)$$dir || exit 1; \ + else \ + test -f $(distdir)/$$file \ + || cp -p $$d/$$file $(distdir)/$$file \ + || exit 1; \ + fi; \ + done +check-am: all-am + $(MAKE) $(AM_MAKEFLAGS) $(check_PROGRAMS) + $(MAKE) $(AM_MAKEFLAGS) check-TESTS +check: check-am +all-am: Makefile $(LTLIBRARIES) $(HEADERS) +installdirs: + for dir in "$(DESTDIR)$(pkgincludedir)"; do \ + test -z "$$dir" || $(mkdir_p) "$$dir"; \ + done +install: install-am +install-exec: install-exec-am +install-data: install-data-am +uninstall: uninstall-am + +install-am: all-am + @$(MAKE) $(AM_MAKEFLAGS) install-exec-am install-data-am + +installcheck: installcheck-am +install-strip: + $(MAKE) $(AM_MAKEFLAGS) INSTALL_PROGRAM="$(INSTALL_STRIP_PROGRAM)" \ + install_sh_PROGRAM="$(INSTALL_STRIP_PROGRAM)" INSTALL_STRIP_FLAG=-s \ + `test -z '$(STRIP)' || \ + echo "INSTALL_PROGRAM_ENV=STRIPPROG='$(STRIP)'"` install +mostlyclean-generic: + +clean-generic: + +distclean-generic: + -test -z "$(CONFIG_CLEAN_FILES)" || rm -f $(CONFIG_CLEAN_FILES) + +maintainer-clean-generic: + @echo "This command is intended for maintainers to use" + @echo "it deletes files that may require special tools to rebuild." +clean: clean-am + +clean-am: clean-checkPROGRAMS clean-generic clean-libtool \ + clean-noinstLTLIBRARIES mostlyclean-am + +distclean: distclean-am + -rm -f Makefile +distclean-am: clean-am distclean-compile distclean-generic \ + distclean-libtool distclean-tags + +dvi: dvi-am + +dvi-am: + +html: html-am + +info: info-am + +info-am: + +install-data-am: install-pkgincludeHEADERS + +install-exec-am: + +install-info: install-info-am + +install-man: + +installcheck-am: + +maintainer-clean: maintainer-clean-am + -rm -f Makefile +maintainer-clean-am: distclean-am maintainer-clean-generic + +mostlyclean: mostlyclean-am + +mostlyclean-am: mostlyclean-compile mostlyclean-generic \ + mostlyclean-libtool + +pdf: pdf-am + +pdf-am: + +ps: ps-am + +ps-am: + +uninstall-am: uninstall-info-am uninstall-pkgincludeHEADERS + +.PHONY: CTAGS GTAGS all all-am check check-TESTS check-am clean \ + clean-checkPROGRAMS clean-generic clean-libtool \ + clean-noinstLTLIBRARIES ctags distclean distclean-compile \ + distclean-generic distclean-libtool distclean-tags distdir dvi \ + dvi-am html html-am info info-am install install-am \ + install-data install-data-am install-exec install-exec-am \ + install-info install-info-am install-man \ + install-pkgincludeHEADERS install-strip installcheck \ + installcheck-am installdirs maintainer-clean \ + maintainer-clean-generic mostlyclean mostlyclean-compile \ + mostlyclean-generic mostlyclean-libtool pdf pdf-am ps ps-am \ + tags uninstall uninstall-am uninstall-info-am \ + uninstall-pkgincludeHEADERS + +# Tell versions [3.59,3.63) of GNU make to not export all variables. +# Otherwise a system limit (for SysV at least) may be exceeded. +.NOEXPORT: diff --git a/gsl-1.9/sum/gsl_sum.h b/gsl-1.9/sum/gsl_sum.h new file mode 100644 index 0000000..4f7a66e --- /dev/null +++ b/gsl-1.9/sum/gsl_sum.h @@ -0,0 +1,162 @@ +/* sum/gsl_sum.h + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman, Brian Gough + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or (at + * your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +/* Author: G. Jungman */ + + +#ifndef __GSL_SUM_H__ +#define __GSL_SUM_H__ + +#include <stdlib.h> + +#undef __BEGIN_DECLS +#undef __END_DECLS +#ifdef __cplusplus +# define __BEGIN_DECLS extern "C" { +# define __END_DECLS } +#else +# define __BEGIN_DECLS /* empty */ +# define __END_DECLS /* empty */ +#endif + +__BEGIN_DECLS + +/* Workspace for Levin U Transform with error estimation, + * + * size = number of terms the workspace can handle + * sum_plain = simple sum of series + * q_num = backward diagonal of numerator; length = size + * q_den = backward diagonal of denominator; length = size + * dq_num = table of numerator derivatives; length = size**2 + * dq_den = table of denominator derivatives; length = size**2 + * dsum = derivative of sum wrt term i; length = size + */ + +typedef struct +{ + size_t size; + size_t i; /* position in array */ + size_t terms_used; /* number of calls */ + double sum_plain; + double *q_num; + double *q_den; + double *dq_num; + double *dq_den; + double *dsum; +} +gsl_sum_levin_u_workspace; + +gsl_sum_levin_u_workspace *gsl_sum_levin_u_alloc (size_t n); +void gsl_sum_levin_u_free (gsl_sum_levin_u_workspace * w); + +/* Basic Levin-u acceleration method. + * + * array = array of series elements + * n = size of array + * sum_accel = result of summation acceleration + * err = estimated error + * + * See [Fessler et al., ACM TOMS 9, 346 (1983) and TOMS-602] + */ + +int gsl_sum_levin_u_accel (const double *array, + const size_t n, + gsl_sum_levin_u_workspace * w, + double *sum_accel, double *abserr); + +/* Basic Levin-u acceleration method with constraints on the terms + * used, + * + * array = array of series elements + * n = size of array + * min_terms = minimum number of terms to sum + * max_terms = maximum number of terms to sum + * sum_accel = result of summation acceleration + * err = estimated error + * + * See [Fessler et al., ACM TOMS 9, 346 (1983) and TOMS-602] + */ + +int gsl_sum_levin_u_minmax (const double *array, + const size_t n, + const size_t min_terms, + const size_t max_terms, + gsl_sum_levin_u_workspace * w, + double *sum_accel, double *abserr); + +/* Basic Levin-u step w/o reference to the array of terms. + * We only need to specify the value of the current term + * to execute the step. See TOMS-745. + * + * sum = t0 + ... + t_{n-1} + term; term = t_{n} + * + * term = value of the series term to be added + * n = position of term in series (starting from 0) + * sum_accel = result of summation acceleration + * sum_plain = simple sum of series + */ + +int +gsl_sum_levin_u_step (const double term, + const size_t n, + const size_t nmax, + gsl_sum_levin_u_workspace * w, + double *sum_accel); + +/* The following functions perform the same calculation without + estimating the errors. They require O(N) storage instead of O(N^2). + This may be useful for summing many similar series where the size + of the error has already been estimated reliably and is not + expected to change. */ + +typedef struct +{ + size_t size; + size_t i; /* position in array */ + size_t terms_used; /* number of calls */ + double sum_plain; + double *q_num; + double *q_den; + double *dsum; +} +gsl_sum_levin_utrunc_workspace; + +gsl_sum_levin_utrunc_workspace *gsl_sum_levin_utrunc_alloc (size_t n); +void gsl_sum_levin_utrunc_free (gsl_sum_levin_utrunc_workspace * w); + +int gsl_sum_levin_utrunc_accel (const double *array, + const size_t n, + gsl_sum_levin_utrunc_workspace * w, + double *sum_accel, double *abserr_trunc); + +int gsl_sum_levin_utrunc_minmax (const double *array, + const size_t n, + const size_t min_terms, + const size_t max_terms, + gsl_sum_levin_utrunc_workspace * w, + double *sum_accel, double *abserr_trunc); + +int gsl_sum_levin_utrunc_step (const double term, + const size_t n, + gsl_sum_levin_utrunc_workspace * w, + double *sum_accel); + +__END_DECLS + +#endif /* __GSL_SUM_H__ */ diff --git a/gsl-1.9/sum/levin_u.c b/gsl-1.9/sum/levin_u.c new file mode 100644 index 0000000..06f188a --- /dev/null +++ b/gsl-1.9/sum/levin_u.c @@ -0,0 +1,255 @@ +/* sum/levin_u.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman, Brian Gough + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or (at + * your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +#include <config.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_test.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_sum.h> + +int +gsl_sum_levin_u_accel (const double *array, const size_t array_size, + gsl_sum_levin_u_workspace * w, + double *sum_accel, double *abserr) +{ + return gsl_sum_levin_u_minmax (array, array_size, + 0, array_size - 1, w, sum_accel, abserr); +} + +int +gsl_sum_levin_u_minmax (const double *array, const size_t array_size, + const size_t min_terms, const size_t max_terms, + gsl_sum_levin_u_workspace * w, + double *sum_accel, double *abserr) +{ + if (array_size == 0) + { + *sum_accel = 0.0; + *abserr = 0.0; + w->sum_plain = 0.0; + w->terms_used = 0; + return GSL_SUCCESS; + } + else if (array_size == 1) + { + *sum_accel = array[0]; + *abserr = 0.0; + w->sum_plain = array[0]; + w->terms_used = 1; + return GSL_SUCCESS; + } + else + { + const double SMALL = 0.01; + const size_t nmax = GSL_MAX (max_terms, array_size) - 1; + double noise_n = 0.0, noise_nm1 = 0.0; + double trunc_n = 0.0, trunc_nm1 = 0.0; + double actual_trunc_n = 0.0, actual_trunc_nm1 = 0.0; + double result_n = 0.0, result_nm1 = 0.0; + double variance = 0; + size_t n; + unsigned int i; + int better = 0; + int before = 0; + int converging = 0; + double least_trunc = GSL_DBL_MAX; + double least_trunc_noise = GSL_DBL_MAX; + double least_trunc_result; + + /* Calculate specified minimum number of terms. No convergence + tests are made, and no truncation information is stored. */ + + for (n = 0; n < min_terms; n++) + { + const double t = array[n]; + result_nm1 = result_n; + gsl_sum_levin_u_step (t, n, nmax, w, &result_n); + } + + least_trunc_result = result_n; + + variance = 0; + for (i = 0; i < n; i++) + { + double dn = w->dsum[i] * GSL_MACH_EPS * array[i]; + variance += dn * dn; + } + noise_n = sqrt (variance); + + /* Calculate up to maximum number of terms. Check truncation + condition. */ + + for (; n <= nmax; n++) + { + const double t = array[n]; + + result_nm1 = result_n; + gsl_sum_levin_u_step (t, n, nmax, w, &result_n); + + /* Compute the truncation error directly */ + + actual_trunc_nm1 = actual_trunc_n; + actual_trunc_n = fabs (result_n - result_nm1); + + /* Average results to make a more reliable estimate of the + real truncation error */ + + trunc_nm1 = trunc_n; + trunc_n = 0.5 * (actual_trunc_n + actual_trunc_nm1); + + noise_nm1 = noise_n; + variance = 0; + + for (i = 0; i <= n; i++) + { + double dn = w->dsum[i] * GSL_MACH_EPS * array[i]; + variance += dn * dn; + } + + noise_n = sqrt (variance); + + /* Determine if we are in the convergence region. */ + + better = (trunc_n < trunc_nm1 || trunc_n < SMALL * fabs (result_n)); + converging = converging || (better && before); + before = better; + + if (converging) + { + if (trunc_n < least_trunc) + { + /* Found a low truncation point in the convergence + region. Save it. */ + + least_trunc_result = result_n; + least_trunc = trunc_n; + least_trunc_noise = noise_n; + } + + if (noise_n > trunc_n / 3.0) + break; + + if (trunc_n < 10.0 * GSL_MACH_EPS * fabs (result_n)) + break; + } + + } + + if (converging) + { + /* Stopped in the convergence region. Return result and + error estimate. */ + + *sum_accel = least_trunc_result; + *abserr = GSL_MAX_DBL (least_trunc, least_trunc_noise); + w->terms_used = n; + return GSL_SUCCESS; + } + else + { + /* Never reached the convergence region. Use the last + calculated values. */ + + *sum_accel = result_n; + *abserr = GSL_MAX_DBL (trunc_n, noise_n); + w->terms_used = n; + return GSL_SUCCESS; + } + } +} + + +int +gsl_sum_levin_u_step (const double term, const size_t n, const size_t nmax, + gsl_sum_levin_u_workspace * w, double *sum_accel) +{ + +#define I(i,j) ((i)*(nmax+1) + (j)) + + if (n == 0) + { + *sum_accel = term; + w->sum_plain = term; + + w->q_den[0] = 1.0 / term; + w->q_num[0] = 1.0; + + w->dq_den[I (0, 0)] = -1.0 / (term * term); + w->dq_num[I (0, 0)] = 0.0; + + w->dsum[0] = 1.0; + + return GSL_SUCCESS; + } + else + { + double result; + double factor = 1.0; + double ratio = (double) n / (n + 1.0); + unsigned int i; + int j; + + w->sum_plain += term; + + w->q_den[n] = 1.0 / (term * (n + 1.0) * (n + 1.0)); + w->q_num[n] = w->sum_plain * w->q_den[n]; + + for (i = 0; i < n; i++) + { + w->dq_den[I (i, n)] = 0; + w->dq_num[I (i, n)] = w->q_den[n]; + } + + w->dq_den[I (n, n)] = -w->q_den[n] / term; + w->dq_num[I (n, n)] = + w->q_den[n] + w->sum_plain * (w->dq_den[I (n, n)]); + + for (j = n - 1; j >= 0; j--) + { + double c = factor * (j + 1) / (n + 1); + factor *= ratio; + w->q_den[j] = w->q_den[j + 1] - c * w->q_den[j]; + w->q_num[j] = w->q_num[j + 1] - c * w->q_num[j]; + + for (i = 0; i < n; i++) + { + w->dq_den[I (i, j)] = + w->dq_den[I (i, j + 1)] - c * w->dq_den[I (i, j)]; + w->dq_num[I (i, j)] = + w->dq_num[I (i, j + 1)] - c * w->dq_num[I (i, j)]; + } + + w->dq_den[I (n, j)] = w->dq_den[I (n, j + 1)]; + w->dq_num[I (n, j)] = w->dq_num[I (n, j + 1)]; + } + + result = w->q_num[0] / w->q_den[0]; + + *sum_accel = result; + + for (i = 0; i <= n; i++) + { + w->dsum[i] = + (w->dq_num[I (i, 0)] - + result * w->dq_den[I (i, 0)]) / w->q_den[0]; + } + + return GSL_SUCCESS; + } +} diff --git a/gsl-1.9/sum/levin_utrunc.c b/gsl-1.9/sum/levin_utrunc.c new file mode 100644 index 0000000..16a3f16 --- /dev/null +++ b/gsl-1.9/sum/levin_utrunc.c @@ -0,0 +1,202 @@ +/* sum/levin_utrunc.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman, Brian Gough + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or (at + * your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +/* Author: G. Jungman */ + +#include <config.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_test.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_sum.h> + +int +gsl_sum_levin_utrunc_accel (const double *array, + const size_t array_size, + gsl_sum_levin_utrunc_workspace * w, + double *sum_accel, double *abserr_trunc) +{ + return gsl_sum_levin_utrunc_minmax (array, array_size, + 0, array_size - 1, + w, sum_accel, abserr_trunc); +} + + +int +gsl_sum_levin_utrunc_minmax (const double *array, + const size_t array_size, + const size_t min_terms, + const size_t max_terms, + gsl_sum_levin_utrunc_workspace * w, + double *sum_accel, double *abserr_trunc) +{ + if (array_size == 0) + { + *sum_accel = 0.0; + *abserr_trunc = 0.0; + w->sum_plain = 0.0; + w->terms_used = 0; + return GSL_SUCCESS; + } + else if (array_size == 1) + { + *sum_accel = array[0]; + *abserr_trunc = GSL_POSINF; + w->sum_plain = array[0]; + w->terms_used = 1; + return GSL_SUCCESS; + } + else + { + const double SMALL = 0.01; + const size_t nmax = GSL_MAX (max_terms, array_size) - 1; + double trunc_n = 0.0, trunc_nm1 = 0.0; + double actual_trunc_n = 0.0, actual_trunc_nm1 = 0.0; + double result_n = 0.0, result_nm1 = 0.0; + size_t n; + int better = 0; + int before = 0; + int converging = 0; + double least_trunc = GSL_DBL_MAX; + double result_least_trunc; + + /* Calculate specified minimum number of terms. No convergence + tests are made, and no truncation information is stored. */ + + for (n = 0; n < min_terms; n++) + { + const double t = array[n]; + + result_nm1 = result_n; + gsl_sum_levin_utrunc_step (t, n, w, &result_n); + } + + /* Assume the result after the minimum calculation is the best. */ + + result_least_trunc = result_n; + + /* Calculate up to maximum number of terms. Check truncation + condition. */ + + for (; n <= nmax; n++) + { + const double t = array[n]; + + result_nm1 = result_n; + gsl_sum_levin_utrunc_step (t, n, w, &result_n); + + /* Compute the truncation error directly */ + + actual_trunc_nm1 = actual_trunc_n; + actual_trunc_n = fabs (result_n - result_nm1); + + /* Average results to make a more reliable estimate of the + real truncation error */ + + trunc_nm1 = trunc_n; + trunc_n = 0.5 * (actual_trunc_n + actual_trunc_nm1); + + /* Determine if we are in the convergence region. */ + + better = (trunc_n < trunc_nm1 || trunc_n < SMALL * fabs (result_n)); + converging = converging || (better && before); + before = better; + + if (converging) + { + if (trunc_n < least_trunc) + { + /* Found a low truncation point in the convergence + region. Save it. */ + + least_trunc = trunc_n; + result_least_trunc = result_n; + } + + if (fabs (trunc_n / result_n) < 10.0 * GSL_MACH_EPS) + break; + } + } + + if (converging) + { + /* Stopped in the convergence region. Return result and + error estimate. */ + + *sum_accel = result_least_trunc; + *abserr_trunc = least_trunc; + w->terms_used = n; + return GSL_SUCCESS; + } + else + { + /* Never reached the convergence region. Use the last + calculated values. */ + + *sum_accel = result_n; + *abserr_trunc = trunc_n; + w->terms_used = n; + return GSL_SUCCESS; + } + } +} + +int +gsl_sum_levin_utrunc_step (const double term, + const size_t n, + gsl_sum_levin_utrunc_workspace * w, double *sum_accel) +{ + if (term == 0.0) + { + /* This is actually harmless when treated in this way. A term + which is exactly zero is simply ignored; the state is not + changed. We return GSL_EZERODIV as an indicator that this + occured. */ + + return GSL_EZERODIV; + } + else if (n == 0) + { + *sum_accel = term; + w->sum_plain = term; + w->q_den[0] = 1.0 / term; + w->q_num[0] = 1.0; + return GSL_SUCCESS; + } + else + { + double factor = 1.0; + double ratio = (double) n / (n + 1.0); + int j; + + w->sum_plain += term; + w->q_den[n] = 1.0 / (term * (n + 1.0) * (n + 1.0)); + w->q_num[n] = w->sum_plain * w->q_den[n]; + + for (j = n - 1; j >= 0; j--) + { + double c = factor * (j + 1) / (n + 1); + factor *= ratio; + w->q_den[j] = w->q_den[j + 1] - c * w->q_den[j]; + w->q_num[j] = w->q_num[j + 1] - c * w->q_num[j]; + } + + *sum_accel = w->q_num[0] / w->q_den[0]; + return GSL_SUCCESS; + } +} diff --git a/gsl-1.9/sum/test.c b/gsl-1.9/sum/test.c new file mode 100644 index 0000000..4ca4f1f --- /dev/null +++ b/gsl-1.9/sum/test.c @@ -0,0 +1,234 @@ +/* sum/test.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman, Brian Gough + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or (at + * your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +/* Author: G. Jungman */ + +#include <config.h> +#include <stdlib.h> +#include <stdio.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_test.h> +#include <gsl/gsl_sum.h> + +#include <gsl/gsl_ieee_utils.h> + +#define N 50 + +void check_trunc (double * t, double expected, const char * desc); +void check_full (double * t, double expected, const char * desc); + +int +main (void) +{ + gsl_ieee_env_setup (); + + { + double t[N]; + int n; + + const double zeta_2 = M_PI * M_PI / 6.0; + + /* terms for zeta(2) */ + + for (n = 0; n < N; n++) + { + double np1 = n + 1.0; + t[n] = 1.0 / (np1 * np1); + } + + check_trunc (t, zeta_2, "zeta(2)"); + check_full (t, zeta_2, "zeta(2)"); + } + + { + double t[N]; + double x, y; + int n; + + /* terms for exp(10.0) */ + x = 10.0; + y = exp(x); + + t[0] = 1.0; + for (n = 1; n < N; n++) + { + t[n] = t[n - 1] * (x / n); + } + + check_trunc (t, y, "exp(10)"); + check_full (t, y, "exp(10)"); + } + + { + double t[N]; + double x, y; + int n; + + /* terms for exp(-10.0) */ + x = -10.0; + y = exp(x); + + t[0] = 1.0; + for (n = 1; n < N; n++) + { + t[n] = t[n - 1] * (x / n); + } + + check_trunc (t, y, "exp(-10)"); + check_full (t, y, "exp(-10)"); + } + + { + double t[N]; + double x, y; + int n; + + /* terms for -log(1-x) */ + x = 0.5; + y = -log(1-x); + t[0] = x; + for (n = 1; n < N; n++) + { + t[n] = t[n - 1] * (x * n) / (n + 1.0); + } + + check_trunc (t, y, "-log(1/2)"); + check_full (t, y, "-log(1/2)"); + } + + { + double t[N]; + double x, y; + int n; + + /* terms for -log(1-x) */ + x = -1.0; + y = -log(1-x); + t[0] = x; + for (n = 1; n < N; n++) + { + t[n] = t[n - 1] * (x * n) / (n + 1.0); + } + + check_trunc (t, y, "-log(2)"); + check_full (t, y, "-log(2)"); + } + + { + double t[N]; + int n; + + double result = 0.192594048773; + + /* terms for an alternating asymptotic series */ + + t[0] = 3.0 / (M_PI * M_PI); + + for (n = 1; n < N; n++) + { + t[n] = -t[n - 1] * (4.0 * (n + 1.0) - 1.0) / (M_PI * M_PI); + } + + check_trunc (t, result, "asymptotic series"); + check_full (t, result, "asymptotic series"); + } + + { + double t[N]; + int n; + + /* Euler's gamma from GNU Calc (precision = 32) */ + + double result = 0.5772156649015328606065120900824; + + /* terms for Euler's gamma */ + + t[0] = 1.0; + + for (n = 1; n < N; n++) + { + t[n] = 1/(n+1.0) + log(n/(n+1.0)); + } + + check_trunc (t, result, "Euler's constant"); + check_full (t, result, "Euler's constant"); + } + + { + double t[N]; + int n; + + /* eta(1/2) = sum_{k=1}^{\infty} (-1)^(k+1) / sqrt(k) + + From Levin, Intern. J. Computer Math. B3:371--388, 1973. + + I=(1-sqrt(2))zeta(1/2) + =(2/sqrt(pi))*integ(1/(exp(x^2)+1),x,0,inf) */ + + double result = 0.6048986434216305; /* approx */ + + /* terms for eta(1/2) */ + + for (n = 0; n < N; n++) + { + t[n] = (n%2 ? -1 : 1) * 1.0 /sqrt(n + 1.0); + } + + check_trunc (t, result, "eta(1/2)"); + check_full (t, result, "eta(1/2)"); + } + + exit (gsl_test_summary ()); +} + +void +check_trunc (double * t, double expected, const char * desc) +{ + double sum_accel, prec; + + gsl_sum_levin_utrunc_workspace * w = gsl_sum_levin_utrunc_alloc (N); + + gsl_sum_levin_utrunc_accel (t, N, w, &sum_accel, &prec); + gsl_test_rel (sum_accel, expected, 1e-8, "trunc result, %s", desc); + + /* No need to check precision for truncated result since this is not + a meaningful number */ + + gsl_sum_levin_utrunc_free (w); +} + +void +check_full (double * t, double expected, const char * desc) +{ + double sum_accel, err_est, sd_actual, sd_est; + + gsl_sum_levin_u_workspace * w = gsl_sum_levin_u_alloc (N); + + gsl_sum_levin_u_accel (t, N, w, &sum_accel, &err_est); + gsl_test_rel (sum_accel, expected, 1e-8, "full result, %s", desc); + + sd_est = -log10 (err_est/fabs(sum_accel)); + sd_actual = -log10 (DBL_EPSILON + fabs ((sum_accel - expected)/expected)); + + /* Allow one digit of slop */ + + gsl_test (sd_est > sd_actual + 1.0, "full significant digits, %s (%g vs %g)", desc, sd_est, sd_actual); + + gsl_sum_levin_u_free (w); +} diff --git a/gsl-1.9/sum/work_u.c b/gsl-1.9/sum/work_u.c new file mode 100644 index 0000000..3f7052b --- /dev/null +++ b/gsl-1.9/sum/work_u.c @@ -0,0 +1,94 @@ +#include <config.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_sum.h> + +gsl_sum_levin_u_workspace * +gsl_sum_levin_u_alloc (size_t n) +{ + gsl_sum_levin_u_workspace * w; + + if (n == 0) + { + GSL_ERROR_VAL ("length n must be positive integer", GSL_EDOM, 0); + } + + w = (gsl_sum_levin_u_workspace *) malloc(sizeof(gsl_sum_levin_u_workspace)); + + if (w == NULL) + { + GSL_ERROR_VAL ("failed to allocate struct", GSL_ENOMEM, 0); + } + + w->q_num = (double *) malloc (n * sizeof (double)); + + if (w->q_num == NULL) + { + free(w) ; /* error in constructor, prevent memory leak */ + + GSL_ERROR_VAL ("failed to allocate space for q_num", GSL_ENOMEM, 0); + } + + w->q_den = (double *) malloc (n * sizeof (double)); + + if (w->q_den == NULL) + { + free (w->q_num); + free (w) ; /* error in constructor, prevent memory leak */ + + GSL_ERROR_VAL ("failed to allocate space for q_den", GSL_ENOMEM, 0); + } + + w->dq_num = (double *) malloc (n * n * sizeof (double)); + + if (w->dq_num == NULL) + { + free (w->q_den); + free (w->q_num); + free(w) ; /* error in constructor, prevent memory leak */ + + GSL_ERROR_VAL ("failed to allocate space for dq_num", GSL_ENOMEM, 0); + } + + w->dq_den = (double *) malloc (n * n * sizeof (double)); + + if (w->dq_den == NULL) + { + free (w->dq_num); + free (w->q_den); + free (w->q_num); + free (w) ; /* error in constructor, prevent memory leak */ + + GSL_ERROR_VAL ("failed to allocate space for dq_den", GSL_ENOMEM, 0); + } + + w->dsum = (double *) malloc (n * sizeof (double)); + + if (w->dsum == NULL) + { + free (w->dq_den); + free (w->dq_num); + free (w->q_den); + free (w->q_num); + free (w) ; /* error in constructor, prevent memory leak */ + + GSL_ERROR_VAL ("failed to allocate space for dsum", GSL_ENOMEM, 0); + } + + w->size = n; + w->terms_used = 0; + w->sum_plain = 0; + + return w; +} + +void +gsl_sum_levin_u_free (gsl_sum_levin_u_workspace * w) +{ + free (w->dsum); + free (w->dq_den); + free (w->dq_num); + free (w->q_den); + free (w->q_num); + free (w); +} diff --git a/gsl-1.9/sum/work_utrunc.c b/gsl-1.9/sum/work_utrunc.c new file mode 100644 index 0000000..67ebb81 --- /dev/null +++ b/gsl-1.9/sum/work_utrunc.c @@ -0,0 +1,67 @@ +#include <config.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_sum.h> + +gsl_sum_levin_utrunc_workspace * +gsl_sum_levin_utrunc_alloc (size_t n) +{ + gsl_sum_levin_utrunc_workspace * w; + + if (n == 0) + { + GSL_ERROR_VAL ("length n must be positive integer", GSL_EDOM, 0); + } + + w = (gsl_sum_levin_utrunc_workspace *) malloc(sizeof(gsl_sum_levin_utrunc_workspace)); + + if (w == NULL) + { + GSL_ERROR_VAL ("failed to allocate struct", GSL_ENOMEM, 0); + } + + w->q_num = (double *) malloc (n * sizeof (double)); + + if (w->q_num == NULL) + { + free(w) ; /* error in constructor, prevent memory leak */ + + GSL_ERROR_VAL ("failed to allocate space for q_num", GSL_ENOMEM, 0); + } + + w->q_den = (double *) malloc (n * sizeof (double)); + + if (w->q_den == NULL) + { + free (w->q_num); + free (w) ; /* error in constructor, prevent memory leak */ + + GSL_ERROR_VAL ("failed to allocate space for q_den", GSL_ENOMEM, 0); + } + + w->dsum = (double *) malloc (n * sizeof (double)); + + if (w->dsum == NULL) + { + free (w->q_den); + free (w->q_num); + free (w) ; /* error in constructor, prevent memory leak */ + + GSL_ERROR_VAL ("failed to allocate space for dsum", GSL_ENOMEM, 0); + } + + w->size = n; + w->terms_used = 0; + w->sum_plain = 0; + + return w; +} + +void +gsl_sum_levin_utrunc_free (gsl_sum_levin_utrunc_workspace * w) +{ + free (w->dsum); + free (w->q_den); + free (w->q_num); + free (w); +} |