From 2b60a12e3d4a3fa43af356e09ab32b82b63b60dc Mon Sep 17 00:00:00 2001 From: Jirka Hladky Date: Wed, 18 Jul 2018 13:56:31 +0200 Subject: Initial commit based on haveged-1.9.2.tar.gz from http://www.issihosts.com/haveged/downloads.html --- nist/Makefile.am | 22 + nist/Makefile.in | 573 ++++++++++++ nist/cephes-protos.h | 180 ++++ nist/cephes.c | 1327 +++++++++++++++++++++++++++ nist/dfft.c | 1381 ++++++++++++++++++++++++++++ nist/matrix.c | 168 ++++ nist/matrix.h | 19 + nist/mconf.h | 169 ++++ nist/nist.c | 57 ++ nist/nist.h | 33 + nist/packtest.c | 2250 ++++++++++++++++++++++++++++++++++++++++++++++ nist/special-functions.c | 32 + nist/special-functions.h | 14 + nist/template9 | 148 +++ 14 files changed, 6373 insertions(+) create mode 100644 nist/Makefile.am create mode 100644 nist/Makefile.in create mode 100644 nist/cephes-protos.h create mode 100644 nist/cephes.c create mode 100644 nist/dfft.c create mode 100644 nist/matrix.c create mode 100644 nist/matrix.h create mode 100644 nist/mconf.h create mode 100644 nist/nist.c create mode 100644 nist/nist.h create mode 100644 nist/packtest.c create mode 100644 nist/special-functions.c create mode 100644 nist/special-functions.h create mode 100644 nist/template9 (limited to 'nist') 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..51fb110 --- /dev/null +++ b/nist/Makefile.in @@ -0,0 +1,573 @@ +# Makefile.in generated by automake 1.15 from Makefile.am. +# @configure_input@ + +# Copyright (C) 1994-2014 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__depfiles_maybe = +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__depfiles_maybe)'; \ + cd $(top_builddir) && $(SHELL) ./config.status $(subdir)/$@ $(am__depfiles_maybe);; \ + esac; + +$(top_builddir)/config.status: $(top_srcdir)/configure $(CONFIG_STATUS_DEPENDENCIES) + cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh + +$(top_srcdir)/configure: $(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: $(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 +#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 + 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 +#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>1; + ipp2=ip; + idp2=ido; + nbd=(ido-1)>>1; + t0=l1*ido; + t10=ip*ido; + + if(ido==1)goto L119; + for(ik=0;ikl1){ + for(j=1;j>1; + np2=n; + + kc=np2; + for(k=1;k>1; + ipp2=ip; + ipph=(ip+1)>>1; + if(idol1)goto L139; + + is= -ido-1; + t1=0; + for(j=1;j>1; + np2=n; + + for(i=2;i +#include +/*#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..0f35ec9 --- /dev/null +++ b/nist/nist.c @@ -0,0 +1,57 @@ +#include +#include +#include +#include +#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 basename[FILENAME_MAX+1]; + +char *GetBaseDir(void) +{ + return basename; +} + +int main(int argc, char **argv) +{ +char *filename = ""; +FILE *fp = stdin; +long result=0; + + if (argc<2 || argc>3) { + printf("Usage sts [