diff options
Diffstat (limited to 'gsl-1.9/poly')
| -rw-r--r-- | gsl-1.9/poly/ChangeLog | 72 | ||||
| -rw-r--r-- | gsl-1.9/poly/Makefile.am | 17 | ||||
| -rw-r--r-- | gsl-1.9/poly/Makefile.in | 546 | ||||
| -rw-r--r-- | gsl-1.9/poly/TODO | 130 | ||||
| -rw-r--r-- | gsl-1.9/poly/balance.c | 127 | ||||
| -rw-r--r-- | gsl-1.9/poly/companion.c | 36 | ||||
| -rw-r--r-- | gsl-1.9/poly/dd.c | 101 | ||||
| -rw-r--r-- | gsl-1.9/poly/eval.c | 36 | ||||
| -rw-r--r-- | gsl-1.9/poly/gsl_poly.h | 132 | ||||
| -rw-r--r-- | gsl-1.9/poly/qr.c | 245 | ||||
| -rw-r--r-- | gsl-1.9/poly/solve_cubic.c | 110 | ||||
| -rw-r--r-- | gsl-1.9/poly/solve_quadratic.c | 85 | ||||
| -rw-r--r-- | gsl-1.9/poly/test.c | 508 | ||||
| -rw-r--r-- | gsl-1.9/poly/zsolve.c | 85 | ||||
| -rw-r--r-- | gsl-1.9/poly/zsolve_cubic.c | 153 | ||||
| -rw-r--r-- | gsl-1.9/poly/zsolve_init.c | 69 | ||||
| -rw-r--r-- | gsl-1.9/poly/zsolve_quadratic.c | 103 |
17 files changed, 2555 insertions, 0 deletions
diff --git a/gsl-1.9/poly/ChangeLog b/gsl-1.9/poly/ChangeLog new file mode 100644 index 0000000..87b227d --- /dev/null +++ b/gsl-1.9/poly/ChangeLog @@ -0,0 +1,72 @@ +2005-07-03 Brian Gough <bjg@network-theory.co.uk> + + * test.c (main): added tests for linear case + + * zsolve_quadratic.c (gsl_poly_complex_solve_quadratic): handle + the linear case + + * solve_quadratic.c (gsl_poly_solve_quadratic): handle the linear + case + +2005-05-19 Brian Gough <bjg@network-theory.co.uk> + + * Makefile.am (noinst_HEADERS): removed norm.c (unused) + +Sun Dec 2 22:02:31 2001 Brian Gough <bjg@network-theory.co.uk> + + * dd.c: added divided differences code from Dan, Ho-Jin + +Wed Oct 31 18:42:10 2001 Brian Gough <bjg@network-theory.co.uk> + + * qr.c (qr_companion): increased maximum number of iterations from + 30 to 60 + +Thu Jun 28 11:24:51 2001 Brian Gough <bjg@network-theory.co.uk> + + * test.c (main): fixed a subtle bug in the test where the array + for the complex results was too small by a factor of two + +Mon Apr 30 12:36:08 2001 Brian Gough <bjg@network-theory.co.uk> + + * eval.c (gsl_poly_eval): added eval function from specfunc + +Tue Aug 24 12:02:47 1999 Brian Gough <bjg@network-theory.co.uk> + + * zsolve.c: added general solver using QR method + +1999-08-22 Mark Galassi <rosalia@lanl.gov> + + * Makefile.am (libgslpoly_a_SOURCES): for now commented out + zsolve.c, since it has some strange errors that make it look like + it was not committed. + +Tue Aug 17 14:23:47 1999 Brian Gough <bjg@network-theory.co.uk> + + * zsolve_cubic.c (gsl_poly_complex_solve_cubic): compute the + discriminant without division so that it is exact for exact inputs + + * solve_cubic.c (gsl_poly_solve_cubic): compute the discriminant + without division so that it is exact for exact inputs + + * changed the way roots are counted to make all the routines + consistent. Conincident roots are no longer a special case, now + they are made up of separate roots which happen to have the same + values. + + * zsolve_quadratic.c (gsl_poly_complex_solve_quadratic): in the + case of a multiple root set all the returned values r0 = r1 = + root, rather than just setting one and leaving the others + potentially unset. + + * solve_quadratic.c (gsl_poly_solve_quadratic): ditto + +Tue Aug 3 19:36:26 1999 Brian Gough <bjg@network-theory.co.uk> + + * gsl_poly.h: changed all functions to take separate variables for + roots (z0, z1, z2) of variable-length arrays + +Sat Feb 20 12:13:46 1999 Brian Gough <bjg@netsci.freeserve.co.uk> + + * split out polynomial root finding algorithms into a new poly/ + directory + diff --git a/gsl-1.9/poly/Makefile.am b/gsl-1.9/poly/Makefile.am new file mode 100644 index 0000000..7de94ac --- /dev/null +++ b/gsl-1.9/poly/Makefile.am @@ -0,0 +1,17 @@ +noinst_LTLIBRARIES = libgslpoly.la + +pkginclude_HEADERS = gsl_poly.h + +INCLUDES= -I$(top_builddir) + +libgslpoly_la_SOURCES = dd.c eval.c solve_quadratic.c solve_cubic.c zsolve_quadratic.c zsolve_cubic.c zsolve.c zsolve_init.c + +noinst_HEADERS = balance.c companion.c qr.c + +TESTS = $(check_PROGRAMS) + +check_PROGRAMS = test + +test_SOURCES = test.c +test_LDADD = libgslpoly.la ../ieee-utils/libgslieeeutils.la ../err/libgslerr.la ../test/libgsltest.la ../sys/libgslsys.la ../utils/libutils.la + diff --git a/gsl-1.9/poly/Makefile.in b/gsl-1.9/poly/Makefile.in new file mode 100644 index 0000000..e0c9645 --- /dev/null +++ b/gsl-1.9/poly/Makefile.in @@ -0,0 +1,546 @@ +# 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 = poly +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) +libgslpoly_la_LIBADD = +am_libgslpoly_la_OBJECTS = dd.lo eval.lo solve_quadratic.lo \ + solve_cubic.lo zsolve_quadratic.lo zsolve_cubic.lo zsolve.lo \ + zsolve_init.lo +libgslpoly_la_OBJECTS = $(am_libgslpoly_la_OBJECTS) +am_test_OBJECTS = test.$(OBJEXT) +test_OBJECTS = $(am_test_OBJECTS) +test_DEPENDENCIES = libgslpoly.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 = $(libgslpoly_la_SOURCES) $(test_SOURCES) +DIST_SOURCES = $(libgslpoly_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 = libgslpoly.la +pkginclude_HEADERS = gsl_poly.h +INCLUDES = -I$(top_builddir) +libgslpoly_la_SOURCES = dd.c eval.c solve_quadratic.c solve_cubic.c zsolve_quadratic.c zsolve_cubic.c zsolve.c zsolve_init.c +noinst_HEADERS = balance.c companion.c qr.c +TESTS = $(check_PROGRAMS) +test_SOURCES = test.c +test_LDADD = libgslpoly.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 poly/Makefile'; \ + cd $(top_srcdir) && \ + $(AUTOMAKE) --gnu --ignore-deps poly/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 +libgslpoly.la: $(libgslpoly_la_OBJECTS) $(libgslpoly_la_DEPENDENCIES) + $(LINK) $(libgslpoly_la_LDFLAGS) $(libgslpoly_la_OBJECTS) $(libgslpoly_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/poly/TODO b/gsl-1.9/poly/TODO new file mode 100644 index 0000000..dbbc4c7 --- /dev/null +++ b/gsl-1.9/poly/TODO @@ -0,0 +1,130 @@ +* Estimate error on general poly roots using Newton method? Allow for +multiple roots and higher derivatives + +* Newton-Maehly (Newton with implicit deflation) + +* Jenkins-Traub + +* Brian Smith's adaptation of Laguerre's method + +* Hirano's method, SIAM J Num Anal 19 (1982) 793-99 by Murota + +* Carstensen, Petkovic, "On iteration methods without derivatives for +the simultaneous determination of polynomial zeros", J. Comput. Appl. +Math., 45 (1993) 251-267 + +* Investigate this, + + > NA Digest Sunday, July 11, 1999 Volume 99 : Issue 28 + > + > From: Murakami Hiroshi <nadigest@tmca.ac.jp> + > Date: Sun, 11 Jul 1999 18:56:54 +0900 (JST) + > Subject: Code for Wilf's Complex Bisection Method + > + > A sample demo source of root finding method for the general complex + > coefficient polynomial is placed to URL <ftp://ftp.tmca.ac.jp/wilf.taz>. + > It is about 8KB in size and is a tar and gnu-zipped file. + > The algorithm is taken from the following reference: + > HERBERT S.WILF,"A Global Bisection Algorithm for Computing the Zeros + > of Polynomials in the Complex Plane",ACM.vol.25,No.3,July 1978,pp.415-420. + > + > The Wilf's method is the complex plane version of the Sturm bi-section + > method for the real polynomial. And theoretically very interesting. + > For the given complex interval (complex rectilinear region and sides are + > parallel to the x-y axis), the procedure gives the count of the complex + > roots of the polynomial inside the interval. Thus, successive bi-section + > or quad-section of the interval can give the sequence of the shrinking + > intervals containing the root inside to attain the required accuracies. + > The source code is written mostly Fortran77 language, however some + > violations are intensionally left as: 1: The DO-ENDDO loop structure. + > 2: The longer than 6 letter names. 3: The use of under-score letter. + > 4: The use of include statement. + > The code will be placed in the public domain. + > + +* Investigate this + +From: "Steven G. Johnson" <stevenj@alum.mit.edu> +To: help-gsl@gnu.org +Subject: [Help-gsl] (in)accuracy of gsl_poly_complex_solve for repeated + roots? +Date: Sun, 05 Jun 2005 16:25:40 -0400 +Precedence: list +Envelope-to: bjg@network-theory.co.uk + +Hi, I noticed that gsl_poly_complex_solve seems to be surprisingly +inaccurate. For example, if you ask it for the roots of 1 + 4x + 6x^2 + +4x^3 + x^4, which should have x = -1 as a four-fold root (note that the +coefficients and solutions are exactly representable), it gives roots: + +-0.999903+9.66605e-05i +-0.999903-9.66605e-05i +-1.0001+9.66834e-05i +-1.0001-9.66834e-05i + +i.e. it is accurate to only 4 significant digits. (On the other hand, +when I have 4 distinct real roots it seems to be accurate to machine +precision.) + +If this kind of catastrophic accuracy loss is intrinsic to the algorithm +when repeated roots are encountered, please note it in the manual. + +However, I suspect that there may be algorithms to obtain higher +accuracy for multiple roots. I found the below references in a +literature search on the topic, which you may want to look into. (The +first reference can be found online at +http://www.neiu.edu/~zzeng/multroot.htm) + +Cordially, +Steven G. Johnson + +--------------------------------------------------------------------- +Algorithm 835: MULTROOT - a Matlab package for computing polynomial +roots and multiplicities +Zeng, Z. (Dept. of Math., Northeastern Illinois Univ., Chicago, IL, USA) +Source: ACM Transactions on Mathematical Software, v 30, n 2, June 2004, +p 218-36 +ISSN: 0098-3500 CODEN: ACMSCU +Publisher: ACM, USA + +Abstract: MULTROOT is a collection of Matlab modules for accurate +computation of polynomial roots, especially roots with nontrivial +multiplicities. As a blackbox-type software, MULTROOT requires the +polynomial coefficients as the only input, and outputs the computed +roots, multiplicities, backward error, estimated forward error, and the +structure-preserving condition number. The most significant features of +MULTROOT are the multiplicity identification capability and high +accuracy on multiple roots without using multiprecision arithmetic, even +if the polynomial coefficients are inexact. A comprehensive test suite +of polynomials that are collected from the literature is included for +numerical experiments and performance comparison (21 refs.) + +--------------------------------------------------------------------- +Ten methods to bound multiple roots of polynomials +Rump, S.M. (Inst. fur Informatik III, Tech. Univ. Hamburg-Harburg, +Hamburg, Germany) Source: Journal of Computational and Applied +Mathematics, v 156, n 2, 15 July 2003, p 403-32 +ISSN: 0377-0427 CODEN: JCAMDI +Publisher: Elsevier, Netherlands + +Abstract: Given a univariate polynomial P with a k-fold multiple root or +a k-fold root cluster near some z, we discuss nine different methods to +compute a disc near z which either contains exactly or contains at least +k roots of P. Many of the presented methods are known and of those some +are new. We are especially interested in the behavior of methods when +implemented in a rigorous way, that is, when taking into account all +possible effects of rounding errors. In other words, every result shall +be mathematically correct. We display extensive test sets comparing the +methods under different circumstances. Based on the results, we present +a tenth, hybrid method combining five of the previous methods which, for +give z, (i) detects the number k of roots near z and (ii) computes an +including disc with in most cases a radius of the order of the numerical +sensitivity of the root cluster. Therefore, the resulting discs are +numerically nearly optimal + + + +_______________________________________________ +Help-gsl mailing list +Help-gsl@gnu.org +http://lists.gnu.org/mailman/listinfo/help-gsl diff --git a/gsl-1.9/poly/balance.c b/gsl-1.9/poly/balance.c new file mode 100644 index 0000000..a1f7916 --- /dev/null +++ b/gsl-1.9/poly/balance.c @@ -0,0 +1,127 @@ +/* poly/balance.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 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. + */ + +static void balance_companion_matrix (double *m, size_t n); + +#define RADIX 2 +#define RADIX2 (RADIX*RADIX) + +static void +balance_companion_matrix (double *m, size_t nc) +{ + int not_converged = 1; + + double row_norm = 0; + double col_norm = 0; + + while (not_converged) + { + size_t i, j; + double g, f, s; + + not_converged = 0; + + for (i = 0; i < nc; i++) + { + /* column norm, excluding the diagonal */ + + if (i != nc - 1) + { + col_norm = fabs (MAT (m, i + 1, i, nc)); + } + else + { + col_norm = 0; + + for (j = 0; j < nc - 1; j++) + { + col_norm += fabs (MAT (m, j, nc - 1, nc)); + } + } + + /* row norm, excluding the diagonal */ + + if (i == 0) + { + row_norm = fabs (MAT (m, 0, nc - 1, nc)); + } + else if (i == nc - 1) + { + row_norm = fabs (MAT (m, i, i - 1, nc)); + } + else + { + row_norm = (fabs (MAT (m, i, i - 1, nc)) + + fabs (MAT (m, i, nc - 1, nc))); + } + + if (col_norm == 0 || row_norm == 0) + { + continue; + } + + g = row_norm / RADIX; + f = 1; + s = col_norm + row_norm; + + while (col_norm < g) + { + f *= RADIX; + col_norm *= RADIX2; + } + + g = row_norm * RADIX; + + while (col_norm > g) + { + f /= RADIX; + col_norm /= RADIX2; + } + + if ((row_norm + col_norm) < 0.95 * s * f) + { + not_converged = 1; + + g = 1 / f; + + if (i == 0) + { + MAT (m, 0, nc - 1, nc) *= g; + } + else + { + MAT (m, i, i - 1, nc) *= g; + MAT (m, i, nc - 1, nc) *= g; + } + + if (i == nc - 1) + { + for (j = 0; j < nc; j++) + { + MAT (m, j, i, nc) *= f; + } + } + else + { + MAT (m, i + 1, i, nc) *= f; + } + } + } + } +} diff --git a/gsl-1.9/poly/companion.c b/gsl-1.9/poly/companion.c new file mode 100644 index 0000000..736483f --- /dev/null +++ b/gsl-1.9/poly/companion.c @@ -0,0 +1,36 @@ +/* poly/companion.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 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. + */ + +static void set_companion_matrix (const double *a, size_t n, double *m); + +static void +set_companion_matrix (const double *a, size_t nc, double *m) +{ + size_t i, j; + + for (i = 0; i < nc; i++) + for (j = 0; j < nc; j++) + MAT (m, i, j, nc) = 0.0; + + for (i = 1; i < nc; i++) + MAT (m, i, i - 1, nc) = 1.0; + + for (i = 0; i < nc; i++) + MAT (m, i, nc - 1, nc) = -a[i] / a[nc]; +} diff --git a/gsl-1.9/poly/dd.c b/gsl-1.9/poly/dd.c new file mode 100644 index 0000000..47bff2c --- /dev/null +++ b/gsl-1.9/poly/dd.c @@ -0,0 +1,101 @@ +/* interpolation/interp_poly.c + * + * Copyright (C) 2001 DAN, HO-JIN + * + * 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. + */ + +/* Modified for standalone use in polynomial directory, B.Gough 2001 */ + +#include <config.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_poly.h> + +int +gsl_poly_dd_init (double dd[], const double xa[], const double ya[], + size_t size) +{ + size_t i, j; + + /* Newton's divided differences */ + + dd[0] = ya[0]; + + for (j = size - 1; j >= 1; j--) + { + dd[j] = (ya[j] - ya[j - 1]) / (xa[j] - xa[j - 1]); + } + + for (i = 2; i < size; i++) + { + for (j = size - 1; j >= i; j--) + { + dd[j] = (dd[j] - dd[j - 1]) / (xa[j] - xa[j - i]); + } + } + + return GSL_SUCCESS; +} + +#ifndef HIDE_INLINE_STATIC +double +gsl_poly_dd_eval (const double dd[], const double xa[], const size_t size, const double x) +{ + size_t i; + double y = dd[size - 1]; + + for (i = size - 1; i--;) + { + y = dd[i] + (x - xa[i]) * y; + } + + return y; +} +#endif + +int +gsl_poly_dd_taylor (double c[], double xp, + const double dd[], const double xa[], size_t size, + double w[]) +{ + size_t i, j; + + for (i = 0; i < size; i++) + { + c[i] = 0.0; + w[i] = 0.0; + } + + w[size - 1] = 1.0; + + c[0] = dd[0]; + + for (i = size - 1; i > 0 && i--;) + { + w[i] = -w[i + 1] * (xa[size - 2 - i] - xp); + + for (j = i + 1; j < size - 1; j++) + { + w[j] = w[j] - w[j + 1] * (xa[size - 2 - i] - xp); + } + + for (j = i; j < size; j++) + { + c[j - i] += w[j] * dd[size - i - 1]; + } + } + + return GSL_SUCCESS; +} diff --git a/gsl-1.9/poly/eval.c b/gsl-1.9/poly/eval.c new file mode 100644 index 0000000..84f6120 --- /dev/null +++ b/gsl-1.9/poly/eval.c @@ -0,0 +1,36 @@ +/* poly/eval.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman + * + * 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_poly.h> + +/*-*-*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*-*/ + +/* checked OK [GJ] Tue May 5 12:19:56 MDT 1998 */ + +#ifndef HIDE_INLINE_STATIC +double +gsl_poly_eval(const double c[], const int len, const double x) +{ + int i; + double ans = c[len-1]; + for(i=len-1; i>0; i--) ans = c[i-1] + x * ans; + return ans; +} +#endif diff --git a/gsl-1.9/poly/gsl_poly.h b/gsl-1.9/poly/gsl_poly.h new file mode 100644 index 0000000..57eda42 --- /dev/null +++ b/gsl-1.9/poly/gsl_poly.h @@ -0,0 +1,132 @@ +/* poly/gsl_poly.h + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004 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_POLY_H__ +#define __GSL_POLY_H__ + +#include <stdlib.h> +#include <gsl/gsl_complex.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 + + +/* Evaluate polynomial + * + * c[0] + c[1] x + c[2] x^2 + ... + c[len-1] x^(len-1) + * + * exceptions: none + */ +double gsl_poly_eval(const double c[], const int len, const double x); + + +#ifdef HAVE_INLINE +extern inline +double gsl_poly_eval(const double c[], const int len, const double x) +{ + int i; + double ans = c[len-1]; + for(i=len-1; i>0; i--) ans = c[i-1] + x * ans; + return ans; +} +#endif /* HAVE_INLINE */ + +/* Work with divided-difference polynomials, Abramowitz & Stegun 25.2.26 */ + +int +gsl_poly_dd_init (double dd[], const double x[], const double y[], + size_t size); + +double +gsl_poly_dd_eval (const double dd[], const double xa[], const size_t size, const double x); + +#ifdef HAVE_INLINE +extern inline +double gsl_poly_dd_eval(const double dd[], const double xa[], const size_t size, const double x) +{ + size_t i; + double y = dd[size - 1]; + for (i = size - 1; i--;) y = dd[i] + (x - xa[i]) * y; + return y; +} +#endif /* HAVE_INLINE */ + + +int +gsl_poly_dd_taylor (double c[], double xp, + const double dd[], const double x[], size_t size, + double w[]); + +/* Solve for real or complex roots of the standard quadratic equation, + * returning the number of real roots. + * + * Roots are returned ordered. + */ +int gsl_poly_solve_quadratic (double a, double b, double c, + double * x0, double * x1); + +int +gsl_poly_complex_solve_quadratic (double a, double b, double c, + gsl_complex * z0, gsl_complex * z1); + + +/* Solve for real roots of the cubic equation + * x^3 + a x^2 + b x + c = 0, returning the + * number of real roots. + * + * Roots are returned ordered. + */ +int gsl_poly_solve_cubic (double a, double b, double c, + double * x0, double * x1, double * x2); + +int +gsl_poly_complex_solve_cubic (double a, double b, double c, + gsl_complex * z0, gsl_complex * z1, + gsl_complex * z2); + + +/* Solve for the complex roots of a general real polynomial */ + +typedef struct +{ + size_t nc ; + double * matrix ; +} +gsl_poly_complex_workspace ; + +gsl_poly_complex_workspace * gsl_poly_complex_workspace_alloc (size_t n); +void gsl_poly_complex_workspace_free (gsl_poly_complex_workspace * w); + +int +gsl_poly_complex_solve (const double * a, size_t n, + gsl_poly_complex_workspace * w, + gsl_complex_packed_ptr z); + +__END_DECLS + +#endif /* __GSL_POLY_H__ */ diff --git a/gsl-1.9/poly/qr.c b/gsl-1.9/poly/qr.c new file mode 100644 index 0000000..23e7bfb --- /dev/null +++ b/gsl-1.9/poly/qr.c @@ -0,0 +1,245 @@ +/* poly/qr.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 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. + */ + +static int qr_companion (double *h, size_t nc, gsl_complex_packed_ptr z); + +static int +qr_companion (double *h, size_t nc, gsl_complex_packed_ptr zroot) +{ + double t = 0.0; + + size_t iterations, e, i, j, k, m; + + double w, x, y, s, z; + + double p = 0, q = 0, r = 0; + + /* FIXME: if p,q,r, are not set to zero then the compiler complains + that they ``might be used uninitialized in this + function''. Looking at the code this does seem possible, so this + should be checked. */ + + int notlast; + + size_t n = nc; + +next_root: + + if (n == 0) + return GSL_SUCCESS ; + + iterations = 0; + +next_iteration: + + for (e = n; e >= 2; e--) + { + double a1 = fabs (FMAT (h, e, e - 1, nc)); + double a2 = fabs (FMAT (h, e - 1, e - 1, nc)); + double a3 = fabs (FMAT (h, e, e, nc)); + + if (a1 <= GSL_DBL_EPSILON * (a2 + a3)) + break; + } + + x = FMAT (h, n, n, nc); + + if (e == n) + { + GSL_SET_COMPLEX_PACKED (zroot, n-1, x + t, 0); /* one real root */ + n--; + goto next_root; + /*continue;*/ + } + + y = FMAT (h, n - 1, n - 1, nc); + w = FMAT (h, n - 1, n, nc) * FMAT (h, n, n - 1, nc); + + if (e == n - 1) + { + p = (y - x) / 2; + q = p * p + w; + y = sqrt (fabs (q)); + + x += t; + + if (q > 0) /* two real roots */ + { + if (p < 0) + y = -y; + y += p; + + GSL_SET_COMPLEX_PACKED (zroot, n-1, x - w / y, 0); + GSL_SET_COMPLEX_PACKED (zroot, n-2, x + y, 0); + } + else + { + GSL_SET_COMPLEX_PACKED (zroot, n-1, x + p, -y); + GSL_SET_COMPLEX_PACKED (zroot, n-2, x + p, y); + } + n -= 2; + + goto next_root; + /*continue;*/ + } + + /* No more roots found yet, do another iteration */ + + if (iterations == 60) /* increased from 30 to 60 */ + { + /* too many iterations - give up! */ + + return GSL_FAILURE ; + } + + if (iterations % 10 == 0 && iterations > 0) + { + /* use an exceptional shift */ + + t += x; + + for (i = 1; i <= n; i++) + { + FMAT (h, i, i, nc) -= x; + } + + s = fabs (FMAT (h, n, n - 1, nc)) + fabs (FMAT (h, n - 1, n - 2, nc)); + y = 0.75 * s; + x = y; + w = -0.4375 * s * s; + } + + iterations++; + + for (m = n - 2; m >= e; m--) + { + double a1, a2, a3; + + z = FMAT (h, m, m, nc); + r = x - z; + s = y - z; + p = FMAT (h, m, m + 1, nc) + (r * s - w) / FMAT (h, m + 1, m, nc); + q = FMAT (h, m + 1, m + 1, nc) - z - r - s; + r = FMAT (h, m + 2, m + 1, nc); + s = fabs (p) + fabs (q) + fabs (r); + p /= s; + q /= s; + r /= s; + + if (m == e) + break; + + a1 = fabs (FMAT (h, m, m - 1, nc)); + a2 = fabs (FMAT (h, m - 1, m - 1, nc)); + a3 = fabs (FMAT (h, m + 1, m + 1, nc)); + + if (a1 * (fabs (q) + fabs (r)) <= GSL_DBL_EPSILON * fabs (p) * (a2 + a3)) + break; + } + + for (i = m + 2; i <= n; i++) + { + FMAT (h, i, i - 2, nc) = 0; + } + + for (i = m + 3; i <= n; i++) + { + FMAT (h, i, i - 3, nc) = 0; + } + + /* double QR step */ + + for (k = m; k <= n - 1; k++) + { + notlast = (k != n - 1); + + if (k != m) + { + p = FMAT (h, k, k - 1, nc); + q = FMAT (h, k + 1, k - 1, nc); + r = notlast ? FMAT (h, k + 2, k - 1, nc) : 0.0; + + x = fabs (p) + fabs (q) + fabs (r); + + if (x == 0) + continue; /* FIXME????? */ + + p /= x; + q /= x; + r /= x; + } + + s = sqrt (p * p + q * q + r * r); + + if (p < 0) + s = -s; + + if (k != m) + { + FMAT (h, k, k - 1, nc) = -s * x; + } + else if (e != m) + { + FMAT (h, k, k - 1, nc) *= -1; + } + + p += s; + x = p / s; + y = q / s; + z = r / s; + q /= p; + r /= p; + + /* do row modifications */ + + for (j = k; j <= n; j++) + { + p = FMAT (h, k, j, nc) + q * FMAT (h, k + 1, j, nc); + + if (notlast) + { + p += r * FMAT (h, k + 2, j, nc); + FMAT (h, k + 2, j, nc) -= p * z; + } + + FMAT (h, k + 1, j, nc) -= p * y; + FMAT (h, k, j, nc) -= p * x; + } + + j = (k + 3 < n) ? (k + 3) : n; + + /* do column modifications */ + + for (i = e; i <= j; i++) + { + p = x * FMAT (h, i, k, nc) + y * FMAT (h, i, k + 1, nc); + + if (notlast) + { + p += z * FMAT (h, i, k + 2, nc); + FMAT (h, i, k + 2, nc) -= p * r; + } + FMAT (h, i, k + 1, nc) -= p * q; + FMAT (h, i, k, nc) -= p; + } + } + + goto next_iteration; +} + diff --git a/gsl-1.9/poly/solve_cubic.c b/gsl-1.9/poly/solve_cubic.c new file mode 100644 index 0000000..979e56a --- /dev/null +++ b/gsl-1.9/poly/solve_cubic.c @@ -0,0 +1,110 @@ +/* poly/solve_cubic.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 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. + */ + +/* solve_cubic.c - finds the real roots of x^3 + a x^2 + b x + c = 0 */ + +#include <config.h> +#include <math.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_poly.h> + +#define SWAP(a,b) do { double tmp = b ; b = a ; a = tmp ; } while(0) + +int +gsl_poly_solve_cubic (double a, double b, double c, + double *x0, double *x1, double *x2) +{ + double q = (a * a - 3 * b); + double r = (2 * a * a * a - 9 * a * b + 27 * c); + + double Q = q / 9; + double R = r / 54; + + double Q3 = Q * Q * Q; + double R2 = R * R; + + double CR2 = 729 * r * r; + double CQ3 = 2916 * q * q * q; + + if (R == 0 && Q == 0) + { + *x0 = - a / 3 ; + *x1 = - a / 3 ; + *x2 = - a / 3 ; + return 3 ; + } + else if (CR2 == CQ3) + { + /* this test is actually R2 == Q3, written in a form suitable + for exact computation with integers */ + + /* Due to finite precision some double roots may be missed, and + considered to be a pair of complex roots z = x +/- epsilon i + close to the real axis. */ + + double sqrtQ = sqrt (Q); + + if (R > 0) + { + *x0 = -2 * sqrtQ - a / 3; + *x1 = sqrtQ - a / 3; + *x2 = sqrtQ - a / 3; + } + else + { + *x0 = - sqrtQ - a / 3; + *x1 = - sqrtQ - a / 3; + *x2 = 2 * sqrtQ - a / 3; + } + return 3 ; + } + else if (CR2 < CQ3) /* equivalent to R2 < Q3 */ + { + double sqrtQ = sqrt (Q); + double sqrtQ3 = sqrtQ * sqrtQ * sqrtQ; + double theta = acos (R / sqrtQ3); + double norm = -2 * sqrtQ; + *x0 = norm * cos (theta / 3) - a / 3; + *x1 = norm * cos ((theta + 2.0 * M_PI) / 3) - a / 3; + *x2 = norm * cos ((theta - 2.0 * M_PI) / 3) - a / 3; + + /* Sort *x0, *x1, *x2 into increasing order */ + + if (*x0 > *x1) + SWAP(*x0, *x1) ; + + if (*x1 > *x2) + { + SWAP(*x1, *x2) ; + + if (*x0 > *x1) + SWAP(*x0, *x1) ; + } + + return 3; + } + else + { + double sgnR = (R >= 0 ? 1 : -1); + double A = -sgnR * pow (fabs (R) + sqrt (R2 - Q3), 1.0/3.0); + double B = Q / A ; + *x0 = A + B - a / 3; + return 1; + } +} diff --git a/gsl-1.9/poly/solve_quadratic.c b/gsl-1.9/poly/solve_quadratic.c new file mode 100644 index 0000000..8af69db --- /dev/null +++ b/gsl-1.9/poly/solve_quadratic.c @@ -0,0 +1,85 @@ +/* poly/solve_quadratic.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 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. + */ + +/* solve_quadratic.c - finds the real roots of a x^2 + b x + c = 0 */ + +#include <config.h> +#include <math.h> + +#include <gsl/gsl_poly.h> + +int +gsl_poly_solve_quadratic (double a, double b, double c, + double *x0, double *x1) +{ + double disc = b * b - 4 * a * c; + + if (a == 0) /* Handle linear case */ + { + if (b == 0) + { + return 0; + } + else + { + *x0 = -c / b; + return 1; + }; + } + + if (disc > 0) + { + if (b == 0) + { + double r = fabs (0.5 * sqrt (disc) / a); + *x0 = -r; + *x1 = r; + } + else + { + double sgnb = (b > 0 ? 1 : -1); + double temp = -0.5 * (b + sgnb * sqrt (disc)); + double r1 = temp / a ; + double r2 = c / temp ; + + if (r1 < r2) + { + *x0 = r1 ; + *x1 = r2 ; + } + else + { + *x0 = r2 ; + *x1 = r1 ; + } + } + return 2; + } + else if (disc == 0) + { + *x0 = -0.5 * b / a ; + *x1 = -0.5 * b / a ; + return 2 ; + } + else + { + return 0; + } +} + diff --git a/gsl-1.9/poly/test.c b/gsl-1.9/poly/test.c new file mode 100644 index 0000000..bd0e141 --- /dev/null +++ b/gsl-1.9/poly/test.c @@ -0,0 +1,508 @@ +/* poly/test.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 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_ieee_utils.h> +#include <gsl/gsl_poly.h> + +int +main (void) +{ + const double eps = 100.0 * GSL_DBL_EPSILON; + + gsl_ieee_env_setup (); + + { + double x, y; + double c[3] = { 1.0, 0.5, 0.3 }; + x = 0.5; + y = gsl_poly_eval (c, 3, x); + gsl_test_rel (y, 1 + 0.5 * x + 0.3 * x * x, eps, + "gsl_poly_eval({1, 0.5, 0.3}, 0.5)"); + } + + { + double x, y; + double d[11] = { 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1 }; + x = 1.0; + y = gsl_poly_eval (d, 11, x); + gsl_test_rel (y, 1.0, eps, + "gsl_poly_eval({1,-1, 1, -1, 1, -1, 1, -1, 1, -1, 1}, 1.0)"); + + } + + /* Quadratic */ + + { + double x0, x1; + + int n = gsl_poly_solve_quadratic (4.0, -20.0, 26.0, &x0, &x1); + + gsl_test (n != 0, "gsl_poly_solve_quadratic, no roots, (2x - 5)^2 = -1"); + } + + { + double x0, x1; + + int n = gsl_poly_solve_quadratic (4.0, -20.0, 25.0, &x0, &x1); + + gsl_test (n != 2, "gsl_poly_solve_quadratic, one root, (2x - 5)^2 = 0"); + gsl_test_rel (x0, 2.5, 1e-9, "x0, (2x - 5)^2 = 0"); + gsl_test_rel (x1, 2.5, 1e-9, "x1, (2x - 5)^2 = 0"); + gsl_test (x0 != x1, "x0 == x1, (2x - 5)^2 = 0"); + } + + { + double x0, x1; + + int n = gsl_poly_solve_quadratic (4.0, -20.0, 21.0, &x0, &x1); + + gsl_test (n != 2, "gsl_poly_solve_quadratic, two roots, (2x - 5)^2 = 4"); + gsl_test_rel (x0, 1.5, 1e-9, "x0, (2x - 5)^2 = 4"); + gsl_test_rel (x1, 3.5, 1e-9, "x1, (2x - 5)^2 = 4"); + } + + { + double x0, x1; + + int n = gsl_poly_solve_quadratic (4.0, 7.0, 0.0, &x0, &x1); + + gsl_test (n != 2, "gsl_poly_solve_quadratic, two roots, x(4x + 7) = 0"); + gsl_test_rel (x0, -1.75, 1e-9, "x0, x(4x + 7) = 0"); + gsl_test_rel (x1, 0.0, 1e-9, "x1, x(4x + 7) = 0"); + } + + { + double x0, x1; + + int n = gsl_poly_solve_quadratic (5.0, 0.0, -20.0, &x0, &x1); + + gsl_test (n != 2, + "gsl_poly_solve_quadratic, two roots b = 0, 5 x^2 = 20"); + gsl_test_rel (x0, -2.0, 1e-9, "x0, 5 x^2 = 20"); + gsl_test_rel (x1, 2.0, 1e-9, "x1, 5 x^2 = 20"); + } + + + { + double x0, x1; + + int n = gsl_poly_solve_quadratic (0.0, 3.0, -21.0, &x0, &x1); + + gsl_test (n != 1, + "gsl_poly_solve_quadratic, one root (linear) 3 x - 21 = 0"); + gsl_test_rel (x0, 7.0, 1e-9, "x0, 3x - 21 = 0"); + } + + + { + double x0, x1; + int n = gsl_poly_solve_quadratic (0.0, 0.0, 1.0, &x0, &x1); + + gsl_test (n != 0, + "gsl_poly_solve_quadratic, no roots 1 = 0"); + } + + + /* Cubic */ + + { + double x0, x1, x2; + + int n = gsl_poly_solve_cubic (0.0, 0.0, -27.0, &x0, &x1, &x2); + + gsl_test (n != 1, "gsl_poly_solve_cubic, one root, x^3 = 27"); + gsl_test_rel (x0, 3.0, 1e-9, "x0, x^3 = 27"); + } + + { + double x0, x1, x2; + + int n = gsl_poly_solve_cubic (-51.0, 867.0, -4913.0, &x0, &x1, &x2); + + gsl_test (n != 3, "gsl_poly_solve_cubic, three roots, (x-17)^3=0"); + gsl_test_rel (x0, 17.0, 1e-9, "x0, (x-17)^3=0"); + gsl_test_rel (x1, 17.0, 1e-9, "x1, (x-17)^3=0"); + gsl_test_rel (x2, 17.0, 1e-9, "x2, (x-17)^3=0"); + } + + { + double x0, x1, x2; + + int n = gsl_poly_solve_cubic (-57.0, 1071.0, -6647.0, &x0, &x1, &x2); + + gsl_test (n != 3, + "gsl_poly_solve_cubic, three roots, (x-17)(x-17)(x-23)=0"); + gsl_test_rel (x0, 17.0, 1e-9, "x0, (x-17)(x-17)(x-23)=0"); + gsl_test_rel (x1, 17.0, 1e-9, "x1, (x-17)(x-17)(x-23)=0"); + gsl_test_rel (x2, 23.0, 1e-9, "x2, (x-17)(x-17)(x-23)=0"); + } + + { + double x0, x1, x2; + + int n = gsl_poly_solve_cubic (-11.0, -493.0, +6647.0, &x0, &x1, &x2); + + gsl_test (n != 3, + "gsl_poly_solve_cubic, three roots, (x+23)(x-17)(x-17)=0"); + gsl_test_rel (x0, -23.0, 1e-9, "x0, (x+23)(x-17)(x-17)=0"); + gsl_test_rel (x1, 17.0, 1e-9, "x1, (x+23)(x-17)(x-17)=0"); + gsl_test_rel (x2, 17.0, 1e-9, "x2, (x+23)(x-17)(x-17)=0"); + } + + { + double x0, x1, x2; + + int n = gsl_poly_solve_cubic (-143.0, 5087.0, -50065.0, &x0, &x1, &x2); + + gsl_test (n != 3, + "gsl_poly_solve_cubic, three roots, (x-17)(x-31)(x-95)=0"); + gsl_test_rel (x0, 17.0, 1e-9, "x0, (x-17)(x-31)(x-95)=0"); + gsl_test_rel (x1, 31.0, 1e-9, "x1, (x-17)(x-31)(x-95)=0"); + gsl_test_rel (x2, 95.0, 1e-9, "x2, (x-17)(x-31)(x-95)=0"); + } + + { + double x0, x1, x2; + + int n = gsl_poly_solve_cubic (-109.0, 803.0, 50065.0, &x0, &x1, &x2); + + gsl_test (n != 3, + "gsl_poly_solve_cubic, three roots, (x+17)(x-31)(x-95)=0"); + gsl_test_rel (x0, -17.0, 1e-9, "x0, (x+17)(x-31)(x-95)=0"); + gsl_test_rel (x1, 31.0, 1e-9, "x1, (x+17)(x-31)(x-95)=0"); + gsl_test_rel (x2, 95.0, 1e-9, "x2, (x+17)(x-31)(x-95)=0"); + } + + /* Quadratic with complex roots */ + + { + gsl_complex z0, z1; + + int n = gsl_poly_complex_solve_quadratic (4.0, -20.0, 26.0, &z0, &z1); + + gsl_test (n != 2, + "gsl_poly_complex_solve_quadratic, 2 roots (2x - 5)^2 = -1"); + gsl_test_rel (GSL_REAL (z0), 2.5, 1e-9, "z0.real, (2x - 5)^2 = -1"); + gsl_test_rel (GSL_IMAG (z0), -0.5, 1e-9, "z0.imag, (2x - 5)^2 = -1"); + + gsl_test_rel (GSL_REAL (z1), 2.5, 1e-9, "z1.real, (2x - 5)^2 = -1"); + gsl_test_rel (GSL_IMAG (z1), 0.5, 1e-9, "z1.imag, (2x - 5)^2 = -1"); + } + + { + gsl_complex z0, z1; + + int n = gsl_poly_complex_solve_quadratic (4.0, -20.0, 25.0, &z0, &z1); + + gsl_test (n != 2, + "gsl_poly_complex_solve_quadratic, one root, (2x - 5)^2 = 0"); + gsl_test_rel (GSL_REAL (z0), 2.5, 1e-9, "z0.real, (2x - 5)^2 = 0"); + gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag (2x - 5)^2 = 0"); + gsl_test_rel (GSL_REAL (z1), 2.5, 1e-9, "z1.real, (2x - 5)^2 = 0"); + gsl_test_rel (GSL_IMAG (z1), 0.0, 1e-9, "z1.imag (2x - 5)^2 = 0"); + gsl_test (GSL_REAL (z0) != GSL_REAL (z1), + "z0.real == z1.real, (2x - 5)^2 = 0"); + gsl_test (GSL_IMAG (z0) != GSL_IMAG (z1), + "z0.imag == z1.imag, (2x - 5)^2 = 0"); + } + + { + gsl_complex z0, z1; + + int n = gsl_poly_complex_solve_quadratic (4.0, -20.0, 21.0, &z0, &z1); + + gsl_test (n != 2, + "gsl_poly_complex_solve_quadratic, two roots, (2x - 5)^2 = 4"); + gsl_test_rel (GSL_REAL (z0), 1.5, 1e-9, "z0.real, (2x - 5)^2 = 4"); + gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag, (2x - 5)^2 = 4"); + gsl_test_rel (GSL_REAL (z1), 3.5, 1e-9, "z1.real, (2x - 5)^2 = 4"); + gsl_test_rel (GSL_IMAG (z1), 0.0, 1e-9, "z1.imag, (2x - 5)^2 = 4"); + } + + { + gsl_complex z0, z1; + + int n = gsl_poly_complex_solve_quadratic (4.0, 7.0, 0.0, &z0, &z1); + + gsl_test (n != 2, + "gsl_poly_complex_solve_quadratic, two roots, x(4x + 7) = 0"); + gsl_test_rel (GSL_REAL (z0), -1.75, 1e-9, "z0.real, x(4x + 7) = 0"); + gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag, x(4x + 7) = 0"); + gsl_test_rel (GSL_REAL (z1), 0.0, 1e-9, "z1.real, x(4x + 7) = 0"); + gsl_test_rel (GSL_IMAG (z1), 0.0, 1e-9, "z1.imag, x(4x + 7) = 0"); + } + + { + gsl_complex z0, z1; + + int n = gsl_poly_complex_solve_quadratic (5.0, 0.0, -20.0, &z0, &z1); + + gsl_test (n != 2, + "gsl_poly_complex_solve_quadratic, two roots b = 0, 5 x^2 = 20"); + gsl_test_rel (GSL_REAL (z0), -2.0, 1e-9, "z0.real, 5 x^2 = 20"); + gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag, 5 x^2 = 20"); + gsl_test_rel (GSL_REAL (z1), 2.0, 1e-9, "z1.real, 5 x^2 = 20"); + gsl_test_rel (GSL_IMAG (z1), 0.0, 1e-9, "z1.imag, 5 x^2 = 20"); + } + + { + gsl_complex z0, z1; + + int n = gsl_poly_complex_solve_quadratic (5.0, 0.0, 20.0, &z0, &z1); + + gsl_test (n != 2, + "gsl_poly_complex_solve_quadratic, two roots b = 0, 5 x^2 = -20"); + gsl_test_rel (GSL_REAL (z0), 0.0, 1e-9, "z0.real, 5 x^2 = -20"); + gsl_test_rel (GSL_IMAG (z0), -2.0, 1e-9, "z0.imag, 5 x^2 = -20"); + gsl_test_rel (GSL_REAL (z1), 0.0, 1e-9, "z1.real, 5 x^2 = -20"); + gsl_test_rel (GSL_IMAG (z1), 2.0, 1e-9, "z1.imag, 5 x^2 = -20"); + } + + + { + gsl_complex z0, z1; + + int n = gsl_poly_complex_solve_quadratic (0.0, 3.0, -21.0, &z0, &z1); + + gsl_test (n != 1, + "gsl_poly_complex_solve_quadratic, one root (linear) 3 x - 21 = 0"); + + gsl_test_rel (GSL_REAL (z0), 7.0, 1e-9, "z0.real, 3x - 21 = 0"); + gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag, 3x - 21 = 0"); + } + + + { + gsl_complex z0, z1; + + int n = gsl_poly_complex_solve_quadratic (0.0, 0.0, 1.0, &z0, &z1); + gsl_test (n != 0, + "gsl_poly_complex_solve_quadratic, no roots 1 = 0"); + } + + + + /* Cubic with complex roots */ + + { + gsl_complex z0, z1, z2; + + int n = gsl_poly_complex_solve_cubic (0.0, 0.0, -27.0, &z0, &z1, &z2); + + gsl_test (n != 3, "gsl_poly_complex_solve_cubic, three root, x^3 = 27"); + gsl_test_rel (GSL_REAL (z0), -1.5, 1e-9, "z0.real, x^3 = 27"); + gsl_test_rel (GSL_IMAG (z0), -1.5 * sqrt (3.0), 1e-9, + "z0.imag, x^3 = 27"); + gsl_test_rel (GSL_REAL (z1), -1.5, 1e-9, "z1.real, x^3 = 27"); + gsl_test_rel (GSL_IMAG (z1), 1.5 * sqrt (3.0), 1e-9, "z1.imag, x^3 = 27"); + gsl_test_rel (GSL_REAL (z2), 3.0, 1e-9, "z2.real, x^3 = 27"); + gsl_test_rel (GSL_IMAG (z2), 0.0, 1e-9, "z2.imag, x^3 = 27"); + } + + { + gsl_complex z0, z1, z2; + + int n = gsl_poly_complex_solve_cubic (-1.0, 1.0, 39.0, &z0, &z1, &z2); + + gsl_test (n != 3, + "gsl_poly_complex_solve_cubic, three root, (x+3)(x^2-4x+13) = 0"); + gsl_test_rel (GSL_REAL (z0), -3.0, 1e-9, "z0.real, (x+3)(x^2+1) = 0"); + gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag, (x+3)(x^2+1) = 0"); + gsl_test_rel (GSL_REAL (z1), 2.0, 1e-9, "z1.real, (x+3)(x^2+1) = 0"); + gsl_test_rel (GSL_IMAG (z1), -3.0, 1e-9, "z1.imag, (x+3)(x^2+1) = 0"); + gsl_test_rel (GSL_REAL (z2), 2.0, 1e-9, "z2.real, (x+3)(x^2+1) = 0"); + gsl_test_rel (GSL_IMAG (z2), 3.0, 1e-9, "z2.imag, (x+3)(x^2+1) = 0"); + } + + { + gsl_complex z0, z1, z2; + + int n = + gsl_poly_complex_solve_cubic (-51.0, 867.0, -4913.0, &z0, &z1, &z2); + + gsl_test (n != 3, + "gsl_poly_complex_solve_cubic, three roots, (x-17)^3=0"); + gsl_test_rel (GSL_REAL (z0), 17.0, 1e-9, "z0.real, (x-17)^3=0"); + gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag, (x-17)^3=0"); + gsl_test_rel (GSL_REAL (z1), 17.0, 1e-9, "z1.real, (x-17)^3=0"); + gsl_test_rel (GSL_IMAG (z1), 0.0, 1e-9, "z1.imag, (x-17)^3=0"); + gsl_test_rel (GSL_REAL (z2), 17.0, 1e-9, "z2.real, (x-17)^3=0"); + gsl_test_rel (GSL_IMAG (z2), 0.0, 1e-9, "z2.imag, (x-17)^3=0"); + } + + { + gsl_complex z0, z1, z2; + + int n = + gsl_poly_complex_solve_cubic (-57.0, 1071.0, -6647.0, &z0, &z1, &z2); + + gsl_test (n != 3, + "gsl_poly_complex_solve_cubic, three roots, (x-17)(x-17)(x-23)=0"); + gsl_test_rel (GSL_REAL (z0), 17.0, 1e-9, "z0.real, (x-17)(x-17)(x-23)=0"); + gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag, (x-17)(x-17)(x-23)=0"); + gsl_test_rel (GSL_REAL (z1), 17.0, 1e-9, "z1.real, (x-17)(x-17)(x-23)=0"); + gsl_test_rel (GSL_IMAG (z1), 0.0, 1e-9, "z1.imag, (x-17)(x-17)(x-23)=0"); + gsl_test_rel (GSL_REAL (z2), 23.0, 1e-9, "z2.real, (x-17)(x-17)(x-23)=0"); + gsl_test_rel (GSL_IMAG (z2), 0.0, 1e-9, "z2.imag, (x-17)(x-17)(x-23)=0"); + } + + { + gsl_complex z0, z1, z2; + + int n = + gsl_poly_complex_solve_cubic (-11.0, -493.0, +6647.0, &z0, &z1, &z2); + + gsl_test (n != 3, + "gsl_poly_complex_solve_cubic, three roots, (x+23)(x-17)(x-17)=0"); + gsl_test_rel (GSL_REAL (z0), -23.0, 1e-9, + "z0.real, (x+23)(x-17)(x-17)=0"); + gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag, (x+23)(x-17)(x-17)=0"); + gsl_test_rel (GSL_REAL (z1), 17.0, 1e-9, "z1.real, (x+23)(x-17)(x-17)=0"); + gsl_test_rel (GSL_IMAG (z1), 0.0, 1e-9, "z1.imag, (x+23)(x-17)(x-17)=0"); + gsl_test_rel (GSL_REAL (z2), 17.0, 1e-9, "z2.real, (x+23)(x-17)(x-17)=0"); + gsl_test_rel (GSL_IMAG (z2), 0.0, 1e-9, "z2.imag, (x+23)(x-17)(x-17)=0"); + } + + + { + gsl_complex z0, z1, z2; + + int n = + gsl_poly_complex_solve_cubic (-143.0, 5087.0, -50065.0, &z0, &z1, &z2); + + gsl_test (n != 3, + "gsl_poly_complex_solve_cubic, three roots, (x-17)(x-31)(x-95)=0"); + gsl_test_rel (GSL_REAL (z0), 17.0, 1e-9, "z0.real, (x-17)(x-31)(x-95)=0"); + gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag, (x-17)(x-31)(x-95)=0"); + gsl_test_rel (GSL_REAL (z1), 31.0, 1e-9, "z1.real, (x-17)(x-31)(x-95)=0"); + gsl_test_rel (GSL_IMAG (z1), 0.0, 1e-9, "z1.imag, (x-17)(x-31)(x-95)=0"); + gsl_test_rel (GSL_REAL (z2), 95.0, 1e-9, "z2.real, (x-17)(x-31)(x-95)=0"); + gsl_test_rel (GSL_IMAG (z2), 0.0, 1e-9, "z2.imag, (x-17)(x-31)(x-95)=0"); + } + + + { + /* Wilkinson polynomial: y = (x-1)(x-2)(x-3)(x-4)(x-5) */ + + double a[6] = { -120, 274, -225, 85, -15, 1 }; + double z[6*2]; + + gsl_poly_complex_workspace *w = gsl_poly_complex_workspace_alloc (6); + + int status = gsl_poly_complex_solve (a, 6, w, z); + + gsl_poly_complex_workspace_free (w); + + gsl_test (status, + "gsl_poly_complex_solve, 5th-order Wilkinson polynomial"); + gsl_test_rel (z[0], 1.0, 1e-9, "z0.real, 5th-order polynomial"); + gsl_test_rel (z[1], 0.0, 1e-9, "z0.imag, 5th-order polynomial"); + gsl_test_rel (z[2], 2.0, 1e-9, "z1.real, 5th-order polynomial"); + gsl_test_rel (z[3], 0.0, 1e-9, "z1.imag, 5th-order polynomial"); + gsl_test_rel (z[4], 3.0, 1e-9, "z2.real, 5th-order polynomial"); + gsl_test_rel (z[5], 0.0, 1e-9, "z2.imag, 5th-order polynomial"); + gsl_test_rel (z[6], 4.0, 1e-9, "z3.real, 5th-order polynomial"); + gsl_test_rel (z[7], 0.0, 1e-9, "z3.imag, 5th-order polynomial"); + gsl_test_rel (z[8], 5.0, 1e-9, "z4.real, 5th-order polynomial"); + gsl_test_rel (z[9], 0.0, 1e-9, "z4.imag, 5th-order polynomial"); + } + + { + /* : 8-th order polynomial y = x^8 + x^4 + 1 */ + + double a[9] = { 1, 0, 0, 0, 1, 0, 0, 0, 1 }; + double z[8*2]; + + double C = 0.5; + double S = sqrt (3.0) / 2.0; + + gsl_poly_complex_workspace *w = gsl_poly_complex_workspace_alloc (9); + + int status = gsl_poly_complex_solve (a, 9, w, z); + + gsl_poly_complex_workspace_free (w); + + gsl_test (status, "gsl_poly_complex_solve, 8th-order polynomial"); + + gsl_test_rel (z[0], -S, 1e-9, "z0.real, 8th-order polynomial"); + gsl_test_rel (z[1], C, 1e-9, "z0.imag, 8th-order polynomial"); + gsl_test_rel (z[2], -S, 1e-9, "z1.real, 8th-order polynomial"); + gsl_test_rel (z[3], -C, 1e-9, "z1.imag, 8th-order polynomial"); + gsl_test_rel (z[4], -C, 1e-9, "z2.real, 8th-order polynomial"); + gsl_test_rel (z[5], S, 1e-9, "z2.imag, 8th-order polynomial"); + gsl_test_rel (z[6], -C, 1e-9, "z3.real, 8th-order polynomial"); + gsl_test_rel (z[7], -S, 1e-9, "z3.imag, 8th-order polynomial"); + gsl_test_rel (z[8], C, 1e-9, "z4.real, 8th-order polynomial"); + gsl_test_rel (z[9], S, 1e-9, "z4.imag, 8th-order polynomial"); + gsl_test_rel (z[10], C, 1e-9, "z5.real, 8th-order polynomial"); + gsl_test_rel (z[11], -S, 1e-9, "z5.imag, 8th-order polynomial"); + gsl_test_rel (z[12], S, 1e-9, "z6.real, 8th-order polynomial"); + gsl_test_rel (z[13], C, 1e-9, "z6.imag, 8th-order polynomial"); + gsl_test_rel (z[14], S, 1e-9, "z7.real, 8th-order polynomial"); + gsl_test_rel (z[15], -C, 1e-9, "z7.imag, 8th-order polynomial"); + + } + + { + int i; + + double xa[7] = {0.16, 0.97, 1.94, 2.74, 3.58, 3.73, 4.70 }; + double ya[7] = {0.73, 1.11, 1.49, 1.84, 2.30, 2.41, 3.07 }; + + double dd_expected[7] = { 7.30000000000000e-01, + 4.69135802469136e-01, + -4.34737219941284e-02, + 2.68681098870099e-02, + -3.22937056934996e-03, + 6.12763259971375e-03, + -6.45402453527083e-03 }; + + double dd[7], coeff[7], work[7]; + + gsl_poly_dd_init (dd, xa, ya, 7); + + for (i = 0; i < 7; i++) + { + gsl_test_rel (dd[i], dd_expected[i], 1e-10, "divided difference dd[%d]", i); + } + + for (i = 0; i < 7; i++) + { + double y = gsl_poly_dd_eval(dd, xa, 7, xa[i]); + gsl_test_rel (y, ya[i], 1e-10, "divided difference y[%d]", i); + } + + gsl_poly_dd_taylor (coeff, 1.5, dd, xa, 7, work); + + for (i = 0; i < 7; i++) + { + double y = gsl_poly_eval(coeff, 7, xa[i] - 1.5); + gsl_test_rel (y, ya[i], 1e-10, "taylor expansion about 1.5 y[%d]", i); + } + } + + + /* now summarize the results */ + + exit (gsl_test_summary ()); +} diff --git a/gsl-1.9/poly/zsolve.c b/gsl-1.9/poly/zsolve.c new file mode 100644 index 0000000..ae2aa51 --- /dev/null +++ b/gsl-1.9/poly/zsolve.c @@ -0,0 +1,85 @@ +/* poly/zsolve.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 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. + */ + +/* zsolve.c - finds the complex roots of = 0 */ + +#include <config.h> +#include <math.h> +#include <stdlib.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_complex.h> +#include <gsl/gsl_poly.h> + +/* C-style matrix elements */ +#define MAT(m,i,j,n) ((m)[(i)*(n) + (j)]) + +/* Fortran-style matrix elements */ +#define FMAT(m,i,j,n) ((m)[((i)-1)*(n) + ((j)-1)]) + +#include "companion.c" +#include "balance.c" +#include "qr.c" + +int +gsl_poly_complex_solve (const double *a, size_t n, + gsl_poly_complex_workspace * w, + gsl_complex_packed_ptr z) +{ + int status; + double *m; + + if (n == 0) + { + GSL_ERROR ("number of terms must be a positive integer", GSL_EINVAL); + } + + if (n == 1) + { + GSL_ERROR ("cannot solve for only one term", GSL_EINVAL); + } + + if (a[n - 1] == 0) + { + GSL_ERROR ("leading term of polynomial must be non-zero", GSL_EINVAL) ; + } + + if (w->nc != n - 1) + { + GSL_ERROR ("size of workspace does not match polynomial", GSL_EINVAL); + } + + m = w->matrix; + + set_companion_matrix (a, n - 1, m); + + balance_companion_matrix (m, n - 1); + + status = qr_companion (m, n - 1, z); + + if (status) + { + GSL_ERROR("root solving qr method failed to converge", GSL_EFAILED); + } + + return GSL_SUCCESS; +} + + + diff --git a/gsl-1.9/poly/zsolve_cubic.c b/gsl-1.9/poly/zsolve_cubic.c new file mode 100644 index 0000000..548423e --- /dev/null +++ b/gsl-1.9/poly/zsolve_cubic.c @@ -0,0 +1,153 @@ +/* poly/zsolve_cubic.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 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. + */ + +/* zsolve_cubic.c - finds the complex roots of x^3 + a x^2 + b x + c = 0 */ + +#include <config.h> +#include <math.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_complex.h> +#include <gsl/gsl_poly.h> + +#define SWAP(a,b) do { double tmp = b ; b = a ; a = tmp ; } while(0) + +int +gsl_poly_complex_solve_cubic (double a, double b, double c, + gsl_complex *z0, gsl_complex *z1, + gsl_complex *z2) +{ + double q = (a * a - 3 * b); + double r = (2 * a * a * a - 9 * a * b + 27 * c); + + double Q = q / 9; + double R = r / 54; + + double Q3 = Q * Q * Q; + double R2 = R * R; + + double CR2 = 729 * r * r; + double CQ3 = 2916 * q * q * q; + + if (R == 0 && Q == 0) + { + GSL_REAL (*z0) = -a / 3; + GSL_IMAG (*z0) = 0; + GSL_REAL (*z1) = -a / 3; + GSL_IMAG (*z1) = 0; + GSL_REAL (*z2) = -a / 3; + GSL_IMAG (*z2) = 0; + return 3; + } + else if (CR2 == CQ3) + { + /* this test is actually R2 == Q3, written in a form suitable + for exact computation with integers */ + + /* Due to finite precision some double roots may be missed, and + will be considered to be a pair of complex roots z = x +/- + epsilon i close to the real axis. */ + + double sqrtQ = sqrt (Q); + + if (R > 0) + { + GSL_REAL (*z0) = -2 * sqrtQ - a / 3; + GSL_IMAG (*z0) = 0; + GSL_REAL (*z1) = sqrtQ - a / 3; + GSL_IMAG (*z1) = 0; + GSL_REAL (*z2) = sqrtQ - a / 3; + GSL_IMAG (*z2) = 0; + } + else + { + GSL_REAL (*z0) = -sqrtQ - a / 3; + GSL_IMAG (*z0) = 0; + GSL_REAL (*z1) = -sqrtQ - a / 3; + GSL_IMAG (*z1) = 0; + GSL_REAL (*z2) = 2 * sqrtQ - a / 3; + GSL_IMAG (*z2) = 0; + } + return 3; + } + else if (CR2 < CQ3) /* equivalent to R2 < Q3 */ + { + double sqrtQ = sqrt (Q); + double sqrtQ3 = sqrtQ * sqrtQ * sqrtQ; + double theta = acos (R / sqrtQ3); + double norm = -2 * sqrtQ; + double r0 = norm * cos (theta / 3) - a / 3; + double r1 = norm * cos ((theta + 2.0 * M_PI) / 3) - a / 3; + double r2 = norm * cos ((theta - 2.0 * M_PI) / 3) - a / 3; + + /* Sort r0, r1, r2 into increasing order */ + + if (r0 > r1) + SWAP (r0, r1); + + if (r1 > r2) + { + SWAP (r1, r2); + + if (r0 > r1) + SWAP (r0, r1); + } + + GSL_REAL (*z0) = r0; + GSL_IMAG (*z0) = 0; + + GSL_REAL (*z1) = r1; + GSL_IMAG (*z1) = 0; + + GSL_REAL (*z2) = r2; + GSL_IMAG (*z2) = 0; + + return 3; + } + else + { + double sgnR = (R >= 0 ? 1 : -1); + double A = -sgnR * pow (fabs (R) + sqrt (R2 - Q3), 1.0 / 3.0); + double B = Q / A; + + if (A + B < 0) + { + GSL_REAL (*z0) = A + B - a / 3; + GSL_IMAG (*z0) = 0; + + GSL_REAL (*z1) = -0.5 * (A + B) - a / 3; + GSL_IMAG (*z1) = -(sqrt (3.0) / 2.0) * fabs(A - B); + + GSL_REAL (*z2) = -0.5 * (A + B) - a / 3; + GSL_IMAG (*z2) = (sqrt (3.0) / 2.0) * fabs(A - B); + } + else + { + GSL_REAL (*z0) = -0.5 * (A + B) - a / 3; + GSL_IMAG (*z0) = -(sqrt (3.0) / 2.0) * fabs(A - B); + + GSL_REAL (*z1) = -0.5 * (A + B) - a / 3; + GSL_IMAG (*z1) = (sqrt (3.0) / 2.0) * fabs(A - B); + + GSL_REAL (*z2) = A + B - a / 3; + GSL_IMAG (*z2) = 0; + } + + return 3; + } +} diff --git a/gsl-1.9/poly/zsolve_init.c b/gsl-1.9/poly/zsolve_init.c new file mode 100644 index 0000000..5c181db --- /dev/null +++ b/gsl-1.9/poly/zsolve_init.c @@ -0,0 +1,69 @@ +/* poly/zsolve_init.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 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 <gsl/gsl_math.h> +#include <gsl/gsl_complex.h> +#include <gsl/gsl_poly.h> +#include <gsl/gsl_errno.h> + +gsl_poly_complex_workspace * +gsl_poly_complex_workspace_alloc (size_t n) +{ + size_t nc ; + + gsl_poly_complex_workspace * w ; + + if (n == 0) + { + GSL_ERROR_VAL ("matrix size n must be positive integer", GSL_EDOM, 0); + } + + w = (gsl_poly_complex_workspace *) + malloc (sizeof(gsl_poly_complex_workspace)); + + if (w == 0) + { + GSL_ERROR_VAL ("failed to allocate space for struct", GSL_ENOMEM, 0); + } + + nc = n - 1; + + w->nc = nc; + + w->matrix = (double *) malloc (nc * nc * sizeof(double)); + + if (w->matrix == 0) + { + free (w) ; /* error in constructor, avoid memory leak */ + + GSL_ERROR_VAL ("failed to allocate space for workspace matrix", + GSL_ENOMEM, 0); + } + + return w ; +} + +void +gsl_poly_complex_workspace_free (gsl_poly_complex_workspace * w) +{ + free(w->matrix) ; + free(w); +} diff --git a/gsl-1.9/poly/zsolve_quadratic.c b/gsl-1.9/poly/zsolve_quadratic.c new file mode 100644 index 0000000..dcdeec0 --- /dev/null +++ b/gsl-1.9/poly/zsolve_quadratic.c @@ -0,0 +1,103 @@ +/* poly/zsolve_quadratic.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 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. + */ + +/* complex_solve_quadratic.c - finds complex roots of a x^2 + b x + c = 0 */ + +#include <config.h> +#include <math.h> +#include <gsl/gsl_complex.h> +#include <gsl/gsl_poly.h> + +int +gsl_poly_complex_solve_quadratic (double a, double b, double c, + gsl_complex *z0, gsl_complex *z1) +{ + double disc = b * b - 4 * a * c; + + if (a == 0) /* Handle linear case */ + { + if (b == 0) + { + return 0; + } + else + { + GSL_REAL(*z0) = -c / b; + GSL_IMAG(*z0) = 0; + return 1; + }; + } + + if (disc > 0) + { + if (b == 0) + { + double s = fabs (0.5 * sqrt (disc) / a); + GSL_REAL (*z0) = -s; + GSL_IMAG (*z0) = 0; + GSL_REAL (*z1) = s; + GSL_IMAG (*z1) = 0; + } + else + { + double sgnb = (b > 0 ? 1 : -1); + double temp = -0.5 * (b + sgnb * sqrt (disc)); + double r1 = temp / a; + double r2 = c / temp; + + if (r1 < r2) + { + GSL_REAL (*z0) = r1; + GSL_IMAG (*z0) = 0; + GSL_REAL (*z1) = r2; + GSL_IMAG (*z1) = 0; + } + else + { + GSL_REAL (*z0) = r2; + GSL_IMAG (*z0) = 0; + GSL_REAL (*z1) = r1; + GSL_IMAG (*z1) = 0; + } + } + return 2; + } + else if (disc == 0) + { + GSL_REAL (*z0) = -0.5 * b / a; + GSL_IMAG (*z0) = 0; + + GSL_REAL (*z1) = -0.5 * b / a; + GSL_IMAG (*z1) = 0; + + return 2; + } + else + { + double s = fabs (0.5 * sqrt (-disc) / a); + + GSL_REAL (*z0) = -0.5 * b / a; + GSL_IMAG (*z0) = -s; + + GSL_REAL (*z1) = -0.5 * b / a; + GSL_IMAG (*z1) = s; + + return 2; + } +} |