diff options
Diffstat (limited to 'gsl-1.9/roots')
| -rw-r--r-- | gsl-1.9/roots/ChangeLog | 159 | ||||
| -rw-r--r-- | gsl-1.9/roots/Makefile.am | 19 | ||||
| -rw-r--r-- | gsl-1.9/roots/Makefile.in | 548 | ||||
| -rw-r--r-- | gsl-1.9/roots/TODO | 9 | ||||
| -rw-r--r-- | gsl-1.9/roots/bisection.c | 134 | ||||
| -rw-r--r-- | gsl-1.9/roots/brent.c | 244 | ||||
| -rw-r--r-- | gsl-1.9/roots/convergence.c | 87 | ||||
| -rw-r--r-- | gsl-1.9/roots/falsepos.c | 178 | ||||
| -rw-r--r-- | gsl-1.9/roots/fdfsolver.c | 88 | ||||
| -rw-r--r-- | gsl-1.9/roots/fsolver.c | 107 | ||||
| -rw-r--r-- | gsl-1.9/roots/gsl_roots.h | 127 | ||||
| -rw-r--r-- | gsl-1.9/roots/newton.c | 106 | ||||
| -rw-r--r-- | gsl-1.9/roots/roots.h | 37 | ||||
| -rw-r--r-- | gsl-1.9/roots/secant.c | 117 | ||||
| -rw-r--r-- | gsl-1.9/roots/steffenson.c | 144 | ||||
| -rw-r--r-- | gsl-1.9/roots/test.c | 315 | ||||
| -rw-r--r-- | gsl-1.9/roots/test.h | 129 | ||||
| -rw-r--r-- | gsl-1.9/roots/test_funcs.c | 229 |
18 files changed, 2777 insertions, 0 deletions
diff --git a/gsl-1.9/roots/ChangeLog b/gsl-1.9/roots/ChangeLog new file mode 100644 index 0000000..aa00c09 --- /dev/null +++ b/gsl-1.9/roots/ChangeLog @@ -0,0 +1,159 @@ +2007-01-04 Brian Gough <bjg@network-theory.co.uk> + + * convergence.c (gsl_root_test_delta): added termination + alternative condition x1==x0 + +2005-03-02 Brian Gough <bjg@network-theory.co.uk> + + * steffenson.c (steffenson_iterate): improved wording of error messages + + * secant.c (secant_iterate): improved wording of error messages + + * roots.h (SAFE_FUNC_CALL): improved wording of error message + + * newton.c (newton_iterate): improved wording of error messages + + * utility.c: removed, not needed any more + +Sun Jul 15 17:53:48 2001 Brian Gough <bjg@network-theory.co.uk> + + * removed interval type + +Sun May 6 14:26:59 2001 Brian Gough <bjg@network-theory.co.uk> + + * test.c: removed tests for macros, which are now in sys/. + +Mon Apr 16 20:17:04 2001 Brian Gough <bjg@network-theory.co.uk> + + * fsolver.c (gsl_root_fsolver_alloc): removed unnecessary status + variable + +Sun Feb 18 15:35:25 2001 Brian Gough <bjg@network-theory.co.uk> + + * fdfsolver.c fsolver.c: changed so that the solver _alloc + function no longer calls _set, the user must do that separately. + +Wed May 17 11:37:15 2000 Brian Gough <bjg@network-theory.co.uk> + + * test_macros.c (test_macros): use GSL_POSINF and GSL_NAN macros + instead of 1/0 and 0/0 + +Mon Feb 14 13:05:30 2000 Brian Gough <bjg@network-theory.co.uk> + + * removed definition of isinf macro (no longer needed) + + * made all internal functions static + +Wed Nov 3 11:59:35 1999 Brian Gough <bjg@network-theory.co.uk> + + * fixed test failures + + * test.c (main): added a call to gsl_ieee_env_setup for testing + + * test_roots.c: increased the maximum number of iterations to 150 + so that the tests still work on the difficult cases. + + * steffenson.c (steffenson_iterate): add a check to avoid division + by zero + +Sat Oct 16 19:43:14 1999 Brian Gough <bjg@network-theory.co.uk> + + * removed GSL_ROOT_EPSILON_BUFFER, not needed anymore + +Wed Jul 21 18:47:01 1999 Brian Gough <bjg@network-theory.co.uk> + + * gsl_roots.h, convergence.c: changed order of relative and + absolute errors to make them the same as quadpack routines + (abs,rel) + +Wed Jul 21 16:30:56 1999 Brian Gough <bjg@network-theory.co.uk> + + * brent.c (brent_iterate): fixed bug where bounding interval could + be incorrect and not include root. + +Mon Mar 1 15:38:06 1999 Brian Gough <bjg@netsci.freeserve.co.uk> + + * moved static class data out of gsl_root_fsolver and + gsl_root_fdfsolver and into gsl_root_fsolver_type and + gsl_root_fdfsolver_type + + +Mon Mar 1 15:38:06 1999 Brian Gough <bjg@netsci.freeserve.co.uk> + + * renamed f_solver to fsolver and fdf_solver to fdfsolver, since + these look neater + +Sun Feb 28 21:11:21 1999 Brian Gough <bjg@netsci.freeserve.co.uk> + + * rewrote the root finding functions in an iterative framework + +Tue Nov 17 16:47:09 1998 Brian Gough <bjg@vvv.lanl.gov> + + * secant.c, falsepos.c newton.c: added gsl_math.h to included + headers to import GSL_MAX and GSL_MIN + +Mon Nov 9 21:21:45 1998 Brian Gough <bjg@vvv.lanl.gov> + + * roots.h: got rid of local MAX(a,b) and MIN(a,b) definitions + since they are now in config.h + +Wed Nov 4 16:08:32 1998 Brian Gough <bjg@vvv.lanl.gov> + + * test.c (test_brent): allow the brent tests to run for more + iterations since they take longer on the pathological cases. + + * brent.c (gsl_root_brent): on each iteration keep track of + current best estimates of the root and the bounds so that they are + returned to the user if the function exits prematurely. + + clean up the brent algorithm based on remarks in the original + paper + +Mon Oct 26 16:31:21 1998 Brian Gough <bjg@vvv.lanl.gov> + + * in all routines with upper and lower bounds if a root is found + exactly then the bracket is collapsed onto the root instead of + being untouched. + +Thu Oct 15 13:59:30 1998 Brian Gough <bjg@vvv.lanl.gov> + + * bisection.c, falsepos.c, secant.c: reordered the tests so that + the minimum number of function evaluations are performed when + there is an early exit due to one of the supplied limits lying on + a root. + +Fri Aug 21 14:48:13 1998 Brian Gough <bjg@vvv.lanl.gov> + + * test.c: clean up of tests to get rid of warnings + +Thu Aug 20 10:21:15 1998 Brian Gough <bjg@vvv.lanl.gov> + + * roots.h (_WITHIN_TOL): added extra parens in macro definition, + for safety + + * falsepos.c (gsl_root_falsepos): removed test for absolute + equality and replaced by a flag indicating which variables + changed. + + * test.c (main): simplified the tests, removed command line + arguments (can use the debugger to select which ones to run) + +Mon Jun 15 22:22:54 1998 Brian Gough <bjg@vvv.lanl.gov> + + * started to eliminate void * arguments for function types (they + are not a good idea and can easily be specified) + +1998-02-09 Mark Galassi <rosalia@cygnus.com> + + * test.c (main): added an extra argument so that the $(srcdir) can + be passed along when "make check" is run in a separate build + directory. + + * test-macros, test-secant, test-bisection, test-newton, + test-falsepos: modified these to use build and source directories + explicitly. Now "make check" in a separate build directory works. + +1998-02-02 Mark Galassi <rosalia@cygnus.com> + + * Makefile.am (TESTS): added $(srcdir) before these scripts, since + the TESTS target picks things from the build directory. diff --git a/gsl-1.9/roots/Makefile.am b/gsl-1.9/roots/Makefile.am new file mode 100644 index 0000000..2ea5015 --- /dev/null +++ b/gsl-1.9/roots/Makefile.am @@ -0,0 +1,19 @@ +# -*-makefile-*- + +noinst_LTLIBRARIES = libgslroots.la + +pkginclude_HEADERS = gsl_roots.h + +noinst_HEADERS = roots.h + +INCLUDES= -I$(top_builddir) + +libgslroots_la_SOURCES = bisection.c brent.c falsepos.c newton.c secant.c steffenson.c convergence.c fsolver.c fdfsolver.c + +check_PROGRAMS = test + +TESTS = $(check_PROGRAMS) + +test_SOURCES = test.c test_funcs.c test.h +test_LDADD = libgslroots.la ../ieee-utils/libgslieeeutils.la ../err/libgslerr.la ../test/libgsltest.la ../sys/libgslsys.la ../utils/libutils.la + diff --git a/gsl-1.9/roots/Makefile.in b/gsl-1.9/roots/Makefile.in new file mode 100644 index 0000000..66a22c7 --- /dev/null +++ b/gsl-1.9/roots/Makefile.in @@ -0,0 +1,548 @@ +# 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@ + +# -*-makefile-*- + + +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 = roots +DIST_COMMON = $(noinst_HEADERS) $(pkginclude_HEADERS) \ + $(srcdir)/Makefile.am $(srcdir)/Makefile.in ChangeLog TODO +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) +libgslroots_la_LIBADD = +am_libgslroots_la_OBJECTS = bisection.lo brent.lo falsepos.lo \ + newton.lo secant.lo steffenson.lo convergence.lo fsolver.lo \ + fdfsolver.lo +libgslroots_la_OBJECTS = $(am_libgslroots_la_OBJECTS) +am_test_OBJECTS = test.$(OBJEXT) test_funcs.$(OBJEXT) +test_OBJECTS = $(am_test_OBJECTS) +test_DEPENDENCIES = libgslroots.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 = $(libgslroots_la_SOURCES) $(test_SOURCES) +DIST_SOURCES = $(libgslroots_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 = $(noinst_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 = libgslroots.la +pkginclude_HEADERS = gsl_roots.h +noinst_HEADERS = roots.h +INCLUDES = -I$(top_builddir) +libgslroots_la_SOURCES = bisection.c brent.c falsepos.c newton.c secant.c steffenson.c convergence.c fsolver.c fdfsolver.c +TESTS = $(check_PROGRAMS) +test_SOURCES = test.c test_funcs.c test.h +test_LDADD = libgslroots.la ../ieee-utils/libgslieeeutils.la ../err/libgslerr.la ../test/libgsltest.la ../sys/libgslsys.la ../utils/libutils.la +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 roots/Makefile'; \ + cd $(top_srcdir) && \ + $(AUTOMAKE) --gnu --ignore-deps roots/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 +libgslroots.la: $(libgslroots_la_OBJECTS) $(libgslroots_la_DEPENDENCIES) + $(LINK) $(libgslroots_la_LDFLAGS) $(libgslroots_la_OBJECTS) $(libgslroots_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/roots/TODO b/gsl-1.9/roots/TODO new file mode 100644 index 0000000..06f17c5 --- /dev/null +++ b/gsl-1.9/roots/TODO @@ -0,0 +1,9 @@ +* Add an inline version of the iterate method for speed? Perhaps not, +the time taken for each iteration surely dominated by the convergence +test. + +* Numerical derivatives? Maybe have a function to manufacture an fdf +from an f and optionally a df. (We'll need to approximate the +derivative if it is not provided; this is something which should be +done outside the root finding package.) + diff --git a/gsl-1.9/roots/bisection.c b/gsl-1.9/roots/bisection.c new file mode 100644 index 0000000..e337ab0 --- /dev/null +++ b/gsl-1.9/roots/bisection.c @@ -0,0 +1,134 @@ +/* roots/bisection.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Reid Priedhorsky, 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. + */ + + +/* bisection.c -- bisection root finding algorithm */ + +#include <config.h> + +#include <stddef.h> +#include <stdlib.h> +#include <stdio.h> +#include <math.h> +#include <float.h> + +#include <gsl/gsl_math.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_roots.h> + +#include "roots.h" + +typedef struct + { + double f_lower, f_upper; + } +bisection_state_t; + +static int bisection_init (void * vstate, gsl_function * f, double * root, double x_lower, double x_upper); +static int bisection_iterate (void * vstate, gsl_function * f, double * root, double * x_lower, double * x_upper); + +static int +bisection_init (void * vstate, gsl_function * f, double * root, double x_lower, double x_upper) +{ + bisection_state_t * state = (bisection_state_t *) vstate; + + double f_lower, f_upper ; + + *root = 0.5 * (x_lower + x_upper) ; + + SAFE_FUNC_CALL (f, x_lower, &f_lower); + SAFE_FUNC_CALL (f, x_upper, &f_upper); + + state->f_lower = f_lower; + state->f_upper = f_upper; + + if ((f_lower < 0.0 && f_upper < 0.0) || (f_lower > 0.0 && f_upper > 0.0)) + { + GSL_ERROR ("endpoints do not straddle y=0", GSL_EINVAL); + } + + return GSL_SUCCESS; + +} + +static int +bisection_iterate (void * vstate, gsl_function * f, double * root, double * x_lower, double * x_upper) +{ + bisection_state_t * state = (bisection_state_t *) vstate; + + double x_bisect, f_bisect; + + const double x_left = *x_lower ; + const double x_right = *x_upper ; + + const double f_lower = state->f_lower; + const double f_upper = state->f_upper; + + if (f_lower == 0.0) + { + *root = x_left ; + *x_upper = x_left; + return GSL_SUCCESS; + } + + if (f_upper == 0.0) + { + *root = x_right ; + *x_lower = x_right; + return GSL_SUCCESS; + } + + x_bisect = (x_left + x_right) / 2.0; + + SAFE_FUNC_CALL (f, x_bisect, &f_bisect); + + if (f_bisect == 0.0) + { + *root = x_bisect; + *x_lower = x_bisect; + *x_upper = x_bisect; + return GSL_SUCCESS; + } + + /* Discard the half of the interval which doesn't contain the root. */ + + if ((f_lower > 0.0 && f_bisect < 0.0) || (f_lower < 0.0 && f_bisect > 0.0)) + { + *root = 0.5 * (x_left + x_bisect) ; + *x_upper = x_bisect; + state->f_upper = f_bisect; + } + else + { + *root = 0.5 * (x_bisect + x_right) ; + *x_lower = x_bisect; + state->f_lower = f_bisect; + } + + return GSL_SUCCESS; +} + + +static const gsl_root_fsolver_type bisection_type = +{"bisection", /* name */ + sizeof (bisection_state_t), + &bisection_init, + &bisection_iterate}; + +const gsl_root_fsolver_type * gsl_root_fsolver_bisection = &bisection_type; diff --git a/gsl-1.9/roots/brent.c b/gsl-1.9/roots/brent.c new file mode 100644 index 0000000..5cd5275 --- /dev/null +++ b/gsl-1.9/roots/brent.c @@ -0,0 +1,244 @@ +/* roots/brent.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Reid Priedhorsky, 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. + */ + +/* brent.c -- brent root finding algorithm */ + +#include <config.h> + +#include <stddef.h> +#include <stdlib.h> +#include <stdio.h> +#include <math.h> +#include <float.h> + +#include <gsl/gsl_math.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_roots.h> + +#include "roots.h" + + +typedef struct + { + double a, b, c, d, e; + double fa, fb, fc; + } +brent_state_t; + +static int brent_init (void * vstate, gsl_function * f, double * root, double x_lower, double x_upper); +static int brent_iterate (void * vstate, gsl_function * f, double * root, double * x_lower, double * x_upper); + + +static int +brent_init (void * vstate, gsl_function * f, double * root, double x_lower, double x_upper) +{ + brent_state_t * state = (brent_state_t *) vstate; + + double f_lower, f_upper ; + + *root = 0.5 * (x_lower + x_upper) ; + + SAFE_FUNC_CALL (f, x_lower, &f_lower); + SAFE_FUNC_CALL (f, x_upper, &f_upper); + + state->a = x_lower; + state->fa = f_lower; + + state->b = x_upper; + state->fb = f_upper; + + state->c = x_upper; + state->fc = f_upper; + + state->d = x_upper - x_lower ; + state->e = x_upper - x_lower ; + + if ((f_lower < 0.0 && f_upper < 0.0) || (f_lower > 0.0 && f_upper > 0.0)) + { + GSL_ERROR ("endpoints do not straddle y=0", GSL_EINVAL); + } + + return GSL_SUCCESS; + +} + +static int +brent_iterate (void * vstate, gsl_function * f, double * root, double * x_lower, double * x_upper) +{ + brent_state_t * state = (brent_state_t *) vstate; + + double tol, m; + + int ac_equal = 0; + + double a = state->a, b = state->b, c = state->c; + double fa = state->fa, fb = state->fb, fc = state->fc; + double d = state->d, e = state->e; + + if ((fb < 0 && fc < 0) || (fb > 0 && fc > 0)) + { + ac_equal = 1; + c = a; + fc = fa; + d = b - a; + e = b - a; + } + + if (fabs (fc) < fabs (fb)) + { + ac_equal = 1; + a = b; + b = c; + c = a; + fa = fb; + fb = fc; + fc = fa; + } + + tol = 0.5 * GSL_DBL_EPSILON * fabs (b); + m = 0.5 * (c - b); + + if (fb == 0) + { + *root = b; + *x_lower = b; + *x_upper = b; + + return GSL_SUCCESS; + } + + if (fabs (m) <= tol) + { + *root = b; + + if (b < c) + { + *x_lower = b; + *x_upper = c; + } + else + { + *x_lower = c; + *x_upper = b; + } + + return GSL_SUCCESS; + } + + if (fabs (e) < tol || fabs (fa) <= fabs (fb)) + { + d = m; /* use bisection */ + e = m; + } + else + { + double p, q, r; /* use inverse cubic interpolation */ + double s = fb / fa; + + if (ac_equal) + { + p = 2 * m * s; + q = 1 - s; + } + else + { + q = fa / fc; + r = fb / fc; + p = s * (2 * m * q * (q - r) - (b - a) * (r - 1)); + q = (q - 1) * (r - 1) * (s - 1); + } + + if (p > 0) + { + q = -q; + } + else + { + p = -p; + } + + if (2 * p < GSL_MIN (3 * m * q - fabs (tol * q), fabs (e * q))) + { + e = d; + d = p / q; + } + else + { + /* interpolation failed, fall back to bisection */ + + d = m; + e = m; + } + } + + a = b; + fa = fb; + + if (fabs (d) > tol) + { + b += d; + } + else + { + b += (m > 0 ? +tol : -tol); + } + + SAFE_FUNC_CALL (f, b, &fb); + + state->a = a ; + state->b = b ; + state->c = c ; + state->d = d ; + state->e = e ; + state->fa = fa ; + state->fb = fb ; + state->fc = fc ; + + /* Update the best estimate of the root and bounds on each + iteration */ + + *root = b; + + if ((fb < 0 && fc < 0) || (fb > 0 && fc > 0)) + { + c = a; + } + + if (b < c) + { + *x_lower = b; + *x_upper = c; + } + else + { + *x_lower = c; + *x_upper = b; + } + + return GSL_SUCCESS ; +} + + +static const gsl_root_fsolver_type brent_type = +{"brent", /* name */ + sizeof (brent_state_t), + &brent_init, + &brent_iterate}; + +const gsl_root_fsolver_type * gsl_root_fsolver_brent = &brent_type; diff --git a/gsl-1.9/roots/convergence.c b/gsl-1.9/roots/convergence.c new file mode 100644 index 0000000..34a78fe --- /dev/null +++ b/gsl-1.9/roots/convergence.c @@ -0,0 +1,87 @@ +/* roots/convergence.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Reid Priedhorsky, 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_errno.h> +#include <gsl/gsl_roots.h> + +int +gsl_root_test_interval (double x_lower, double x_upper, double epsabs, double epsrel) +{ + const double abs_lower = fabs(x_lower) ; + const double abs_upper = fabs(x_upper) ; + + double min_abs, tolerance; + + if (epsrel < 0.0) + GSL_ERROR ("relative tolerance is negative", GSL_EBADTOL); + + if (epsabs < 0.0) + GSL_ERROR ("absolute tolerance is negative", GSL_EBADTOL); + + if (x_lower > x_upper) + GSL_ERROR ("lower bound larger than upper bound", GSL_EINVAL); + + if ((x_lower > 0.0 && x_upper > 0.0) || (x_lower < 0.0 && x_upper < 0.0)) + { + min_abs = GSL_MIN_DBL(abs_lower, abs_upper) ; + } + else + { + min_abs = 0; + } + + tolerance = epsabs + epsrel * min_abs ; + + if (fabs(x_upper - x_lower) < tolerance) + return GSL_SUCCESS; + + return GSL_CONTINUE ; +} + +int +gsl_root_test_delta (double x1, double x0, double epsabs, double epsrel) +{ + const double tolerance = epsabs + epsrel * fabs(x1) ; + + if (epsrel < 0.0) + GSL_ERROR ("relative tolerance is negative", GSL_EBADTOL); + + if (epsabs < 0.0) + GSL_ERROR ("absolute tolerance is negative", GSL_EBADTOL); + + if (fabs(x1 - x0) < tolerance || x1 == x0) + return GSL_SUCCESS; + + return GSL_CONTINUE ; +} + +int +gsl_root_test_residual (double f, double epsabs) +{ + if (epsabs < 0.0) + GSL_ERROR ("absolute tolerance is negative", GSL_EBADTOL); + + if (fabs(f) < epsabs) + return GSL_SUCCESS; + + return GSL_CONTINUE ; +} + diff --git a/gsl-1.9/roots/falsepos.c b/gsl-1.9/roots/falsepos.c new file mode 100644 index 0000000..d264f69 --- /dev/null +++ b/gsl-1.9/roots/falsepos.c @@ -0,0 +1,178 @@ +/* roots/falsepos.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Reid Priedhorsky, 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. + */ + +/* falsepos.c -- falsepos root finding algorithm + + The false position algorithm uses bracketing by linear interpolation. + + If a linear interpolation step would decrease the size of the + bracket by less than a bisection step would then the algorithm + takes a bisection step instead. + + The last linear interpolation estimate of the root is used. If a + bisection step causes it to fall outside the brackets then it is + replaced by the bisection estimate (x_upper + x_lower)/2. + +*/ + +#include <config.h> + +#include <stddef.h> +#include <stdlib.h> +#include <stdio.h> +#include <math.h> +#include <float.h> + +#include <gsl/gsl_math.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_roots.h> + +#include "roots.h" + +typedef struct + { + double f_lower, f_upper; + } +falsepos_state_t; + +static int falsepos_init (void * vstate, gsl_function * f, double * root, double x_lower, double x_upper); +static int falsepos_iterate (void * vstate, gsl_function * f, double * root, double * x_lower, double * x_upper); + +static int +falsepos_init (void * vstate, gsl_function * f, double * root, double x_lower, double x_upper) +{ + falsepos_state_t * state = (falsepos_state_t *) vstate; + + double f_lower, f_upper ; + + *root = 0.5 * (x_lower + x_upper); + + SAFE_FUNC_CALL (f, x_lower, &f_lower); + SAFE_FUNC_CALL (f, x_upper, &f_upper); + + state->f_lower = f_lower; + state->f_upper = f_upper; + + if ((f_lower < 0.0 && f_upper < 0.0) || (f_lower > 0.0 && f_upper > 0.0)) + { + GSL_ERROR ("endpoints do not straddle y=0", GSL_EINVAL); + } + + return GSL_SUCCESS; + +} + +static int +falsepos_iterate (void * vstate, gsl_function * f, double * root, double * x_lower, double * x_upper) +{ + falsepos_state_t * state = (falsepos_state_t *) vstate; + + double x_linear, f_linear; + double x_bisect, f_bisect; + + double x_left = *x_lower ; + double x_right = *x_upper ; + + double f_lower = state->f_lower; + double f_upper = state->f_upper; + + double w ; + + if (f_lower == 0.0) + { + *root = x_left ; + *x_upper = x_left; + return GSL_SUCCESS; + } + + if (f_upper == 0.0) + { + *root = x_right ; + *x_lower = x_right; + return GSL_SUCCESS; + } + + /* Draw a line between f(*lower_bound) and f(*upper_bound) and + note where it crosses the X axis; that's where we will + split the interval. */ + + x_linear = x_right - (f_upper * (x_left - x_right) / (f_lower - f_upper)); + + SAFE_FUNC_CALL (f, x_linear, &f_linear); + + if (f_linear == 0.0) + { + *root = x_linear; + *x_lower = x_linear; + *x_upper = x_linear; + return GSL_SUCCESS; + } + + /* Discard the half of the interval which doesn't contain the root. */ + + if ((f_lower > 0.0 && f_linear < 0.0) || (f_lower < 0.0 && f_linear > 0.0)) + { + *root = x_linear ; + *x_upper = x_linear; + state->f_upper = f_linear; + w = x_linear - x_left ; + } + else + { + *root = x_linear ; + *x_lower = x_linear; + state->f_lower = f_linear; + w = x_right - x_linear; + } + + if (w < 0.5 * (x_right - x_left)) + { + return GSL_SUCCESS ; + } + + x_bisect = 0.5 * (x_left + x_right); + + SAFE_FUNC_CALL (f, x_bisect, &f_bisect); + + if ((f_lower > 0.0 && f_bisect < 0.0) || (f_lower < 0.0 && f_bisect > 0.0)) + { + *x_upper = x_bisect; + state->f_upper = f_bisect; + if (*root > x_bisect) + *root = 0.5 * (x_left + x_bisect) ; + } + else + { + *x_lower = x_bisect; + state->f_lower = f_bisect; + if (*root < x_bisect) + *root = 0.5 * (x_bisect + x_right) ; + } + + return GSL_SUCCESS; +} + + +static const gsl_root_fsolver_type falsepos_type = +{"falsepos", /* name */ + sizeof (falsepos_state_t), + &falsepos_init, + &falsepos_iterate}; + +const gsl_root_fsolver_type * gsl_root_fsolver_falsepos = &falsepos_type; diff --git a/gsl-1.9/roots/fdfsolver.c b/gsl-1.9/roots/fdfsolver.c new file mode 100644 index 0000000..07d0b27 --- /dev/null +++ b/gsl-1.9/roots/fdfsolver.c @@ -0,0 +1,88 @@ +/* roots/fdfsolver.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Reid Priedhorsky, 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 <stdlib.h> +#include <string.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_roots.h> + +gsl_root_fdfsolver * +gsl_root_fdfsolver_alloc (const gsl_root_fdfsolver_type * T) +{ + + gsl_root_fdfsolver * s = (gsl_root_fdfsolver *) malloc (sizeof (gsl_root_fdfsolver)); + + if (s == 0) + { + GSL_ERROR_VAL ("failed to allocate space for root solver struct", + GSL_ENOMEM, 0); + }; + + s->state = malloc (T->size); + + if (s->state == 0) + { + free (s); /* exception in constructor, avoid memory leak */ + + GSL_ERROR_VAL ("failed to allocate space for root solver state", + GSL_ENOMEM, 0); + }; + + s->type = T ; + s->fdf = NULL; + + return s; +} + +int +gsl_root_fdfsolver_set (gsl_root_fdfsolver * s, gsl_function_fdf * f, double root) +{ + s->fdf = f; + s->root = root; + + return (s->type->set) (s->state, s->fdf, &(s->root)); +} + +int +gsl_root_fdfsolver_iterate (gsl_root_fdfsolver * s) +{ + return (s->type->iterate) (s->state, s->fdf, &(s->root)); +} + +void +gsl_root_fdfsolver_free (gsl_root_fdfsolver * s) +{ + free (s->state); + free (s); +} + +const char * +gsl_root_fdfsolver_name (const gsl_root_fdfsolver * s) +{ + return s->type->name; +} + +double +gsl_root_fdfsolver_root (const gsl_root_fdfsolver * s) +{ + return s->root; +} + + diff --git a/gsl-1.9/roots/fsolver.c b/gsl-1.9/roots/fsolver.c new file mode 100644 index 0000000..039c554 --- /dev/null +++ b/gsl-1.9/roots/fsolver.c @@ -0,0 +1,107 @@ +/* roots/fsolver.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Reid Priedhorsky, 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 <stdlib.h> +#include <string.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_roots.h> + +gsl_root_fsolver * +gsl_root_fsolver_alloc (const gsl_root_fsolver_type * T) +{ + gsl_root_fsolver * s = (gsl_root_fsolver *) malloc (sizeof (gsl_root_fsolver)); + + if (s == 0) + { + GSL_ERROR_VAL ("failed to allocate space for root solver struct", + GSL_ENOMEM, 0); + }; + + s->state = malloc (T->size); + + if (s->state == 0) + { + free (s); /* exception in constructor, avoid memory leak */ + + GSL_ERROR_VAL ("failed to allocate space for root solver state", + GSL_ENOMEM, 0); + }; + + s->type = T ; + s->function = NULL ; + + return s; +} + +int +gsl_root_fsolver_set (gsl_root_fsolver * s, gsl_function * f, double x_lower, double x_upper) +{ + if (x_lower > x_upper) + { + GSL_ERROR ("invalid interval (lower > upper)", GSL_EINVAL); + } + + s->function = f; + s->root = 0.5 * (x_lower + x_upper); /* initial estimate */ + s->x_lower = x_lower; + s->x_upper = x_upper; + + return (s->type->set) (s->state, s->function, &(s->root), x_lower, x_upper); +} + +int +gsl_root_fsolver_iterate (gsl_root_fsolver * s) +{ + return (s->type->iterate) (s->state, + s->function, &(s->root), + &(s->x_lower), &(s->x_upper)); +} + +void +gsl_root_fsolver_free (gsl_root_fsolver * s) +{ + free (s->state); + free (s); +} + +const char * +gsl_root_fsolver_name (const gsl_root_fsolver * s) +{ + return s->type->name; +} + +double +gsl_root_fsolver_root (const gsl_root_fsolver * s) +{ + return s->root; +} + +double +gsl_root_fsolver_x_lower (const gsl_root_fsolver * s) +{ + return s->x_lower; +} + +double +gsl_root_fsolver_x_upper (const gsl_root_fsolver * s) +{ + return s->x_upper; +} + diff --git a/gsl-1.9/roots/gsl_roots.h b/gsl-1.9/roots/gsl_roots.h new file mode 100644 index 0000000..7103418 --- /dev/null +++ b/gsl-1.9/roots/gsl_roots.h @@ -0,0 +1,127 @@ +/* roots/gsl_roots.h + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Reid Priedhorsky, 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. + */ + +#ifndef __GSL_ROOTS_H__ +#define __GSL_ROOTS_H__ + +#include <stdlib.h> +#include <gsl/gsl_types.h> +#include <gsl/gsl_math.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 + +typedef struct + { + const char *name; + size_t size; + int (*set) (void *state, gsl_function * f, double * root, double x_lower, double x_upper); + int (*iterate) (void *state, gsl_function * f, double * root, double * x_lower, double * x_upper); + } +gsl_root_fsolver_type; + +typedef struct + { + const gsl_root_fsolver_type * type; + gsl_function * function ; + double root ; + double x_lower; + double x_upper; + void *state; + } +gsl_root_fsolver; + +typedef struct + { + const char *name; + size_t size; + int (*set) (void *state, gsl_function_fdf * f, double * root); + int (*iterate) (void *state, gsl_function_fdf * f, double * root); + } +gsl_root_fdfsolver_type; + +typedef struct + { + const gsl_root_fdfsolver_type * type; + gsl_function_fdf * fdf ; + double root ; + void *state; + } +gsl_root_fdfsolver; + +gsl_root_fsolver * +gsl_root_fsolver_alloc (const gsl_root_fsolver_type * T); +void gsl_root_fsolver_free (gsl_root_fsolver * s); + +int gsl_root_fsolver_set (gsl_root_fsolver * s, + gsl_function * f, + double x_lower, double x_upper); + +int gsl_root_fsolver_iterate (gsl_root_fsolver * s); + +const char * gsl_root_fsolver_name (const gsl_root_fsolver * s); +double gsl_root_fsolver_root (const gsl_root_fsolver * s); +double gsl_root_fsolver_x_lower (const gsl_root_fsolver * s); +double gsl_root_fsolver_x_upper (const gsl_root_fsolver * s); + + +gsl_root_fdfsolver * +gsl_root_fdfsolver_alloc (const gsl_root_fdfsolver_type * T); + +int +gsl_root_fdfsolver_set (gsl_root_fdfsolver * s, + gsl_function_fdf * fdf, double root); + +int +gsl_root_fdfsolver_iterate (gsl_root_fdfsolver * s); + +void +gsl_root_fdfsolver_free (gsl_root_fdfsolver * s); + +const char * gsl_root_fdfsolver_name (const gsl_root_fdfsolver * s); +double gsl_root_fdfsolver_root (const gsl_root_fdfsolver * s); + +int +gsl_root_test_interval (double x_lower, double x_upper, double epsabs, double epsrel); + +int +gsl_root_test_residual (double f, double epsabs); + +int +gsl_root_test_delta (double x1, double x0, double epsabs, double epsrel); + +GSL_VAR const gsl_root_fsolver_type * gsl_root_fsolver_bisection; +GSL_VAR const gsl_root_fsolver_type * gsl_root_fsolver_brent; +GSL_VAR const gsl_root_fsolver_type * gsl_root_fsolver_falsepos; +GSL_VAR const gsl_root_fdfsolver_type * gsl_root_fdfsolver_newton; +GSL_VAR const gsl_root_fdfsolver_type * gsl_root_fdfsolver_secant; +GSL_VAR const gsl_root_fdfsolver_type * gsl_root_fdfsolver_steffenson; + +__END_DECLS + +#endif /* __GSL_ROOTS_H__ */ diff --git a/gsl-1.9/roots/newton.c b/gsl-1.9/roots/newton.c new file mode 100644 index 0000000..09a4351 --- /dev/null +++ b/gsl-1.9/roots/newton.c @@ -0,0 +1,106 @@ +/* roots/newton.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Reid Priedhorsky, 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. + */ + +/* newton.c -- newton root finding algorithm + + This is the classical Newton-Raphson iteration. + + x[i+1] = x[i] - f(x[i])/f'(x[i]) + +*/ + +#include <config.h> + +#include <stddef.h> +#include <stdlib.h> +#include <stdio.h> +#include <math.h> +#include <float.h> + +#include <gsl/gsl_math.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_roots.h> + +#include "roots.h" + +typedef struct + { + double f, df; + } +newton_state_t; + +static int newton_init (void * vstate, gsl_function_fdf * fdf, double * root); +static int newton_iterate (void * vstate, gsl_function_fdf * fdf, double * root); + +static int +newton_init (void * vstate, gsl_function_fdf * fdf, double * root) +{ + newton_state_t * state = (newton_state_t *) vstate; + + const double x = *root ; + + state->f = GSL_FN_FDF_EVAL_F (fdf, x); + state->df = GSL_FN_FDF_EVAL_DF (fdf, x) ; + + return GSL_SUCCESS; + +} + +static int +newton_iterate (void * vstate, gsl_function_fdf * fdf, double * root) +{ + newton_state_t * state = (newton_state_t *) vstate; + + double root_new, f_new, df_new; + + if (state->df == 0.0) + { + GSL_ERROR("derivative is zero", GSL_EZERODIV); + } + + root_new = *root - (state->f / state->df); + + *root = root_new ; + + GSL_FN_FDF_EVAL_F_DF(fdf, root_new, &f_new, &df_new); + + state->f = f_new ; + state->df = df_new ; + + if (!finite(f_new)) + { + GSL_ERROR ("function value is not finite", GSL_EBADFUNC); + } + + if (!finite (df_new)) + { + GSL_ERROR ("derivative value is not finite", GSL_EBADFUNC); + } + + return GSL_SUCCESS; +} + + +static const gsl_root_fdfsolver_type newton_type = +{"newton", /* name */ + sizeof (newton_state_t), + &newton_init, + &newton_iterate}; + +const gsl_root_fdfsolver_type * gsl_root_fdfsolver_newton = &newton_type; diff --git a/gsl-1.9/roots/roots.h b/gsl-1.9/roots/roots.h new file mode 100644 index 0000000..7b6245e --- /dev/null +++ b/gsl-1.9/roots/roots.h @@ -0,0 +1,37 @@ +/* roots/roots.h + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Reid Priedhorsky, 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. + */ + +/* roots.h -- declarations for internal root finding and RF support stuff. */ + +#ifndef __ROOTS_H__ +#define __ROOTS_H__ + +/* Call the pointed-to function with argument x, put its result in y, and + return an error if the function value is Inf/Nan. */ + +#define SAFE_FUNC_CALL(f, x, yp) \ +do { \ + *yp = GSL_FN_EVAL(f,x); \ + if (!finite(*yp)) \ + GSL_ERROR("function value is not finite", GSL_EBADFUNC); \ +} while (0) + +#endif /* __ROOTS_H__ */ + + diff --git a/gsl-1.9/roots/secant.c b/gsl-1.9/roots/secant.c new file mode 100644 index 0000000..c44eea5 --- /dev/null +++ b/gsl-1.9/roots/secant.c @@ -0,0 +1,117 @@ +/* roots/secant.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Reid Priedhorsky, 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. + */ + +/* secant.c -- secant root finding algorithm + + The secant algorithm is a variant of the Newton algorithm with the + derivative term replaced by a numerical estimate from the last two + function evaluations. + + x[i+1] = x[i] - f(x[i]) / f'_est + + where f'_est = (f(x[i]) - f(x[i-1])) / (x[i] - x[i-1]) + + The exact derivative is used for the initial value of f'_est. + +*/ + +#include <config.h> + +#include <stddef.h> +#include <stdlib.h> +#include <stdio.h> +#include <math.h> +#include <float.h> + +#include <gsl/gsl_math.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_roots.h> + +#include "roots.h" + +typedef struct + { + double f; + double df; + } +secant_state_t; + +static int secant_init (void * vstate, gsl_function_fdf * fdf, double * root); +static int secant_iterate (void * vstate, gsl_function_fdf * fdf, double * root); + +static int +secant_init (void * vstate, gsl_function_fdf * fdf, double * root) +{ + secant_state_t * state = (secant_state_t *) vstate; + + const double x = *root; + + GSL_FN_FDF_EVAL_F_DF (fdf, x, &(state->f), &(state->df)); + + return GSL_SUCCESS; + +} + +static int +secant_iterate (void * vstate, gsl_function_fdf * fdf, double * root) +{ + secant_state_t * state = (secant_state_t *) vstate; + + const double x = *root ; + const double f = state->f; + const double df = state->df; + + double x_new, f_new, df_new; + + if (state->df == 0.0) + { + GSL_ERROR("derivative is zero", GSL_EZERODIV); + } + + x_new = x - (f / df); + + f_new = GSL_FN_FDF_EVAL_F(fdf, x_new) ; + df_new = (f_new - f) / (x_new - x) ; + + *root = x_new ; + + state->f = f_new ; + state->df = df_new ; + + if (!finite (f_new)) + { + GSL_ERROR ("function value is not finite", GSL_EBADFUNC); + } + + if (!finite (df_new)) + { + GSL_ERROR ("derivative value is not finite", GSL_EBADFUNC); + } + + return GSL_SUCCESS; +} + + +static const gsl_root_fdfsolver_type secant_type = +{"secant", /* name */ + sizeof (secant_state_t), + &secant_init, + &secant_iterate}; + +const gsl_root_fdfsolver_type * gsl_root_fdfsolver_secant = &secant_type; diff --git a/gsl-1.9/roots/steffenson.c b/gsl-1.9/roots/steffenson.c new file mode 100644 index 0000000..e4169f6 --- /dev/null +++ b/gsl-1.9/roots/steffenson.c @@ -0,0 +1,144 @@ +/* roots/steffenson.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Reid Priedhorsky, 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. + */ + +/* steffenson.c -- steffenson root finding algorithm + + This is Newton's method with an Aitken "delta-squared" + acceleration of the iterates. This can improve the convergence on + multiple roots where the ordinary Newton algorithm is slow. + + x[i+1] = x[i] - f(x[i]) / f'(x[i]) + + x_accelerated[i] = x[i] - (x[i+1] - x[i])**2 / (x[i+2] - 2*x[i+1] + x[i]) + + We can only use the accelerated estimate after three iterations, + and use the unaccelerated value until then. + + */ + +#include <config.h> + +#include <stddef.h> +#include <stdlib.h> +#include <stdio.h> +#include <math.h> +#include <float.h> + +#include <gsl/gsl_math.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_roots.h> + +#include "roots.h" + +typedef struct + { + double f, df; + double x; + double x_1; + double x_2; + int count; + } +steffenson_state_t; + +static int steffenson_init (void * vstate, gsl_function_fdf * fdf, double * root); +static int steffenson_iterate (void * vstate, gsl_function_fdf * fdf, double * root); + +static int +steffenson_init (void * vstate, gsl_function_fdf * fdf, double * root) +{ + steffenson_state_t * state = (steffenson_state_t *) vstate; + + const double x = *root ; + + state->f = GSL_FN_FDF_EVAL_F (fdf, x); + state->df = GSL_FN_FDF_EVAL_DF (fdf, x) ; + + state->x = x; + state->x_1 = 0.0; + state->x_2 = 0.0; + + state->count = 1; + + return GSL_SUCCESS; + +} + +static int +steffenson_iterate (void * vstate, gsl_function_fdf * fdf, double * root) +{ + steffenson_state_t * state = (steffenson_state_t *) vstate; + + double x_new, f_new, df_new; + + double x_1 = state->x_1 ; + double x = state->x ; + + if (state->df == 0.0) + { + GSL_ERROR("derivative is zero", GSL_EZERODIV); + } + + x_new = x - (state->f / state->df); + + GSL_FN_FDF_EVAL_F_DF(fdf, x_new, &f_new, &df_new); + + state->x_2 = x_1 ; + state->x_1 = x ; + state->x = x_new; + + state->f = f_new ; + state->df = df_new ; + + if (!finite (f_new)) + { + GSL_ERROR ("function value is not finite", GSL_EBADFUNC); + } + + if (state->count < 3) + { + *root = x_new ; + state->count++ ; + } + else + { + double u = (x - x_1) ; + double v = (x_new - 2 * x + x_1); + + if (v == 0) + *root = x_new; /* avoid division by zero */ + else + *root = x_1 - u * u / v ; /* accelerated value */ + } + + if (!finite (df_new)) + { + GSL_ERROR ("derivative value is not finite", GSL_EBADFUNC); + } + + return GSL_SUCCESS; +} + + +static const gsl_root_fdfsolver_type steffenson_type = +{"steffenson", /* name */ + sizeof (steffenson_state_t), + &steffenson_init, + &steffenson_iterate}; + +const gsl_root_fdfsolver_type * gsl_root_fdfsolver_steffenson = &steffenson_type; diff --git a/gsl-1.9/roots/test.c b/gsl-1.9/roots/test.c new file mode 100644 index 0000000..4e9b94b --- /dev/null +++ b/gsl-1.9/roots/test.c @@ -0,0 +1,315 @@ +/* roots/test.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Reid Priedhorsky, 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 <stdlib.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_test.h> +#include <gsl/gsl_roots.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_ieee_utils.h> + +#include "roots.h" +#include "test.h" + +/* stopping parameters */ +const double EPSREL = (10 * GSL_DBL_EPSILON); +const double EPSABS = (10 * GSL_DBL_EPSILON); +const unsigned int MAX_ITERATIONS = 150; + +void my_error_handler (const char *reason, const char *file, + int line, int err); + +#define WITHIN_TOL(a, b, epsrel, epsabs) \ + ((fabs((a) - (b)) < (epsrel) * GSL_MIN(fabs(a), fabs(b)) + (epsabs))) + +int +main (void) +{ + gsl_function F_sin, F_cos, F_func1, F_func2, F_func3, F_func4, + F_func5, F_func6; + + gsl_function_fdf FDF_sin, FDF_cos, FDF_func1, FDF_func2, FDF_func3, FDF_func4, + FDF_func5, FDF_func6; + + const gsl_root_fsolver_type * fsolver[4] ; + const gsl_root_fdfsolver_type * fdfsolver[4] ; + + const gsl_root_fsolver_type ** T; + const gsl_root_fdfsolver_type ** S; + + gsl_ieee_env_setup(); + + fsolver[0] = gsl_root_fsolver_bisection; + fsolver[1] = gsl_root_fsolver_brent; + fsolver[2] = gsl_root_fsolver_falsepos; + fsolver[3] = 0; + + fdfsolver[0] = gsl_root_fdfsolver_newton; + fdfsolver[1] = gsl_root_fdfsolver_secant; + fdfsolver[2] = gsl_root_fdfsolver_steffenson; + fdfsolver[3] = 0; + + F_sin = create_function (sin_f) ; + F_cos = create_function (cos_f) ; + F_func1 = create_function (func1) ; + F_func2 = create_function (func2) ; + F_func3 = create_function (func3) ; + F_func4 = create_function (func4) ; + F_func5 = create_function (func5) ; + F_func6 = create_function (func6) ; + + FDF_sin = create_fdf (sin_f, sin_df, sin_fdf) ; + FDF_cos = create_fdf (cos_f, cos_df, cos_fdf) ; + FDF_func1 = create_fdf (func1, func1_df, func1_fdf) ; + FDF_func2 = create_fdf (func2, func2_df, func2_fdf) ; + FDF_func3 = create_fdf (func3, func3_df, func3_fdf) ; + FDF_func4 = create_fdf (func4, func4_df, func4_fdf) ; + FDF_func5 = create_fdf (func5, func5_df, func5_fdf) ; + FDF_func6 = create_fdf (func6, func6_df, func6_fdf) ; + + gsl_set_error_handler (&my_error_handler); + + for (T = fsolver ; *T != 0 ; T++) + { + test_f (*T, "sin(x) [3, 4]", &F_sin, 3.0, 4.0, M_PI); + test_f (*T, "sin(x) [-4, -3]", &F_sin, -4.0, -3.0, -M_PI); + test_f (*T, "sin(x) [-1/3, 1]", &F_sin, -1.0 / 3.0, 1.0, 0.0); + test_f (*T, "cos(x) [0, 3]", &F_cos, 0.0, 3.0, M_PI / 2.0); + test_f (*T, "cos(x) [-3, 0]", &F_cos, -3.0, 0.0, -M_PI / 2.0); + test_f (*T, "x^20 - 1 [0.1, 2]", &F_func1, 0.1, 2.0, 1.0); + test_f (*T, "sqrt(|x|)*sgn(x)", &F_func2, -1.0 / 3.0, 1.0, 0.0); + test_f (*T, "x^2 - 1e-8 [0, 1]", &F_func3, 0.0, 1.0, sqrt (1e-8)); + test_f (*T, "x exp(-x) [-1/3, 2]", &F_func4, -1.0 / 3.0, 2.0, 0.0); + test_f (*T, "(x - 1)^7 [0.9995, 1.0002]", &F_func6, 0.9995, 1.0002, 1.0); + + test_f_e (*T, "invalid range check [4, 0]", &F_sin, 4.0, 0.0, M_PI); + test_f_e (*T, "invalid range check [1, 1]", &F_sin, 1.0, 1.0, M_PI); + test_f_e (*T, "invalid range check [0.1, 0.2]", &F_sin, 0.1, 0.2, M_PI); + } + + for (S = fdfsolver ; *S != 0 ; S++) + { + test_fdf (*S,"sin(x) {3.4}", &FDF_sin, 3.4, M_PI); + test_fdf (*S,"sin(x) {-3.3}", &FDF_sin, -3.3, -M_PI); + test_fdf (*S,"sin(x) {0.5}", &FDF_sin, 0.5, 0.0); + test_fdf (*S,"cos(x) {0.6}", &FDF_cos, 0.6, M_PI / 2.0); + test_fdf (*S,"cos(x) {-2.5}", &FDF_cos, -2.5, -M_PI / 2.0); + test_fdf (*S,"x^{20} - 1 {0.9}", &FDF_func1, 0.9, 1.0); + test_fdf (*S,"x^{20} - 1 {1.1}", &FDF_func1, 1.1, 1.0); + test_fdf (*S,"sqrt(|x|)*sgn(x) {1.001}", &FDF_func2, 0.001, 0.0); + test_fdf (*S,"x^2 - 1e-8 {1}", &FDF_func3, 1.0, sqrt (1e-8)); + test_fdf (*S,"x exp(-x) {-2}", &FDF_func4, -2.0, 0.0); + test_fdf_e (*S,"max iterations x -> +Inf, x exp(-x) {2}", &FDF_func4, 2.0, 0.0); + test_fdf_e (*S,"max iterations x -> -Inf, 1/(1 + exp(-x)) {0}", &FDF_func5, 0.0, 0.0); + } + + test_fdf (gsl_root_fdfsolver_steffenson, + "(x - 1)^7 {0.9}", &FDF_func6, 0.9, 1.0); + + /* now summarize the results */ + + exit (gsl_test_summary ()); +} + + +/* Using gsl_root_bisection, find the root of the function pointed to by f, + using the interval [lower_bound, upper_bound]. Check if f succeeded and + that it was accurate enough. */ + +void +test_f (const gsl_root_fsolver_type * T, const char * description, gsl_function *f, + double lower_bound, double upper_bound, double correct_root) +{ + int status; + size_t iterations = 0; + double r, a, b; + double x_lower, x_upper; + gsl_root_fsolver * s; + + x_lower = lower_bound; + x_upper = upper_bound; + + s = gsl_root_fsolver_alloc(T); + gsl_root_fsolver_set(s, f, x_lower, x_upper) ; + + do + { + iterations++ ; + + gsl_root_fsolver_iterate (s); + + r = gsl_root_fsolver_root(s); + + a = gsl_root_fsolver_x_lower(s); + b = gsl_root_fsolver_x_upper(s); + + if (a > b) + gsl_test (GSL_FAILURE, "interval is invalid (%g,%g)", a, b); + + if (r < a || r > b) + gsl_test (GSL_FAILURE, "r lies outside interval %g (%g,%g)", r, a, b); + + status = gsl_root_test_interval (a,b, EPSABS, EPSREL); + } + while (status == GSL_CONTINUE && iterations < MAX_ITERATIONS); + + + gsl_test (status, "%s, %s (%g obs vs %g expected) ", + gsl_root_fsolver_name(s), description, + gsl_root_fsolver_root(s), correct_root); + + if (iterations == MAX_ITERATIONS) + { + gsl_test (GSL_FAILURE, "exceeded maximum number of iterations"); + } + + /* check the validity of the returned result */ + + if (!WITHIN_TOL (r, correct_root, EPSREL, EPSABS)) + { + gsl_test (GSL_FAILURE, "incorrect precision (%g obs vs %g expected)", + r, correct_root); + + } + + gsl_root_fsolver_free(s); +} + +void +test_f_e (const gsl_root_fsolver_type * T, + const char * description, gsl_function *f, + double lower_bound, double upper_bound, double correct_root) +{ + int status; + size_t iterations = 0; + double x_lower, x_upper; + gsl_root_fsolver * s; + + x_lower = lower_bound; + x_upper = upper_bound; + + s = gsl_root_fsolver_alloc(T); + status = gsl_root_fsolver_set(s, f, x_lower, x_upper) ; + + gsl_test (status != GSL_EINVAL, "%s (set), %s", T->name, description); + + if (status == GSL_EINVAL) + { + gsl_root_fsolver_free(s); + return ; + } + + do + { + iterations++ ; + gsl_root_fsolver_iterate (s); + x_lower = gsl_root_fsolver_x_lower(s); + x_upper = gsl_root_fsolver_x_lower(s); + status = gsl_root_test_interval (x_lower, x_upper, + EPSABS, EPSREL); + } + while (status == GSL_CONTINUE && iterations < MAX_ITERATIONS); + + gsl_test (!status, "%s, %s", gsl_root_fsolver_name(s), description, + gsl_root_fsolver_root(s) - correct_root); + + gsl_root_fsolver_free(s); +} + +void +test_fdf (const gsl_root_fdfsolver_type * T, const char * description, + gsl_function_fdf *fdf, double root, double correct_root) +{ + int status; + size_t iterations = 0; + double prev = 0 ; + + gsl_root_fdfsolver * s = gsl_root_fdfsolver_alloc(T); + gsl_root_fdfsolver_set (s, fdf, root) ; + + do + { + iterations++ ; + prev = gsl_root_fdfsolver_root(s); + gsl_root_fdfsolver_iterate (s); + status = gsl_root_test_delta(gsl_root_fdfsolver_root(s), prev, + EPSABS, EPSREL); + } + while (status == GSL_CONTINUE && iterations < MAX_ITERATIONS); + + gsl_test (status, "%s, %s (%g obs vs %g expected) ", + gsl_root_fdfsolver_name(s), description, + gsl_root_fdfsolver_root(s), correct_root); + + if (iterations == MAX_ITERATIONS) + { + gsl_test (GSL_FAILURE, "exceeded maximum number of iterations"); + } + + /* check the validity of the returned result */ + + if (!WITHIN_TOL (gsl_root_fdfsolver_root(s), correct_root, + EPSREL, EPSABS)) + { + gsl_test (GSL_FAILURE, "incorrect precision (%g obs vs %g expected)", + gsl_root_fdfsolver_root(s), correct_root); + + } + gsl_root_fdfsolver_free(s); +} + +void +test_fdf_e (const gsl_root_fdfsolver_type * T, + const char * description, gsl_function_fdf *fdf, + double root, double correct_root) +{ + int status; + size_t iterations = 0; + double prev = 0 ; + + gsl_root_fdfsolver * s = gsl_root_fdfsolver_alloc(T); + status = gsl_root_fdfsolver_set (s, fdf, root) ; + + gsl_test (status, "%s (set), %s", T->name, description); + + do + { + iterations++ ; + prev = gsl_root_fdfsolver_root(s); + gsl_root_fdfsolver_iterate (s); + status = gsl_root_test_delta(gsl_root_fdfsolver_root(s), prev, + EPSABS, EPSREL); + } + while (status == GSL_CONTINUE && iterations < MAX_ITERATIONS); + + gsl_test (!status, "%s, %s", gsl_root_fdfsolver_name(s), + description, gsl_root_fdfsolver_root(s) - correct_root); + gsl_root_fdfsolver_free(s); +} + +void +my_error_handler (const char *reason, const char *file, int line, int err) +{ + if (0) + printf ("(caught [%s:%d: %s (%d)])\n", file, line, reason, err); +} + + + diff --git a/gsl-1.9/roots/test.h b/gsl-1.9/roots/test.h new file mode 100644 index 0000000..c00febf --- /dev/null +++ b/gsl-1.9/roots/test.h @@ -0,0 +1,129 @@ +/* roots/test.h + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Reid Priedhorsky, 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. + */ + +gsl_function create_function (double (*f)(double, void *)) ; +gsl_function_fdf create_fdf (double (*f)(double, void *), + double (*df)(double, void *), + void (*fdf)(double, void *, double *, double *)); + +void + test_macros (void); + +void + test_roots (void); + +void + test_poly (void); + +void +test_f (const gsl_root_fsolver_type * T, + const char * description, gsl_function *f, + double lower_bound, double upper_bound, double correct_root); + +void +test_f_e (const gsl_root_fsolver_type * T, const char * description, + gsl_function *f, + double lower_bound, double upper_bound, double correct_root); + +void +test_fdf (const gsl_root_fdfsolver_type * T, const char * description, + gsl_function_fdf *fdf, double root, double correct_root); + +void +test_fdf_e (const gsl_root_fdfsolver_type * T, const char * description, + gsl_function_fdf *fdf, double root, double correct_root); + + +void + usage (void); + +void + error_handler (const char *reason, const char *file, int line); + +double + func1 (double x, void * p); + +double + func1_df (double x, void * p); + +void + func1_fdf (double x, void * p, double *y, double *yprime); + +double + func2 (double x, void * p); + +double + func2_df (double x, void * p); + +void + func2_fdf (double x, void * p, double *y, double *yprime); + +double + func3 (double x, void * p); + +double + func3_df (double x, void * p); + +void + func3_fdf (double x, void * p, double *y, double *yprime); + +double + func4 (double x, void * p); + +double + func4_df (double x, void * p); + +void + func4_fdf (double x, void * p, double *y, double *yprime); + +double + func5 (double x, void * p); + +double + func5_df (double x, void * p); + +void + func5_fdf (double x, void * p, double *y, double *yprime); + +double + func6 (double x, void * p); + +double + func6_df (double x, void * p); + +void + func6_fdf (double x, void * p, double *y, double *yprime); + +double + sin_f (double x, void * p); + +double + sin_df (double x, void * p); + +void + sin_fdf (double x, void * p, double *y, double *yprime); + +double + cos_f (double x, void * p); + +double + cos_df (double x, void * p); + +void + cos_fdf (double x, void * p, double *y, double *yprime); diff --git a/gsl-1.9/roots/test_funcs.c b/gsl-1.9/roots/test_funcs.c new file mode 100644 index 0000000..b9773a3 --- /dev/null +++ b/gsl-1.9/roots/test_funcs.c @@ -0,0 +1,229 @@ +/* roots/test_funcs.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Reid Priedhorsky, 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 <math.h> +#include <stdlib.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_roots.h> + +#include "test.h" + +gsl_function create_function (double (*f)(double, void *)) +{ + gsl_function F ; + F.function = f; + F.params = 0; + return F ; +} + +gsl_function_fdf create_fdf (double (*f)(double, void *), + double (*df)(double, void *), + void (*fdf)(double, void *, double *, double *)) +{ + gsl_function_fdf FDF ; + FDF.f = f ; + FDF.df = df ; + FDF.fdf = fdf ; + FDF.params = 0 ; + return FDF ; +} + +/* f(x) = x^{20} - 1 */ +/* f'(x) = 20x^{19} */ +/* zero at x = 1 or -1 */ + +double +func1 (double x, void *p) +{ + return pow (x, 20.0) - 1; +} + +double +func1_df (double x, void * p) +{ + return 20.0 * pow (x, 19.0); +} + +void +func1_fdf (double x, void * p, double *y, double *yprime) +{ + *y = func1 (x, p); + *yprime = 20.0 * pow (x, 19.0); +} + +/* f(x) = sqrt(abs(x))*sgn(x) */ +/* f'(x) = 1 / sqrt(abs(x) */ +/* zero at x = 0 */ +double +func2 (double x, void * p) +{ + double delta; + + if (x > 0) + delta = 1.0; + else if (x < 0) + delta = -1.0; + else + delta = 0.0; + + return sqrt (fabs (x)) * delta; +} + +double +func2_df (double x, void * p) +{ + return 1 / sqrt (fabs (x)); +} + +void +func2_fdf (double x, void * p, double *y, double *yprime) +{ + *y = func2 (x, p); + *yprime = 1 / sqrt (fabs (x)); +} + + +/* f(x) = x^2 - 1e-8 */ +/* f'(x) = 2x */ +/* zero at x = sqrt(1e-8) or -sqrt(1e-8) */ +double +func3 (double x, void * p) +{ + return pow (x, 2.0) - 1e-8; +} + +double +func3_df (double x, void * p) +{ + return 2 * x; +} + +void +func3_fdf (double x, void * p, double *y, double *yprime) +{ + *y = func3 (x, p); + *yprime = 2 * x; +} + +/* f(x) = x exp(-x) */ +/* f'(x) = exp(-x) - x exp(-x) */ +/* zero at x = 0 */ +double +func4 (double x, void * p) +{ + return x * exp (-x); +} + +double +func4_df (double x, void * p) +{ + return exp (-x) - x * exp (-x); +} + +void +func4_fdf (double x, void * p, double *y, double *yprime) +{ + *y = func4 (x, p); + *yprime = exp (-x) - x * exp (-x); +} + +/* f(x) = 1/(1+exp(x)) */ +/* f'(x) = -exp(x) / (1 + exp(x))^2 */ +/* no roots! */ +double +func5 (double x, void * p) +{ + return 1 / (1 + exp (x)); +} + +double +func5_df (double x, void * p) +{ + return -exp (x) / pow (1 + exp (x), 2.0); +} + +void +func5_fdf (double x, void * p, double *y, double *yprime) +{ + *y = func5 (x, p); + *yprime = -exp (x) / pow (1 + exp (x), 2.0); +} + +/* f(x) = (x - 1)^7 */ +/* f'(x) = 7 * (x - 1)^6 */ +/* zero at x = 1 */ +double +func6 (double x, void * p) +{ + return pow (x - 1, 7.0); +} + +double +func6_df (double x, void * p) +{ + return 7.0 * pow (x - 1, 6.0); +} + +void +func6_fdf (double x, void * p, double *y, double *yprime) +{ + *y = func6 (x, p); + *yprime = 7.0 * pow (x - 1, 6.0); +} + +/* sin(x) packaged up nicely. */ +double +sin_f (double x, void * p) +{ + return sin (x); +} + +double +sin_df (double x, void * p) +{ + return cos (x); +} + +void +sin_fdf (double x, void * p, double *y, double *yprime) +{ + *y = sin (x); + *yprime = cos (x); +} + +/* cos(x) packaged up nicely. */ +double +cos_f (double x, void * p) +{ + return cos (x); +} + +double +cos_df (double x, void * p) +{ + return -sin (x); +} + +void +cos_fdf (double x, void * p, double *y, double *yprime) +{ + *y = cos (x); + *yprime = -sin (x); +} |