diff options
Diffstat (limited to 'gsl-1.9/qrng')
| -rw-r--r-- | gsl-1.9/qrng/ChangeLog | 23 | ||||
| -rw-r--r-- | gsl-1.9/qrng/Makefile.am | 14 | ||||
| -rw-r--r-- | gsl-1.9/qrng/Makefile.in | 543 | ||||
| -rw-r--r-- | gsl-1.9/qrng/TODO | 10 | ||||
| -rw-r--r-- | gsl-1.9/qrng/gsl_qrng.h | 110 | ||||
| -rw-r--r-- | gsl-1.9/qrng/niederreiter-2.c | 353 | ||||
| -rw-r--r-- | gsl-1.9/qrng/qrng.c | 127 | ||||
| -rw-r--r-- | gsl-1.9/qrng/sobol.c | 231 | ||||
| -rw-r--r-- | gsl-1.9/qrng/test.c | 122 |
9 files changed, 1533 insertions, 0 deletions
diff --git a/gsl-1.9/qrng/ChangeLog b/gsl-1.9/qrng/ChangeLog new file mode 100644 index 0000000..e46c87f --- /dev/null +++ b/gsl-1.9/qrng/ChangeLog @@ -0,0 +1,23 @@ +2006-04-21 Brian Gough <bjg@network-theory.co.uk> + + * qrng.c gsl_qrng.h: use q instead of r in all prototypes + +2004-05-26 Brian Gough <bjg@network-theory.co.uk> + + * test.c: added <stdlib.h> for exit() + +2002-11-16 Brian Gough <bjg@network-theory.co.uk> + + * niederreiter-2.c (calculate_cj): clear uninitialized memory + +Mon Apr 22 19:20:05 2002 Brian Gough <bjg@network-theory.co.uk> + + * test.c: output results for each individual test + + * qrng.c (gsl_qrng_init): added missing init function + +Mon Oct 23 19:59:46 2000 Brian Gough <bjg@network-theory.co.uk> + + * qrng.c (gsl_qrng_get): use an array type for the array + arguments, rather than a pointer + diff --git a/gsl-1.9/qrng/Makefile.am b/gsl-1.9/qrng/Makefile.am new file mode 100644 index 0000000..b33787a --- /dev/null +++ b/gsl-1.9/qrng/Makefile.am @@ -0,0 +1,14 @@ +noinst_LTLIBRARIES = libgslqrng.la + +pkginclude_HEADERS = gsl_qrng.h + +INCLUDES= -I$(top_builddir) + +libgslqrng_la_SOURCES = gsl_qrng.h qrng.c niederreiter-2.c sobol.c + +TESTS = $(check_PROGRAMS) +check_PROGRAMS = test + +test_SOURCES = test.c +test_LDADD = libgslqrng.la ../ieee-utils/libgslieeeutils.la ../err/libgslerr.la ../test/libgsltest.la ../sys/libgslsys.la ../utils/libutils.la + diff --git a/gsl-1.9/qrng/Makefile.in b/gsl-1.9/qrng/Makefile.in new file mode 100644 index 0000000..7847ef7 --- /dev/null +++ b/gsl-1.9/qrng/Makefile.in @@ -0,0 +1,543 @@ +# 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 = qrng +DIST_COMMON = $(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) +libgslqrng_la_LIBADD = +am_libgslqrng_la_OBJECTS = qrng.lo niederreiter-2.lo sobol.lo +libgslqrng_la_OBJECTS = $(am_libgslqrng_la_OBJECTS) +am_test_OBJECTS = test.$(OBJEXT) +test_OBJECTS = $(am_test_OBJECTS) +test_DEPENDENCIES = libgslqrng.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 = $(libgslqrng_la_SOURCES) $(test_SOURCES) +DIST_SOURCES = $(libgslqrng_la_SOURCES) $(test_SOURCES) +am__vpath_adj_setup = srcdirstrip=`echo "$(srcdir)" | sed 's|.|.|g'`; +am__vpath_adj = case $$p in \ + $(srcdir)/*) f=`echo "$$p" | sed "s|^$$srcdirstrip/||"`;; \ + *) f=$$p;; \ + esac; +am__strip_dir = `echo $$p | sed -e 's|^.*/||'`; +am__installdirs = "$(DESTDIR)$(pkgincludedir)" +pkgincludeHEADERS_INSTALL = $(INSTALL_HEADER) +HEADERS = $(pkginclude_HEADERS) +ETAGS = etags +CTAGS = ctags +DISTFILES = $(DIST_COMMON) $(DIST_SOURCES) $(TEXINFOS) $(EXTRA_DIST) +ACLOCAL = @ACLOCAL@ +AMTAR = @AMTAR@ +AR = @AR@ +AUTOCONF = @AUTOCONF@ +AUTOHEADER = @AUTOHEADER@ +AUTOMAKE = @AUTOMAKE@ +AWK = @AWK@ +CC = @CC@ +CFLAGS = @CFLAGS@ +CPP = @CPP@ +CPPFLAGS = @CPPFLAGS@ +CYGPATH_W = @CYGPATH_W@ +DEFS = @DEFS@ +ECHO = @ECHO@ +ECHO_C = @ECHO_C@ +ECHO_N = @ECHO_N@ +ECHO_T = @ECHO_T@ +EGREP = @EGREP@ +EXEEXT = @EXEEXT@ +GSL_CFLAGS = @GSL_CFLAGS@ +GSL_LIBS = @GSL_LIBS@ +GSL_LT_CBLAS_VERSION = @GSL_LT_CBLAS_VERSION@ +GSL_LT_VERSION = @GSL_LT_VERSION@ +HAVE_AIX_IEEE_INTERFACE = @HAVE_AIX_IEEE_INTERFACE@ +HAVE_DARWIN86_IEEE_INTERFACE = @HAVE_DARWIN86_IEEE_INTERFACE@ +HAVE_DARWIN_IEEE_INTERFACE = @HAVE_DARWIN_IEEE_INTERFACE@ +HAVE_EXTENDED_PRECISION_REGISTERS = @HAVE_EXTENDED_PRECISION_REGISTERS@ +HAVE_FREEBSD_IEEE_INTERFACE = @HAVE_FREEBSD_IEEE_INTERFACE@ +HAVE_GNUM68K_IEEE_INTERFACE = @HAVE_GNUM68K_IEEE_INTERFACE@ +HAVE_GNUPPC_IEEE_INTERFACE = @HAVE_GNUPPC_IEEE_INTERFACE@ +HAVE_GNUSPARC_IEEE_INTERFACE = @HAVE_GNUSPARC_IEEE_INTERFACE@ +HAVE_GNUX86_IEEE_INTERFACE = @HAVE_GNUX86_IEEE_INTERFACE@ +HAVE_HPUX11_IEEE_INTERFACE = @HAVE_HPUX11_IEEE_INTERFACE@ +HAVE_HPUX_IEEE_INTERFACE = @HAVE_HPUX_IEEE_INTERFACE@ +HAVE_IEEE_COMPARISONS = @HAVE_IEEE_COMPARISONS@ +HAVE_IEEE_DENORMALS = @HAVE_IEEE_DENORMALS@ +HAVE_INLINE = @HAVE_INLINE@ +HAVE_IRIX_IEEE_INTERFACE = @HAVE_IRIX_IEEE_INTERFACE@ +HAVE_NETBSD_IEEE_INTERFACE = @HAVE_NETBSD_IEEE_INTERFACE@ +HAVE_OPENBSD_IEEE_INTERFACE = @HAVE_OPENBSD_IEEE_INTERFACE@ +HAVE_OS2EMX_IEEE_INTERFACE = @HAVE_OS2EMX_IEEE_INTERFACE@ +HAVE_PRINTF_LONGDOUBLE = @HAVE_PRINTF_LONGDOUBLE@ +HAVE_SOLARIS_IEEE_INTERFACE = @HAVE_SOLARIS_IEEE_INTERFACE@ +HAVE_SUNOS4_IEEE_INTERFACE = @HAVE_SUNOS4_IEEE_INTERFACE@ +HAVE_TRU64_IEEE_INTERFACE = @HAVE_TRU64_IEEE_INTERFACE@ +INSTALL_DATA = @INSTALL_DATA@ +INSTALL_PROGRAM = @INSTALL_PROGRAM@ +INSTALL_SCRIPT = @INSTALL_SCRIPT@ +INSTALL_STRIP_PROGRAM = @INSTALL_STRIP_PROGRAM@ +LDFLAGS = @LDFLAGS@ +LIBOBJS = @LIBOBJS@ +LIBS = @LIBS@ +LIBTOOL = @LIBTOOL@ +LN_S = @LN_S@ +LTLIBOBJS = @LTLIBOBJS@ +MAINT = @MAINT@ +MAINTAINER_MODE_FALSE = @MAINTAINER_MODE_FALSE@ +MAINTAINER_MODE_TRUE = @MAINTAINER_MODE_TRUE@ +MAKEINFO = @MAKEINFO@ +OBJEXT = @OBJEXT@ +PACKAGE = @PACKAGE@ +PACKAGE_BUGREPORT = @PACKAGE_BUGREPORT@ +PACKAGE_NAME = @PACKAGE_NAME@ +PACKAGE_STRING = @PACKAGE_STRING@ +PACKAGE_TARNAME = @PACKAGE_TARNAME@ +PACKAGE_VERSION = @PACKAGE_VERSION@ +PATH_SEPARATOR = @PATH_SEPARATOR@ +RANLIB = @RANLIB@ +RELEASED = @RELEASED@ +SET_MAKE = @SET_MAKE@ +SHELL = @SHELL@ +STRIP = @STRIP@ +VERSION = @VERSION@ +ac_ct_AR = @ac_ct_AR@ +ac_ct_CC = @ac_ct_CC@ +ac_ct_RANLIB = @ac_ct_RANLIB@ +ac_ct_STRIP = @ac_ct_STRIP@ +am__leading_dot = @am__leading_dot@ +am__tar = @am__tar@ +am__untar = @am__untar@ +bindir = @bindir@ +build = @build@ +build_alias = @build_alias@ +build_cpu = @build_cpu@ +build_os = @build_os@ +build_vendor = @build_vendor@ +datadir = @datadir@ +exec_prefix = @exec_prefix@ +host = @host@ +host_alias = @host_alias@ +host_cpu = @host_cpu@ +host_os = @host_os@ +host_vendor = @host_vendor@ +includedir = @includedir@ +infodir = @infodir@ +install_sh = @install_sh@ +libdir = @libdir@ +libexecdir = @libexecdir@ +localstatedir = @localstatedir@ +mandir = @mandir@ +mkdir_p = @mkdir_p@ +oldincludedir = @oldincludedir@ +prefix = @prefix@ +program_transform_name = @program_transform_name@ +sbindir = @sbindir@ +sharedstatedir = @sharedstatedir@ +sysconfdir = @sysconfdir@ +target_alias = @target_alias@ +noinst_LTLIBRARIES = libgslqrng.la +pkginclude_HEADERS = gsl_qrng.h +INCLUDES = -I$(top_builddir) +libgslqrng_la_SOURCES = gsl_qrng.h qrng.c niederreiter-2.c sobol.c +TESTS = $(check_PROGRAMS) +test_SOURCES = test.c +test_LDADD = libgslqrng.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 qrng/Makefile'; \ + cd $(top_srcdir) && \ + $(AUTOMAKE) --gnu --ignore-deps qrng/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 +libgslqrng.la: $(libgslqrng_la_OBJECTS) $(libgslqrng_la_DEPENDENCIES) + $(LINK) $(libgslqrng_la_LDFLAGS) $(libgslqrng_la_OBJECTS) $(libgslqrng_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/qrng/TODO b/gsl-1.9/qrng/TODO new file mode 100644 index 0000000..50ca918 --- /dev/null +++ b/gsl-1.9/qrng/TODO @@ -0,0 +1,10 @@ +* Sort out "const" in prototypes, it looks odd since there are consts +for many functions which modify the state. (Applies to both qrng and rng) + +* It would be useful to include functions to generate quasi-random + number distributions (at least the Gaussian distribution)? obviously + the Box-Mueller algorithm cannot be used to maintain the low + discrepancy characteristics of the sequence, there are several + methods suggested in the literature but it isn't at all my field so + I'm not sure which is to be preferred. (DENISON Francis + <francis.denison@irsn.fr>) diff --git a/gsl-1.9/qrng/gsl_qrng.h b/gsl-1.9/qrng/gsl_qrng.h new file mode 100644 index 0000000..889799a --- /dev/null +++ b/gsl-1.9/qrng/gsl_qrng.h @@ -0,0 +1,110 @@ +/* Author: G. Jungman + */ +#ifndef __GSL_QRNG_H__ +#define __GSL_QRNG_H__ + +#include <stdlib.h> +#include <gsl/gsl_types.h> +#include <gsl/gsl_errno.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 + + +/* Once again, more inane C-style OOP... kill me now. */ + +/* Structure describing a type of generator. + */ +typedef struct +{ + const char * name; + unsigned int max_dimension; + size_t (*state_size) (unsigned int dimension); + int (*init_state) (void * state, unsigned int dimension); + int (*get) (void * state, unsigned int dimension, double x[]); +} +gsl_qrng_type; + +/* Structure describing a generator instance of a + * specified type, with generator-specific state info + * and dimension-specific info. + */ +typedef struct +{ + const gsl_qrng_type * type; + unsigned int dimension; + size_t state_size; + void * state; +} +gsl_qrng; + + +/* Supported generator types. + */ +GSL_VAR const gsl_qrng_type * gsl_qrng_niederreiter_2; +GSL_VAR const gsl_qrng_type * gsl_qrng_sobol; + + +/* Allocate and initialize a generator + * of the specified type, in the given + * space dimension. + */ +gsl_qrng * gsl_qrng_alloc (const gsl_qrng_type * T, unsigned int dimension); + + +/* Copy a generator. */ +int gsl_qrng_memcpy (gsl_qrng * dest, const gsl_qrng * src); + + +/* Clone a generator. */ +gsl_qrng * gsl_qrng_clone (const gsl_qrng * q); + + +/* Free a generator. */ +void gsl_qrng_free (gsl_qrng * q); + + +/* Intialize a generator. */ +void gsl_qrng_init (gsl_qrng * q); + + +/* Get the standardized name of the generator. */ +const char * gsl_qrng_name (const gsl_qrng * q); + + +/* ISN'T THIS CONFUSING FOR PEOPLE? + WHAT IF SOMEBODY TRIES TO COPY WITH THIS ??? + */ +size_t gsl_qrng_size (const gsl_qrng * q); + + +void * gsl_qrng_state (const gsl_qrng * q); + + +/* Retrieve next vector in sequence. */ +int gsl_qrng_get (const gsl_qrng * q, double x[]); + + +#ifdef HAVE_INLINE +extern inline int gsl_qrng_get (const gsl_qrng * q, double x[]); +extern inline int gsl_qrng_get (const gsl_qrng * q, double x[]) +{ + return (q->type->get) (q->state, q->dimension, x); +} + +#endif /* HAVE_INLINE */ + + +__END_DECLS + + +#endif /* !__GSL_QRNG_H__ */ diff --git a/gsl-1.9/qrng/niederreiter-2.c b/gsl-1.9/qrng/niederreiter-2.c new file mode 100644 index 0000000..1a28009 --- /dev/null +++ b/gsl-1.9/qrng/niederreiter-2.c @@ -0,0 +1,353 @@ +/* Author: G. Jungman + */ +/* Implement Niederreiter base 2 generator. + * See: + * Bratley, Fox, Niederreiter, ACM Trans. Model. Comp. Sim. 2, 195 (1992) + */ +#include <config.h> +#include <gsl/gsl_qrng.h> + + +#define NIED2_CHARACTERISTIC 2 +#define NIED2_BASE 2 + +#define NIED2_MAX_DIMENSION 12 +#define NIED2_MAX_PRIM_DEGREE 5 +#define NIED2_MAX_DEGREE 50 + +#define NIED2_BIT_COUNT 30 +#define NIED2_NBITS (NIED2_BIT_COUNT+1) + +#define MAXV (NIED2_NBITS + NIED2_MAX_DEGREE) + +/* Z_2 field operations */ +#define NIED2_ADD(x,y) (((x)+(y))%2) +#define NIED2_MUL(x,y) (((x)*(y))%2) +#define NIED2_SUB(x,y) NIED2_ADD((x),(y)) + + +static size_t nied2_state_size(unsigned int dimension); +static int nied2_init(void * state, unsigned int dimension); +static int nied2_get(void * state, unsigned int dimension, double * v); + + +static const gsl_qrng_type nied2_type = +{ + "niederreiter-base-2", + NIED2_MAX_DIMENSION, + nied2_state_size, + nied2_init, + nied2_get +}; + +const gsl_qrng_type * gsl_qrng_niederreiter_2 = &nied2_type; + +/* primitive polynomials in binary encoding */ +static const int primitive_poly[NIED2_MAX_DIMENSION+1][NIED2_MAX_PRIM_DEGREE+1] = +{ + { 1, 0, 0, 0, 0, 0 }, /* 1 */ + { 0, 1, 0, 0, 0, 0 }, /* x */ + { 1, 1, 0, 0, 0, 0 }, /* 1 + x */ + { 1, 1, 1, 0, 0, 0 }, /* 1 + x + x^2 */ + { 1, 1, 0, 1, 0, 0 }, /* 1 + x + x^3 */ + { 1, 0, 1, 1, 0, 0 }, /* 1 + x^2 + x^3 */ + { 1, 1, 0, 0, 1, 0 }, /* 1 + x + x^4 */ + { 1, 0, 0, 1, 1, 0 }, /* 1 + x^3 + x^4 */ + { 1, 1, 1, 1, 1, 0 }, /* 1 + x + x^2 + x^3 + x^4 */ + { 1, 0, 1, 0, 0, 1 }, /* 1 + x^2 + x^5 */ + { 1, 0, 0, 1, 0, 1 }, /* 1 + x^3 + x^5 */ + { 1, 1, 1, 1, 0, 1 }, /* 1 + x + x^2 + x^3 + x^5 */ + { 1, 1, 1, 0, 1, 1 } /* 1 + x + x^2 + x^4 + x^5 */ +}; + +/* degrees of primitive polynomials */ +static const int poly_degree[NIED2_MAX_DIMENSION+1] = +{ + 0, 1, 1, 2, 3, 3, 4, 4, 4, 5, 5, 5, 5 +}; + + +typedef struct +{ + unsigned int sequence_count; + int cj[NIED2_NBITS][NIED2_MAX_DIMENSION]; + int nextq[NIED2_MAX_DIMENSION]; +} nied2_state_t; + + +static size_t nied2_state_size(unsigned int dimension) +{ + return sizeof(nied2_state_t); +} + + +/* Multiply polynomials over Z_2. + * Notice use of a temporary vector, + * side-stepping aliasing issues when + * one of inputs is the same as the output + * [especially important in the original fortran version, I guess]. + */ +static void poly_multiply( + const int pa[], int pa_degree, + const int pb[], int pb_degree, + int pc[], int * pc_degree + ) +{ + int j, k; + int pt[NIED2_MAX_DEGREE+1]; + int pt_degree = pa_degree + pb_degree; + + for(k=0; k<=pt_degree; k++) { + int term = 0; + for(j=0; j<=k; j++) { + const int conv_term = NIED2_MUL(pa[k-j], pb[j]); + term = NIED2_ADD(term, conv_term); + } + pt[k] = term; + } + + for(k=0; k<=pt_degree; k++) { + pc[k] = pt[k]; + } + for(k=pt_degree+1; k<=NIED2_MAX_DEGREE; k++) { + pc[k] = 0; + } + + *pc_degree = pt_degree; +} + + +/* Calculate the values of the constants V(J,R) as + * described in BFN section 3.3. + * + * px = appropriate irreducible polynomial for current dimension + * pb = polynomial defined in section 2.3 of BFN. + * pb is modified + */ +static void calculate_v( + const int px[], int px_degree, + int pb[], int * pb_degree, + int v[], int maxv + ) +{ + const int nonzero_element = 1; /* nonzero element of Z_2 */ + const int arbitrary_element = 1; /* arbitray element of Z_2 */ + + /* The polynomial ph is px**(J-1), which is the value of B on arrival. + * In section 3.3, the values of Hi are defined with a minus sign: + * don't forget this if you use them later ! + */ + int ph[NIED2_MAX_DEGREE+1]; + /* int ph_degree = *pb_degree; */ + int bigm = *pb_degree; /* m from section 3.3 */ + int m; /* m from section 2.3 */ + int r, k, kj; + + for(k=0; k<=NIED2_MAX_DEGREE; k++) { + ph[k] = pb[k]; + } + + /* Now multiply B by PX so B becomes PX**J. + * In section 2.3, the values of Bi are defined with a minus sign : + * don't forget this if you use them later ! + */ + poly_multiply(px, px_degree, pb, *pb_degree, pb, pb_degree); + m = *pb_degree; + + /* Now choose a value of Kj as defined in section 3.3. + * We must have 0 <= Kj < E*J = M. + * The limit condition on Kj does not seem very relevant + * in this program. + */ + /* Quoting from BFN: "Our program currently sets each K_q + * equal to eq. This has the effect of setting all unrestricted + * values of v to 1." + * Actually, it sets them to the arbitrary chosen value. + * Whatever. + */ + kj = bigm; + + /* Now choose values of V in accordance with + * the conditions in section 3.3. + */ + for(r=0; r<kj; r++) { + v[r] = 0; + } + v[kj] = 1; + + + if(kj >= bigm) { + for(r=kj+1; r<m; r++) { + v[r] = arbitrary_element; + } + } + else { + /* This block is never reached. */ + + int term = NIED2_SUB(0, ph[kj]); + + for(r=kj+1; r<bigm; r++) { + v[r] = arbitrary_element; + + /* Check the condition of section 3.3, + * remembering that the H's have the opposite sign. [????????] + */ + term = NIED2_SUB(term, NIED2_MUL(ph[r], v[r])); + } + + /* Now v[bigm] != term. */ + v[bigm] = NIED2_ADD(nonzero_element, term); + + for(r=bigm+1; r<m; r++) { + v[r] = arbitrary_element; + } + } + + /* Calculate the remaining V's using the recursion of section 2.3, + * remembering that the B's have the opposite sign. + */ + for(r=0; r<=maxv-m; r++) { + int term = 0; + for(k=0; k<m; k++) { + term = NIED2_SUB(term, NIED2_MUL(pb[k], v[r+k])); + } + v[r+m] = term; + } +} + + +static void calculate_cj(nied2_state_t * ns, unsigned int dimension) +{ + int ci[NIED2_NBITS][NIED2_NBITS]; + int v[MAXV+1]; + int r; + unsigned int i_dim; + + for(i_dim=0; i_dim<dimension; i_dim++) { + + const int poly_index = i_dim + 1; + int j, k; + + /* Niederreiter (page 56, after equation (7), defines two + * variables Q and U. We do not need Q explicitly, but we + * do need U. + */ + int u = 0; + + /* For each dimension, we need to calculate powers of an + * appropriate irreducible polynomial, see Niederreiter + * page 65, just below equation (19). + * Copy the appropriate irreducible polynomial into PX, + * and its degree into E. Set polynomial B = PX ** 0 = 1. + * M is the degree of B. Subsequently B will hold higher + * powers of PX. + */ + int pb[NIED2_MAX_DEGREE+1]; + int px[NIED2_MAX_DEGREE+1]; + int px_degree = poly_degree[poly_index]; + int pb_degree = 0; + + for(k=0; k<=px_degree; k++) { + px[k] = primitive_poly[poly_index][k]; + pb[k] = 0; + } + + for (;k<NIED2_MAX_DEGREE+1;k++) { + px[k] = 0; + pb[k] = 0; + } + + pb[0] = 1; + + for(j=0; j<NIED2_NBITS; j++) { + + /* If U = 0, we need to set B to the next power of PX + * and recalculate V. + */ + if(u == 0) calculate_v(px, px_degree, pb, &pb_degree, v, MAXV); + + /* Now C is obtained from V. Niederreiter + * obtains A from V (page 65, near the bottom), and then gets + * C from A (page 56, equation (7)). However this can be done + * in one step. Here CI(J,R) corresponds to + * Niederreiter's C(I,J,R). + */ + for(r=0; r<NIED2_NBITS; r++) { + ci[r][j] = v[r+u]; + } + + /* Advance Niederreiter's state variables. */ + ++u; + if(u == px_degree) u = 0; + } + + /* The array CI now holds the values of C(I,J,R) for this value + * of I. We pack them into array CJ so that CJ(I,R) holds all + * the values of C(I,J,R) for J from 1 to NBITS. + */ + for(r=0; r<NIED2_NBITS; r++) { + int term = 0; + for(j=0; j<NIED2_NBITS; j++) { + term = 2*term + ci[r][j]; + } + ns->cj[r][i_dim] = term; + } + + } +} + + +static int nied2_init(void * state, unsigned int dimension) +{ + nied2_state_t * n_state = (nied2_state_t *) state; + unsigned int i_dim; + + if(dimension < 1 || dimension > NIED2_MAX_DIMENSION) return GSL_EINVAL; + + calculate_cj(n_state, dimension); + + for(i_dim=0; i_dim<dimension; i_dim++) n_state->nextq[i_dim] = 0; + n_state->sequence_count = 0; + + return GSL_SUCCESS; +} + + +static int nied2_get(void * state, unsigned int dimension, double * v) +{ + static const double recip = 1.0/(double)(1U << NIED2_NBITS); /* 2^(-nbits) */ + nied2_state_t * n_state = (nied2_state_t *) state; + int r; + int c; + unsigned int i_dim; + + /* Load the result from the saved state. */ + for(i_dim=0; i_dim<dimension; i_dim++) { + v[i_dim] = n_state->nextq[i_dim] * recip; + } + + /* Find the position of the least-significant zero in sequence_count. + * This is the bit that changes in the Gray-code representation as + * the count is advanced. + */ + r = 0; + c = n_state->sequence_count; + while(1) { + if((c % 2) == 1) { + ++r; + c /= 2; + } + else break; + } + + if(r >= NIED2_NBITS) return GSL_EFAILED; /* FIXME: better error code here */ + + /* Calculate the next state. */ + for(i_dim=0; i_dim<dimension; i_dim++) { + n_state->nextq[i_dim] ^= n_state->cj[r][i_dim]; + } + + n_state->sequence_count++; + + return GSL_SUCCESS; +} diff --git a/gsl-1.9/qrng/qrng.c b/gsl-1.9/qrng/qrng.c new file mode 100644 index 0000000..d85eccf --- /dev/null +++ b/gsl-1.9/qrng/qrng.c @@ -0,0 +1,127 @@ +/* Author: G. Jungman + */ +#include <config.h> +#include <stdlib.h> +#include <string.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_qrng.h> + + +gsl_qrng * +gsl_qrng_alloc (const gsl_qrng_type * T, unsigned int dimension) +{ + + gsl_qrng * q = (gsl_qrng *) malloc (sizeof (gsl_qrng)); + + if (q == 0) + { + GSL_ERROR_VAL ("allocation failed for qrng struct", + GSL_ENOMEM, 0); + }; + + q->dimension = dimension; + q->state_size = T->state_size(dimension); + q->state = malloc (q->state_size); + + if (q->state == 0) + { + free (q); + GSL_ERROR_VAL ("allocation failed for qrng state", + GSL_ENOMEM, 0); + }; + + q->type = T; + + T->init_state(q->state, q->dimension); + + return q; +} + +void +gsl_qrng_init (gsl_qrng * q) +{ + (q->type->init_state) (q->state, q->dimension); +} + +int +gsl_qrng_memcpy (gsl_qrng * dest, const gsl_qrng * src) +{ + if (dest->type != src->type) + { + GSL_ERROR ("generators must be of the same type", GSL_EINVAL); + } + + dest->dimension = src->dimension; + dest->state_size = src->state_size; + memcpy (dest->state, src->state, src->state_size); + + return GSL_SUCCESS; +} + + +gsl_qrng * +gsl_qrng_clone (const gsl_qrng * q) +{ + gsl_qrng * r = (gsl_qrng *) malloc (sizeof (gsl_qrng)); + + if (r == 0) + { + GSL_ERROR_VAL ("failed to allocate space for rng struct", + GSL_ENOMEM, 0); + }; + + r->dimension = q->dimension; + r->state_size = q->state_size; + r->state = malloc (r->state_size); + + if (r->state == 0) + { + free (r); + GSL_ERROR_VAL ("failed to allocate space for rng state", + GSL_ENOMEM, 0); + }; + + r->type = q->type; + + memcpy (r->state, q->state, q->state_size); + + return r; +} + +#ifndef HIDE_INLINE_STATIC +int +gsl_qrng_get (const gsl_qrng * q, double x[]) +{ + return (q->type->get) (q->state, q->dimension, x); +} +#endif + +const char * +gsl_qrng_name (const gsl_qrng * q) +{ + return q->type->name; +} + + +size_t +gsl_qrng_size (const gsl_qrng * q) +{ + return q->state_size; +} + + +void * +gsl_qrng_state (const gsl_qrng * q) +{ + return q->state; +} + + +void +gsl_qrng_free (gsl_qrng * q) +{ + if(q != 0) { + if(q->state != 0) free (q->state); + free (q); + } +} diff --git a/gsl-1.9/qrng/sobol.c b/gsl-1.9/qrng/sobol.c new file mode 100644 index 0000000..51c9baa --- /dev/null +++ b/gsl-1.9/qrng/sobol.c @@ -0,0 +1,231 @@ +/* Author: G. Jungman + */ +/* Implementation for Sobol generator. + * See + * [Bratley+Fox, TOMS 14, 88 (1988)] + * [Antonov+Saleev, USSR Comput. Maths. Math. Phys. 19, 252 (1980)] + */ +#include <config.h> +#include <gsl/gsl_qrng.h> + + +/* maximum allowed space dimension */ +#define SOBOL_MAX_DIMENSION 40 + +/* bit count; assumes sizeof(int) >= 32-bit */ +#define SOBOL_BIT_COUNT 30 + +/* prototypes for generator type functions */ +static size_t sobol_state_size(unsigned int dimension); +static int sobol_init(void * state, unsigned int dimension); +static int sobol_get(void * state, unsigned int dimension, double * v); + +/* global Sobol generator type object */ +static const gsl_qrng_type sobol_type = +{ + "sobol", + SOBOL_MAX_DIMENSION, + sobol_state_size, + sobol_init, + sobol_get +}; +const gsl_qrng_type * gsl_qrng_sobol = &sobol_type; + + +/* primitive polynomials in binary encoding + */ +static const int primitive_polynomials[SOBOL_MAX_DIMENSION] = +{ + 1, 3, 7, 11, 13, 19, 25, 37, 59, 47, + 61, 55, 41, 67, 97, 91, 109, 103, 115, 131, + 193, 137, 145, 143, 241, 157, 185, 167, 229, 171, + 213, 191, 253, 203, 211, 239, 247, 285, 369, 299 +}; + +/* degrees of the primitive polynomials */ +static const int degree_table[SOBOL_MAX_DIMENSION] = +{ + 0, 1, 2, 3, 3, 4, 4, 5, 5, 5, + 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, + 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, + 7, 7, 7, 7, 7, 7, 7, 8, 8, 8 +}; + + +/* initial values for direction tables, following + * Bratley+Fox, taken from [Sobol+Levitan, preprint 1976] + */ +static const int v_init[8][SOBOL_MAX_DIMENSION] = +{ + { + 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 + }, + { + 0, 0, 1, 3, 1, 3, 1, 3, 3, 1, + 3, 1, 3, 1, 3, 1, 1, 3, 1, 3, + 1, 3, 1, 3, 3, 1, 3, 1, 3, 1, + 3, 1, 1, 3, 1, 3, 1, 3, 1, 3 + }, + { + 0, 0, 0, 7, 5, 1, 3, 3, 7, 5, + 5, 7, 7, 1, 3, 3, 7, 5, 1, 1, + 5, 3, 3, 1, 7, 5, 1, 3, 3, 7, + 5, 1, 1, 5, 7, 7, 5, 1, 3, 3 + }, + { + 0, 0, 0, 0, 0, 1, 7, 9, 13, 11, + 1, 3, 7, 9, 5, 13, 13, 11, 3, 15, + 5, 3, 15, 7, 9, 13, 9, 1, 11, 7, + 5, 15, 1, 15, 11, 5, 3, 1, 7, 9 + }, + { + 0, 0, 0, 0, 0, 0, 0, 9, 3, 27, + 15, 29, 21, 23, 19, 11, 25, 7, 13, 17, + 1, 25, 29, 3, 31, 11, 5, 23, 27, 19, + 21, 5, 1, 17, 13, 7, 15, 9, 31, 9 + }, + { + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 37, 33, 7, 5, 11, 39, 63, + 27, 17, 15, 23, 29, 3, 21, 13, 31, 25, + 9, 49, 33, 19, 29, 11, 19, 27, 15, 25 + }, + { + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 13, + 33, 115, 41, 79, 17, 29, 119, 75, 73, 105, + 7, 59, 65, 21, 3, 113, 61, 89, 45, 107 + }, + { + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 7, 23, 39 + } +}; + + +/* Sobol generator state. + * sequence_count = number of calls with this generator + * last_numerator_vec = last generated numerator vector + * last_denominator_inv = 1/denominator for last numerator vector + * v_direction = direction number table + */ +typedef struct +{ + unsigned int sequence_count; + double last_denominator_inv; + int last_numerator_vec[SOBOL_MAX_DIMENSION]; + int v_direction[SOBOL_BIT_COUNT][SOBOL_MAX_DIMENSION]; +} sobol_state_t; + +/* NOTE: I fixed the size for the preliminary implementation. + This could/should be relaxed, as an optimization. + */ + +static size_t sobol_state_size(unsigned int dimension) +{ + return sizeof(sobol_state_t); +} + + +static int sobol_init(void * state, unsigned int dimension) +{ + sobol_state_t * s_state = (sobol_state_t *) state; + unsigned int i_dim; + int j, k; + int ell; + + if(dimension < 1 || dimension > SOBOL_MAX_DIMENSION) { + return GSL_EINVAL; + } + + /* Initialize direction table in dimension 0. */ + for(k=0; k<SOBOL_BIT_COUNT; k++) s_state->v_direction[k][0] = 1; + + /* Initialize in remaining dimensions. */ + for(i_dim=1; i_dim<dimension; i_dim++) { + + const int poly_index = i_dim; + const int degree_i = degree_table[poly_index]; + int includ[8]; + + /* Expand the polynomial bit pattern to separate + * components of the logical array includ[]. + */ + int p_i = primitive_polynomials[poly_index]; + for(k = degree_i-1; k >= 0; k--) { + includ[k] = ((p_i % 2) == 1); + p_i /= 2; + } + + /* Leading elements for dimension i come from v_init[][]. */ + for(j=0; j<degree_i; j++) s_state->v_direction[j][i_dim] = v_init[j][i_dim]; + + /* Calculate remaining elements for this dimension, + * as explained in Bratley+Fox, section 2. + */ + for(j=degree_i; j<SOBOL_BIT_COUNT; j++) { + int newv = s_state->v_direction[j-degree_i][i_dim]; + ell = 1; + for(k=0; k<degree_i; k++) { + ell *= 2; + if(includ[k]) newv ^= (ell * s_state->v_direction[j-k-1][i_dim]); + } + s_state->v_direction[j][i_dim] = newv; + } + } + + /* Multiply columns of v by appropriate power of 2. */ + ell = 1; + for(j=SOBOL_BIT_COUNT-1-1; j>=0; j--) { + ell *= 2; + for(i_dim=0; i_dim<dimension; i_dim++) { + s_state->v_direction[j][i_dim] *= ell; + } + } + + /* 1/(common denominator of the elements in v_direction) */ + s_state->last_denominator_inv = 1.0 /(2.0 * ell); + + /* final setup */ + s_state->sequence_count = 0; + for(i_dim=0; i_dim<dimension; i_dim++) s_state->last_numerator_vec[i_dim] = 0; + + return GSL_SUCCESS; +} + + +static int sobol_get(void * state, unsigned int dimension, double * v) +{ + sobol_state_t * s_state = (sobol_state_t *) state; + + unsigned int i_dimension; + + /* Find the position of the least-significant zero in count. */ + int ell = 0; + int c = s_state->sequence_count; + while(1) { + ++ell; + if((c % 2) == 1) c /= 2; + else break; + } + + /* Check for exhaustion. */ + if(ell > SOBOL_BIT_COUNT) return GSL_EFAILED; /* FIXME: good return code here */ + + for(i_dimension=0; i_dimension<dimension; i_dimension++) { + const int direction_i = s_state->v_direction[ell-1][i_dimension]; + const int old_numerator_i = s_state->last_numerator_vec[i_dimension]; + const int new_numerator_i = old_numerator_i ^ direction_i; + s_state->last_numerator_vec[i_dimension] = new_numerator_i; + v[i_dimension] = new_numerator_i * s_state->last_denominator_inv; + } + + s_state->sequence_count++; + + return GSL_SUCCESS; +} diff --git a/gsl-1.9/qrng/test.c b/gsl-1.9/qrng/test.c new file mode 100644 index 0000000..4b7d08f --- /dev/null +++ b/gsl-1.9/qrng/test.c @@ -0,0 +1,122 @@ +/* Author: G. Jungman + */ +#include <config.h> +#include <stdlib.h> +#include <gsl/gsl_ieee_utils.h> + +#include <gsl/gsl_qrng.h> +#include <gsl/gsl_test.h> + +void test_sobol(void) +{ + int status = 0; + double v[3]; + /* int i; */ + + /* test in dimension 2 */ + gsl_qrng * g = gsl_qrng_alloc(gsl_qrng_sobol, 2); + gsl_qrng_get(g, v); + gsl_qrng_get(g, v); + gsl_qrng_get(g, v); + status += ( v[0] != 0.25 || v[1] != 0.75 ); + gsl_qrng_get(g, v); + status += ( v[0] != 0.375 || v[1] != 0.375 ); + gsl_qrng_free(g); + + gsl_test (status, "Sobol d=2"); + + status = 0; + /* test in dimension 3 */ + g = gsl_qrng_alloc(gsl_qrng_sobol, 3); + gsl_qrng_get(g, v); + gsl_qrng_get(g, v); + gsl_qrng_get(g, v); + status += ( v[0] != 0.25 || v[1] != 0.75 || v[2] != 0.25 ); + gsl_qrng_get(g, v); + status += ( v[0] != 0.375 || v[1] != 0.375 || v[2] != 0.625 ); + + gsl_test (status, "Sobol d=3"); + + status = 0; + gsl_qrng_init(g); + gsl_qrng_get(g, v); + gsl_qrng_get(g, v); + gsl_qrng_get(g, v); + status += ( v[0] != 0.25 || v[1] != 0.75 || v[2] != 0.25 ); + gsl_qrng_get(g, v); + status += ( v[0] != 0.375 || v[1] != 0.375 || v[2] != 0.625 ); + gsl_qrng_free(g); + + gsl_test (status, "Sobol d=3 (reinitialized)"); +} + + +void test_nied2(void) +{ + int status = 0; + double v[3]; + /* int i; */ + + /* test in dimension 2 */ + gsl_qrng * g = gsl_qrng_alloc(gsl_qrng_niederreiter_2, 2); + gsl_qrng_get(g, v); + gsl_qrng_get(g, v); + gsl_qrng_get(g, v); + status += ( v[0] != 0.75 || v[1] != 0.25 ); + gsl_qrng_get(g, v); + status += ( v[0] != 0.25 || v[1] != 0.75 ); + gsl_qrng_get(g, v); + gsl_qrng_get(g, v); + gsl_qrng_get(g, v); + status += ( v[0] != 0.625 || v[1] != 0.125 ); + gsl_qrng_free(g); + + gsl_test (status, "Niederreiter d=2"); + + status = 0; + + /* test in dimension 3 */ + g = gsl_qrng_alloc(gsl_qrng_niederreiter_2, 3); + gsl_qrng_get(g, v); + gsl_qrng_get(g, v); + gsl_qrng_get(g, v); + status += ( v[0] != 0.75 || v[1] != 0.25 || v[2] != 0.3125 ); + gsl_qrng_get(g, v); + status += ( v[0] != 0.25 || v[1] != 0.75 || v[2] != 0.5625 ); + gsl_qrng_get(g, v); + gsl_qrng_get(g, v); + gsl_qrng_get(g, v); + status += ( v[0] != 0.625 || v[1] != 0.125 || v[2] != 0.6875 ); + + gsl_test (status, "Niederreiter d=3"); + + status = 0; + + gsl_qrng_init(g); + gsl_qrng_get(g, v); + gsl_qrng_get(g, v); + gsl_qrng_get(g, v); + status += ( v[0] != 0.75 || v[1] != 0.25 || v[2] != 0.3125 ); + gsl_qrng_get(g, v); + status += ( v[0] != 0.25 || v[1] != 0.75 || v[2] != 0.5625 ); + gsl_qrng_get(g, v); + gsl_qrng_get(g, v); + gsl_qrng_get(g, v); + status += ( v[0] != 0.625 || v[1] != 0.125 || v[2] != 0.6875 ); + gsl_qrng_free(g); + + + gsl_test (status, "Niederreiter d=3 (reinitialized)"); +} + + +int main() +{ + + gsl_ieee_env_setup (); + + test_sobol(); + test_nied2(); + + exit (gsl_test_summary ()); +} |