diff options
Diffstat (limited to 'nist')
| -rw-r--r-- | nist/Makefile.am | 22 | ||||
| -rw-r--r-- | nist/Makefile.in | 576 | ||||
| -rw-r--r-- | nist/cephes-protos.h | 180 | ||||
| -rw-r--r-- | nist/cephes.c | 1327 | ||||
| -rw-r--r-- | nist/dfft.c | 1381 | ||||
| -rw-r--r-- | nist/matrix.c | 168 | ||||
| -rw-r--r-- | nist/matrix.h | 19 | ||||
| -rw-r--r-- | nist/mconf.h | 169 | ||||
| -rw-r--r-- | nist/nist.c | 57 | ||||
| -rw-r--r-- | nist/nist.h | 33 | ||||
| -rw-r--r-- | nist/packtest.c | 2252 | ||||
| -rw-r--r-- | nist/special-functions.c | 32 | ||||
| -rw-r--r-- | nist/special-functions.h | 14 | ||||
| -rw-r--r-- | nist/template9 | 148 |
14 files changed, 6378 insertions, 0 deletions
diff --git a/nist/Makefile.am b/nist/Makefile.am new file mode 100644 index 0000000..bbaaf17 --- /dev/null +++ b/nist/Makefile.am @@ -0,0 +1,22 @@ +## Process this file with automake to produce Makefile.in. + + +if ENABLE_NIST_TEST +check_PROGRAMS = nist +endif + +AM_CFLAGS=-Wall + +nist_LDADD = -lm + +nist_SOURCES = cephes.c dfft.c matrix.c nist.c packtest.c special-functions.c cephes-protos.h matrix.h mconf.h special-functions.h nist.h + +CLEANFILES = sample nist.out + +MAINTAINERCLEANFILES = Makefile.in + +if ENABLE_NIST_TEST +check-local: + ../src/haveged -n 16M $* + ./nist sample ${srcdir} +endif
\ No newline at end of file diff --git a/nist/Makefile.in b/nist/Makefile.in new file mode 100644 index 0000000..82999a1 --- /dev/null +++ b/nist/Makefile.in @@ -0,0 +1,576 @@ +# Makefile.in generated by automake 1.16.1 from Makefile.am. +# @configure_input@ + +# Copyright (C) 1994-2018 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@ +VPATH = @srcdir@ +am__is_gnu_make = { \ + if test -z '$(MAKELEVEL)'; then \ + false; \ + elif test -n '$(MAKE_HOST)'; then \ + true; \ + elif test -n '$(MAKE_VERSION)' && test -n '$(CURDIR)'; then \ + true; \ + else \ + false; \ + fi; \ +} +am__make_running_with_option = \ + case $${target_option-} in \ + ?) ;; \ + *) echo "am__make_running_with_option: internal error: invalid" \ + "target option '$${target_option-}' specified" >&2; \ + exit 1;; \ + esac; \ + has_opt=no; \ + sane_makeflags=$$MAKEFLAGS; \ + if $(am__is_gnu_make); then \ + sane_makeflags=$$MFLAGS; \ + else \ + case $$MAKEFLAGS in \ + *\\[\ \ ]*) \ + bs=\\; \ + sane_makeflags=`printf '%s\n' "$$MAKEFLAGS" \ + | sed "s/$$bs$$bs[$$bs $$bs ]*//g"`;; \ + esac; \ + fi; \ + skip_next=no; \ + strip_trailopt () \ + { \ + flg=`printf '%s\n' "$$flg" | sed "s/$$1.*$$//"`; \ + }; \ + for flg in $$sane_makeflags; do \ + test $$skip_next = yes && { skip_next=no; continue; }; \ + case $$flg in \ + *=*|--*) continue;; \ + -*I) strip_trailopt 'I'; skip_next=yes;; \ + -*I?*) strip_trailopt 'I';; \ + -*O) strip_trailopt 'O'; skip_next=yes;; \ + -*O?*) strip_trailopt 'O';; \ + -*l) strip_trailopt 'l'; skip_next=yes;; \ + -*l?*) strip_trailopt 'l';; \ + -[dEDm]) skip_next=yes;; \ + -[JT]) skip_next=yes;; \ + esac; \ + case $$flg in \ + *$$target_option*) has_opt=yes; break;; \ + esac; \ + done; \ + test $$has_opt = yes +am__make_dryrun = (target_option=n; $(am__make_running_with_option)) +am__make_keepgoing = (target_option=k; $(am__make_running_with_option)) +pkgdatadir = $(datadir)/@PACKAGE@ +pkgincludedir = $(includedir)/@PACKAGE@ +pkglibdir = $(libdir)/@PACKAGE@ +pkglibexecdir = $(libexecdir)/@PACKAGE@ +am__cd = CDPATH="$${ZSH_VERSION+.}$(PATH_SEPARATOR)" && cd +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@ +@ENABLE_NIST_TEST_TRUE@check_PROGRAMS = nist$(EXEEXT) +subdir = nist +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) +DIST_COMMON = $(srcdir)/Makefile.am $(am__DIST_COMMON) +mkinstalldirs = $(install_sh) -d +CONFIG_HEADER = $(top_builddir)/config.h +CONFIG_CLEAN_FILES = +CONFIG_CLEAN_VPATH_FILES = +am_nist_OBJECTS = cephes.$(OBJEXT) dfft.$(OBJEXT) matrix.$(OBJEXT) \ + nist.$(OBJEXT) packtest.$(OBJEXT) special-functions.$(OBJEXT) +nist_OBJECTS = $(am_nist_OBJECTS) +nist_DEPENDENCIES = +AM_V_lt = $(am__v_lt_@AM_V@) +am__v_lt_ = $(am__v_lt_@AM_DEFAULT_V@) +am__v_lt_0 = --silent +am__v_lt_1 = +AM_V_P = $(am__v_P_@AM_V@) +am__v_P_ = $(am__v_P_@AM_DEFAULT_V@) +am__v_P_0 = false +am__v_P_1 = : +AM_V_GEN = $(am__v_GEN_@AM_V@) +am__v_GEN_ = $(am__v_GEN_@AM_DEFAULT_V@) +am__v_GEN_0 = @echo " GEN " $@; +am__v_GEN_1 = +AM_V_at = $(am__v_at_@AM_V@) +am__v_at_ = $(am__v_at_@AM_DEFAULT_V@) +am__v_at_0 = @ +am__v_at_1 = +DEFAULT_INCLUDES = -I.@am__isrc@ -I$(top_builddir) +depcomp = +am__maybe_remake_depfiles = +COMPILE = $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) \ + $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) +LTCOMPILE = $(LIBTOOL) $(AM_V_lt) --tag=CC $(AM_LIBTOOLFLAGS) \ + $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) \ + $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) \ + $(AM_CFLAGS) $(CFLAGS) +AM_V_CC = $(am__v_CC_@AM_V@) +am__v_CC_ = $(am__v_CC_@AM_DEFAULT_V@) +am__v_CC_0 = @echo " CC " $@; +am__v_CC_1 = +CCLD = $(CC) +LINK = $(LIBTOOL) $(AM_V_lt) --tag=CC $(AM_LIBTOOLFLAGS) \ + $(LIBTOOLFLAGS) --mode=link $(CCLD) $(AM_CFLAGS) $(CFLAGS) \ + $(AM_LDFLAGS) $(LDFLAGS) -o $@ +AM_V_CCLD = $(am__v_CCLD_@AM_V@) +am__v_CCLD_ = $(am__v_CCLD_@AM_DEFAULT_V@) +am__v_CCLD_0 = @echo " CCLD " $@; +am__v_CCLD_1 = +SOURCES = $(nist_SOURCES) +DIST_SOURCES = $(nist_SOURCES) +am__can_run_installinfo = \ + case $$AM_UPDATE_INFO_DIR in \ + n|no|NO) false;; \ + *) (install-info --version) >/dev/null 2>&1;; \ + esac +am__tagged_files = $(HEADERS) $(SOURCES) $(TAGS_FILES) $(LISP) +# Read a list of newline-separated strings from the standard input, +# and print each of them once, without duplicates. Input order is +# *not* preserved. +am__uniquify_input = $(AWK) '\ + BEGIN { nonempty = 0; } \ + { items[$$0] = 1; nonempty = 1; } \ + END { if (nonempty) { for (i in items) print i; }; } \ +' +# Make sure the list of sources is unique. This is necessary because, +# e.g., the same source file might be shared among _SOURCES variables +# for different programs/libraries. +am__define_uniq_tagged_files = \ + list='$(am__tagged_files)'; \ + unique=`for i in $$list; do \ + if test -f "$$i"; then echo $$i; else echo $(srcdir)/$$i; fi; \ + done | $(am__uniquify_input)` +ETAGS = etags +CTAGS = ctags +am__DIST_COMMON = $(srcdir)/Makefile.in +DISTFILES = $(DIST_COMMON) $(DIST_SOURCES) $(TEXINFOS) $(EXTRA_DIST) +ACLOCAL = @ACLOCAL@ +AMTAR = @AMTAR@ +AM_DEFAULT_VERBOSITY = @AM_DEFAULT_VERBOSITY@ +AR = @AR@ +AUTOCONF = @AUTOCONF@ +AUTOHEADER = @AUTOHEADER@ +AUTOMAKE = @AUTOMAKE@ +AWK = @AWK@ +CC = @CC@ +CFLAGS = @CFLAGS@ +CPP = @CPP@ +CPPFLAGS = @CPPFLAGS@ +CYGPATH_W = @CYGPATH_W@ +DEFS = @DEFS@ +DLLTOOL = @DLLTOOL@ +DSYMUTIL = @DSYMUTIL@ +DUMPBIN = @DUMPBIN@ +ECHO_C = @ECHO_C@ +ECHO_N = @ECHO_N@ +ECHO_T = @ECHO_T@ +EGREP = @EGREP@ +EXEEXT = @EXEEXT@ +FGREP = @FGREP@ +GREP = @GREP@ +HAVEGE_LT_VERSION = @HAVEGE_LT_VERSION@ +HA_DISTRO = @HA_DISTRO@ +HA_LDFLAGS = @HA_LDFLAGS@ +HA_UNITD = @HA_UNITD@ +INSTALL = @INSTALL@ +INSTALL_DATA = @INSTALL_DATA@ +INSTALL_PROGRAM = @INSTALL_PROGRAM@ +INSTALL_SCRIPT = @INSTALL_SCRIPT@ +INSTALL_STRIP_PROGRAM = @INSTALL_STRIP_PROGRAM@ +LD = @LD@ +LDFLAGS = @LDFLAGS@ +LIBOBJS = @LIBOBJS@ +LIBS = @LIBS@ +LIBTOOL = @LIBTOOL@ +LIPO = @LIPO@ +LN_S = @LN_S@ +LTLIBOBJS = @LTLIBOBJS@ +LT_SYS_LIBRARY_PATH = @LT_SYS_LIBRARY_PATH@ +MAKEINFO = @MAKEINFO@ +MANIFEST_TOOL = @MANIFEST_TOOL@ +MKDIR_P = @MKDIR_P@ +NM = @NM@ +NMEDIT = @NMEDIT@ +OBJDUMP = @OBJDUMP@ +OBJEXT = @OBJEXT@ +OTOOL = @OTOOL@ +OTOOL64 = @OTOOL64@ +PACKAGE = @PACKAGE@ +PACKAGE_BUGREPORT = @PACKAGE_BUGREPORT@ +PACKAGE_NAME = @PACKAGE_NAME@ +PACKAGE_STRING = @PACKAGE_STRING@ +PACKAGE_TARNAME = @PACKAGE_TARNAME@ +PACKAGE_URL = @PACKAGE_URL@ +PACKAGE_VERSION = @PACKAGE_VERSION@ +PATH_SEPARATOR = @PATH_SEPARATOR@ +RANLIB = @RANLIB@ +SED = @SED@ +SET_MAKE = @SET_MAKE@ +SHELL = @SHELL@ +STRIP = @STRIP@ +VERSION = @VERSION@ +abs_builddir = @abs_builddir@ +abs_srcdir = @abs_srcdir@ +abs_top_builddir = @abs_top_builddir@ +abs_top_srcdir = @abs_top_srcdir@ +ac_ct_AR = @ac_ct_AR@ +ac_ct_CC = @ac_ct_CC@ +ac_ct_DUMPBIN = @ac_ct_DUMPBIN@ +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@ +builddir = @builddir@ +datadir = @datadir@ +datarootdir = @datarootdir@ +docdir = @docdir@ +dvidir = @dvidir@ +exec_prefix = @exec_prefix@ +host = @host@ +host_alias = @host_alias@ +host_cpu = @host_cpu@ +host_os = @host_os@ +host_vendor = @host_vendor@ +htmldir = @htmldir@ +includedir = @includedir@ +infodir = @infodir@ +install_sh = @install_sh@ +libdir = @libdir@ +libexecdir = @libexecdir@ +localedir = @localedir@ +localstatedir = @localstatedir@ +mandir = @mandir@ +mkdir_p = @mkdir_p@ +oldincludedir = @oldincludedir@ +pdfdir = @pdfdir@ +prefix = @prefix@ +program_transform_name = @program_transform_name@ +psdir = @psdir@ +sbindir = @sbindir@ +sharedstatedir = @sharedstatedir@ +srcdir = @srcdir@ +sysconfdir = @sysconfdir@ +target_alias = @target_alias@ +top_build_prefix = @top_build_prefix@ +top_builddir = @top_builddir@ +top_srcdir = @top_srcdir@ +AM_CFLAGS = -Wall +nist_LDADD = -lm +nist_SOURCES = cephes.c dfft.c matrix.c nist.c packtest.c special-functions.c cephes-protos.h matrix.h mconf.h special-functions.h nist.h +CLEANFILES = sample nist.out +MAINTAINERCLEANFILES = Makefile.in +all: all-am + +.SUFFIXES: +.SUFFIXES: .c .lo .o .obj +$(srcdir)/Makefile.in: $(srcdir)/Makefile.am $(am__configure_deps) + @for dep in $?; do \ + case '$(am__configure_deps)' in \ + *$$dep*) \ + ( cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh ) \ + && { if test -f $@; then exit 0; else break; fi; }; \ + exit 1;; \ + esac; \ + done; \ + echo ' cd $(top_srcdir) && $(AUTOMAKE) --gnu --ignore-deps nist/Makefile'; \ + $(am__cd) $(top_srcdir) && \ + $(AUTOMAKE) --gnu --ignore-deps nist/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__maybe_remake_depfiles)'; \ + cd $(top_builddir) && $(SHELL) ./config.status $(subdir)/$@ $(am__maybe_remake_depfiles);; \ + esac; + +$(top_builddir)/config.status: $(top_srcdir)/configure $(CONFIG_STATUS_DEPENDENCIES) + cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh + +$(top_srcdir)/configure: $(am__configure_deps) + cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh +$(ACLOCAL_M4): $(am__aclocal_m4_deps) + cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh +$(am__aclocal_m4_deps): + +clean-checkPROGRAMS: + @list='$(check_PROGRAMS)'; test -n "$$list" || exit 0; \ + echo " rm -f" $$list; \ + rm -f $$list || exit $$?; \ + test -n "$(EXEEXT)" || exit 0; \ + list=`for p in $$list; do echo "$$p"; done | sed 's/$(EXEEXT)$$//'`; \ + echo " rm -f" $$list; \ + rm -f $$list + +nist$(EXEEXT): $(nist_OBJECTS) $(nist_DEPENDENCIES) $(EXTRA_nist_DEPENDENCIES) + @rm -f nist$(EXEEXT) + $(AM_V_CCLD)$(LINK) $(nist_OBJECTS) $(nist_LDADD) $(LIBS) + +mostlyclean-compile: + -rm -f *.$(OBJEXT) + +distclean-compile: + -rm -f *.tab.c + +.c.o: + $(AM_V_CC)$(COMPILE) -c -o $@ $< + +.c.obj: + $(AM_V_CC)$(COMPILE) -c -o $@ `$(CYGPATH_W) '$<'` + +.c.lo: + $(AM_V_CC)$(LTCOMPILE) -c -o $@ $< + +mostlyclean-libtool: + -rm -f *.lo + +clean-libtool: + -rm -rf .libs _libs + +ID: $(am__tagged_files) + $(am__define_uniq_tagged_files); mkid -fID $$unique +tags: tags-am +TAGS: tags + +tags-am: $(TAGS_DEPENDENCIES) $(am__tagged_files) + set x; \ + here=`pwd`; \ + $(am__define_uniq_tagged_files); \ + shift; \ + if test -z "$(ETAGS_ARGS)$$*$$unique"; then :; else \ + test -n "$$unique" || unique=$$empty_fix; \ + if test $$# -gt 0; then \ + $(ETAGS) $(ETAGSFLAGS) $(AM_ETAGSFLAGS) $(ETAGS_ARGS) \ + "$$@" $$unique; \ + else \ + $(ETAGS) $(ETAGSFLAGS) $(AM_ETAGSFLAGS) $(ETAGS_ARGS) \ + $$unique; \ + fi; \ + fi +ctags: ctags-am + +CTAGS: ctags +ctags-am: $(TAGS_DEPENDENCIES) $(am__tagged_files) + $(am__define_uniq_tagged_files); \ + test -z "$(CTAGS_ARGS)$$unique" \ + || $(CTAGS) $(CTAGSFLAGS) $(AM_CTAGSFLAGS) $(CTAGS_ARGS) \ + $$unique + +GTAGS: + here=`$(am__cd) $(top_builddir) && pwd` \ + && $(am__cd) $(top_srcdir) \ + && gtags -i $(GTAGS_ARGS) "$$here" +cscopelist: cscopelist-am + +cscopelist-am: $(am__tagged_files) + list='$(am__tagged_files)'; \ + case "$(srcdir)" in \ + [\\/]* | ?:[\\/]*) sdir="$(srcdir)" ;; \ + *) sdir=$(subdir)/$(srcdir) ;; \ + esac; \ + for i in $$list; do \ + if test -f "$$i"; then \ + echo "$(subdir)/$$i"; \ + else \ + echo "$$sdir/$$i"; \ + fi; \ + done >> $(top_builddir)/cscope.files + +distclean-tags: + -rm -f TAGS ID GTAGS GRTAGS GSYMS GPATH tags + +distdir: $(BUILT_SOURCES) + $(MAKE) $(AM_MAKEFLAGS) distdir-am + +distdir-am: $(DISTFILES) + @srcdirstrip=`echo "$(srcdir)" | sed 's/[].[^$$\\*]/\\\\&/g'`; \ + topsrcdirstrip=`echo "$(top_srcdir)" | sed 's/[].[^$$\\*]/\\\\&/g'`; \ + list='$(DISTFILES)'; \ + dist_files=`for file in $$list; do echo $$file; done | \ + sed -e "s|^$$srcdirstrip/||;t" \ + -e "s|^$$topsrcdirstrip/|$(top_builddir)/|;t"`; \ + case $$dist_files in \ + */*) $(MKDIR_P) `echo "$$dist_files" | \ + sed '/\//!d;s|^|$(distdir)/|;s,/[^/]*$$,,' | \ + sort -u` ;; \ + esac; \ + for file in $$dist_files; do \ + if test -f $$file || test -d $$file; then d=.; else d=$(srcdir); fi; \ + if test -d $$d/$$file; then \ + dir=`echo "/$$file" | sed -e 's,/[^/]*$$,,'`; \ + if test -d "$(distdir)/$$file"; then \ + find "$(distdir)/$$file" -type d ! -perm -700 -exec chmod u+rwx {} \;; \ + fi; \ + if test -d $(srcdir)/$$file && test $$d != $(srcdir); then \ + cp -fpR $(srcdir)/$$file "$(distdir)$$dir" || exit 1; \ + find "$(distdir)/$$file" -type d ! -perm -700 -exec chmod u+rwx {} \;; \ + fi; \ + cp -fpR $$d/$$file "$(distdir)$$dir" || exit 1; \ + else \ + test -f "$(distdir)/$$file" \ + || cp -p $$d/$$file "$(distdir)/$$file" \ + || exit 1; \ + fi; \ + done +@ENABLE_NIST_TEST_FALSE@check-local: +check-am: all-am + $(MAKE) $(AM_MAKEFLAGS) $(check_PROGRAMS) + $(MAKE) $(AM_MAKEFLAGS) check-local +check: check-am +all-am: Makefile +installdirs: +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: + if test -z '$(STRIP)'; then \ + $(MAKE) $(AM_MAKEFLAGS) INSTALL_PROGRAM="$(INSTALL_STRIP_PROGRAM)" \ + install_sh_PROGRAM="$(INSTALL_STRIP_PROGRAM)" INSTALL_STRIP_FLAG=-s \ + install; \ + else \ + $(MAKE) $(AM_MAKEFLAGS) INSTALL_PROGRAM="$(INSTALL_STRIP_PROGRAM)" \ + install_sh_PROGRAM="$(INSTALL_STRIP_PROGRAM)" INSTALL_STRIP_FLAG=-s \ + "INSTALL_PROGRAM_ENV=STRIPPROG='$(STRIP)'" install; \ + fi +mostlyclean-generic: + +clean-generic: + -test -z "$(CLEANFILES)" || rm -f $(CLEANFILES) + +distclean-generic: + -test -z "$(CONFIG_CLEAN_FILES)" || rm -f $(CONFIG_CLEAN_FILES) + -test . = "$(srcdir)" || test -z "$(CONFIG_CLEAN_VPATH_FILES)" || rm -f $(CONFIG_CLEAN_VPATH_FILES) + +maintainer-clean-generic: + @echo "This command is intended for maintainers to use" + @echo "it deletes files that may require special tools to rebuild." + -test -z "$(MAINTAINERCLEANFILES)" || rm -f $(MAINTAINERCLEANFILES) +clean: clean-am + +clean-am: clean-checkPROGRAMS clean-generic clean-libtool \ + mostlyclean-am + +distclean: distclean-am + -rm -f Makefile +distclean-am: clean-am distclean-compile distclean-generic \ + distclean-tags + +dvi: dvi-am + +dvi-am: + +html: html-am + +html-am: + +info: info-am + +info-am: + +install-data-am: + +install-dvi: install-dvi-am + +install-dvi-am: + +install-exec-am: + +install-html: install-html-am + +install-html-am: + +install-info: install-info-am + +install-info-am: + +install-man: + +install-pdf: install-pdf-am + +install-pdf-am: + +install-ps: install-ps-am + +install-ps-am: + +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: + +.MAKE: check-am install-am install-strip + +.PHONY: CTAGS GTAGS TAGS all all-am check check-am check-local clean \ + clean-checkPROGRAMS clean-generic clean-libtool cscopelist-am \ + ctags ctags-am 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-dvi install-dvi-am install-exec \ + install-exec-am install-html install-html-am install-info \ + install-info-am install-man install-pdf install-pdf-am \ + install-ps install-ps-am 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 tags-am uninstall uninstall-am + +.PRECIOUS: Makefile + + +@ENABLE_NIST_TEST_TRUE@check-local: +@ENABLE_NIST_TEST_TRUE@ ../src/haveged -n 16M $* +@ENABLE_NIST_TEST_TRUE@ ./nist sample ${srcdir} + +# 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/nist/cephes-protos.h b/nist/cephes-protos.h new file mode 100644 index 0000000..1041936 --- /dev/null +++ b/nist/cephes-protos.h @@ -0,0 +1,180 @@ +/* + * This file was automatically generated by version 1.7 of cextract. + * Manual editing not recommended. + * + * Created: Fri Mar 31 19:17:33 1995 + */ +extern double acosh ( double x ); +extern int airy ( double x, double *ai, double *aip, double *bi, double *bip ); +extern double asin ( double x ); +extern double acos ( double x ); +extern double asinh ( double x ); +extern double atan ( double x ); +extern double atan2 ( double y, double x ); +extern double atanh ( double x ); +extern double bdtrc ( int k, int n, double p ); +extern double bdtr ( int k, int n, double p ); +extern double bdtri ( int k, int n, double y ); +extern double beta ( double a, double b ); +extern double lbeta ( double a, double b ); +extern double btdtr ( double a, double b, double x ); +extern double cbrt ( double x ); +extern double chbevl ( double x, void *P, int n ); +extern double chdtrc ( double df, double x ); +extern double chdtr ( double df, double x ); +extern double chdtri ( double df, double y ); +// extern void clog ( cmplx *z, cmplx *w ); +extern void _cexp ( cmplx *z, cmplx *w ); +extern void _csin ( cmplx *z, cmplx *w ); +extern void _ccos ( cmplx *z, cmplx *w ); +extern void _ctan ( cmplx *z, cmplx *w ); +extern void _ccot ( cmplx *z, cmplx *w ); +extern void _casin ( cmplx *z, cmplx *w ); +extern void _cacos ( cmplx *z, cmplx *w ); +extern void _catan ( cmplx *z, cmplx *w ); +extern void cadd ( cmplx *a, cmplx *b, cmplx *c ); +extern void csub ( cmplx *a, cmplx *b, cmplx *c ); +extern void cmul ( cmplx *a, cmplx *b, cmplx *c ); +extern void cdiv ( cmplx *a, cmplx *b, cmplx *c ); +extern void cmov ( void *a, void *b ); +extern void cneg ( cmplx *a ); +//extern double _cabs ( cmplx *z ); +extern void _csqrt ( cmplx *z, cmplx *w ); +extern double hypot ( double x, double y ); +extern double cosh ( double x ); +extern double dawsn ( double xx ); +extern void eigens ( double A[], double RR[], double E[], int N ); +extern double ellie ( double phi, double m ); +extern double ellik ( double phi, double m ); +extern double ellpe ( double x ); +extern int ellpj ( double u, double m, double *sn, double *cn, double *dn, double *ph ); +extern double ellpk ( double x ); +extern double exp ( double x ); +extern double exp10 ( double x ); +extern double exp1m ( double x ); +extern double exp2 ( double x ); +extern double expn ( int n, double x ); +extern double fabs ( double x ); +extern double fac ( int i ); +extern double fdtrc ( int ia, int ib, double x ); +extern double fdtr ( int ia, int ib, double x ); +extern double fdtri ( int ia, int ib, double y ); +extern int fftr ( double x[], int m0, double sine[] ); +extern double ceil ( double x ); +extern double floor ( double x ); +extern double frexp ( double x, int *pw2 ); +extern double ldexp ( double x, int pw2 ); +extern int _signbit ( double x ); +/* extern int isnan ( double x ); */ +/*extern int isfinite ( double x );*/ +extern int fresnl ( double xxa, double *ssa, double *cca ); +extern double gamma ( double x ); +extern double lgam ( double x ); +extern double gdtr ( double a, double b, double x ); +extern double gdtrc ( double a, double b, double x ); +extern int gels ( double A[], double R[], int M, double EPS, double AUX[] ); +extern double hyp2f1 ( double a, double b, double c, double x ); +extern double hyperg ( double a, double b, double x ); +extern double hyp2f0 ( double a, double b, double x, int type, double *err ); +extern double i0 ( double x ); +extern double i0e ( double x ); +extern double i1 ( double x ); +extern double i1e ( double x ); +extern double igamc ( double a, double x ); +extern double igam ( double a, double x ); +extern double igami ( double a, double y0 ); +extern double incbet ( double aa, double bb, double xx ); +extern double incbi ( double aa, double bb, double yy0 ); +extern double iv ( double v, double x ); +extern double j0 ( double x ); +extern double y0 ( double x ); +extern double j1 ( double x ); +extern double y1 ( double x ); +extern double jn ( int n, double x ); +extern double jv ( double n, double x ); +extern double k0 ( double x ); +extern double k0e ( double x ); +extern double k1 ( double x ); +extern double k1e ( double x ); +extern double kn ( int nn, double x ); +extern int levnsn ( int n, double r[], double a[], double e[], double refl[] ); +extern double log ( double x ); +extern double log10 ( double x ); +extern double _log2 ( double x ); +extern long lrand ( void ); +extern long lsqrt ( long x ); +extern int minv ( double A[], double X[], int n, double B[], int IPS[] ); +extern int mmmpy ( int r, int c, double *A, double *B, double *Y ); +extern int mtherr ( char *name, int code ); +/*extern double polevl ( double x, void *P, int N ); +extern double p1evl ( double x, void *P, int N );*/ +extern int mtransp ( int n, double *A, double *T ); +extern int mvmpy ( int r, int c, double *A, double *V, double *Y ); +extern double nbdtrc ( int k, int n, double p ); +extern double nbdtr ( int k, int n, double p ); +extern double nbdtri ( int k, int n, double p ); +extern double ndtr ( double a ); +extern double erfc ( double a ); +extern double erf ( double x ); +extern double ndtri ( double y0 ); +extern double pdtrc ( int k, double m ); +extern double pdtr ( int k, double m ); +extern double pdtri ( int k, double y ); +extern double pow ( double x, double y ); +extern double powi ( double x, int nn ); +extern double psi ( double x ); +extern void revers ( double y[], double x[], int n ); +extern double rgamma ( double x ); +extern double round ( double x ); +extern int sprec ( void ); +extern int dprec ( void ); +extern int ldprec ( void ); +extern int shichi ( double x, double *si, double *ci ); +extern int sici ( double x, double *si, double *ci ); +extern double simpsn ( double f[], double delta ); +extern int simq ( double A[], double B[], double X[], int n, int flag, int IPS[] ); +extern double sin ( double x ); +extern double cos ( double x ); +extern double radian ( double d, double m, double s ); +/* +extern int sincos ( double x, double *s, double *c, int flg ); +*/ +extern double sindg ( double x ); +extern double cosdg ( double x ); +extern double sinh ( double x ); +extern double spence ( double x ); +extern double sqrt ( double x ); +extern double stdtr ( int k, double t ); +extern double stdtri ( int k, double p ); +extern double onef2 ( double a, double b, double c, double x, double *err ); +extern double threef0 ( double a, double b, double c, double x, double *err ); +extern double struve ( double v, double x ); +extern double tan ( double x ); +extern double cot ( double x ); +extern double tandg ( double x ); +extern double cotdg ( double x ); +extern double tanh ( double x ); +extern double log1p ( double x ); +extern double expm1 ( double x ); +extern double cosm1 ( double x ); +extern double yn ( int n, double x ); +extern double zeta ( double x, double q ); +extern double zetac ( double x ); +extern int drand ( double *a ); + +/* polyn.c */ +extern void polini ( int maxdeg ); +extern void polprt ( double a[], int na, int d ); +extern void polclr ( double *a, int n ); +extern void polmov ( double *a, int na, double *b ); +extern void polmul ( double a[], int na, double b[], int nb, double c[] ); +extern void poladd ( double a[], int na, double b[], int nb, double c[] ); +extern void polsub ( double a[], int na, double b[], int nb, double c[] ); +extern int poldiv ( double a[], int na, double b[], int nb, double c[] ); +extern void polsbt ( double a[], int na, double b[], int nb, double c[] ); +extern double poleva ( double a[], int na, double x ); +/* polmisc.c */ +extern void polatn ( double num[], double den[], double ans[], int nn ); +extern void polsqt ( double pol[], double ans[], int nn ); +extern void polsin ( double x[], double y[], int nn ); +extern void polcos ( double x[], double y[], int nn ); diff --git a/nist/cephes.c b/nist/cephes.c new file mode 100644 index 0000000..6fc1500 --- /dev/null +++ b/nist/cephes.c @@ -0,0 +1,1327 @@ +#include <stdio.h> +#include "mconf.h" + +double p1evl(double, double [], int); +double polevl(double, double [], int); + +#ifndef ANSIPROT +double lgam(), exp(), log(), fabs(), igam(), igamc(); +#endif + +int merror = 0; + +/* Notice: the order of appearance of the following + * messages is bound to the error codes defined + * in mconf.h. + */ +static char *ermsg[7] = { +"unknown", /* error code 0 */ +"domain", /* error code 1 */ +"singularity", /* et seq. */ +"overflow", +"underflow", +"total loss of precision", +"partial loss of precision" +}; + + + +/* const.c + * + * Globally declared constants + * + * + * + * SYNOPSIS: + * + * extern double nameofconstant; + * + * + * + * + * DESCRIPTION: + * + * This file contains a number of mathematical constants and + * also some needed size parameters of the computer arithmetic. + * The values are supplied as arrays of hexadecimal integers + * for IEEE arithmetic; arrays of octal constants for DEC + * arithmetic; and in a normal decimal scientific notation for + * other machines. The particular notation used is determined + * by a symbol (DEC, IBMPC, or UNK) defined in the include file + * mconf.h. + * + * The default size parameters are as follows. + * + * For DEC and UNK modes: + * MACHEP = 1.38777878078144567553E-17 2**-56 + * MAXLOG = 8.8029691931113054295988E1 log(2**127) + * MINLOG = -8.872283911167299960540E1 log(2**-128) + * MAXNUM = 1.701411834604692317316873e38 2**127 + * + * For IEEE arithmetic (IBMPC): + * MACHEP = 1.11022302462515654042E-16 2**-53 + * MAXLOG = 7.09782712893383996843E2 log(2**1024) + * MINLOG = -7.08396418532264106224E2 log(2**-1022) + * MAXNUM = 1.7976931348623158E308 2**1024 + * + * The global symbols for mathematical constants are + * PI = 3.14159265358979323846 pi + * PIO2 = 1.57079632679489661923 pi/2 + * PIO4 = 7.85398163397448309616E-1 pi/4 + * SQRT2 = 1.41421356237309504880 sqrt(2) + * SQRTH = 7.07106781186547524401E-1 sqrt(2)/2 + * LOG2E = 1.4426950408889634073599 1/log(2) + * SQ2OPI = 7.9788456080286535587989E-1 sqrt( 2/pi ) + * LOGE2 = 6.93147180559945309417E-1 log(2) + * LOGSQ2 = 3.46573590279972654709E-1 log(2)/2 + * THPIO4 = 2.35619449019234492885 3*pi/4 + * TWOOPI = 6.36619772367581343075535E-1 2/pi + * + * These lists are subject to change. + */ + +/* const.c */ + +/* +Cephes Math Library Release 2.3: March, 1995 +Copyright 1984, 1995 by Stephen L. Moshier +*/ + +#ifdef UNK +#if 1 +double MACHEP = 1.11022302462515654042E-16; /* 2**-53 */ +#else +double MACHEP = 1.38777878078144567553E-17; /* 2**-56 */ +#endif +double UFLOWTHRESH = 2.22507385850720138309E-308; /* 2**-1022 */ +#ifdef DENORMAL +double MAXLOG = 7.09782712893383996732E2; /* log(MAXNUM) */ +/* double MINLOG = -7.44440071921381262314E2; */ /* log(2**-1074) */ +double MINLOG = -7.451332191019412076235E2; /* log(2**-1075) */ +#else +double MAXLOG = 7.08396418532264106224E2; /* log 2**1022 */ +double MINLOG = -7.08396418532264106224E2; /* log 2**-1022 */ +#endif +double MAXNUM = 1.79769313486231570815E308; /* 2**1024*(1-MACHEP) */ +double PI = 3.14159265358979323846; /* pi */ +double PIO2 = 1.57079632679489661923; /* pi/2 */ +double PIO4 = 7.85398163397448309616E-1; /* pi/4 */ +double SQRT2 = 1.41421356237309504880; /* sqrt(2) */ +double SQRTH = 7.07106781186547524401E-1; /* sqrt(2)/2 */ +double LOG2E = 1.4426950408889634073599; /* 1/log(2) */ +double SQ2OPI = 7.9788456080286535587989E-1; /* sqrt( 2/pi ) */ +double LOGE2 = 6.93147180559945309417E-1; /* log(2) */ +double LOGSQ2 = 3.46573590279972654709E-1; /* log(2)/2 */ +double THPIO4 = 2.35619449019234492885; /* 3*pi/4 */ +double TWOOPI = 6.36619772367581343075535E-1; /* 2/pi */ +#ifdef INFINITIES +double INFINITY = 1.0/0.0; /* 99e999; */ +#else +double INFINITY = 1.79769313486231570815E308; /* 2**1024*(1-MACHEP) */ +#endif +#ifdef NANS +double NAN = 1.0/0.0 - 1.0/0.0; +#else +double NAN = 0.0; +#endif +#ifdef MINUSZERO +double NEGZERO = -0.0; +#else +double NEGZERO = 0.0; +#endif +#endif + +#ifdef IBMPC + /* 2**-53 = 1.11022302462515654042E-16 */ +unsigned short MACHEP[4] = {0x0000,0x0000,0x0000,0x3ca0}; +unsigned short UFLOWTHRESH[4] = {0x0000,0x0000,0x0000,0x0010}; +#ifdef DENORMAL + /* log(MAXNUM) = 7.09782712893383996732224E2 */ +unsigned short MAXLOG[4] = {0x39ef,0xfefa,0x2e42,0x4086}; + /* log(2**-1074) = - -7.44440071921381262314E2 */ +/*unsigned short MINLOG[4] = {0x71c3,0x446d,0x4385,0xc087};*/ +unsigned short MINLOG[4] = {0x3052,0xd52d,0x4910,0xc087}; +#else + /* log(2**1022) = 7.08396418532264106224E2 */ +unsigned short MAXLOG[4] = {0xbcd2,0xdd7a,0x232b,0x4086}; + /* log(2**-1022) = - 7.08396418532264106224E2 */ +unsigned short MINLOG[4] = {0xbcd2,0xdd7a,0x232b,0xc086}; +#endif + /* 2**1024*(1-MACHEP) = 1.7976931348623158E308 */ +unsigned short MAXNUM[4] = {0xffff,0xffff,0xffff,0x7fef}; +unsigned short PI[4] = {0x2d18,0x5444,0x21fb,0x4009}; +unsigned short PIO2[4] = {0x2d18,0x5444,0x21fb,0x3ff9}; +unsigned short PIO4[4] = {0x2d18,0x5444,0x21fb,0x3fe9}; +unsigned short SQRT2[4] = {0x3bcd,0x667f,0xa09e,0x3ff6}; +unsigned short SQRTH[4] = {0x3bcd,0x667f,0xa09e,0x3fe6}; +unsigned short LOG2E[4] = {0x82fe,0x652b,0x1547,0x3ff7}; +unsigned short SQ2OPI[4] = {0x3651,0x33d4,0x8845,0x3fe9}; +unsigned short LOGE2[4] = {0x39ef,0xfefa,0x2e42,0x3fe6}; +unsigned short LOGSQ2[4] = {0x39ef,0xfefa,0x2e42,0x3fd6}; +unsigned short THPIO4[4] = {0x21d2,0x7f33,0xd97c,0x4002}; +unsigned short TWOOPI[4] = {0xc883,0x6dc9,0x5f30,0x3fe4}; +#ifdef INFINITIES +unsigned short INFINITY[4] = {0x0000,0x0000,0x0000,0x7ff0}; +#else +unsigned short INFINITY[4] = {0xffff,0xffff,0xffff,0x7fef}; +#endif +#ifdef NANS +unsigned short NAN[4] = {0x0000,0x0000,0x0000,0x7ffc}; +#else +unsigned short NAN[4] = {0x0000,0x0000,0x0000,0x0000}; +#endif +#ifdef MINUSZERO +unsigned short NEGZERO[4] = {0x0000,0x0000,0x0000,0x8000}; +#else +unsigned short NEGZERO[4] = {0x0000,0x0000,0x0000,0x0000}; +#endif +#endif + +#ifdef MIEEE + /* 2**-53 = 1.11022302462515654042E-16 */ +unsigned short MACHEP[4] = {0x3ca0,0x0000,0x0000,0x0000}; +unsigned short UFLOWTHRESH[4] = {0x0010,0x0000,0x0000,0x0000}; +#ifdef DENORMAL + /* log(2**1024) = 7.09782712893383996843E2 */ +unsigned short MAXLOG[4] = {0x4086,0x2e42,0xfefa,0x39ef}; + /* log(2**-1074) = - -7.44440071921381262314E2 */ +/* unsigned short MINLOG[4] = {0xc087,0x4385,0x446d,0x71c3}; */ +unsigned short MINLOG[4] = {0xc087,0x4910,0xd52d,0x3052}; +#else + /* log(2**1022) = 7.08396418532264106224E2 */ +unsigned short MAXLOG[4] = {0x4086,0x232b,0xdd7a,0xbcd2}; + /* log(2**-1022) = - 7.08396418532264106224E2 */ +unsigned short MINLOG[4] = {0xc086,0x232b,0xdd7a,0xbcd2}; +#endif + /* 2**1024*(1-MACHEP) = 1.7976931348623158E308 */ +unsigned short MAXNUM[4] = {0x7fef,0xffff,0xffff,0xffff}; +unsigned short PI[4] = {0x4009,0x21fb,0x5444,0x2d18}; +unsigned short PIO2[4] = {0x3ff9,0x21fb,0x5444,0x2d18}; +unsigned short PIO4[4] = {0x3fe9,0x21fb,0x5444,0x2d18}; +unsigned short SQRT2[4] = {0x3ff6,0xa09e,0x667f,0x3bcd}; +unsigned short SQRTH[4] = {0x3fe6,0xa09e,0x667f,0x3bcd}; +unsigned short LOG2E[4] = {0x3ff7,0x1547,0x652b,0x82fe}; +unsigned short SQ2OPI[4] = {0x3fe9,0x8845,0x33d4,0x3651}; +unsigned short LOGE2[4] = {0x3fe6,0x2e42,0xfefa,0x39ef}; +unsigned short LOGSQ2[4] = {0x3fd6,0x2e42,0xfefa,0x39ef}; +unsigned short THPIO4[4] = {0x4002,0xd97c,0x7f33,0x21d2}; +unsigned short TWOOPI[4] = {0x3fe4,0x5f30,0x6dc9,0xc883}; +#ifdef INFINITIES +unsigned short INFINITY[4] = {0x7ff0,0x0000,0x0000,0x0000}; +#else +unsigned short INFINITY[4] = {0x7fef,0xffff,0xffff,0xffff}; +#endif +#ifdef NANS +unsigned short NAN[4] = {0x7ff8,0x0000,0x0000,0x0000}; +#else +unsigned short NAN[4] = {0x0000,0x0000,0x0000,0x0000}; +#endif +#ifdef MINUSZERO +unsigned short NEGZERO[4] = {0x8000,0x0000,0x0000,0x0000}; +#else +unsigned short NEGZERO[4] = {0x0000,0x0000,0x0000,0x0000}; +#endif +#endif + +#ifdef DEC + /* 2**-56 = 1.38777878078144567553E-17 */ +unsigned short MACHEP[4] = {0022200,0000000,0000000,0000000}; +unsigned short UFLOWTHRESH[4] = {0x0080,0x0000,0x0000,0x0000}; + /* log 2**127 = 88.029691931113054295988 */ +unsigned short MAXLOG[4] = {041660,007463,0143742,025733,}; + /* log 2**-128 = -88.72283911167299960540 */ +unsigned short MINLOG[4] = {0141661,071027,0173721,0147572,}; + /* 2**127 = 1.701411834604692317316873e38 */ +unsigned short MAXNUM[4] = {077777,0177777,0177777,0177777,}; +unsigned short PI[4] = {040511,007732,0121041,064302,}; +unsigned short PIO2[4] = {040311,007732,0121041,064302,}; +unsigned short PIO4[4] = {040111,007732,0121041,064302,}; +unsigned short SQRT2[4] = {040265,002363,031771,0157145,}; +unsigned short SQRTH[4] = {040065,002363,031771,0157144,}; +unsigned short LOG2E[4] = {040270,0125073,024534,013761,}; +unsigned short SQ2OPI[4] = {040114,041051,0117241,0131204,}; +unsigned short LOGE2[4] = {040061,071027,0173721,0147572,}; +unsigned short LOGSQ2[4] = {037661,071027,0173721,0147572,}; +unsigned short THPIO4[4] = {040426,0145743,0174631,007222,}; +unsigned short TWOOPI[4] = {040042,0174603,067116,042025,}; +/* Approximate infinity by MAXNUM. */ +unsigned short INFINITY[4] = {077777,0177777,0177777,0177777,}; +unsigned short NAN[4] = {0000000,0000000,0000000,0000000}; +#ifdef MINUSZERO +unsigned short NEGZERO[4] = {0000000,0000000,0000000,0100000}; +#else +unsigned short NEGZERO[4] = {0000000,0000000,0000000,0000000}; +#endif +#endif + +#ifndef UNK +extern unsigned short MACHEP[]; +extern unsigned short UFLOWTHRESH[]; +extern unsigned short MAXLOG[]; +extern unsigned short UNDLOG[]; +extern unsigned short MINLOG[]; +extern unsigned short MAXNUM[]; +extern unsigned short PI[]; +extern unsigned short PIO2[]; +extern unsigned short PIO4[]; +extern unsigned short SQRT2[]; +extern unsigned short SQRTH[]; +extern unsigned short LOG2E[]; +extern unsigned short SQ2OPI[]; +extern unsigned short LOGE2[]; +extern unsigned short LOGSQ2[]; +extern unsigned short THPIO4[]; +extern unsigned short TWOOPI[]; +extern unsigned short INFINITY[]; +extern unsigned short NAN[]; +extern unsigned short NEGZERO[]; +#endif + +extern double MACHEP, MAXLOG; +static double big = 4.503599627370496e15; +static double biginv = 2.22044604925031308085e-16; + +#ifdef UNK +static double P[] = { + 1.60119522476751861407E-4, + 1.19135147006586384913E-3, + 1.04213797561761569935E-2, + 4.76367800457137231464E-2, + 2.07448227648435975150E-1, + 4.94214826801497100753E-1, + 9.99999999999999996796E-1 +}; +static double Q[] = { +-2.31581873324120129819E-5, + 5.39605580493303397842E-4, +-4.45641913851797240494E-3, + 1.18139785222060435552E-2, + 3.58236398605498653373E-2, +-2.34591795718243348568E-1, + 7.14304917030273074085E-2, + 1.00000000000000000320E0 +}; +#define MAXGAM 171.624376956302725 +static double LOGPI = 1.14472988584940017414; +#endif + +#ifdef DEC +static unsigned short P[] = { +0035047,0162701,0146301,0005234, +0035634,0023437,0032065,0176530, +0036452,0137157,0047330,0122574, +0037103,0017310,0143041,0017232, +0037524,0066516,0162563,0164605, +0037775,0004671,0146237,0014222, +0040200,0000000,0000000,0000000 +}; +static unsigned short Q[] = { +0134302,0041724,0020006,0116565, +0035415,0072121,0044251,0025634, +0136222,0003447,0035205,0121114, +0036501,0107552,0154335,0104271, +0037022,0135717,0014776,0171471, +0137560,0034324,0165024,0037021, +0037222,0045046,0047151,0161213, +0040200,0000000,0000000,0000000 +}; +#define MAXGAM 34.84425627277176174 +static unsigned short LPI[4] = { +0040222,0103202,0043475,0006750, +}; +#define LOGPI *(double *)LPI +#endif + +#ifdef IBMPC +static unsigned short P[] = { +0x2153,0x3998,0xfcb8,0x3f24, +0xbfab,0xe686,0x84e3,0x3f53, +0x14b0,0xe9db,0x57cd,0x3f85, +0x23d3,0x18c4,0x63d9,0x3fa8, +0x7d31,0xdcae,0x8da9,0x3fca, +0xe312,0x3993,0xa137,0x3fdf, +0x0000,0x0000,0x0000,0x3ff0 +}; +static unsigned short Q[] = { +0xd3af,0x8400,0x487a,0xbef8, +0x2573,0x2915,0xae8a,0x3f41, +0xb44a,0xe750,0x40e4,0xbf72, +0xb117,0x5b1b,0x31ed,0x3f88, +0xde67,0xe33f,0x5779,0x3fa2, +0x87c2,0x9d42,0x071a,0xbfce, +0x3c51,0xc9cd,0x4944,0x3fb2, +0x0000,0x0000,0x0000,0x3ff0 +}; +#define MAXGAM 171.624376956302725 +static unsigned short LPI[4] = { +0xa1bd,0x48e7,0x50d0,0x3ff2, +}; +#define LOGPI *(double *)LPI +#endif + +#ifdef MIEEE +static unsigned short P[] = { +0x3f24,0xfcb8,0x3998,0x2153, +0x3f53,0x84e3,0xe686,0xbfab, +0x3f85,0x57cd,0xe9db,0x14b0, +0x3fa8,0x63d9,0x18c4,0x23d3, +0x3fca,0x8da9,0xdcae,0x7d31, +0x3fdf,0xa137,0x3993,0xe312, +0x3ff0,0x0000,0x0000,0x0000 +}; +static unsigned short Q[] = { +0xbef8,0x487a,0x8400,0xd3af, +0x3f41,0xae8a,0x2915,0x2573, +0xbf72,0x40e4,0xe750,0xb44a, +0x3f88,0x31ed,0x5b1b,0xb117, +0x3fa2,0x5779,0xe33f,0xde67, +0xbfce,0x071a,0x9d42,0x87c2, +0x3fb2,0x4944,0xc9cd,0x3c51, +0x3ff0,0x0000,0x0000,0x0000 +}; +#define MAXGAM 171.624376956302725 +static unsigned short LPI[4] = { +0x3ff2,0x50d0,0x48e7,0xa1bd, +}; +#define LOGPI *(double *)LPI +#endif + +/* Stirling's formula for the gamma function */ +#if UNK +static double STIR[5] = { + 7.87311395793093628397E-4, +-2.29549961613378126380E-4, +-2.68132617805781232825E-3, + 3.47222221605458667310E-3, + 8.33333333333482257126E-2, +}; +#define MAXSTIR 143.01608 +static double SQTPI = 2.50662827463100050242E0; +#endif +#if DEC +static unsigned short STIR[20] = { +0035516,0061622,0144553,0112224, +0135160,0131531,0037460,0165740, +0136057,0134460,0037242,0077270, +0036143,0107070,0156306,0027751, +0037252,0125252,0125252,0146064, +}; +#define MAXSTIR 26.77 +static unsigned short SQT[4] = { +0040440,0066230,0177661,0034055, +}; +#define SQTPI *(double *)SQT +#endif +#if IBMPC +static unsigned short STIR[20] = { +0x7293,0x592d,0xcc72,0x3f49, +0x1d7c,0x27e6,0x166b,0xbf2e, +0x4fd7,0x07d4,0xf726,0xbf65, +0xc5fd,0x1b98,0x71c7,0x3f6c, +0x5986,0x5555,0x5555,0x3fb5, +}; +#define MAXSTIR 143.01608 +static unsigned short SQT[4] = { +0x2706,0x1ff6,0x0d93,0x4004, +}; +#define SQTPI *(double *)SQT +#endif +#if MIEEE +static unsigned short STIR[20] = { +0x3f49,0xcc72,0x592d,0x7293, +0xbf2e,0x166b,0x27e6,0x1d7c, +0xbf65,0xf726,0x07d4,0x4fd7, +0x3f6c,0x71c7,0x1b98,0xc5fd, +0x3fb5,0x5555,0x5555,0x5986, +}; +#define MAXSTIR 143.01608 +static unsigned short SQT[4] = { +0x4004,0x0d93,0x1ff6,0x2706, +}; +#define SQTPI *(double *)SQT +#endif + +int sgngam = 0; +extern int sgngam; +extern double MAXLOG, MAXNUM, PI; +#ifndef ANSIPROT +double pow(), log(), exp(), sin(), polevl(), p1evl(), floor(), fabs(); +int isnan(), isfinite(); +#endif +#ifdef INFINITIES +extern double INFINITY; +#endif +#ifdef NANS +extern double NAN; +#endif + +/* igam.c + * + * Incomplete gamma integral + * + * + * + * SYNOPSIS: + * + * double a, x, y, igam(); + * + * y = igam( a, x ); + * + * DESCRIPTION: + * + * The function is defined by + * + * x + * - + * 1 | | -t a-1 + * igam(a,x) = ----- | e t dt. + * - | | + * | (a) - + * 0 + * + * + * In this implementation both arguments must be positive. + * The integral is evaluated by either a power series or + * continued fraction expansion, depending on the relative + * values of a and x. + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE 0,30 200000 3.6e-14 2.9e-15 + * IEEE 0,100 300000 9.9e-14 1.5e-14 + */ +/* igamc() + * + * Complemented incomplete gamma integral + * + * + * + * SYNOPSIS: + * + * double a, x, y, igamc(); + * + * y = igamc( a, x ); + * + * DESCRIPTION: + * + * The function is defined by + * + * + * igamc(a,x) = 1 - igam(a,x) + * + * inf. + * - + * 1 | | -t a-1 + * = ----- | e t dt. + * - | | + * | (a) - + * x + * + * + * In this implementation both arguments must be positive. + * The integral is evaluated by either a power series or + * continued fraction expansion, depending on the relative + * values of a and x. + * + * ACCURACY: + * + * Tested at random a, x. + * a x Relative error: + * arithmetic domain domain # trials peak rms + * IEEE 0.5,100 0,100 200000 1.9e-14 1.7e-15 + * IEEE 0.01,0.5 0,100 200000 1.4e-13 1.6e-15 + */ + +/* +Cephes Math Library Release 2.0: April, 1987 +Copyright 1985, 1987 by Stephen L. Moshier +Direct inquiries to 30 Frost Street, Cambridge, MA 02140 +*/ + +double igamc( a, x ) +double a, x; +{ +double ans, ax, c, yc, r, t, y, z; +double pk, pkm1, pkm2, qk, qkm1, qkm2; + +if( (x <= 0) || ( a <= 0) ) + return( 1.0 ); + +if( (x < 1.0) || (x < a) ) + return( 1.e0 - igam(a,x) ); + +ax = a * log(x) - x - lgam(a); +if( ax < -MAXLOG ) + { + mtherr( "igamc", UNDERFLOW ); + return( 0.0 ); + } +ax = exp(ax); + +/* continued fraction */ +y = 1.0 - a; +z = x + y + 1.0; +c = 0.0; +pkm2 = 1.0; +qkm2 = x; +pkm1 = x + 1.0; +qkm1 = z * x; +ans = pkm1/qkm1; + +do + { + c += 1.0; + y += 1.0; + z += 2.0; + yc = y * c; + pk = pkm1 * z - pkm2 * yc; + qk = qkm1 * z - qkm2 * yc; + if( qk != 0 ) + { + r = pk/qk; + t = fabs( (ans - r)/r ); + ans = r; + } + else + t = 1.0; + pkm2 = pkm1; + pkm1 = pk; + qkm2 = qkm1; + qkm1 = qk; + if( fabs(pk) > big ) + { + pkm2 *= biginv; + pkm1 *= biginv; + qkm2 *= biginv; + qkm1 *= biginv; + } + } +while( t > MACHEP ); + +return( ans * ax ); +} + + + +/* left tail of incomplete gamma function: + * + * inf. k + * a -x - x + * x e > ---------- + * - - + * k=0 | (a+k+1) + * + */ + +double igam( a, x ) +double a, x; +{ +double ans, ax, c, r; + +if( (x <= 0) || ( a <= 0) ) + return( 0.0 ); + +if( (x > 1.0) && (x > a ) ) + return( 1.e0 - igamc(a,x) ); + +/* Compute x**a * exp(-x) / gamma(a) */ +ax = a * log(x) - x - lgam(a); +if( ax < -MAXLOG ) + { + mtherr( "igam", UNDERFLOW ); + return( 0.0 ); + } +ax = exp(ax); + +/* power series */ +r = a; +c = 1.0; +ans = 1.0; + +do + { + r += 1.0; + c *= x/r; + ans += c; + } +while( c/ans > MACHEP ); + +return( ans * ax/a ); +} + +/* gamma.c + * + * Gamma function + * + * + * + * SYNOPSIS: + * + * double x, y, gamma(); + * extern int sgngam; + * + * y = gamma( x ); + * + * + * + * DESCRIPTION: + * + * Returns gamma function of the argument. The result is + * correctly signed, and the sign (+1 or -1) is also + * returned in a global (extern) variable named sgngam. + * This variable is also filled in by the logarithmic gamma + * function lgam(). + * + * Arguments |x| <= 34 are reduced by recurrence and the function + * approximated by a rational function of degree 6/7 in the + * interval (2,3). Large arguments are handled by Stirling's + * formula. Large negative arguments are made positive using + * a reflection formula. + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC -34, 34 10000 1.3e-16 2.5e-17 + * IEEE -170,-33 20000 2.3e-15 3.3e-16 + * IEEE -33, 33 20000 9.4e-16 2.2e-16 + * IEEE 33, 171.6 20000 2.3e-15 3.2e-16 + * + * Error for arguments outside the test range will be larger + * owing to error amplification by the exponential function. + * + */ +/* lgam() + * + * Natural logarithm of gamma function + * + * + * + * SYNOPSIS: + * + * double x, y, lgam(); + * extern int sgngam; + * + * y = lgam( x ); + * + * + * + * DESCRIPTION: + * + * Returns the base e (2.718...) logarithm of the absolute + * value of the gamma function of the argument. + * The sign (+1 or -1) of the gamma function is returned in a + * global (extern) variable named sgngam. + * + * For arguments greater than 13, the logarithm of the gamma + * function is approximated by the logarithmic version of + * Stirling's formula using a polynomial approximation of + * degree 4. Arguments between -33 and +33 are reduced by + * recurrence to the interval [2,3] of a rational approximation. + * The cosecant reflection formula is employed for arguments + * less than -33. + * + * Arguments greater than MAXLGM return MAXNUM and an error + * message. MAXLGM = 2.035093e36 for DEC + * arithmetic or 2.556348e305 for IEEE arithmetic. + * + * + * + * ACCURACY: + * + * + * arithmetic domain # trials peak rms + * DEC 0, 3 7000 5.2e-17 1.3e-17 + * DEC 2.718, 2.035e36 5000 3.9e-17 9.9e-18 + * IEEE 0, 3 28000 5.4e-16 1.1e-16 + * IEEE 2.718, 2.556e305 40000 3.5e-16 8.3e-17 + * The error criterion was relative when the function magnitude + * was greater than one but absolute when it was less than one. + * + * The following test used the relative error criterion, though + * at certain points the relative error could be much higher than + * indicated. + * IEEE -200, -4 10000 4.8e-16 1.3e-16 + * + */ + +/* gamma.c */ +/* gamma function */ + +/* +Cephes Math Library Release 2.2: July, 1992 +Copyright 1984, 1987, 1989, 1992 by Stephen L. Moshier +Direct inquiries to 30 Frost Street, Cambridge, MA 02140 +*/ + +/* Gamma function computed by Stirling's formula. + * The polynomial STIR is valid for 33 <= x <= 172. + */ +static double stirf(x) +double x; +{ +double y, w, v; + +w = 1.0/x; +w = 1.0 + w * polevl( w, STIR, 4 ); +y = exp(x); +if( x > MAXSTIR ) + { /* Avoid overflow in pow() */ + v = pow( x, 0.5 * x - 0.25 ); + y = v * (v / y); + } +else + { + y = pow( x, x - 0.5 ) / y; + } +y = SQTPI * y * w; +return( y ); +} + + + +double gamma(x) +double x; +{ +double p, q, z; +int i; + +sgngam = 1; +#ifdef NANS +if( isnan(x) ) + return(x); +#endif +#ifdef INFINITIES +#ifdef NANS +if( x == INFINITY ) + return(x); +if( x == -INFINITY ) + return(NAN); +#else +if( !isfinite(x) ) + return(x); +#endif +#endif +q = fabs(x); + +if( q > 33.0 ) + { + if( x < 0.0 ) + { + p = floor(q); + if( p == q ) + { +#ifdef NANS +gamnan: + mtherr( "gamma", DOMAIN ); + return (NAN); +#else + goto goverf; +#endif + } + i = p; + if( (i & 1) == 0 ) + sgngam = -1; + z = q - p; + if( z > 0.5 ) + { + p += 1.0; + z = q - p; + } + z = q * sin( PI * z ); + if( z == 0.0 ) + { +#ifdef INFINITIES + return( sgngam * INFINITY); +#else +goverf: + mtherr( "gamma", OVERFLOW ); + return( sgngam * MAXNUM); +#endif + } + z = fabs(z); + z = PI/(z * stirf(q) ); + } + else + { + z = stirf(x); + } + return( sgngam * z ); + } + +z = 1.0; +while( x >= 3.0 ) + { + x -= 1.0; + z *= x; + } + +while( x < 0.0 ) + { + if( x > -1.E-9 ) + goto small; + z /= x; + x += 1.0; + } + +while( x < 2.0 ) + { + if( x < 1.e-9 ) + goto small; + z /= x; + x += 1.0; + } + +if( x == 2.0 ) + return(z); + +x -= 2.0; +p = polevl( x, P, 6 ); +q = polevl( x, Q, 7 ); +return( z * p / q ); + +small: +if( x == 0.0 ) + { +#ifdef INFINITIES +#ifdef NANS + goto gamnan; +#else + return( INFINITY ); +#endif +#else + mtherr( "gamma", SING ); + return( MAXNUM ); +#endif + } +else + return( z/((1.0 + 0.5772156649015329 * x) * x) ); +} + + + +/* A[]: Stirling's formula expansion of log gamma + * B[], C[]: log gamma function between 2 and 3 + */ +#ifdef UNK +static double A[] = { + 8.11614167470508450300E-4, +-5.95061904284301438324E-4, + 7.93650340457716943945E-4, +-2.77777777730099687205E-3, + 8.33333333333331927722E-2 +}; +static double B[] = { +-1.37825152569120859100E3, +-3.88016315134637840924E4, +-3.31612992738871184744E5, +-1.16237097492762307383E6, +-1.72173700820839662146E6, +-8.53555664245765465627E5 +}; +static double C[] = { +/* 1.00000000000000000000E0, */ +-3.51815701436523470549E2, +-1.70642106651881159223E4, +-2.20528590553854454839E5, +-1.13933444367982507207E6, +-2.53252307177582951285E6, +-2.01889141433532773231E6 +}; +/* log( sqrt( 2*pi ) ) */ +static double LS2PI = 0.91893853320467274178; +#define MAXLGM 2.556348e305 +#endif + +#ifdef DEC +static unsigned short A[] = { +0035524,0141201,0034633,0031405, +0135433,0176755,0126007,0045030, +0035520,0006371,0003342,0172730, +0136066,0005540,0132605,0026407, +0037252,0125252,0125252,0125132 +}; +static unsigned short B[] = { +0142654,0044014,0077633,0035410, +0144027,0110641,0125335,0144760, +0144641,0165637,0142204,0047447, +0145215,0162027,0146246,0155211, +0145322,0026110,0010317,0110130, +0145120,0061472,0120300,0025363 +}; +static unsigned short C[] = { +/*0040200,0000000,0000000,0000000*/ +0142257,0164150,0163630,0112622, +0143605,0050153,0156116,0135272, +0144527,0056045,0145642,0062332, +0145213,0012063,0106250,0001025, +0145432,0111254,0044577,0115142, +0145366,0071133,0050217,0005122 +}; +/* log( sqrt( 2*pi ) ) */ +static unsigned short LS2P[] = {040153,037616,041445,0172645,}; +#define LS2PI *(double *)LS2P +#define MAXLGM 2.035093e36 +#endif + +#ifdef IBMPC +static unsigned short A[] = { +0x6661,0x2733,0x9850,0x3f4a, +0xe943,0xb580,0x7fbd,0xbf43, +0x5ebb,0x20dc,0x019f,0x3f4a, +0xa5a1,0x16b0,0xc16c,0xbf66, +0x554b,0x5555,0x5555,0x3fb5 +}; +static unsigned short B[] = { +0x6761,0x8ff3,0x8901,0xc095, +0xb93e,0x355b,0xf234,0xc0e2, +0x89e5,0xf890,0x3d73,0xc114, +0xdb51,0xf994,0xbc82,0xc131, +0xf20b,0x0219,0x4589,0xc13a, +0x055e,0x5418,0x0c67,0xc12a +}; +static unsigned short C[] = { +/*0x0000,0x0000,0x0000,0x3ff0,*/ +0x12b2,0x1cf3,0xfd0d,0xc075, +0xd757,0x7b89,0xaa0d,0xc0d0, +0x4c9b,0xb974,0xeb84,0xc10a, +0x0043,0x7195,0x6286,0xc131, +0xf34c,0x892f,0x5255,0xc143, +0xe14a,0x6a11,0xce4b,0xc13e +}; +/* log( sqrt( 2*pi ) ) */ +static unsigned short LS2P[] = { +0xbeb5,0xc864,0x67f1,0x3fed +}; +#define LS2PI *(double *)LS2P +#define MAXLGM 2.556348e305 +#endif + +#ifdef MIEEE +static unsigned short A[] = { +0x3f4a,0x9850,0x2733,0x6661, +0xbf43,0x7fbd,0xb580,0xe943, +0x3f4a,0x019f,0x20dc,0x5ebb, +0xbf66,0xc16c,0x16b0,0xa5a1, +0x3fb5,0x5555,0x5555,0x554b +}; +static unsigned short B[] = { +0xc095,0x8901,0x8ff3,0x6761, +0xc0e2,0xf234,0x355b,0xb93e, +0xc114,0x3d73,0xf890,0x89e5, +0xc131,0xbc82,0xf994,0xdb51, +0xc13a,0x4589,0x0219,0xf20b, +0xc12a,0x0c67,0x5418,0x055e +}; +static unsigned short C[] = { +0xc075,0xfd0d,0x1cf3,0x12b2, +0xc0d0,0xaa0d,0x7b89,0xd757, +0xc10a,0xeb84,0xb974,0x4c9b, +0xc131,0x6286,0x7195,0x0043, +0xc143,0x5255,0x892f,0xf34c, +0xc13e,0xce4b,0x6a11,0xe14a +}; +/* log( sqrt( 2*pi ) ) */ +static unsigned short LS2P[] = { +0x3fed,0x67f1,0xc864,0xbeb5 +}; +#define LS2PI *(double *)LS2P +#define MAXLGM 2.556348e305 +#endif + + +/* Logarithm of gamma function */ + + +double lgam(x) +double x; +{ +double p, q, u, w, z; +int i; + +sgngam = 1; +#ifdef NANS +if( isnan(x) ) + return(x); +#endif + +#ifdef INFINITIES +if( !isfinite(x) ) + return(INFINITY); +#endif + +if( x < -34.0 ) + { + q = -x; + w = lgam(q); /* note this modifies sgngam! */ + p = floor(q); + if( p == q ) + { +lgsing: +#ifdef INFINITIES + mtherr( "lgam", SING ); + return (INFINITY); +#else + goto loverf; +#endif + } + i = p; + if( (i & 1) == 0 ) + sgngam = -1; + else + sgngam = 1; + z = q - p; + if( z > 0.5 ) + { + p += 1.0; + z = p - q; + } + z = q * sin( PI * z ); + if( z == 0.0 ) + goto lgsing; +/* z = log(PI) - log( z ) - w;*/ + z = LOGPI - log( z ) - w; + return( z ); + } + +if( x < 13.0 ) + { + z = 1.0; + p = 0.0; + u = x; + while( u >= 3.0 ) + { + p -= 1.0; + u = x + p; + z *= u; + } + while( u < 2.0 ) + { + if( u == 0.0 ) + goto lgsing; + z /= u; + p += 1.0; + u = x + p; + } + if( z < 0.0 ) + { + sgngam = -1; + z = -z; + } + else + sgngam = 1; + if( u == 2.0 ) + return( log(z) ); + p -= 2.0; + x = x + p; + p = x * polevl( x, B, 5 ) / p1evl( x, C, 6); + return( log(z) + p ); + } + +if( x > MAXLGM ) + { +#ifdef INFINITIES + return( sgngam * INFINITY ); +#else +loverf: + mtherr( "lgam", OVERFLOW ); + return( sgngam * MAXNUM ); +#endif + } + +q = ( x - 0.5 ) * log(x) - x + LS2PI; +if( x > 1.0e8 ) + return( q ); + +p = 1.0/(x*x); +if( x >= 1000.0 ) + q += (( 7.9365079365079365079365e-4 * p + - 2.7777777777777777777778e-3) *p + + 0.0833333333333333333333) / x; +else + q += polevl( p, A, 4 ) / x; +return( q ); +} + +/* mtherr.c + * + * Library common error handling routine + * + * + * + * SYNOPSIS: + * + * char *fctnam; + * int code; + * int mtherr(); + * + * mtherr( fctnam, code ); + * + * + * + * DESCRIPTION: + * + * This routine may be called to report one of the following + * error conditions (in the include file mconf.h). + * + * Mnemonic Value Significance + * + * DOMAIN 1 argument domain error + * SING 2 function singularity + * OVERFLOW 3 overflow range error + * UNDERFLOW 4 underflow range error + * TLOSS 5 total loss of precision + * PLOSS 6 partial loss of precision + * EDOM 33 Unix domain error code + * ERANGE 34 Unix range error code + * + * The default version of the file prints the function name, + * passed to it by the pointer fctnam, followed by the + * error condition. The display is directed to the standard + * output device. The routine then returns to the calling + * program. Users may wish to modify the program to abort by + * calling exit() under severe error conditions such as domain + * errors. + * + * Since all error conditions pass control to this function, + * the display may be easily changed, eliminated, or directed + * to an error logging device. + * + * SEE ALSO: + * + * mconf.h + * + */ + +/* +Cephes Math Library Release 2.0: April, 1987 +Copyright 1984, 1987 by Stephen L. Moshier +Direct inquiries to 30 Frost Street, Cambridge, MA 02140 +*/ + +int mtherr( name, code ) +char *name; +int code; +{ + +/* Display string passed by calling program, + * which is supposed to be the name of the + * function in which the error occurred: + */ +printf("\n%s ", name); + +/* Set global error message word */ +merror = code; + +/* Display error message defined + * by the code argument. + */ +if( (code <= 0) || (code >= 7) ) + code = 0; +printf( "%s error\n", ermsg[code] ); + +/* Return to calling + * program + */ +return( 0 ); +} + +/* polevl.c + * p1evl.c + * + * Evaluate polynomial + * + * + * + * SYNOPSIS: + * + * int N; + * double x, y, coef[N+1], polevl[]; + * + * y = polevl( x, coef, N ); + * + * + * + * DESCRIPTION: + * + * Evaluates polynomial of degree N: + * + * 2 N + * y = C + C x + C x +...+ C x + * 0 1 2 N + * + * Coefficients are stored in reverse order: + * + * coef[0] = C , ..., coef[N] = C . + * N 0 + * + * The function p1evl() assumes that coef[N] = 1.0 and is + * omitted from the array. Its calling arguments are + * otherwise the same as polevl(). + * + * + * SPEED: + * + * In the interest of speed, there are no checks for out + * of bounds arithmetic. This routine is used by most of + * the functions in the library. Depending on available + * equipment features, the user may wish to rewrite the + * program in microcode or assembly language. + * + */ + + +/* +Cephes Math Library Release 2.1: December, 1988 +Copyright 1984, 1987, 1988 by Stephen L. Moshier +Direct inquiries to 30 Frost Street, Cambridge, MA 02140 +*/ + +double polevl( x, coef, N ) +double x; +double coef[]; +int N; +{ +double ans; +int i; +double *p; + +p = coef; +ans = *p++; +i = N; + +do + ans = ans * x + *p++; +while( --i ); + +return( ans ); +} + +/* p1evl() */ +/* N + * Evaluate polynomial when coefficient of x is 1.0. + * Otherwise same as polevl. + */ + +double p1evl( x, coef, N ) +double x; +double coef[]; +int N; +{ +double ans; +double *p; +int i; + +p = coef; +ans = x + *p++; +i = N-1; + +do + ans = ans * x + *p++; +while( --i ); + +return( ans ); +} diff --git a/nist/dfft.c b/nist/dfft.c new file mode 100644 index 0000000..482a868 --- /dev/null +++ b/nist/dfft.c @@ -0,0 +1,1381 @@ +/* Notes from RFB: + + Looks like the user-level routines are: + + Real FFT + + void __ogg_fdrffti(int n, double *wsave, int *ifac) + void __ogg_fdrfftf(int n,double *r,double *wsave,int *ifac) + void __ogg_fdrfftb(int n, double *r, double *wsave, int *ifac) + + __ogg_fdrffti == initialization + __ogg_fdrfftf == forward transform + __ogg_fdrfftb == backward transform + + Parameters are + n == length of sequence + r == sequence to be transformed (input) + == transformed sequence (output) + wsave == work array of length 2n (allocated by caller) + ifac == work array of length 15 (allocated by caller) + + Cosine quarter-wave FFT + + void __ogg_fdcosqi(int n, double *wsave, int *ifac) + void __ogg_fdcosqf(int n,double *x,double *wsave,int *ifac) + void __ogg_fdcosqb(int n,double *x,double *wsave,int *ifac) +*/ + +/******************************************************************** + * * + * THIS FILE IS PART OF THE OggSQUISH SOFTWARE CODEC SOURCE CODE. * + * * + ******************************************************************** + + file: fft.c + function: Fast discrete Fourier and cosine transforms and inverses + author: Monty <xiphmont@mit.edu> + modifications by: Monty + last modification date: Jul 1 1996 + + ********************************************************************/ + +/* These Fourier routines were originally based on the Fourier + routines of the same names from the NETLIB bihar and fftpack + fortran libraries developed by Paul N. Swarztrauber at the National + Center for Atmospheric Research in Boulder, CO USA. They have been + reimplemented in C and optimized in a few ways for OggSquish. */ + +/* As the original fortran libraries are public domain, the C Fourier + routines in this file are hereby released to the public domain as + well. The C routines here produce output exactly equivalent to the + original fortran routines. Of particular interest are the facts + that (like the original fortran), these routines can work on + arbitrary length vectors that need not be powers of two in + length. */ + +#include <math.h> +#define STIN static + +static void drfti1(int n, double *wa, int *ifac){ + static int ntryh[4] = { 4,2,3,5 }; + static double tpi = 6.28318530717958647692528676655900577; + double arg,argh,argld,fi; + int ntry=0,i,j=-1; + int k1, l1, l2, ib; + int ld, ii, ip, is, nq, nr; + int ido, ipm, nfm1; + int nl=n; + int nf=0; + + L101: + j++; + if (j < 4) + ntry=ntryh[j]; + else + ntry+=2; + + L104: + nq=nl/ntry; + nr=nl-ntry*nq; + if (nr!=0) goto L101; + + nf++; + ifac[nf+1]=ntry; + nl=nq; + if(ntry!=2)goto L107; + if(nf==1)goto L107; + + for (i=1;i<nf;i++){ + ib=nf-i+1; + ifac[ib+1]=ifac[ib]; + } + ifac[2] = 2; + + L107: + if(nl!=1)goto L104; + ifac[0]=n; + ifac[1]=nf; + argh=tpi/n; + is=0; + nfm1=nf-1; + l1=1; + + if(nfm1==0)return; + + for (k1=0;k1<nfm1;k1++){ + ip=ifac[k1+2]; + ld=0; + l2=l1*ip; + ido=n/l2; + ipm=ip-1; + + for (j=0;j<ipm;j++){ + ld+=l1; + i=is; + argld=(double)ld*argh; + fi=0.; + for (ii=2;ii<ido;ii+=2){ + fi+=1.; + arg=fi*argld; + wa[i++]=cos(arg); + wa[i++]=sin(arg); + } + is+=ido; + } + l1=l2; + } +} + +void __ogg_fdrffti(int n, double *wsave, int *ifac){ + + if (n == 1) return; + drfti1(n, wsave+n, ifac); +} + +void __ogg_fdcosqi(int n, double *wsave, int *ifac){ + static double pih = 1.57079632679489661923132169163975; + static int k; + static double fk, dt; + + dt=pih/n; + fk=0.; + for(k=0;k<n;k++){ + fk+=1.; + wsave[k] = cos(fk*dt); + } + + __ogg_fdrffti(n, wsave+n,ifac); +} + +STIN void dradf2(int ido,int l1,double *cc,double *ch,double *wa1){ + int i,k; + double ti2,tr2; + int t0,t1,t2,t3,t4,t5,t6; + + t1=0; + t0=(t2=l1*ido); + t3=ido<<1; + for(k=0;k<l1;k++){ + ch[t1<<1]=cc[t1]+cc[t2]; + ch[(t1<<1)+t3-1]=cc[t1]-cc[t2]; + t1+=ido; + t2+=ido; + } + + if(ido<2)return; + if(ido==2)goto L105; + + t1=0; + t2=t0; + for(k=0;k<l1;k++){ + t3=t2; + t4=(t1<<1)+(ido<<1); + t5=t1; + t6=t1+t1; + for(i=2;i<ido;i+=2){ + t3+=2; + t4-=2; + t5+=2; + t6+=2; + tr2=wa1[i-2]*cc[t3-1]+wa1[i-1]*cc[t3]; + ti2=wa1[i-2]*cc[t3]-wa1[i-1]*cc[t3-1]; + ch[t6]=cc[t5]+ti2; + ch[t4]=ti2-cc[t5]; + ch[t6-1]=cc[t5-1]+tr2; + ch[t4-1]=cc[t5-1]-tr2; + } + t1+=ido; + t2+=ido; + } + + if(ido%2==1)return; + + L105: + t3=(t2=(t1=ido)-1); + t2+=t0; + for(k=0;k<l1;k++){ + ch[t1]=-cc[t2]; + ch[t1-1]=cc[t3]; + t1+=ido<<1; + t2+=ido; + t3+=ido; + } +} + +STIN void dradf4(int ido,int l1,double *cc,double *ch,double *wa1, + double *wa2,double *wa3){ + static double hsqt2 = .70710678118654752440084436210485; + int i,k,t0,t1,t2,t3,t4,t5,t6; + double ci2,ci3,ci4,cr2,cr3,cr4,ti1,ti2,ti3,ti4,tr1,tr2,tr3,tr4; + t0=l1*ido; + + t1=t0; + t4=t1<<1; + t2=t1+(t1<<1); + t3=0; + + for(k=0;k<l1;k++){ + tr1=cc[t1]+cc[t2]; + tr2=cc[t3]+cc[t4]; + ch[t5=t3<<2]=tr1+tr2; + ch[(ido<<2)+t5-1]=tr2-tr1; + ch[(t5+=(ido<<1))-1]=cc[t3]-cc[t4]; + ch[t5]=cc[t2]-cc[t1]; + + t1+=ido; + t2+=ido; + t3+=ido; + t4+=ido; + } + + if(ido<2)return; + if(ido==2)goto L105; + + t1=0; + for(k=0;k<l1;k++){ + t2=t1; + t4=t1<<2; + t5=(t6=ido<<1)+t4; + for(i=2;i<ido;i+=2){ + t3=(t2+=2); + t4+=2; + t5-=2; + + t3+=t0; + cr2=wa1[i-2]*cc[t3-1]+wa1[i-1]*cc[t3]; + ci2=wa1[i-2]*cc[t3]-wa1[i-1]*cc[t3-1]; + t3+=t0; + cr3=wa2[i-2]*cc[t3-1]+wa2[i-1]*cc[t3]; + ci3=wa2[i-2]*cc[t3]-wa2[i-1]*cc[t3-1]; + t3+=t0; + cr4=wa3[i-2]*cc[t3-1]+wa3[i-1]*cc[t3]; + ci4=wa3[i-2]*cc[t3]-wa3[i-1]*cc[t3-1]; + + tr1=cr2+cr4; + tr4=cr4-cr2; + ti1=ci2+ci4; + ti4=ci2-ci4; + ti2=cc[t2]+ci3; + ti3=cc[t2]-ci3; + tr2=cc[t2-1]+cr3; + tr3=cc[t2-1]-cr3; + + + ch[t4-1]=tr1+tr2; + ch[t4]=ti1+ti2; + + ch[t5-1]=tr3-ti4; + ch[t5]=tr4-ti3; + + ch[t4+t6-1]=ti4+tr3; + ch[t4+t6]=tr4+ti3; + + ch[t5+t6-1]=tr2-tr1; + ch[t5+t6]=ti1-ti2; + } + t1+=ido; + } + if(ido%2==1)return; + + L105: + + t2=(t1=t0+ido-1)+(t0<<1); + t3=ido<<2; + t4=ido; + t5=ido<<1; + t6=ido; + + for(k=0;k<l1;k++){ + ti1=-hsqt2*(cc[t1]+cc[t2]); + tr1=hsqt2*(cc[t1]-cc[t2]); + ch[t4-1]=tr1+cc[t6-1]; + ch[t4+t5-1]=cc[t6-1]-tr1; + ch[t4]=ti1-cc[t1+t0]; + ch[t4+t5]=ti1+cc[t1+t0]; + t1+=ido; + t2+=ido; + t4+=t3; + t6+=ido; + } +} + +STIN void dradfg(int ido,int ip,int l1,int idl1,double *cc,double *c1, + double *c2,double *ch,double *ch2,double *wa){ + + static double tpi=6.28318530717958647692528676655900577; + int idij,ipph,i,j,k,l,ic,ik,is; + int t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10; + double dc2,ai1,ai2,ar1,ar2,ds2; + int nbd; + double dcp,arg,dsp,ar1h,ar2h; + int idp2,ipp2; + + arg=tpi/(double)ip; + dcp=cos(arg); + dsp=sin(arg); + ipph=(ip+1)>>1; + ipp2=ip; + idp2=ido; + nbd=(ido-1)>>1; + t0=l1*ido; + t10=ip*ido; + + if(ido==1)goto L119; + for(ik=0;ik<idl1;ik++)ch2[ik]=c2[ik]; + + t1=0; + for(j=1;j<ip;j++){ + t1+=t0; + t2=t1; + for(k=0;k<l1;k++){ + ch[t2]=c1[t2]; + t2+=ido; + } + } + + is=-ido; + t1=0; + if(nbd>l1){ + for(j=1;j<ip;j++){ + t1+=t0; + is+=ido; + t2= -ido+t1; + for(k=0;k<l1;k++){ + idij=is-1; + t2+=ido; + t3=t2; + for(i=2;i<ido;i+=2){ + idij+=2; + t3+=2; + ch[t3-1]=wa[idij-1]*c1[t3-1]+wa[idij]*c1[t3]; + ch[t3]=wa[idij-1]*c1[t3]-wa[idij]*c1[t3-1]; + } + } + } + }else{ + + for(j=1;j<ip;j++){ + is+=ido; + idij=is-1; + t1+=t0; + t2=t1; + for(i=2;i<ido;i+=2){ + idij+=2; + t2+=2; + t3=t2; + for(k=0;k<l1;k++){ + ch[t3-1]=wa[idij-1]*c1[t3-1]+wa[idij]*c1[t3]; + ch[t3]=wa[idij-1]*c1[t3]-wa[idij]*c1[t3-1]; + t3+=ido; + } + } + } + } + + t1=0; + t2=ipp2*t0; + if(nbd<l1){ + for(j=1;j<ipph;j++){ + t1+=t0; + t2-=t0; + t3=t1; + t4=t2; + for(i=2;i<ido;i+=2){ + t3+=2; + t4+=2; + t5=t3-ido; + t6=t4-ido; + for(k=0;k<l1;k++){ + t5+=ido; + t6+=ido; + c1[t5-1]=ch[t5-1]+ch[t6-1]; + c1[t6-1]=ch[t5]-ch[t6]; + c1[t5]=ch[t5]+ch[t6]; + c1[t6]=ch[t6-1]-ch[t5-1]; + } + } + } + }else{ + for(j=1;j<ipph;j++){ + t1+=t0; + t2-=t0; + t3=t1; + t4=t2; + for(k=0;k<l1;k++){ + t5=t3; + t6=t4; + for(i=2;i<ido;i+=2){ + t5+=2; + t6+=2; + c1[t5-1]=ch[t5-1]+ch[t6-1]; + c1[t6-1]=ch[t5]-ch[t6]; + c1[t5]=ch[t5]+ch[t6]; + c1[t6]=ch[t6-1]-ch[t5-1]; + } + t3+=ido; + t4+=ido; + } + } + } + +L119: + for(ik=0;ik<idl1;ik++)c2[ik]=ch2[ik]; + + t1=0; + t2=ipp2*idl1; + for(j=1;j<ipph;j++){ + t1+=t0; + t2-=t0; + t3=t1-ido; + t4=t2-ido; + for(k=0;k<l1;k++){ + t3+=ido; + t4+=ido; + c1[t3]=ch[t3]+ch[t4]; + c1[t4]=ch[t4]-ch[t3]; + } + } + + ar1=1.; + ai1=0.; + t1=0; + t2=ipp2*idl1; + t3=(ip-1)*idl1; + for(l=1;l<ipph;l++){ + t1+=idl1; + t2-=idl1; + ar1h=dcp*ar1-dsp*ai1; + ai1=dcp*ai1+dsp*ar1; + ar1=ar1h; + t4=t1; + t5=t2; + t6=t3; + t7=idl1; + + for(ik=0;ik<idl1;ik++){ + ch2[t4++]=c2[ik]+ar1*c2[t7++]; + ch2[t5++]=ai1*c2[t6++]; + } + + dc2=ar1; + ds2=ai1; + ar2=ar1; + ai2=ai1; + + t4=idl1; + t5=(ipp2-1)*idl1; + for(j=2;j<ipph;j++){ + t4+=idl1; + t5-=idl1; + + ar2h=dc2*ar2-ds2*ai2; + ai2=dc2*ai2+ds2*ar2; + ar2=ar2h; + + t6=t1; + t7=t2; + t8=t4; + t9=t5; + for(ik=0;ik<idl1;ik++){ + ch2[t6++]+=ar2*c2[t8++]; + ch2[t7++]+=ai2*c2[t9++]; + } + } + } + + t1=0; + for(j=1;j<ipph;j++){ + t1+=idl1; + t2=t1; + for(ik=0;ik<idl1;ik++)ch2[ik]+=c2[t2++]; + } + + if(ido<l1)goto L132; + + t1=0; + t2=0; + for(k=0;k<l1;k++){ + t3=t1; + t4=t2; + for(i=0;i<ido;i++)cc[t4++]=ch[t3++]; + t1+=ido; + t2+=t10; + } + + goto L135; + + L132: + for(i=0;i<ido;i++){ + t1=i; + t2=i; + for(k=0;k<l1;k++){ + cc[t2]=ch[t1]; + t1+=ido; + t2+=t10; + } + } + + L135: + t1=0; + t2=ido<<1; + t3=0; + t4=ipp2*t0; + for(j=1;j<ipph;j++){ + + t1+=t2; + t3+=t0; + t4-=t0; + + t5=t1; + t6=t3; + t7=t4; + + for(k=0;k<l1;k++){ + cc[t5-1]=ch[t6]; + cc[t5]=ch[t7]; + t5+=t10; + t6+=ido; + t7+=ido; + } + } + + if(ido==1)return; + if(nbd<l1)goto L141; + + t1=-ido; + t3=0; + t4=0; + t5=ipp2*t0; + for(j=1;j<ipph;j++){ + t1+=t2; + t3+=t2; + t4+=t0; + t5-=t0; + t6=t1; + t7=t3; + t8=t4; + t9=t5; + for(k=0;k<l1;k++){ + for(i=2;i<ido;i+=2){ + ic=idp2-i; + cc[i+t7-1]=ch[i+t8-1]+ch[i+t9-1]; + cc[ic+t6-1]=ch[i+t8-1]-ch[i+t9-1]; + cc[i+t7]=ch[i+t8]+ch[i+t9]; + cc[ic+t6]=ch[i+t9]-ch[i+t8]; + } + t6+=t10; + t7+=t10; + t8+=ido; + t9+=ido; + } + } + return; + + L141: + + t1=-ido; + t3=0; + t4=0; + t5=ipp2*t0; + for(j=1;j<ipph;j++){ + t1+=t2; + t3+=t2; + t4+=t0; + t5-=t0; + for(i=2;i<ido;i+=2){ + t6=idp2+t1-i; + t7=i+t3; + t8=i+t4; + t9=i+t5; + for(k=0;k<l1;k++){ + cc[t7-1]=ch[t8-1]+ch[t9-1]; + cc[t6-1]=ch[t8-1]-ch[t9-1]; + cc[t7]=ch[t8]+ch[t9]; + cc[t6]=ch[t9]-ch[t8]; + t6+=t10; + t7+=t10; + t8+=ido; + t9+=ido; + } + } + } +} + +STIN void drftf1(int n,double *c,double *ch,double *wa,int *ifac){ + int i,k1,l1,l2; + int na,kh,nf; + int ip,iw,ido,idl1,ix2,ix3; + + nf=ifac[1]; + na=1; + l2=n; + iw=n; + + for(k1=0;k1<nf;k1++){ + kh=nf-k1; + ip=ifac[kh+1]; + l1=l2/ip; + ido=n/l2; + idl1=ido*l1; + iw-=(ip-1)*ido; + na=1-na; + + if(ip!=4)goto L102; + + ix2=iw+ido; + ix3=ix2+ido; + if(na!=0) + dradf4(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1); + else + dradf4(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1); + goto L110; + + L102: + if(ip!=2)goto L104; + if(na!=0)goto L103; + + dradf2(ido,l1,c,ch,wa+iw-1); + goto L110; + + L103: + dradf2(ido,l1,ch,c,wa+iw-1); + goto L110; + + L104: + if(ido==1)na=1-na; + if(na!=0)goto L109; + + dradfg(ido,ip,l1,idl1,c,c,c,ch,ch,wa+iw-1); + na=1; + goto L110; + + L109: + dradfg(ido,ip,l1,idl1,ch,ch,ch,c,c,wa+iw-1); + na=0; + + L110: + l2=l1; + } + + if(na==1)return; + + for(i=0;i<n;i++)c[i]=ch[i]; +} + +void __ogg_fdrfftf(int n,double *r,double *wsave,int *ifac){ + if(n==1)return; + drftf1(n,r,wsave,wsave+n,ifac); +} + +STIN void dcsqf1(int n,double *x,double *w,double *xh,int *ifac){ + int modn,i,k,kc; + int np2,ns2; + double xim1; + + ns2=(n+1)>>1; + np2=n; + + kc=np2; + for(k=1;k<ns2;k++){ + kc--; + xh[k]=x[k]+x[kc]; + xh[kc]=x[k]-x[kc]; + } + + modn=n%2; + if(modn==0)xh[ns2]=x[ns2]+x[ns2]; + + for(k=1;k<ns2;k++){ + kc=np2-k; + x[k]=w[k-1]*xh[kc]+w[kc-1]*xh[k]; + x[kc]=w[k-1]*xh[k]-w[kc-1]*xh[kc]; + } + + if(modn==0)x[ns2]=w[ns2-1]*xh[ns2]; + + __ogg_fdrfftf(n,x,xh,ifac); + + for(i=2;i<n;i+=2){ + xim1=x[i-1]-x[i]; + x[i]=x[i-1]+x[i]; + x[i-1]=xim1; + } +} + +void __ogg_fdcosqf(int n,double *x,double *wsave,int *ifac){ + static double sqrt2=1.4142135623730950488016887242097; + double tsqx; + + switch(n){ + case 0:case 1: + return; + case 2: + tsqx=sqrt2*x[1]; + x[1]=x[0]-tsqx; + x[0]+=tsqx; + return; + default: + dcsqf1(n,x,wsave,wsave+n,ifac); + return; + } +} + +STIN void dradb2(int ido,int l1,double *cc,double *ch,double *wa1){ + int i,k,t0,t1,t2,t3,t4,t5,t6; + double ti2,tr2; + + t0=l1*ido; + + t1=0; + t2=0; + t3=(ido<<1)-1; + for(k=0;k<l1;k++){ + ch[t1]=cc[t2]+cc[t3+t2]; + ch[t1+t0]=cc[t2]-cc[t3+t2]; + t2=(t1+=ido)<<1; + } + + if(ido<2)return; + if(ido==2)goto L105; + + t1=0; + t2=0; + for(k=0;k<l1;k++){ + t3=t1; + t5=(t4=t2)+(ido<<1); + t6=t0+t1; + for(i=2;i<ido;i+=2){ + t3+=2; + t4+=2; + t5-=2; + t6+=2; + ch[t3-1]=cc[t4-1]+cc[t5-1]; + tr2=cc[t4-1]-cc[t5-1]; + ch[t3]=cc[t4]-cc[t5]; + ti2=cc[t4]+cc[t5]; + ch[t6-1]=wa1[i-2]*tr2-wa1[i-1]*ti2; + ch[t6]=wa1[i-2]*ti2+wa1[i-1]*tr2; + } + t2=(t1+=ido)<<1; + } + + if(ido%2==1)return; + +L105: + t1=ido-1; + t2=ido-1; + for(k=0;k<l1;k++){ + ch[t1]=cc[t2]+cc[t2]; + ch[t1+t0]=-(cc[t2+1]+cc[t2+1]); + t1+=ido; + t2+=ido<<1; + } +} + +STIN void dradb3(int ido,int l1,double *cc,double *ch,double *wa1, + double *wa2){ + static double taur = -.5; + static double taui = .86602540378443864676372317075293618; + int i,k,t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10; + double ci2,ci3,di2,di3,cr2,cr3,dr2,dr3,ti2,tr2; + t0=l1*ido; + + t1=0; + t2=t0<<1; + t3=ido<<1; + t4=ido+(ido<<1); + t5=0; + for(k=0;k<l1;k++){ + tr2=cc[t3-1]+cc[t3-1]; + cr2=cc[t5]+(taur*tr2); + ch[t1]=cc[t5]+tr2; + ci3=taui*(cc[t3]+cc[t3]); + ch[t1+t0]=cr2-ci3; + ch[t1+t2]=cr2+ci3; + t1+=ido; + t3+=t4; + t5+=t4; + } + + if(ido==1)return; + + t1=0; + t3=ido<<1; + for(k=0;k<l1;k++){ + t7=t1+(t1<<1); + t6=(t5=t7+t3); + t8=t1; + t10=(t9=t1+t0)+t0; + + for(i=2;i<ido;i+=2){ + t5+=2; + t6-=2; + t7+=2; + t8+=2; + t9+=2; + t10+=2; + tr2=cc[t5-1]+cc[t6-1]; + cr2=cc[t7-1]+(taur*tr2); + ch[t8-1]=cc[t7-1]+tr2; + ti2=cc[t5]-cc[t6]; + ci2=cc[t7]+(taur*ti2); + ch[t8]=cc[t7]+ti2; + cr3=taui*(cc[t5-1]-cc[t6-1]); + ci3=taui*(cc[t5]+cc[t6]); + dr2=cr2-ci3; + dr3=cr2+ci3; + di2=ci2+cr3; + di3=ci2-cr3; + ch[t9-1]=wa1[i-2]*dr2-wa1[i-1]*di2; + ch[t9]=wa1[i-2]*di2+wa1[i-1]*dr2; + ch[t10-1]=wa2[i-2]*dr3-wa2[i-1]*di3; + ch[t10]=wa2[i-2]*di3+wa2[i-1]*dr3; + } + t1+=ido; + } +} + +STIN void dradb4(int ido,int l1,double *cc,double *ch,double *wa1, + double *wa2,double *wa3){ + static double sqrt2=1.4142135623730950488016887242097; + int i,k,t0,t1,t2,t3,t4,t5,t6,t7,t8; + double ci2,ci3,ci4,cr2,cr3,cr4,ti1,ti2,ti3,ti4,tr1,tr2,tr3,tr4; + t0=l1*ido; + + t1=0; + t2=ido<<2; + t3=0; + t6=ido<<1; + for(k=0;k<l1;k++){ + t4=t3+t6; + t5=t1; + tr3=cc[t4-1]+cc[t4-1]; + tr4=cc[t4]+cc[t4]; + tr1=cc[t3]-cc[(t4+=t6)-1]; + tr2=cc[t3]+cc[t4-1]; + ch[t5]=tr2+tr3; + ch[t5+=t0]=tr1-tr4; + ch[t5+=t0]=tr2-tr3; + ch[t5+=t0]=tr1+tr4; + t1+=ido; + t3+=t2; + } + + if(ido<2)return; + if(ido==2)goto L105; + + t1=0; + for(k=0;k<l1;k++){ + t5=(t4=(t3=(t2=t1<<2)+t6))+t6; + t7=t1; + for(i=2;i<ido;i+=2){ + t2+=2; + t3+=2; + t4-=2; + t5-=2; + t7+=2; + ti1=cc[t2]+cc[t5]; + ti2=cc[t2]-cc[t5]; + ti3=cc[t3]-cc[t4]; + tr4=cc[t3]+cc[t4]; + tr1=cc[t2-1]-cc[t5-1]; + tr2=cc[t2-1]+cc[t5-1]; + ti4=cc[t3-1]-cc[t4-1]; + tr3=cc[t3-1]+cc[t4-1]; + ch[t7-1]=tr2+tr3; + cr3=tr2-tr3; + ch[t7]=ti2+ti3; + ci3=ti2-ti3; + cr2=tr1-tr4; + cr4=tr1+tr4; + ci2=ti1+ti4; + ci4=ti1-ti4; + + ch[(t8=t7+t0)-1]=wa1[i-2]*cr2-wa1[i-1]*ci2; + ch[t8]=wa1[i-2]*ci2+wa1[i-1]*cr2; + ch[(t8+=t0)-1]=wa2[i-2]*cr3-wa2[i-1]*ci3; + ch[t8]=wa2[i-2]*ci3+wa2[i-1]*cr3; + ch[(t8+=t0)-1]=wa3[i-2]*cr4-wa3[i-1]*ci4; + ch[t8]=wa3[i-2]*ci4+wa3[i-1]*cr4; + } + t1+=ido; + } + + if(ido%2 == 1)return; + + L105: + + t1=ido; + t2=ido<<2; + t3=ido-1; + t4=ido+(ido<<1); + for(k=0;k<l1;k++){ + t5=t3; + ti1=cc[t1]+cc[t4]; + ti2=cc[t4]-cc[t1]; + tr1=cc[t1-1]-cc[t4-1]; + tr2=cc[t1-1]+cc[t4-1]; + ch[t5]=tr2+tr2; + ch[t5+=t0]=sqrt2*(tr1-ti1); + ch[t5+=t0]=ti2+ti2; + ch[t5+=t0]=-sqrt2*(tr1+ti1); + + t3+=ido; + t1+=t2; + t4+=t2; + } +} + +STIN void dradbg(int ido,int ip,int l1,int idl1,double *cc,double *c1, + double *c2,double *ch,double *ch2,double *wa){ + static double tpi=6.28318530717958647692528676655900577; + int idij,ipph,i,j,k,l,ik,is,t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10, + t11,t12; + double dc2,ai1,ai2,ar1,ar2,ds2; + int nbd; + double dcp,arg,dsp,ar1h,ar2h; + int ipp2; + + t10=ip*ido; + t0=l1*ido; + arg=tpi/(double)ip; + dcp=cos(arg); + dsp=sin(arg); + nbd=(ido-1)>>1; + ipp2=ip; + ipph=(ip+1)>>1; + if(ido<l1)goto L103; + + t1=0; + t2=0; + for(k=0;k<l1;k++){ + t3=t1; + t4=t2; + for(i=0;i<ido;i++){ + ch[t3]=cc[t4]; + t3++; + t4++; + } + t1+=ido; + t2+=t10; + } + goto L106; + + L103: + t1=0; + for(i=0;i<ido;i++){ + t2=t1; + t3=t1; + for(k=0;k<l1;k++){ + ch[t2]=cc[t3]; + t2+=ido; + t3+=t10; + } + t1++; + } + + L106: + t1=0; + t2=ipp2*t0; + t7=(t5=ido<<1); + for(j=1;j<ipph;j++){ + t1+=t0; + t2-=t0; + t3=t1; + t4=t2; + t6=t5; + for(k=0;k<l1;k++){ + ch[t3]=cc[t6-1]+cc[t6-1]; + ch[t4]=cc[t6]+cc[t6]; + t3+=ido; + t4+=ido; + t6+=t10; + } + t5+=t7; + } + + if (ido == 1)goto L116; + if(nbd<l1)goto L112; + + t1=0; + t2=ipp2*t0; + t7=0; + for(j=1;j<ipph;j++){ + t1+=t0; + t2-=t0; + t3=t1; + t4=t2; + + t7+=(ido<<1); + t8=t7; + for(k=0;k<l1;k++){ + t5=t3; + t6=t4; + t9=t8; + t11=t8; + for(i=2;i<ido;i+=2){ + t5+=2; + t6+=2; + t9+=2; + t11-=2; + ch[t5-1]=cc[t9-1]+cc[t11-1]; + ch[t6-1]=cc[t9-1]-cc[t11-1]; + ch[t5]=cc[t9]-cc[t11]; + ch[t6]=cc[t9]+cc[t11]; + } + t3+=ido; + t4+=ido; + t8+=t10; + } + } + goto L116; + + L112: + t1=0; + t2=ipp2*t0; + t7=0; + for(j=1;j<ipph;j++){ + t1+=t0; + t2-=t0; + t3=t1; + t4=t2; + t7+=(ido<<1); + t8=t7; + t9=t7; + for(i=2;i<ido;i+=2){ + t3+=2; + t4+=2; + t8+=2; + t9-=2; + t5=t3; + t6=t4; + t11=t8; + t12=t9; + for(k=0;k<l1;k++){ + ch[t5-1]=cc[t11-1]+cc[t12-1]; + ch[t6-1]=cc[t11-1]-cc[t12-1]; + ch[t5]=cc[t11]-cc[t12]; + ch[t6]=cc[t11]+cc[t12]; + t5+=ido; + t6+=ido; + t11+=t10; + t12+=t10; + } + } + } + +L116: + ar1=1.; + ai1=0.; + t1=0; + t9=(t2=ipp2*idl1); + t3=(ip-1)*idl1; + for(l=1;l<ipph;l++){ + t1+=idl1; + t2-=idl1; + + ar1h=dcp*ar1-dsp*ai1; + ai1=dcp*ai1+dsp*ar1; + ar1=ar1h; + t4=t1; + t5=t2; + t6=0; + t7=idl1; + t8=t3; + for(ik=0;ik<idl1;ik++){ + c2[t4++]=ch2[t6++]+ar1*ch2[t7++]; + c2[t5++]=ai1*ch2[t8++]; + } + dc2=ar1; + ds2=ai1; + ar2=ar1; + ai2=ai1; + + t6=idl1; + t7=t9-idl1; + for(j=2;j<ipph;j++){ + t6+=idl1; + t7-=idl1; + ar2h=dc2*ar2-ds2*ai2; + ai2=dc2*ai2+ds2*ar2; + ar2=ar2h; + t4=t1; + t5=t2; + t11=t6; + t12=t7; + for(ik=0;ik<idl1;ik++){ + c2[t4++]+=ar2*ch2[t11++]; + c2[t5++]+=ai2*ch2[t12++]; + } + } + } + + t1=0; + for(j=1;j<ipph;j++){ + t1+=idl1; + t2=t1; + for(ik=0;ik<idl1;ik++)ch2[ik]+=ch2[t2++]; + } + + t1=0; + t2=ipp2*t0; + for(j=1;j<ipph;j++){ + t1+=t0; + t2-=t0; + t3=t1; + t4=t2; + for(k=0;k<l1;k++){ + ch[t3]=c1[t3]-c1[t4]; + ch[t4]=c1[t3]+c1[t4]; + t3+=ido; + t4+=ido; + } + } + + if(ido==1)goto L132; + if(nbd<l1)goto L128; + + t1=0; + t2=ipp2*t0; + for(j=1;j<ipph;j++){ + t1+=t0; + t2-=t0; + t3=t1; + t4=t2; + for(k=0;k<l1;k++){ + t5=t3; + t6=t4; + for(i=2;i<ido;i+=2){ + t5+=2; + t6+=2; + ch[t5-1]=c1[t5-1]-c1[t6]; + ch[t6-1]=c1[t5-1]+c1[t6]; + ch[t5]=c1[t5]+c1[t6-1]; + ch[t6]=c1[t5]-c1[t6-1]; + } + t3+=ido; + t4+=ido; + } + } + goto L132; + + L128: + t1=0; + t2=ipp2*t0; + for(j=1;j<ipph;j++){ + t1+=t0; + t2-=t0; + t3=t1; + t4=t2; + for(i=2;i<ido;i+=2){ + t3+=2; + t4+=2; + t5=t3; + t6=t4; + for(k=0;k<l1;k++){ + ch[t5-1]=c1[t5-1]-c1[t6]; + ch[t6-1]=c1[t5-1]+c1[t6]; + ch[t5]=c1[t5]+c1[t6-1]; + ch[t6]=c1[t5]-c1[t6-1]; + t5+=ido; + t6+=ido; + } + } + } + +L132: + if(ido==1)return; + + for(ik=0;ik<idl1;ik++)c2[ik]=ch2[ik]; + + t1=0; + for(j=1;j<ip;j++){ + t2=(t1+=t0); + for(k=0;k<l1;k++){ + c1[t2]=ch[t2]; + t2+=ido; + } + } + + if(nbd>l1)goto L139; + + is= -ido-1; + t1=0; + for(j=1;j<ip;j++){ + is+=ido; + t1+=t0; + idij=is; + t2=t1; + for(i=2;i<ido;i+=2){ + t2+=2; + idij+=2; + t3=t2; + for(k=0;k<l1;k++){ + c1[t3-1]=wa[idij-1]*ch[t3-1]-wa[idij]*ch[t3]; + c1[t3]=wa[idij-1]*ch[t3]+wa[idij]*ch[t3-1]; + t3+=ido; + } + } + } + return; + + L139: + is= -ido-1; + t1=0; + for(j=1;j<ip;j++){ + is+=ido; + t1+=t0; + t2=t1; + for(k=0;k<l1;k++){ + idij=is; + t3=t2; + for(i=2;i<ido;i+=2){ + idij+=2; + t3+=2; + c1[t3-1]=wa[idij-1]*ch[t3-1]-wa[idij]*ch[t3]; + c1[t3]=wa[idij-1]*ch[t3]+wa[idij]*ch[t3-1]; + } + t2+=ido; + } + } +} + +STIN void drftb1(int n, double *c, double *ch, double *wa, int *ifac){ + int i,k1,l1,l2; + int na; + int nf,ip,iw,ix2,ix3,ido,idl1; + + nf=ifac[1]; + na=0; + l1=1; + iw=1; + + for(k1=0;k1<nf;k1++){ + ip=ifac[k1 + 2]; + l2=ip*l1; + ido=n/l2; + idl1=ido*l1; + if(ip!=4)goto L103; + ix2=iw+ido; + ix3=ix2+ido; + + if(na!=0) + dradb4(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1); + else + dradb4(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1); + na=1-na; + goto L115; + + L103: + if(ip!=2)goto L106; + + if(na!=0) + dradb2(ido,l1,ch,c,wa+iw-1); + else + dradb2(ido,l1,c,ch,wa+iw-1); + na=1-na; + goto L115; + + L106: + if(ip!=3)goto L109; + + ix2=iw+ido; + if(na!=0) + dradb3(ido,l1,ch,c,wa+iw-1,wa+ix2-1); + else + dradb3(ido,l1,c,ch,wa+iw-1,wa+ix2-1); + na=1-na; + goto L115; + + L109: +/* The radix five case can be translated later..... */ +/* if(ip!=5)goto L112; + + ix2=iw+ido; + ix3=ix2+ido; + ix4=ix3+ido; + if(na!=0) + dradb5(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1,wa+ix4-1); + else + dradb5(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1,wa+ix4-1); + na=1-na; + goto L115; + + L112:*/ + if(na!=0) + dradbg(ido,ip,l1,idl1,ch,ch,ch,c,c,wa+iw-1); + else + dradbg(ido,ip,l1,idl1,c,c,c,ch,ch,wa+iw-1); + if(ido==1)na=1-na; + + L115: + l1=l2; + iw+=(ip-1)*ido; + } + + if(na==0)return; + + for(i=0;i<n;i++)c[i]=ch[i]; +} + +void __ogg_fdrfftb(int n, double *r, double *wsave, int *ifac){ + if (n == 1)return; + drftb1(n, r, wsave, wsave+n, ifac); +} + +STIN void dcsqb1(int n,double *x,double *w,double *xh,int *ifac){ + int modn,i,k,kc; + int np2,ns2; + double xim1; + + ns2=(n+1)>>1; + np2=n; + + for(i=2;i<n;i+=2){ + xim1=x[i-1]+x[i]; + x[i]-=x[i-1]; + x[i-1]=xim1; + } + + x[0]+=x[0]; + modn=n%2; + if(modn==0)x[n-1]+=x[n-1]; + + __ogg_fdrfftb(n,x,xh,ifac); + + kc=np2; + for(k=1;k<ns2;k++){ + kc--; + xh[k]=w[k-1]*x[kc]+w[kc-1]*x[k]; + xh[kc]=w[k-1]*x[k]-w[kc-1]*x[kc]; + } + + if(modn==0)x[ns2]=w[ns2-1]*(x[ns2]+x[ns2]); + + kc=np2; + for(k=1;k<ns2;k++){ + kc--; + x[k]=xh[k]+xh[kc]; + x[kc]=xh[k]-xh[kc]; + } + x[0]+=x[0]; +} + +void __ogg_fdcosqb(int n,double *x,double *wsave,int *ifac){ + static double tsqrt2 = 2.8284271247461900976033774484194; + double x1; + + if(n<2){ + x[0]*=4; + return; + } + if(n==2){ + x1=(x[0]+x[1])*4; + x[1]=tsqrt2*(x[0]-x[1]); + x[0]=x1; + return; + } + + dcsqb1(n,x,wsave,wsave+n,ifac); +} diff --git a/nist/matrix.c b/nist/matrix.c new file mode 100644 index 0000000..50bf007 --- /dev/null +++ b/nist/matrix.c @@ -0,0 +1,168 @@ + + +#include <stdio.h> +#include <stdlib.h> +/*#include "defs.h" +#include "proto.h"*/ +#define MIN(x,y) ((x) > (y) ? (y) : (x)) + +#include "matrix.h" + +/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * + R A N K A L G O R I T H M R O U T I N E S + * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ + +int computeRank(int M, int Q, BitField** matrix) +{ + int i; + int rank; + int m = MIN(M,Q); + + /* FORWARD APPLICATION OF ELEMENTARY ROW OPERATIONS */ + + for(i = 0; i < m-1; i++) { + if (matrix[i][i].b == 1) + perform_elementary_row_operations(0, i, M, Q, matrix); + else { /* matrix[i][i] = 0 */ + if (find_unit_element_and_swap(0, i, M, Q, matrix) == 1) + perform_elementary_row_operations(0, i, M, Q, matrix); + } + } + + /* BACKWARD APPLICATION OF ELEMENTARY ROW OPERATIONS */ + for(i = m-1; i > 0; i--) { + if (matrix[i][i].b == 1) + perform_elementary_row_operations(1, i, M, Q, matrix); + else { /* matrix[i][i] = 0 */ + if (find_unit_element_and_swap(1, i, M, Q, matrix) == 1) + perform_elementary_row_operations(1, i, M, Q, matrix); + } + } + rank = determine_rank(m,M,Q,matrix); + + return rank; +} + +void perform_elementary_row_operations(int flag, int i, int M, int Q, + BitField** A) +{ + int j,k; + + switch(flag) + { + case 0: for(j = i+1; j < M; j++) + if (A[j][i].b == 1) + for(k = i; k < Q; k++) + A[j][k].b = (A[j][k].b + A[i][k].b) % 2; + break; + + case 1: for(j = i-1; j >= 0; j--) + if (A[j][i].b == 1) + for(k = 0; k < Q; k++) + A[j][k].b = (A[j][k].b + A[i][k].b) % 2; + break; + default: fprintf(stderr,"ERROR IN CALL TO perform_elementary_row_operations\n"); + break; + } + return; +} + +int find_unit_element_and_swap(int flag, int i, int M, int Q, BitField** A) +{ + int index; + int row_op = 0; + + switch(flag) + { + case 0: index = i+1; + while ((index < M) && (A[index][i].b == 0)) + index++; + if (index < M) + row_op = swap_rows(i,index,Q,A); + break; + case 1: + index = i-1; + while ((index >= 0) && (A[index][i].b == 0)) + index--; + if (index >= 0) + row_op = swap_rows(i,index,Q,A); + break; + default: fprintf(stderr,"ERROW IN CALL TO find_unit_element_and_swap\n"); + break; + } + return row_op; +} + +int swap_rows(int i, int index, int Q, BitField** A) +{ + int p; + BitField temp; + + for(p = 0; p < Q; p++) { + temp.b = A[i][p].b; + A[i][p].b = A[index][p].b; + A[index][p].b = temp.b; + } + return 1; +} + +int determine_rank(int m, int M, int Q, BitField** A) +{ + int i, j, rank, allZeroes; + + /* DETERMINE RANK, THAT IS, COUNT THE NUMBER OF NONZERO ROWS */ + + rank = m; + for(i = 0; i < M; i++) { + allZeroes = 1; + for(j=0; j < Q; j++) { + if (A[i][j].b == 1) { + allZeroes = 0; + break; + } + } + if (allZeroes == 1) rank--; + } + return rank; +} + + +void display_matrix(int M, int Q, BitField** m) +{ + int i,j; + + for (i = 0; i < M; i++) { + for (j = 0; j < Q; j++) + fprintf(stderr,"%d ", m[i][j].b); + fprintf(stderr,"\n"); + } + return; +} + +void def_matrix(int M, int Q, BitField** m,int k, int *pt, int *PT, int*DATA, int *ARRAY) +{ + int i,j; + for (i = 0; i < M; i++) + for (j = 0; j < Q; j++) { + m[i][j].b = *DATA & 1; +(*PT)++; + if ((*PT) == 32) + { + *PT = 0; + (*pt)++; + *DATA = ARRAY[*pt]; + } + else + *DATA = (*DATA) >> 1; + + } + return; +} + +void delete_matrix(int M, BitField** matrix) +{ + int i; + for (i = 0; i < M; i++) + free(matrix[i]); + free(matrix); +} diff --git a/nist/matrix.h b/nist/matrix.h new file mode 100644 index 0000000..d66e9a2 --- /dev/null +++ b/nist/matrix.h @@ -0,0 +1,19 @@ +/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * + R A N K A L G O R I T H M F U N C T I O N P R O T O T Y P E S + * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ +struct bit { + unsigned char b:1; +}; +typedef struct bit BitField; + + +void perform_elementary_row_operations(int,int,int,int,BitField**); +int find_unit_element_and_swap(int,int,int,int, BitField**); +int swap_rows(int, int, int, BitField**); +int determine_rank(int, int, int, BitField**); +int computeRank(int,int,BitField**); +void define_matrix(int,int,int,BitField**); +BitField** create_matrix(int, int); +void display_matrix(int, int, BitField**); +void def_matrix(int M, int Q, BitField** m,int k, int *pt, int *PT, int*DATA, int *ARRAY); +void delete_matrix(int, BitField**); diff --git a/nist/mconf.h b/nist/mconf.h new file mode 100644 index 0000000..84c15d7 --- /dev/null +++ b/nist/mconf.h @@ -0,0 +1,169 @@ +/* mconf.h + * + * Common include file for math routines + * + * + * + * SYNOPSIS: + * + * #include "mconf.h" + * + * + * + * DESCRIPTION: + * + * This file contains definitions for error codes that are + * passed to the common error handling routine mtherr() + * (which see). + * + * The file also includes a conditional assembly definition + * for the type of computer arithmetic (IEEE, DEC, Motorola + * IEEE, or UNKnown). + * + * For Digital Equipment PDP-11 and VAX computers, certain + * IBM systems, and others that use numbers with a 56-bit + * significand, the symbol DEC should be defined. In this + * mode, most floating point constants are given as arrays + * of octal integers to eliminate decimal to binary conversion + * errors that might be introduced by the compiler. + * + * For little-endian computers, such as IBM PC, that follow the + * IEEE Standard for Binary Floating Point Arithmetic (ANSI/IEEE + * Std 754-1985), the symbol IBMPC should be defined. These + * numbers have 53-bit significands. In this mode, constants + * are provided as arrays of hexadecimal 16 bit integers. + * + * Big-endian IEEE format is denoted MIEEE. On some RISC + * systems such as Sun SPARC, double precision constants + * must be stored on 8-byte address boundaries. Since integer + * arrays may be aligned differently, the MIEEE configuration + * may fail on such machines. + * + * To accommodate other types of computer arithmetic, all + * constants are also provided in a normal decimal radix + * which one can hope are correctly converted to a suitable + * format by the available C language compiler. To invoke + * this mode, define the symbol UNK. + * + * An important difference among these modes is a predefined + * set of machine arithmetic constants for each. The numbers + * MACHEP (the machine roundoff error), MAXNUM (largest number + * represented), and several other parameters are preset by + * the configuration symbol. Check the file const.c to + * ensure that these values are correct for your computer. + * + * Configurations NANS, INFINITIES, MINUSZERO, and DENORMAL + * may fail on many systems. Verify that they are supposed + * to work on your computer. + */ + +/* +Cephes Math Library Release 2.3: June, 1995 +Copyright 1984, 1987, 1989, 1995 by Stephen L. Moshier +*/ + + +/* Constant definitions for math error conditions + */ +#ifndef DOMAIN +#define DOMAIN 1 /* argument domain error */ +#define SING 2 /* argument singularity */ +#define OVERFLOW 3 /* overflow range error */ +#define UNDERFLOW 4 /* underflow range error */ +#define TLOSS 5 /* total loss of precision */ +#define PLOSS 6 /* partial loss of precision */ + +#define EDOM 33 +#define ERANGE 34 +#endif + +/* Complex numeral. */ +typedef struct + { + double r; + double i; + } cmplx; + +/* Long double complex numeral. */ +/* +typedef struct + { + long double r; + long double i; + } cmplxl; +*/ + +/* Type of computer arithmetic */ + +/* PDP-11, Pro350, VAX: + */ +/* #define DEC 1 */ + +/* Intel IEEE, low order words come first: + */ +/*#define IBMPC 1 */ + +/* Motorola IEEE, high order words come first + * (Sun 680x0 workstation): + */ +/* #define MIEEE 1 */ + +/* UNKnown arithmetic, invokes coefficients given in + * normal decimal format. Beware of range boundary + * problems (MACHEP, MAXLOG, etc. in const.c) and + * roundoff problems in pow.c: + * (Sun SPARCstation) + */ +#define UNK 1 + +/* If you define UNK, then be sure to set BIGENDIAN properly. */ +#define BIGENDIAN 0 + +/* Define this `volatile' if your compiler thinks + * that floating point arithmetic obeys the associative + * and distributive laws. It will defeat some optimizations + * (but probably not enough of them). + * + * #define VOLATILE volatile + */ +#define VOLATILE + +/* For 12-byte long doubles on an i386, pad a 16-bit short 0 + * to the end of real constants initialized by integer arrays. + * + * #define XPD 0, + * + * Otherwise, the type is 10 bytes long and XPD should be + * defined blank (e.g., Microsoft C). + * + * #define XPD + */ +#define XPD 0, + +/* Define to support tiny denormal numbers, else undefine. */ +#define DENORMAL 1 + +/* Define to ask for infinity support, else undefine. */ +/* #define INFINITIES 1 */ + +/* Define to ask for support of numbers that are Not-a-Number, + else undefine. This may automatically define INFINITIES in some files. */ +/* #define NANS 1 */ + +/* Define to distinguish between -0.0 and +0.0. */ +/* #define MINUSZERO 1 */ + +/* Define 1 for ANSI C atan2() function + See atan.c and clog.c. */ +#define ANSIC 1 + +/* Get ANSI function prototypes, if you want them. */ +#ifdef __STDC__ +#define ANSIPROT +#include "cephes-protos.h" +#else +int mtherr(); +#endif + +/* Variable for error reporting. See mtherr.c. */ +extern int merror; diff --git a/nist/nist.c b/nist/nist.c new file mode 100644 index 0000000..e68c27b --- /dev/null +++ b/nist/nist.c @@ -0,0 +1,57 @@ +#include <stdio.h> +#include <stdlib.h> +#include <unistd.h> +#include <string.h> +#include "nist.h" + +#define _32MB (32 * 1024 * 1024) +#define _08MB (8 * 1024 * 1024) +#define _04MB (4 * 1024 * 1024) +#define _02MB (2 * 1024 * 1024) + + +static int random_pool1 [_32MB]; +char basedirname[FILENAME_MAX+1]; + +char *GetBaseDir(void) +{ + return basedirname; +} + +int main(int argc, char **argv) +{ +char *filename = ""; +FILE *fp = stdin; +long result=0; + + if (argc<2 || argc>3) { + printf("Usage sts <file> [<template directory>]\n"); + return 1; + } + /** + * get optional directory name + */ + basedirname[0] = 0; + if (argc>2) { + strcat(basedirname, argv[2]); + strcat(basedirname, "/"); + } + +filename = argv[1]; +if ((fp = fopen(filename, "rb")) == NULL) { + printf("Cannot open file %s\n", filename); + return 2; + } +result = fread(random_pool1,sizeof(int),_04MB,fp); +fclose(fp); +if (result!=_04MB) { + printf("16MB sample required %ld != %d\n",result, _04MB); + return 3; + } +if (PackTestF (random_pool1, _04MB, "nist.out") < 8) { + if (PackTestL (random_pool1, _04MB, "nist.out") < 8) + return 0; + return 5; + } +return 4; +} diff --git a/nist/nist.h b/nist/nist.h new file mode 100644 index 0000000..b5962cd --- /dev/null +++ b/nist/nist.h @@ -0,0 +1,33 @@ +/** + * Supply function prototypes for the test suite. + */ +int PackTestF (int *ARRAY, int ArraySize, char *C); +int PackTestL (int *ARRAY, int ArraySize, char *C); + +int ApproximateEntropy (int mmin, int mmax, int n, int *ARRAY); +int BlockFrequency (int ArraySize, int m, int *ARRAY); +int CumulativeSums (int n, int *ARRAY); +int DiscreteFourierTransform (int N, int *ARRAY); +int Frequency (int n, int *ARRAY); +int LempelZivCompression (int n, int *ARRAY, int *DATA, int *pt, int *PT); +int LinearComplexity (int M, int N, int *ARRAY, int PT); +int LongestRunOfOnes (int n, int *ARRAY); +int NonOverlappingTemplateMatchings (int m, int n, int *ARRAY); +int OverlappingTemplateMatchings (int m, int n, int *ARRAY); +int RandomExcursions (int n, int *ARRAY); +int RandomExcursionsVariant (int n, int *ARRAY); +int Rank (int n, int *ARRAY); +int Runs (int n, int *ARRAY); +int Serial (int m, int n, int *ARRAY, int PT); +int UNIVERSAL (int n, int *ARRAY); +int Universal (int n, int *ARRAY); + +double psi2 (int m, int n, int *ARRAY, int PT); +/** + * From dfft.c + */ +void __ogg_fdrffti(int n, double *wsave, int *ifac); +void __ogg_fdrfftf(int n,double *r,double *wsave,int *ifac); + +char *GetBaseDir(void); + diff --git a/nist/packtest.c b/nist/packtest.c new file mode 100644 index 0000000..879cd73 --- /dev/null +++ b/nist/packtest.c @@ -0,0 +1,2252 @@ +/* -------------------------------------------------------------------------- + Title : The NIST Statistical Test Suite + + Date : December 1999 + + Programmer : Juan Soto + + Summary : For use in the evaluation of the randomness of bitstreams + produced by cryptographic random number generators. + + Package : Version 1.0 + + Copyright : (c) 1999 by the National Institute Of Standards & Technology + + History : Version 1.0 by J. Soto, October 1999 + Revised by J. Soto, November 1999 + + Keywords : Pseudorandom Number Generator (PRNG), Randomness, Statistical + Tests, Complementary Error functions, Incomplete Gamma + Function, Random Walks, Rank, Fast Fourier Transform, + Template, Cryptographically Secure PRNG (CSPRNG), + Approximate Entropy (ApEn), Secure Hash Algorithm (SHA-1), + Blum-Blum-Shub (BBS) CSPRNG, Micali-Schnorr (MS) CSPRNG, + + Source : David Banks, Elaine Barker, James Dray, Allen Heckert, + Stefan Leigh, Mark Levenson, James Nechvatal, Andrew Rukhin, + Miles Smid, Juan Soto, Mark Vangel, and San Vo. + + Technical + Assistance : Lawrence Bassham, Ron Boisvert, James Filliben, Sharon Keller, + Daniel Lozier, and Bert Rust. + + Warning : Portability Issues. + + Limitation : Amount of memory allocated for workspace. + + Restrictions: Permission to use, copy, and modify this software without + fee is hereby granted, provided that this entire notice is + included in all copies of any software which is or includes + a copy or modification of this software and in all copies + of the supporting documentation for such software. + -------------------------------------------------------------------------- */ + +/* -------------------------------------------------------------------------- + * Derived by Andre Seznec (IRISA/INRIA - www.irisa.fr/caps) to overcome the + * limitations on the sequences size. It should now work on sequences larger + * than 1MB and smaller than 256MB. + * + * Modified on March 2002, last update on June 2002 + * + * Modified Aug 2008 for clean compilation under gcc 4.4 + * -------------------------------------------------------------------------- + */ +#ifndef NOIO +#define IO +#endif +#define LinuxorUnix +#ifdef WIN +#ifndef CYGWIN +#undef LinuxorUnix +/* same libraries are available*/ +#endif +#endif + +#ifdef LinuxorUnix + +#define MAX(x,y) ((x) < (y) ? (y) : (x)) +#define MOD(x,y) (((x) % (y)) < 0 ? ((x) % (y)) + (y) : ((x) % (y))) + + +#include <stdio.h> +#include <math.h> +#include <stdlib.h> +#include "special-functions.h" +#include "mconf.h" +#include "matrix.h" +#include "nist.h" + +#define THRESHOLDFAIL 0.001 +static int TotalTEST = 0; + +/* + * on Windows systems, we encountered difficulties with the memory allocation + * libraries, therefor we used our own allocation technique in a big array + */ + +#define BIGSIZE 32768*1024 +char BigArray[BIGSIZE]; +int PtBigArray = 0; +char * +MYCALLOC (int n, int S) +{ + int i, OldPt; + OldPt = PtBigArray; + PtBigArray += ((n * S + 16) & 0xfffffff0); + if (PtBigArray >= BIGSIZE) + { + printf ("Pb memory allocation in PackTest\n"); + exit (1); + } + for (i = OldPt; i < PtBigArray; i++) + BigArray[i] = 0; + + return (&BigArray[OldPt]); +} + +void +MYFREE (char *X) +{ + PtBigArray = (int) (X - BigArray); +} + +BitField **create_matrix (int M, int Q); + +FILE *fileout; + +int +PackTestF (int *ARRAY, int ArraySize, char *C) +{ + int i; + int TEST = 0; + int failure = 0; + +#ifdef IO + fileout = fopen (C, "w"); +#endif + if (ArraySize >= (1 << 26)) + { + fprintf (fileout, + "This test does not work for an array larger than 256Mbytes\n"); + return (failure); + } +#ifdef IO + { + printf + ("16 NIST tests to be executed: results will be written in file %s\n\n", + C); + printf + ("To maintain response time in a reasonable range,\nsome of the tests are executed on only subsequences\n\n"); + + TEST++; + fprintf (stderr, "Test %d, ", TEST); + fprintf (stderr, "\t Frequency test \n\n"); + fprintf (fileout, " \n\n Frequency test \n\n"); + } +#endif + for (i = 4; i < ArraySize; i = (i << 1)) + { + int inter; + inter = MOD (ARRAY[ArraySize - i], ArraySize - i); + if (failure >= 8) + { +#ifdef IO + fprintf (stderr, + "%d failures at %f threshold,\nsequence should be considered as not random\n", + failure, THRESHOLDFAIL); + +#endif + return (failure); + } + + + failure += Frequency (i, &ARRAY[inter]); + } + + if (failure >= 8) + { +#ifdef IO + fprintf (stderr, + "%d failures at %f threshold,\nsequence should be considered as not random\n", + failure, THRESHOLDFAIL); + +#endif + return (failure); + } + + failure += Frequency (ArraySize, ARRAY); +#ifdef IO + { + TEST++; + fprintf (stderr, "Test %d, ", TEST); + fprintf (stderr, "\t Block Frequency test \n\n"); + fprintf (fileout, " \n\n Block Frequency test \n\n"); + } +#endif + for (i = 4; i < ArraySize; i = (i << 1)) + { + int inter; + inter = MOD (ARRAY[ArraySize - i], ArraySize - i); + failure += BlockFrequency (i, (32 * i / 99), &ARRAY[inter]); + } + + if (failure >= 8) + { +#ifdef IO + fprintf (stderr, + "%d failures at %f threshold,\nsequence should be considered as not random\n", + failure, THRESHOLDFAIL); + +#endif + return (failure); + } + + failure += BlockFrequency (ArraySize, (32 * ArraySize / 99), ARRAY); +#ifdef IO + { + TEST++; + fprintf (stderr, "Test %d, ", TEST); + fprintf (stderr, "\t Runs test \n\n"); + fprintf (fileout, " \n\n Runs test \n\n"); + } +#endif + for (i = 4; i < ArraySize; i = (i << 1)) + { + int inter; + inter = MOD (ARRAY[ArraySize - i], ArraySize - i); + if (failure >= 8) + { +#ifdef IO + + + fprintf (stderr, + "%d failures at %f threshold,\nsequence should be considered as not random\n", + failure, THRESHOLDFAIL); + +#endif + return (failure); + } + + failure += Runs (i, &ARRAY[inter]); + } + Runs (ArraySize, ARRAY); +#ifdef IO + { + TEST++; + fprintf (stderr, "Test %d, ", TEST); + fprintf (stderr, "\t LongestRunOfOnes \n\n"); + fprintf (fileout, " \n\n LongestRunOfOnes \n\n"); + } +#endif + for (i = 0; i < 8; i++) + { + int index; + index = (ArraySize / 8) * i; + if (ArraySize > 262144) + index = index + (MOD (ARRAY[index], (ArraySize / 8) - 32768)); + failure += LongestRunOfOnes (((750000) >> 5) + 1, &ARRAY[index]); +if (failure >= 8) + { +#ifdef IO + + + fprintf (stderr, + "%d failures at %f threshold,\nsequence should be considered as not random\n", + failure, THRESHOLDFAIL); + +#endif + return (failure); + } + } + +fclose(fileout); +return(failure); +} +int +PackTestL (int *ARRAY, int ArraySize, char *C) +{ + int i; + int TEST = 4; + int DATA, pt, PT; + int failure = 0; + +#ifdef IO + fileout = fopen (C, "a"); +#endif +#ifdef IO + { + TEST++; + fprintf (stderr, "Test %d, ", TEST); + fprintf (stderr, + "\t Binary Matrix Rank test on 8 random 38,912 bit slices \n\n"); + fprintf (fileout, + " \n\n Binary Matrix Rank test on 8 random 38,912 bit slices \n\n"); + } +#endif + if (ArraySize >= (1 << 26)) + { + fprintf (fileout, + "This test does not work for an array larger than 256Mbytes\n"); + return (failure); + } + + + for (i = 0; i < 8; i++) + { + int index; + index = (ArraySize / 8) * i; + if (ArraySize > 262144) + index = index + (MOD (ARRAY[index], (ArraySize / 8) - 32768)); + failure += Rank (38192 >> 5, &ARRAY[index]); + if (failure >= 8) + { +#ifdef IO + + + fprintf (stderr, + "%d failures at %f threshold,\nsequence should be considered as not random\n", + failure, THRESHOLDFAIL); + +#endif + return (failure); + } + + + } +#ifdef IO + { + TEST++; + fprintf (stderr, "Test %d, ", TEST); + fprintf (stderr, "\t Discrete Fourier test by slice of 1 Mbits\n\n"); + fprintf (fileout, " \n\n Discrete Fourier test by slice of 1 Mbits\n\n"); + fprintf (fileout, " \n\n 8 random slices are picked \n\n"); + } +#endif + for (i = 0; i < 8; i++) + { + int index; + index = (ArraySize / 8) * i; + if (ArraySize > 262144) + index = index + (MOD (ARRAY[index], (ArraySize / 8) - 32768)); + if (failure >= 8) + { +#ifdef IO + + + fprintf (stderr, + "%d failures at %f threshold,\nsequence should be considered as not random\n", + failure, THRESHOLDFAIL); + +#endif + return (failure); + } + + failure += DiscreteFourierTransform (32768, &ARRAY[index]); + } +#ifdef IO + { + TEST++; + fprintf (stderr, "Test %d, ", TEST); + fprintf (stderr, "\t NON OVERLAPPING TEMPLATE MATCHING TEST \n"); + fprintf (fileout, "\n\n\t NON OVERLAPPING TEMPLATE MATCHING TEST \n"); + fprintf (fileout, "\t 1 random 1Mbit slices \n\n"); + } +#endif + for (i = 0; i < 1; i++) + { + int index; + index = (ArraySize) * i; + if (ArraySize > 262144) + index = index + MOD (ARRAY[index], ((ArraySize) - 32768)); + if (failure >= 8) + { +#ifdef IO + + + fprintf (stderr, + "%d failures at %f threshold,\nsequence should be considered as not random\n", + failure, THRESHOLDFAIL); + +#endif + return (failure); + } + + failure += + NonOverlappingTemplateMatchings (9, 1024 * 1024, &ARRAY[index]); + } +#ifdef IO + { + TEST++; + fprintf (stderr, "Test %d, ", TEST); + fprintf (stderr, "\t OVERLAPPING TEMPLATE MATCHING TEST \n"); + fprintf (fileout, "\n\n\t OVERLAPPING TEMPLATE MATCHING TEST \n"); + fprintf (fileout, "\t 8 random 1000000 bits slices \n\n"); + } +#endif + for (i = 0; i < 8; i++) + { + int index; + index = (ArraySize / 8) * i; + if (ArraySize > 262144) + index = index + (MOD (ARRAY[index], (ArraySize / 8) - 32768)); + if (failure >= 8) + { +#ifdef IO + + + fprintf (stderr, + "%d failures at %f threshold,\nsequence should be considered as not random\n", + failure, THRESHOLDFAIL); + +#endif + return (failure); + } + + + failure += + OverlappingTemplateMatchings (9, 1000000 >> 5, &ARRAY[index]); + } +#ifdef IO + { + TEST++; + fprintf (stderr, "Test %d, ", TEST); + fprintf (stderr, "\t Maurer's Universal test\n"); + fprintf (fileout, "\n\n Maurer's Universal test\n"); + fprintf + (fileout, + "\n For each of the L parameters, we test the beginning of the sequence\n\n"); + } +#endif + + if (failure >= 8) + { +#ifdef IO + fprintf (stderr, + "%d failures at %f threshold,\nsequence should be considered as not random\n", + failure, THRESHOLDFAIL); + +#endif + return (failure); + } + + + failure += Universal (ArraySize, ARRAY); + + TEST++; +#ifdef IO + { + fprintf (stderr, "Test %d, ", TEST); + fprintf (stderr, + "\t LEMPEL-ZIV COMPRESSION TEST: 8 random slices of 1000000 bits\n\n"); + fprintf (fileout, + "\n\n\t LEMPEL-ZIV COMPRESSION TEST: 8 random slices of 1000000 bits\n\n"); + } +#endif + for (i = 0; i < 8; i++) + { + int index; + index = (ArraySize / 8) * i; + if (ArraySize > 262144) + index = index + (MOD (ARRAY[index], (ArraySize / 8) - 32768)); + DATA = ARRAY[index]; + pt = 0; + PT = 0; + if (failure >= 8) + { +#ifdef IO + + + fprintf (stderr, + "%d failures at %f threshold,\nsequence should be considered as not random\n", + failure, THRESHOLDFAIL); + +#endif + return (failure); + } + + failure += + LempelZivCompression (1000000, &ARRAY[index], &DATA, &pt, &PT); + } + if (failure >= 8) + { +/* in order to enable fast detection of strong failure for random sequences*/ +#ifdef IO + fprintf (fileout, + "%d failed individual tests with THRESHOLD %f on %d individual tests before LINEAR COMPLEXITY TEST\ndon't waste your CPU time\n", + failure, THRESHOLDFAIL, TotalTEST); + fprintf (stderr, + "%d failed individual tests with THRESHOLD %f on %d individual tests before LINEAR COMPLEXITY TEST\ndon't waste your CPU time\n", + failure, THRESHOLDFAIL, TotalTEST); + +#endif + fclose (fileout); + return (failure); + } +#ifdef IO + { + TEST++; + fprintf (stderr, "Test %d, ", TEST); + fprintf (stderr, + "\t LINEAR COMPLEXITY TEST: 1 SLICE OF 4731 * 371 bits\n\n"); + fprintf (fileout, + "\n\n\t LINEAR COMPLEXITY TEST: 1 SLICE OF 4731 * 371 bits\n\n"); + } +#endif + +/* these parameters were chosen arbitraly, NIST recommendation are 500<=M<=5000, N>=200, NM >=1000000*/ + for (i = 5; i <= 5; i++) + { + int index, inter; + inter = (371 * 4731 + 32) >> 5; + index = (ArraySize / 8) * i; + if (ArraySize > 262144) + index = index + (MOD (ARRAY[index], (ArraySize / 8) - 32768)); + if (index + inter > ArraySize) + index = ArraySize - inter; + + + failure += LinearComplexity (4731, 371, &ARRAY[index], 0); + } +#ifdef IO + { + TEST++; + fprintf (stderr, "Test %d, ", TEST); + fprintf (stderr, "\t SERIAL TEST BY SLICE OF 1000000 bits\n\n"); + fprintf (fileout, "\n\n\t SERIAL TEST BY SLICE OF 1000000 bits\n\n"); + } +#endif + for (i = 0; i < 8; i++) + { + int index; + index = (ArraySize / 8) * i; + if (ArraySize > 262144) + index = index + (MOD (ARRAY[index], (ArraySize / 8) - 32768)); + + if (failure >= 8) + { +#ifdef IO + + + fprintf (stderr, + "%d failures at %f threshold,\nsequence should be considered as not random\n", + failure, THRESHOLDFAIL); + +#endif + return (failure); + } + + failure += Serial (14, 1000000, &ARRAY[index], 0); + } +#ifdef IO + { + TEST++; + fprintf (stderr, "Test %d, ", TEST); + fprintf (stderr, "\t APPROXIMATE ENTROPY TEST\n\n"); + fprintf (fileout, "\n\n\t APPROXIMATE ENTROPY TEST\n\n"); + } +#endif + failure += ApproximateEntropy (3, 17, 32 * ArraySize, ARRAY); +#ifdef IO + { + TEST++; + fprintf (stderr, "Test %d, ", TEST); + fprintf (stderr, "\t CUSUM TEST\n\n"); + fprintf (fileout, "\n\n\t CUSUM TEST\n\n"); + } +#endif + for (i = 4; i < ArraySize; i = (i << 1)) + { + int inter; + fprintf (stderr, "\t\t ArraySize: %d\n", i); + inter = MOD (ARRAY[ArraySize - i], ArraySize - i); + if (failure >= 8) + { +#ifdef IO + + + fprintf (stderr, + "%d failures at %f threshold,\nsequence should be considered as not random\n", + failure, THRESHOLDFAIL); + +#endif + return (failure); + } + + failure += CumulativeSums (32 * i, &ARRAY[inter]); + } + + if (failure >= 8) + { +#ifdef IO + fprintf (stderr, + "%d failures at %f threshold,\nsequence should be considered as not random\n", + failure, THRESHOLDFAIL); + +#endif + return (failure); + } + + failure += CumulativeSums (32 * ArraySize, ARRAY); + +#ifdef IO + { + TEST++; + fprintf (stderr, "Test %d, ", TEST); + fprintf (stderr, + "\t RANDOM EXCURSION TEST: 8 SLICE OF 1000000 bits\n\n"); + fprintf (fileout, + "\n\n\t RANDOM EXCURSION TEST: 8 SLICE OF 1000000 bits\n\n"); + } +#endif + for (i = 0; i < 8; i++) + { + int index; + fprintf (stderr, "\t\t Slice number: %d\n", i); + index = (ArraySize / 8) * i; + if (ArraySize > 262144) + index = index + (MOD (ARRAY[index], (ArraySize / 8) - 32768)); + + if (failure >= 8) + { +#ifdef IO + + + fprintf (stderr, + "%d failures at %f threshold,\nsequence should be considered as not random\n", + failure, THRESHOLDFAIL); + +#endif + return (failure); + } + + failure += RandomExcursions (1000000, &ARRAY[index >> 5]); + } +#ifdef IO + { + TEST++; + fprintf (stderr, "Test %d, ", TEST); + fprintf (stderr, + "\t RANDOM EXCURSION TEST VARIANT: 8 SLICES OF 1000000 bits\n\n"); + fprintf (fileout, + "\n\n\t RANDOM EXCURSION TEST VARIANT: 8 SLICES OF 1000000 bits\n\n"); + } +#endif + for (i = 0; i < 8; i++) + { + int index; + index = (ArraySize / 8) * i; + if (ArraySize > 262144) + index = index + (MOD (ARRAY[index], (ArraySize / 8) - 32768)); + + if (failure >= 8) + { +#ifdef IO + + + fprintf (stderr, + "%d failures at %f threshold,\nsequence should be considered as not random\n", + failure, THRESHOLDFAIL); + +#endif + return (failure); + } + + failure += RandomExcursionsVariant (1000000, &ARRAY[index >> 5]); + } + + +#ifdef IO + { + fprintf (stderr, + "%d failed individual tests with THRESHOLD %f on %d individual tests\n", + failure, THRESHOLDFAIL, TotalTEST); + fprintf (fileout, + "%d failed individual tests with THRESHOLD %f on %d individual tests\n", + failure, THRESHOLDFAIL, TotalTEST); + } +#endif + fclose (fileout); + return (failure); +} + +#define MAXNUMOFTEMPLATES 148 +int TEMPLATE[MAXNUMOFTEMPLATES]; +int Skip[MAXNUMOFTEMPLATES]; +int nu[6][MAXNUMOFTEMPLATES]; +int Wj[8][MAXNUMOFTEMPLATES]; +int W_obs[MAXNUMOFTEMPLATES]; +/* was adapted to m =9 by A. Seznec 03/04/2002*/ +int +NonOverlappingTemplateMatchings (int m, int n, int *ARRAY) +{ + FILE *fp; + double sum, chi2, p_value, lambda; + int i, j, k; + char directory[FILENAME_MAX]; + int M, N, K = 5; + int fail = 0; + double pi[6], varWj; + int DATA, PT, pt; + N = 8; + M = floor (n / N); + + + lambda = (M - m + 1) / pow (2, m); + varWj = M * (1. / pow (2., m) - (2. * m - 1.) / pow (2., 2. * m)); + sprintf (directory, "%stemplate%d", GetBaseDir(), m); + if ((fp = fopen (directory, "r")) == NULL) + { +#ifdef IO + { + fprintf (fileout, "\tTemplate file %s not existing\n", directory); + fprintf (stderr, "\tTemplate file %s not existing\n", directory); + } +#endif + exit (1); + + } + else + { + + sum = 0.0; + for (i = 0; i < 2; i++) + { /* Compute Probabilities */ + pi[i] = exp (-lambda + i * log (lambda) - lgam (i + 1)); + sum += pi[i]; + } + pi[0] = sum; + for (i = 2; i <= K; i++) + { /* Compute Probabilities */ + pi[i - 1] = exp (-lambda + i * log (lambda) - lgam (i + 1)); + sum += pi[i - 1]; + } + pi[K] = 1 - sum; + for (i = 0; i < MAXNUMOFTEMPLATES; i++) + { + int inter; + TEMPLATE[i] = 0; + for (j = 0; j < m; j++) + { + if (fscanf (fp, "%d", &inter)<1) inter=0; + TEMPLATE[i] += (inter << j); + } + } + DATA = ARRAY[0] & ((1 << m) - 1); + pt = 0; + PT = m; + + for (i = 0; i < MAXNUMOFTEMPLATES; i++) + for (k = 0; k <= K; k++) + nu[k][i] = 0; + + for (i = 0; i < N; i++) + { + for (k = 0; k < MAXNUMOFTEMPLATES; k++) + W_obs[k] = 0; + for (j = 0; j < M - m + 1; j++) + { + for (k = 0; k < MAXNUMOFTEMPLATES; k++) + { + if (Skip[k] == 0) + { + if (DATA == TEMPLATE[k]) + { + W_obs[k]++; + Skip[k] = m - 1; + } + } + else + { + Skip[k]--; + } + } + PT++; + if (PT == 32) + { + PT = 0; + pt++; + } + DATA = ((DATA << 1) + ((ARRAY[pt] >> PT) & 1)) & ((1 << m) - 1); + } +/* skipping the final values in the slice*/ + for (j = M - m + 1; j < M; j++) + { + PT++; + if (PT == 32) + { + PT = 0; + pt++; + } + DATA = ((DATA << 1) + ((ARRAY[pt] >> PT) & 1)) & ((1 << m) - 1); + } + + for (k = 0; k < MAXNUMOFTEMPLATES; k++) + { + Wj[i][k] = W_obs[k]; + } + } + } + for (k = 0; k < MAXNUMOFTEMPLATES; k++) + { + sum = 0; + chi2 = 0.0; /* Compute Chi Square */ + for (i = 0; i < N; i++) + { + chi2 += pow (((double) Wj[i][k] - lambda) / pow (varWj, 0.5), 2); + } + p_value = igamc (N / 2.0, chi2 / 2.0); + if (p_value < THRESHOLDFAIL) + fail++; + else + TotalTEST++; +#ifdef IO + { + fprintf (fileout, "p_value= %f\n", p_value); + } +#endif + } + fclose (fp); + return (fail); +} + +int +OverlappingTemplateMatchings (int m, int n, int *ARRAY) +{ +/* A. Seznec 03/03/2002: +I got some troubles with porting the function from the NIST package: +computation of array pi did not work satisfactory. +For m=9, I picked back values from page 34 in NIST report. +I do not have the values for m=10. + +*/ + + int i; + double sum, chi2; + int W_obs; + double p_value; + int DATA; + int M, N, j, K = 5; + unsigned int nu[6] = { 0, 0, 0, 0, 0, 0 }; + int fail = 0; + double pi[6] = + { 0.367879, 0.183940, 0.137955, 0.099634, 0.069935, 0.140657 }; + +/* double pi[6] = + { 0.143783, 0.139430, 0.137319, 0.124314, 0.106209, 0.348945 };*/ + + + int PT, pt; + pt = 0; + PT = 0; + + M = 1032; +/* N = (int) floor (n / M);*/ + N = 968; + + +/* + lambda = (double) (M - m + 1) / pow (2, m); + eta = lambda / 2.0; + sum = 0.0; + for (i = 0; i < K; i++) + { + + pi[i] = Pr (i, eta); + sum += pi[i]; + } + pi[K] = 1 - sum;*/ + DATA = ARRAY[0] & ((1 << m) - 1); + pt = 0; + PT = m; + for (i = 0; i < N; i++) + { + W_obs = 0; + for (j = 0; j < M - m + 1; j++) + { + if (DATA == ((1 << m) - 1)) + W_obs++; + PT++; + if (PT == 32) + { + PT = 0; + pt++; + } + DATA = ((DATA << 1) + ((ARRAY[pt] >> PT) & 1)) & ((1 << m) - 1); + } +/* skipping the final values in the slice*/ + for (j = M - m + 1; j < M; j++) + { + PT++; + if (PT == 32) + { + PT = 0; + pt++; + } + DATA = ((DATA << 1) + ((ARRAY[pt] >> PT) & 1)) & ((1 << m) - 1); + } + + if (W_obs <= 4) + nu[W_obs]++; + else + nu[K]++; + } + + sum = 0; + chi2 = 0.0; /* Compute Chi Square */ + for (i = 0; i < K + 1; i++) + { + chi2 += + pow ((double) nu[i] - ((double) N * pi[i]), 2) / ((double) N * pi[i]); + sum += nu[i]; + + } + p_value = igamc (K / 2., chi2 / 2.); + if (p_value < THRESHOLDFAIL) + fail++; + else + TotalTEST++; +#ifdef IO + { + fprintf (fileout, "m= %d, p_value= %f\n", m, p_value); + } +#endif + return (fail); +} + +int +RandomExcursionsVariant (int n, int *ARRAY) +{ + int i, p, J, x, constraint; + double p_value; + int stateX[18] = + { -9, -8, -7, -6, -5, -4, -3, -2, -1, 1, 2, 3, 4, 5, 6, 7, 8, 9 }; + int count; + int fail = 0; + int *S_k; + int pt, PT; + if ((S_k = (int *) MYCALLOC (n, sizeof (int))) == NULL) + { + + } + else + { + J = 0; + S_k[0] = 2 * (ARRAY[0] & 1) - 1; + pt = 0; + PT = 1; + for (i = 1; i < n; i++) + { + S_k[i] = S_k[i - 1] + 2 * ((ARRAY[pt] >> PT) & 1) - 1; + PT++; + if (PT == 32) + { + pt++; + PT = 0; + } + if (S_k[i] == 0) + J++; + } + if (S_k[n - 1] != 0) + J++; + + constraint = MAX (0.005 * pow (n, 0.5), 500); + if (J < constraint) + { + +#ifdef IO + { + fprintf (fileout, + "\n\t\tWARNING: TEST NOT APPLICABLE. THERE ARE "); + fprintf (fileout, "AN\n\t INSUFFICIENT NUMBER OF CYCLES.\n"); + fprintf (fileout, + "\t\t---------------------------------------------"); + fprintf (fileout, "\n"); + } +#endif + } + + else + { + for (p = 0; p <= 17; p++) + { + x = stateX[p]; + count = 0; + for (i = 0; i < n; i++) + if (S_k[i] == x) + count++; + p_value = + erfc (fabs (count - J) / + (sqrt (2. * J * (4. * fabs (x) - 2)))); + if (p_value < THRESHOLDFAIL) + fail++; + else + TotalTEST++; + +#ifdef IO + { + fprintf (fileout, "p_value= %f \n", p_value); + } +#endif + + } + } + } +#ifdef IO + { + fprintf (fileout, "\n"); + } +#endif + MYFREE ((char*)S_k); + return (fail); +} + +int +RandomExcursions (int n, int *ARRAY) +{ + int b, i, j, k, J, x; + int cycleStart, cycleStop, *cycle, *S_k; + int stateX[8] = { -4, -3, -2, -1, 1, 2, 3, 4 }; + int counter[8] = { 0, 0, 0, 0, 0, 0, 0, 0 }; + int PT, pt; + int fail = 0; + double p_value, sum, constraint, nu[6][8]; + double pi[5][6] = { + {0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000}, + {0.5000000000, 0.2500000000, 0.1250000000, 0.06250000000, 0.03125000000, 0.0312500000}, + {0.7500000000, 0.06250000000, 0.04687500000, 0.03515625000, 0.02636718750, 0.0791015625}, + {0.8333333333, 0.02777777778, 0.02314814815, 0.01929012346, 0.01607510288, 0.0803755143}, + {0.8750000000, 0.01562500000, 0.01367187500, 0.01196289063, 0.01046752930, 0.0732727051} + }; + + + S_k = (int *) MYCALLOC (n, sizeof (int)); + cycle = (int *) MYCALLOC (MAX (1000, n / 200), sizeof (int)); + { + J = 0; /* DETERMINE CYCLES */ + + S_k[0] = 2 * (ARRAY[0] & 1) - 1; + pt = 0; + PT = 1; + for (i = 1; i < n; i++) + { + S_k[i] = S_k[i - 1] + 2 * ((ARRAY[pt] >> PT) & 1) - 1; + PT++; + if (PT == 32) + { + pt++; + PT = 0; + } + if (S_k[i] == 0) + { + J++; + if (J > MAX (1000, n / 128)) + { +#ifdef IO + { + fprintf + (fileout, + "ERROR IN FUNCTION randomExcursions: EXCEEDING THE MAX"); + fprintf (fileout, " NUMBER OF CYCLES EXPECTED\n."); + } +#endif + MYFREE ((char*)cycle); + MYFREE ((char*)S_k); + + return (fail); + } + cycle[J] = i; + } + } + if (S_k[n - 1] != 0) + { + J++; + + } + + + constraint = MAX (0.005 * pow (n, 0.5), 500); + if (J < constraint) + { + +#ifdef IO + { + fprintf (fileout, + "\t\t---------------------------------------------"); + fprintf (fileout, + "\n\t\tWARNING: TEST NOT APPLICABLE. THERE ARE "); + fprintf (fileout, "AN\n\t INSUFFICIENT NUMBER OF CYCLES.\n"); + fprintf (fileout, + "\t\t---------------------------------------------"); + fprintf (fileout, "\n"); + } +#endif + + } + else + { + cycleStart = 0; + cycleStop = cycle[1]; + for (k = 0; k < 6; k++) + for (i = 0; i < 8; i++) + nu[k][i] = 0.; + for (j = 1; j <= J; j++) + { /* FOR EACH CYCLE */ + for (i = 0; i < 8; i++) + counter[i] = 0; + for (i = cycleStart; i < cycleStop; i++) + { + if ((S_k[i] >= 1 && S_k[i] <= 4) + || (S_k[i] >= -4 && S_k[i] <= -1)) + { + if (S_k[i] < 0) + b = 4; + else + b = 3; + counter[S_k[i] + b]++; + } + } + cycleStart = cycle[j] + 1; + if (j < J) + cycleStop = cycle[j + 1]; + else + cycleStop = n; + + for (i = 0; i < 8; i++) + { + if (counter[i] >= 0 && counter[i] <= 4) + nu[counter[i]][i]++; + else if (counter[i] >= 5) + nu[5][i]++; + } + } + for (i = 0; i < 8; i++) + { + x = stateX[i]; + sum = 0.; + for (k = 0; k < 6; k++) + { + sum += pow (nu[k][i] - J * pi[(int) fabs (x)][k], 2) / + (J * pi[(int) fabs (x)][k]); + } + p_value = igamc (2.5, sum / 2.); + + if (p_value < THRESHOLDFAIL) + fail++; + else + TotalTEST++; + +#ifdef IO + { + fprintf (fileout, "p_value= %f\n", p_value); + } +#endif + + + + } + } + +#ifdef IO + { + fprintf (fileout, "\n"); + } +#endif + + MYFREE ((char*)cycle); + MYFREE ((char*)S_k); + } + return (fail); +} + +int +CumulativeSums (int n, int *ARRAY) +{ + int i, k, start, finish, mode; + double p_value, cusum, sum, sum1, sum2; + int z; + int pt, PT; + int fail = 0; + for (mode = 0; mode < 2; mode++) + { /* mode = {0,1} => {forward,reverse} */ + sum = 0.0; + cusum = 1.0; + if (mode == 0) + { + pt = 0; + PT = 0; + + for (i = 0; i < n; i++) + { + sum += (double) (2 * ((ARRAY[pt] >> PT) & 1) - 1); + PT++; + if (PT == 32) + { + PT = 0; + pt++; + } + cusum = MAX (cusum, fabs (sum)); + } + } + else if (mode == 1) + { + pt = (n >> 5); + PT = 31; + for (i = n - 1; i >= 0; i--) + { + sum += (double) (2 * ((ARRAY[pt] >> PT) & 1) - 1); + PT--; + if (PT == -1) + { + PT = 31; + pt--; + } + cusum = MAX (cusum, fabs (sum)); + } + } + z = (int) cusum; + + sum1 = 0.0; + start = (-n / z + 1) / 4; + finish = (n / z - 1) / 4; + for (k = start; k <= finish; k++) + sum1 += + (normal ((4 * k + 1) * z / sqrt (n)) - + normal ((4 * k - 1) * z / sqrt (n))); + + sum2 = 0.0; + start = (-n / z - 3) / 4; + finish = (n / z - 1) / 4; + for (k = start; k <= finish; k++) + sum2 += + (normal ((double) ((4 * k + 3) * z) / sqrt ((double) n)) - + normal ((double) ((4 * k + 1) * z) / sqrt (n))); + p_value = 1.0 - sum1 + sum2; + if (mode == 1) +{ if (p_value < THRESHOLDFAIL) + fail++; + else + TotalTEST++; + +#ifdef IO + { + fprintf (fileout, "%d bits sequence, reverse p_value= %f \n", n, + p_value); + } +#endif + } + if (mode == 0){ + if (p_value < THRESHOLDFAIL) + fail++; + else + TotalTEST++; + +#ifdef IO + { + fprintf (fileout, "%d bits sequence, forward p_value= %f \n", n, + p_value); + } +#endif +} + } + return (fail); + +} + + + +int +ApproximateEntropy (int mmin, int mmax, int n, int *ARRAY) +{ + int i, blockSize, seqLength; + int powLen; + double sum, numOfBlocks, ApEn[25], apen, chi_squared, p_value; + unsigned int *P[25]; + int pt, PT, DATA; + int fail = 0; + int MaskEnt[25] = + { 0, 1, 3, 7, 0xf, 0x1f, 0x3f, 0x7f, 0xff, 0x1ff, 0x3ff, 0x7ff, 0xfff, + 0x1fff, 0x3fff, 0x7fff, 0xffff, 0x1ffff, 0x3ffff, 0x7ffff, 0xfffff, + 0x1fffff, + 0x3fffff, 0x7fffff, 0xffffff + }; + + seqLength = n - (mmax); + numOfBlocks = (double) seqLength; + for (blockSize = mmin; blockSize <= mmax; blockSize++) + { + powLen = (int) pow (2, blockSize + 1); + if ((P[blockSize] = + (unsigned int *) MYCALLOC (powLen, sizeof (unsigned int))) == NULL) + { + +#ifdef IO + { + fprintf (fileout, "ApEn: Insufficient memory available.\n"); + } +#endif + exit (1); + return (fail); + } + for (i = 0; i < powLen; i++) + P[blockSize][i] = 0; + } + + for (blockSize = mmin; blockSize <= mmax; blockSize++) + { + DATA = ARRAY[0]; + pt = 0; + PT = mmax; + for (i = 0; i < seqLength; i++) + { /* COMPUTE FREQUENCY */ + + (P[blockSize][DATA & MaskEnt[blockSize]])++; + PT++; + if (PT == 32) + { + PT = 0; + pt++; + }; + DATA = (DATA << 1) + ((ARRAY[pt] >> PT) & 1); + } + } + + for (blockSize = mmin; blockSize <= mmax; blockSize++) + { + sum = 0.0; + + for (i = 0; i < (int) pow (2, blockSize); i++) + { + if (P[blockSize][i] > 0) + sum += + (((double) P[blockSize][i]) / (numOfBlocks)) * + log ((double) P[blockSize][i] / numOfBlocks); + } + + ApEn[blockSize] = sum; + } + for (blockSize = mmax; blockSize >= mmin; blockSize--) + MYFREE ((char*)P[blockSize]); + + for (i = mmin; i < mmax; i++) + { + apen = ApEn[i] - ApEn[i + 1]; + chi_squared = 2.0 * ((double) seqLength) * (log (2) - apen); + p_value = igamc (pow (2, i - 1), chi_squared / 2.); + if (p_value < THRESHOLDFAIL) + fail++; + else + TotalTEST++; +#ifdef IO + { + fprintf (fileout, "m= %d,\t p_value= %f, Entropy per %d bits %f \n", + i, p_value, i, -ApEn[i] / log (2)); + } +#endif + } + + return (fail); +} + +int +Serial (int m, int n, int *ARRAY, int PT) +{ + double p_value1, p_value2, psim0, psim1, psim2, del1, del2; + int fail = 0; + psim0 = psi2 (m, n, &ARRAY[PT >> 5], (PT & 31)); + psim1 = psi2 (m - 1, n, &ARRAY[PT >> 5], (PT & 31)); + psim2 = psi2 (m - 2, n, &ARRAY[PT >> 5], (PT & 31)); + del1 = psim0 - psim1; + del2 = psim0 - 2.0 * psim1 + psim2; + p_value1 = igamc (pow (2, m - 1) / 2, del1 / 2.0); + p_value2 = igamc (pow (2, m - 2) / 2, del2 / 2.0); + + if (p_value1 < THRESHOLDFAIL) + fail++; + else + TotalTEST++; +#ifdef IO + { + fprintf (fileout, "p_value1= %f\n", p_value1); + } +#endif + if (p_value1 < THRESHOLDFAIL) + fail++; + else + TotalTEST++; +#ifdef IO + { + fprintf (fileout, "p_value2= %f\n", p_value2); + } +#endif + return (fail); +} + + +double +psi2 (int m, int n, int *ARRAY, int PT) +{ + int i, j, k; + double sum, numOfBlocks; + unsigned int *P; + + if ((m == 0) || (m == -1)) + return 0.0; + numOfBlocks = n; + + P = (unsigned int*) MYCALLOC (65536, sizeof (unsigned int)); + for (i = 0; i < numOfBlocks; i++) + { /* COMPUTE FREQUENCY */ + k = 1; + for (j = 0; j < m; j++) + { + + if (((ARRAY[(MOD (PT + i + j, n) >> 5)] >> + (MOD (PT + i + j, n) & 31)) & 1) == 0) + k *= 2; + else + k = 2 * k + 1; + } + P[k - 1]++; + } + sum = 0.0; + for (i = (int) pow (2, m) - 1; i < (int) pow (2, m + 1) - 1; i++) + sum += pow (P[i], 2); + sum = (sum * pow (2, m) / (double) n) - (double) n; + MYFREE ((char*)P); + return sum; +} + + +int +LinearComplexity (int M, int N, int *ARRAY, int PT) +{ + int i, ii, j, d; + int L, m, N_, parity, sign; + double p_value, T_, mean; + int fail = 0; + int K = 6; + double pi[7] = + { 0.01047, 0.03125, 0.12500, 0.50000, 0.25000, 0.06250, 0.020833 }; + double nu[7], chi2; + int *T, *P, *B_, *C; + + B_ = (int *) MYCALLOC (M, sizeof (int)); + C = (int *) MYCALLOC (M, sizeof (int)); + P = (int *) MYCALLOC (M, sizeof (int)); + T = (int *) MYCALLOC (M, sizeof (int)); + + for (i = 0; i < K + 1; i++) + nu[i] = 0.00; + for (ii = 0; ii < N; ii++) + { + for (i = 0; i < M; i++) + { + B_[i] = 0; + C[i] = 0; + T[i] = 0; + P[i] = 0; + } + L = 0; + m = -1; + d = 0; + C[0] = 1; + B_[0] = 1; + /* DETERMINE LINEAR COMPLEXITY */ + N_ = 0; + + while (N_ < M) + { + d = + ((ARRAY[(ii * M + N_ + PT) >> 5]) >> + ((ii * M + N_ + PT) & 31)) & 1; + + for (i = 1; i <= L; i++) + + d += + (C[i] & + (((ARRAY[(ii * M + N_ - i + PT) >> 5]) >> + ((ii * M + N_ - i + PT) & 31)) & 1)); + d = d & 1; + if (d == 1) + { + for (i = 0; i < M; i++) + { + T[i] = C[i]; + P[i] = 0; + } + for (j = 0; j < M; j++) + if (B_[j] == 1) + P[j + N_ - m] = 1; + for (i = 0; i < M; i++) + C[i] = (C[i] + P[i]) & 1; + if (L <= N_ / 2) + { + L = N_ + 1 - L; + m = N_; + for (i = 0; i < M; i++) + B_[i] = T[i]; + } + } + N_++; + } + if ((parity = (M + 1) % 2) == 0) + sign = -1; + else + sign = 1; + mean = + M / 2. + (9. + sign) / 36. - 1. / pow (2, M) * (M / 3. + 2. / 9.); + if ((parity = M % 2) == 0) + sign = 1; + else + sign = -1; + T_ = sign * (L - mean) + 2. / 9.; + + if (T_ <= -2.5) + nu[0]++; + else if (T_ > -2.5 && T_ <= -1.5) + nu[1]++; + else if (T_ > -1.5 && T_ <= -0.5) + nu[2]++; + else if (T_ > -0.5 && T_ <= 0.5) + nu[3]++; + else if (T_ > 0.5 && T_ <= 1.5) + nu[4]++; + else if (T_ > 1.5 && T_ <= 2.5) + nu[5]++; + else + nu[6]++; + + + } + chi2 = 0.00; + for (i = 0; i < K + 1; i++) + chi2 += pow (nu[i] - N * pi[i], 2) / (N * pi[i]); + p_value = igamc (K / 2.0, chi2 / 2.0); + if (p_value < THRESHOLDFAIL) + fail++; + else + TotalTEST++; +#ifdef IO + { + fprintf (fileout, "p_value= %f\n", p_value); + } +#endif + + MYFREE ((char*)B_); + return (fail); +} + + +int +LempelZivCompression (int n, int *ARRAY, int *DATA, int *pt, int *PT) +{ + int W; /* Number of words */ + int i, j, k, prev_I, powLen, max_len; + int done = 0; + double p_value, mean=0.0, variance=0.0; + BitField *P; + int fail = 0; + W = 0; + k = (int) log (n) / log (2) + 6; + powLen = pow (2, k); + + if ((P = (BitField *) MYCALLOC (powLen, sizeof (BitField))) == NULL) + { + +#ifdef IO + { + fprintf (fileout, "\t\tLempel-Ziv Compression has been aborted,\n"); + fprintf (fileout, "\t\tdue to insufficient workspace!\n"); + } +#endif + } + else + { + for (i = 0; i < powLen; i++) + P[i].b = 0; + i = 0; + max_len = 0; + while (i <= n - 1) + { + done = 0; + j = 1; + prev_I = i; + while (!done) + { + if (2 * j + 1 <= powLen) + { + if ((*DATA & 1) == 0) + { + if (P[2 * j].b == 1) + { + j *= 2; + } + else + { + P[2 * j].b = 1; + done = 1; + } + } + else + { + if (P[2 * j + 1].b == 1) + { + j = 2 * j + 1; + } + else + { + P[2 * j + 1].b = 1; + done = 1; + } + } + (*PT)++; + if (*PT == 32) + { + (*pt)++; + *DATA = ARRAY[*pt]; + *PT = 0; + } + else + *DATA = *DATA >> 1; + i++; + if (i > n - 1) + { + done = 1; + } + if (done) + { + W++; + max_len = MAX (max_len, i - prev_I); + } + } + else + { + +#ifdef IO + { + fprintf (fileout, + "\t\tWARNING: Segmentation Violation Imminent."); + fprintf (fileout, + "\n\t Lempel-Ziv Compression Terminated.\n"); + fprintf (fileout, + "\t\t-----------------------------------------\n"); + fflush (fileout); + } +#endif + done = 1; + i = n; + } + } + } + } + switch (n) + { + case 100000: + mean = 8782.500000; + variance = 11.589744; + break; + case 200000: + mean = 16292.1000; + variance = 21.4632; + break; + case 400000: + mean = 30361.9500; + variance = 58.7868; + break; + case 600000: + mean = 43787.5000; + variance = 70.1579; + break; + case 800000: + mean = 56821.9500; + variance = 67.4184; + break; + case 1000000: /* Updated Mean and Variance 10/26/99 */ + mean = 69588.20190000; + variance = 73.23726011; + /* Previous Mean and Variance + mean = 69586.250000; + variance = 70.448718; + */ + break; + default: + break; + } + p_value = 0.5 * erfc ((mean - W) / pow (2.0 * variance, 0.5)); + if (p_value < THRESHOLDFAIL) + fail++; + else + TotalTEST++; +#ifdef IO + { + fprintf (fileout, "p_value= %f\n", p_value); + } +#endif + MYFREE ((char*)P); + return (fail); +} + + +int +DiscreteFourierTransform (int N, int *ARRAY) +{ + double p_value, upperBound, *m, *X; + int i, count; + double N_l, N_o, d; + double *wsave; + int n, j, J, *ifac; + int fail = 0; + n = 32 * N; + X = (double *) MYCALLOC (n, sizeof (double)); + wsave = (double *) MYCALLOC (2 * n + 15, sizeof (double)); + ifac = (int *) MYCALLOC (15, sizeof (double)); + m = (double *) MYCALLOC (n / 2 + 1, sizeof (double)); + { + J = 0; + for (i = 0; i < N; i++) + for (j = 0; j < 32; j++) + { + X[J] = (2 * ((ARRAY[i] >> j) & 1)) - 1; + J++; + } + __ogg_fdrffti (n, wsave, ifac); /* INITIALIZE WORK ARRAYS */ + __ogg_fdrfftf (n, X, wsave, ifac); /* APPLY FORWARD FFT */ + + m[0] = sqrt (X[0] * X[0]); /* COMPUTE MAGNITUDE */ + + + for (i = 0; i < n / 2; i++) + { /* DISPLAY FOURIER POINTS */ + m[i + 1] = sqrt (pow (X[2 * i + 1], 2) + pow (X[2 * i + 2], 2)); + } + count = 0; /* CONFIDENCE INTERVAL */ + upperBound = sqrt (3 * n); + for (i = 0; i < n / 2; i++) + if (m[i] < upperBound) + count++; + N_l = (double) count; /* number of peaks less than h = sqrt(3*n) */ + N_o = (double) 0.95 *n / 2.; + d = (N_l - N_o) / sqrt (n / 2. * 0.95 * 0.05); + p_value = erfc (fabs (d) / sqrt (2.)); + + if (p_value < THRESHOLDFAIL) + fail++; + else + TotalTEST++; + +#ifdef IO + { + fprintf (fileout, "p_value= %f\n", p_value); + } +#endif + } + MYFREE ((char*)m); + MYFREE ((char*)ifac); + MYFREE ((char*)wsave); + MYFREE ((char*)X); + return (fail); +} + +int +LongestRunOfOnes (int n, int *ARRAY) +{ + double p_value, sum, chi2; + int N, i, j; + int run, v_n_obs; + int fail = 0; +/* since we are not interested in short sequences we used only 10000*/ + + + double pi[7] = { 0.0882, 0.2092, 0.2483, 0.1933, 0.1208, 0.0675, 0.0727 }; + double K = 6; + unsigned int nu[7] = { 0, 0, 0, 0, 0, 0, 0 }; + int M = 10000; + int PT; + int pt; + int DATA; + pt = 0; + PT = 0; + DATA = ARRAY[0]; + N = (int) floor ((32 * n) / M); + for (i = 0; i < N; i++) + { + v_n_obs = 0; + run = 0; + for (j = i * M; j < (i + 1) * M; j++) + { + if (DATA & 1) + { + run++; + v_n_obs = MAX (v_n_obs, run); + } + else + run = 0; + PT++; + if (PT == 32) + { + PT = 0; + pt++; + DATA = ARRAY[pt]; + } + else + DATA = DATA >> 1; + } + if (v_n_obs <= 10) + nu[0]++; + else if (v_n_obs == 11) + nu[1]++; + else if (v_n_obs == 12) + nu[2]++; + else if (v_n_obs == 13) + nu[3]++; + else if (v_n_obs == 14) + nu[4]++; + else if (v_n_obs == 15) + nu[5]++; + else if (v_n_obs >= 16) + nu[6]++; + } + chi2 = 0.0; /* Compute Chi Square */ + sum = 0; + for (i = 0; i < ((int) K) + 1; i++) + { + chi2 += + pow ((double) nu[i] - ((double) N * pi[i]), 2) / ((double) N * pi[i]); + sum += nu[i]; + } + p_value = igamc (K / 2., chi2 / 2.); + if (p_value < THRESHOLDFAIL) + fail++; + else + TotalTEST++; +#ifdef IO + { + fprintf (fileout, "p_value= %f\n", p_value); + } +#endif + return (fail); +} + +int +Runs (int n, int *ARRAY) +{ + int i, j; + double argument, pi, V_n_obs, tau; + double p_value, product; + int SUM; + double N; + int fail = 0; + int V_N_obs; + N = (double) (32 * (n - 1) + 1); + + SUM = 0; + for (i = 0; i < n - 1; i++) + for (j = 0; j < 32; j++) + { + SUM += (ARRAY[i] >> j) & 1; + } + SUM += (ARRAY[n] & 1); + + pi = ((double) SUM) / N; + tau = 2.0 / sqrt (N); + + if (fabs (pi - 0.5) < tau) + { + V_N_obs = 0; + for (i = 0; i < n - 1; i++) + { + for (j = 0; j < 31; j++) + V_N_obs += (((ARRAY[i] >> j) ^ (ARRAY[i] >> (j + 1))) & 1); + V_N_obs += ((ARRAY[i] >> 31) ^ (ARRAY[i + 1])) & 1; + } + + V_N_obs++; + V_n_obs = (double) V_N_obs; + product = pi * (1.e0 - pi); + argument = + fabs (V_n_obs - + 2.e0 * N * product) / (2.e0 * sqrt (2.e0 * N) * product); + p_value = erfc (argument); + } + else + { + p_value = 0.0; + } + if (p_value < THRESHOLDFAIL) + fail++; + else + TotalTEST++; +#ifdef IO + { + fprintf (fileout, "p_value=%f\n", p_value); + } +#endif + return (fail); +} + +int +Frequency (int n, int *ARRAY) +{ + int i, j; + double f, s_obs, p_value; + double sqrt2 = 1.41421356237309504880; + int SUM; + int fail = 0; + SUM = 0; + + for (i = 0; i < n; i++) + { + for (j = 0; j < 32; j++) + SUM += (2 * ((ARRAY[i] >> j) & 1)) - 1; + } + s_obs = fabs ((double) SUM) / sqrt (32.0 * ((double) n)); + f = s_obs / sqrt2; + p_value = erfc (f); + if (p_value < THRESHOLDFAIL) + fail++; + else + TotalTEST++; +#ifdef IO + { + fprintf (fileout, "%d bits sequence, p_value= %f\n", 32 * n, p_value); + } +#endif + return (fail); +} + +int +BlockFrequency (int ArraySize, int m, int *ARRAY) +{ + int i, j, N, n; + double arg1, arg2, p_value; + double sum, pi, v, chi_squared; + int BlockSum; + int PT; + int pt; + int DATA; + int fail = 0; + n = ArraySize * 32; + pt = 0; + PT = 0; + DATA = ARRAY[0]; + N = (int) floor ((double) n / (double) m); + sum = 0.0; + for (i = 0; i < N; i++) + { + pi = 0.0; + BlockSum = 0; + for (j = 0; j < m; j++) + { + BlockSum += (DATA & 1); + PT++; + if (PT == 32) + { + PT = 0; + pt++; + DATA = ARRAY[pt]; + } + else + DATA = DATA >> 1; + } + + pi = (double) BlockSum / (double) m; + v = pi - 0.5; + sum += v * v; + } + chi_squared = 4.0 * m * sum; + arg1 = (double) N / 2.e0; + arg2 = chi_squared / 2.e0; + p_value = igamc (arg1, arg2); + + if (p_value < THRESHOLDFAIL) + fail++; + else + TotalTEST++; +#ifdef IO + { + fprintf (fileout, " %d bits sequence, p_value= %f\n", n, p_value); + } +#endif + return (fail); +} + + + + + + +int SizeMin[18] = + { 0, 0, 0, 0, 0, 0, 387840, 904960, 2068480, 4654080, 10342400, 22753280, + 49643520, 107560960, 231669760, 496435200, 1059061760, 0x7fffffff +}; + + +int +Universal (int n, int *ARRAY) +{ + int I, param; + int fail = 0; + if (32 * n < SizeMin[6]) + { + +#ifdef IO + { + fprintf (fileout, "Too small\n"); + } +#endif + return (fail); + } + I = 6; + while ((32 * n >= SizeMin[I]) & (I < 17)) + { + if (32 * n >= SizeMin[I + 1]) + param = SizeMin[I + 1] - 1; + else + param = 32 * n; + fail += UNIVERSAL (param, ARRAY); + I++; + } + return (fail); +} + + + + +int +UNIVERSAL (int n, int *ARRAY) +{ + int i, j, p, K, L, Q; + double arg, sqrt2, sigma, phi, sum, p_value, c; + long *T, decRep; + double expected_value[17] = { + 0, 0, 0, 0, 0, 0, 5.2177052, 6.1962507, 7.1836656, + 8.1764248, 9.1723243, 10.170032, 11.168765, + 12.168070, 13.167693, 14.167488, 15.167379 + }; + double variance[17] = { + 0, 0, 0, 0, 0, 0, 2.954, 3.125, 3.238, 3.311, 3.356, 3.384, + 3.401, 3.410, 3.416, 3.419, 3.421 + }; + int PT, DATA, pt; + int fail = 0; + double Pow[20]; + double POW; + if (n >= 387840) + L = 6; + if (n >= 904960) + L = 7; + if (n >= 2068480) + L = 8; + if (n >= 4654080) + L = 9; + if (n >= 10342400) + L = 10; + if (n >= 22753280) + L = 11; + if (n >= 49643520) + L = 12; + if (n >= 107560960) + L = 13; + if (n >= 231669760) + L = 14; + if (n >= 496435200) + L = 15; + if (n >= 1059061760) + L = 16; + PT = 0; + pt = 0; + DATA = ARRAY[pt]; + POW = pow (2, L); + for (i = 0; i <= L; i++) + Pow[i] = pow (2, i); + Q = (int) 10 *POW; + K = (int) (floor (n / L) - (double) Q); + c = + 0.7 - 0.8 / (double) L + (4 + 32 / (double) L) * pow (K, + -3 / + (double) L) / 15; + sigma = c * sqrt (variance[L] / (double) K); + sqrt2 = sqrt (2); + sum = 0.0; + p = (int) pow (2, L); + T = (long *) MYCALLOC (p, sizeof (long)); + for (i = 0; i < p; i++) + T[i] = 0; + for (i = 1; i <= Q; i++) + { + decRep = 0; + for (j = 0; j < L; j++) + { + decRep += (DATA & 1) * Pow[L - 1 - j]; + PT++; + if (PT == 32) + { + pt++; + DATA = ARRAY[pt]; + PT = 0; + } + else + DATA = DATA >> 1; + } + T[decRep] = i; + } + for (i = Q + 1; i <= Q + K; i++) + { + decRep = 0; + for (j = 0; j < L; j++) + { + decRep += (DATA & 1) * Pow[L - 1 - j]; + PT++; + if (PT == 32) + { + pt++; + DATA = ARRAY[pt]; + PT = 0; + } + else + DATA = DATA >> 1; + } + sum += log (i - T[decRep]) / log (2); + T[decRep] = i; + } + phi = (double) (sum / (double) K); + arg = fabs (phi - expected_value[L]) / (sqrt2 * sigma); + p_value = erfc (arg); + MYFREE ((char*)T); + if (p_value < THRESHOLDFAIL) + fail++; + else + TotalTEST++; +#ifdef IO + { + fprintf (fileout, "L= %d\tp_value %f\texp_value= %f\tphi = %f\n", L, + p_value, expected_value[L], phi); + } +#endif + return (fail); +} + + +BitField ** +create_matrix (int M, int Q) +{ + int i; + BitField **matrix; + + if ((matrix = (BitField **) MYCALLOC (M, sizeof (BitField *))) == NULL) + { + printf + ("ERROR IN FUNCTION create_matrix: Insufficient memory available."); + printf ("\n"); + fprintf (stderr, + "CREATE_MATRIX: Insufficient memory for %dx%d matrix.\n", M, + M); + return matrix; + } + else + { + for (i = 0; i < M; i++) + { + if ((matrix[i] = (BitField *)MYCALLOC (Q, sizeof (BitField))) == NULL) + { + fprintf (stderr, + "CREATE_MATRIX: Insufficient memory for %dx%d ", M, + M); + fprintf (stderr, "matrix.\n"); + printf + ("ERROR IN FUNCTION create_matrix: Insufficient memory for "); + printf ("%dx%d matrix.\n", M, M); + return NULL; + } + } + return matrix; + } +} + + +int +Rank (int n, int *ARRAY) +{ + int N = (int) floor (n / (32)); /* NUMBER OF MATRICES */ + int r; + double p_value, product; + int i, k; + double chi_squared, arg1; + double p_32, p_31, p_30; /* PROBABILITIES */ + double R; /* RANK */ + double F_32, F_31, F_30; /* FREQUENCIES */ + BitField **matrix = create_matrix (32, 32); + int pt, PT, DATA; + int fail = 0; + pt = 0; + PT = 0; + DATA = ARRAY[0]; + r = 32; /* COMPUTE PROBABILITIES */ + product = 1; + for (i = 0; i <= r - 1; i++) + product *= + ((1.e0 - pow (2, i - 32)) * (1.e0 - pow (2, i - 32))) / (1.e0 - + pow (2, + i - r)); + p_32 = pow (2, r * (32 + 32 - r) - 32 * 32) * product; + + r = 31; + product = 1; + for (i = 0; i <= r - 1; i++) + product *= + ((1.e0 - pow (2, i - 32)) * (1.e0 - pow (2, i - 32))) / (1.e0 - + pow (2, + i - r)); + p_31 = pow (2, r * (32 + 32 - r) - 32 * 32) * product; + + p_30 = 1 - (p_32 + p_31); + + F_32 = 0; + F_31 = 0; + for (k = 0; k < N; k++) + { /* FOR EACH 32x32 MATRIX */ + def_matrix (32, 32, matrix, k, &pt, &PT, &DATA, ARRAY); + R = computeRank (32, 32, matrix); + if (R == 32) + F_32++; /* DETERMINE FREQUENCIES */ + if (R == 31) + F_31++; + } + F_30 = (double) N - (F_32 + F_31); + + chi_squared = (pow (F_32 - N * p_32, 2) / (double) (N * p_32) + + pow (F_31 - N * p_31, 2) / (double) (N * p_31) + + pow (F_30 - N * p_30, 2) / (double) (N * p_30)); + + arg1 = -chi_squared / 2.e0; + + p_value = exp (arg1); + + MYFREE ((char*)matrix); + + if (p_value < THRESHOLDFAIL) + fail++; + else + TotalTEST++; +#ifdef IO + { + fprintf (fileout, "p_value= %f\n", p_value); + } +#endif + return (fail); +} +#else +#include <stdio.h> +int +PackTestF (int *ARRAY, int ArraySize, char *C) +{ + fprintf (stderr, + "Sorry, on-line analysis is implemented using the CygWin math libraries,\nyou must recompile the application under the CygWin environment to allow online analysis\n"); +} +int +PackTestL (int *ARRAY, int ArraySize, char *C) +{ + fprintf (stderr, + "Sorry, on-line analysis is implemented using the CygWin math libraries,\nyou must recompile the application under the CygWin environment to allow online analysis\n"); +} +#endif diff --git a/nist/special-functions.c b/nist/special-functions.c new file mode 100644 index 0000000..c15b7e1 --- /dev/null +++ b/nist/special-functions.c @@ -0,0 +1,32 @@ +#include <stdio.h> +#include <math.h> +/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * + S P E C I A L F U N C T I O N S + * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ +#define LinuxorUnix +#ifdef WIN +#ifndef CYGWIN +#undef LinuxorUnix +/* same libraries are available*/ +#endif +#endif + +#ifdef LinuxorUnix +double normal(double x) +{ + double arg, result, sqrt2=1.414213562373095048801688724209698078569672; + + if (x > 0) { + arg = x/sqrt2; + result = 0.5 * ( 1 + erf(arg) ); + } + else { + arg = -x/sqrt2; + result = 0.5 * ( 1 - erf(arg) ); + } + return( result); +} + +double normal2(double a) +{ return (1.0-0.5*erfc(a/sqrt(2.0))); } +#endif diff --git a/nist/special-functions.h b/nist/special-functions.h new file mode 100644 index 0000000..61e4839 --- /dev/null +++ b/nist/special-functions.h @@ -0,0 +1,14 @@ +/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * + S P E C I A L F U N C T I O N P R O T O T Y P E S +* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ + +double derf_(); /* ERROR FUNCTION */ +double derfc(); /* COMPLEMENTARY ERROR */ +double normal(double); /* NORMAL DISTRIBUTION */ +double normal2(double); /* NORMAL DISTRIBUTION */ +/* +double gser(double,double,double,double); +double gcf(double,double,double,double); +double gammln(double); LOG GAMMA FUNCTION +double gammq(double,double); INCOMPLETE GAMMA FUNCTION +*/ diff --git a/nist/template9 b/nist/template9 new file mode 100644 index 0000000..604bc53 --- /dev/null +++ b/nist/template9 @@ -0,0 +1,148 @@ +0 0 0 0 0 0 0 0 1 +0 0 0 0 0 0 0 1 1 +0 0 0 0 0 0 1 0 1 +0 0 0 0 0 0 1 1 1 +0 0 0 0 0 1 0 0 1 +0 0 0 0 0 1 0 1 1 +0 0 0 0 0 1 1 0 1 +0 0 0 0 0 1 1 1 1 +0 0 0 0 1 0 0 0 1 +0 0 0 0 1 0 0 1 1 +0 0 0 0 1 0 1 0 1 +0 0 0 0 1 0 1 1 1 +0 0 0 0 1 1 0 0 1 +0 0 0 0 1 1 0 1 1 +0 0 0 0 1 1 1 0 1 +0 0 0 0 1 1 1 1 1 +0 0 0 1 0 0 0 1 1 +0 0 0 1 0 0 1 0 1 +0 0 0 1 0 0 1 1 1 +0 0 0 1 0 1 0 0 1 +0 0 0 1 0 1 0 1 1 +0 0 0 1 0 1 1 0 1 +0 0 0 1 0 1 1 1 1 +0 0 0 1 1 0 0 1 1 +0 0 0 1 1 0 1 0 1 +0 0 0 1 1 0 1 1 1 +0 0 0 1 1 1 0 0 1 +0 0 0 1 1 1 0 1 1 +0 0 0 1 1 1 1 0 1 +0 0 0 1 1 1 1 1 1 +0 0 1 0 0 0 0 1 1 +0 0 1 0 0 0 1 0 1 +0 0 1 0 0 0 1 1 1 +0 0 1 0 0 1 0 1 1 +0 0 1 0 0 1 1 0 1 +0 0 1 0 0 1 1 1 1 +0 0 1 0 1 0 0 1 1 +0 0 1 0 1 0 1 0 1 +0 0 1 0 1 0 1 1 1 +0 0 1 0 1 1 0 1 1 +0 0 1 0 1 1 1 0 1 +0 0 1 0 1 1 1 1 1 +0 0 1 1 0 0 1 0 1 +0 0 1 1 0 0 1 1 1 +0 0 1 1 0 1 0 1 1 +0 0 1 1 0 1 1 0 1 +0 0 1 1 0 1 1 1 1 +0 0 1 1 1 0 1 0 1 +0 0 1 1 1 0 1 1 1 +0 0 1 1 1 1 0 1 1 +0 0 1 1 1 1 1 0 1 +0 0 1 1 1 1 1 1 1 +0 1 0 0 0 0 0 1 1 +0 1 0 0 0 0 1 1 1 +0 1 0 0 0 1 0 1 1 +0 1 0 0 0 1 1 1 1 +0 1 0 0 1 0 0 1 1 +0 1 0 0 1 0 1 1 1 +0 1 0 0 1 1 0 1 1 +0 1 0 0 1 1 1 1 1 +0 1 0 1 0 0 0 1 1 +0 1 0 1 0 0 1 1 1 +0 1 0 1 0 1 0 1 1 +0 1 0 1 0 1 1 1 1 +0 1 0 1 1 0 0 1 1 +0 1 0 1 1 0 1 1 1 +0 1 0 1 1 1 0 1 1 +0 1 0 1 1 1 1 1 1 +0 1 1 0 0 0 1 1 1 +0 1 1 0 0 1 1 1 1 +0 1 1 0 1 0 1 1 1 +0 1 1 0 1 1 1 1 1 +0 1 1 1 0 1 1 1 1 +0 1 1 1 1 1 1 1 1 +1 0 0 0 0 0 0 0 0 +1 0 0 0 1 0 0 0 0 +1 0 0 1 0 0 0 0 0 +1 0 0 1 0 1 0 0 0 +1 0 0 1 1 0 0 0 0 +1 0 0 1 1 1 0 0 0 +1 0 1 0 0 0 0 0 0 +1 0 1 0 0 0 1 0 0 +1 0 1 0 0 1 0 0 0 +1 0 1 0 0 1 1 0 0 +1 0 1 0 1 0 0 0 0 +1 0 1 0 1 0 1 0 0 +1 0 1 0 1 1 0 0 0 +1 0 1 0 1 1 1 0 0 +1 0 1 1 0 0 0 0 0 +1 0 1 1 0 0 1 0 0 +1 0 1 1 0 1 0 0 0 +1 0 1 1 0 1 1 0 0 +1 0 1 1 1 0 0 0 0 +1 0 1 1 1 0 1 0 0 +1 0 1 1 1 1 0 0 0 +1 0 1 1 1 1 1 0 0 +1 1 0 0 0 0 0 0 0 +1 1 0 0 0 0 0 1 0 +1 1 0 0 0 0 1 0 0 +1 1 0 0 0 1 0 0 0 +1 1 0 0 0 1 0 1 0 +1 1 0 0 1 0 0 0 0 +1 1 0 0 1 0 0 1 0 +1 1 0 0 1 0 1 0 0 +1 1 0 0 1 1 0 0 0 +1 1 0 0 1 1 0 1 0 +1 1 0 1 0 0 0 0 0 +1 1 0 1 0 0 0 1 0 +1 1 0 1 0 0 1 0 0 +1 1 0 1 0 1 0 0 0 +1 1 0 1 0 1 0 1 0 +1 1 0 1 0 1 1 0 0 +1 1 0 1 1 0 0 0 0 +1 1 0 1 1 0 0 1 0 +1 1 0 1 1 0 1 0 0 +1 1 0 1 1 1 0 0 0 +1 1 0 1 1 1 0 1 0 +1 1 0 1 1 1 1 0 0 +1 1 1 0 0 0 0 0 0 +1 1 1 0 0 0 0 1 0 +1 1 1 0 0 0 1 0 0 +1 1 1 0 0 0 1 1 0 +1 1 1 0 0 1 0 0 0 +1 1 1 0 0 1 0 1 0 +1 1 1 0 0 1 1 0 0 +1 1 1 0 1 0 0 0 0 +1 1 1 0 1 0 0 1 0 +1 1 1 0 1 0 1 0 0 +1 1 1 0 1 0 1 1 0 +1 1 1 0 1 1 0 0 0 +1 1 1 0 1 1 0 1 0 +1 1 1 0 1 1 1 0 0 +1 1 1 1 0 0 0 0 0 +1 1 1 1 0 0 0 1 0 +1 1 1 1 0 0 1 0 0 +1 1 1 1 0 0 1 1 0 +1 1 1 1 0 1 0 0 0 +1 1 1 1 0 1 0 1 0 +1 1 1 1 0 1 1 0 0 +1 1 1 1 0 1 1 1 0 +1 1 1 1 1 0 0 0 0 +1 1 1 1 1 0 0 1 0 +1 1 1 1 1 0 1 0 0 +1 1 1 1 1 0 1 1 0 +1 1 1 1 1 1 0 0 0 +1 1 1 1 1 1 0 1 0 +1 1 1 1 1 1 1 0 0 +1 1 1 1 1 1 1 1 0 |