diff --git a/.github/workflows/codeql-analysis.yml b/.github/workflows/codeql-analysis.yml new file mode 100644 index 0000000..34cd7a0 --- /dev/null +++ b/.github/workflows/codeql-analysis.yml @@ -0,0 +1,74 @@ +# For most projects, this workflow file will not need changing; you simply need +# to commit it to your repository. +# +# You may wish to alter this file to override the set of languages analyzed, +# or to provide custom queries or build logic. +# +# ******** NOTE ******** +# We have attempted to detect the languages in your repository. Please check +# the `language` matrix defined below to confirm you have the correct set of +# supported CodeQL languages. +# +name: "CodeQL" + +on: + push: + branches: [ "master" ] + pull_request: + # The branches below must be a subset of the branches above + branches: [ "master" ] + schedule: + - cron: '36 14 4 * 1-5' + +jobs: + analyze: + name: Analyze + runs-on: ubuntu-latest + permissions: + actions: read + contents: read + security-events: write + + strategy: + fail-fast: false + matrix: + language: [ 'cpp' ] + # CodeQL supports [ 'cpp', 'csharp', 'go', 'java', 'javascript', 'python', 'ruby' ] + # Learn more about CodeQL language support at https://aka.ms/codeql-docs/language-support + + steps: + - name: Checkout repository + uses: actions/checkout@v3 + + # Initializes the CodeQL tools for scanning. + - name: Initialize CodeQL + uses: github/codeql-action/init@v2 + with: + languages: ${{ matrix.language }} + # If you wish to specify custom queries, you can do so here or in a config file. + # By default, queries listed here will override any specified in a config file. + # Prefix the list here with "+" to use these queries and those in the config file. + + # Details on CodeQL's query packs refer to : https://docs.github.com/en/code-security/code-scanning/automatically-scanning-your-code-for-vulnerabilities-and-errors/configuring-code-scanning#using-queries-in-ql-packs + # queries: security-extended,security-and-quality + + + # Autobuild attempts to build any compiled languages (C/C++, C#, or Java). + # If this step fails, then you should remove it and run the build manually (see below) + - name: Autobuild + uses: github/codeql-action/autobuild@v2 + + # ℹī¸ Command-line programs to run using the OS shell. + # 📚 See https://docs.github.com/en/actions/using-workflows/workflow-syntax-for-github-actions#jobsjob_idstepsrun + + # If the Autobuild fails above, remove it and uncomment the following three lines. + # modify them (or add more) to build your code if your project, please refer to the EXAMPLE below for guidance. + + # - run: | + # echo "Run, Build Application using script" + # ./location_of_script_within_repo/buildscript.sh + + - name: Perform CodeQL Analysis + uses: github/codeql-action/analyze@v2 + with: + category: "/language:${{matrix.language}}" diff --git a/ArchLinux/PKGBUILD b/ArchLinux/PKGBUILD index ab6c208..d020c0d 100644 --- a/ArchLinux/PKGBUILD +++ b/ArchLinux/PKGBUILD @@ -31,7 +31,7 @@ md5sums=('SKIP') # a description of each element in the source array. pkgver() { - cd "$srcdir/${pkgname%-git}" + cd "$srcdir/${pkgname%-git}" # The examples below are not absolute and need to be adapted to each repo. The # primary goal is to generate version numbers that will increase according to @@ -40,43 +40,43 @@ pkgver() { # are not available, is recommended. # Bazaar -# printf "r%s" "$(bzr revno)" +# printf "r%s" "$(bzr revno)" # Git, tags available -# printf "%s" "$(git describe --long --tags | sed 's/\([^-]*-\)g/r\1/;s/-/./g')" +# printf "%s" "$(git describe --long --tags | sed 's/\([^-]*-\)g/r\1/;s/-/./g')" # Git, no tags available -# printf "r%s.%s" "$(git rev-list --count HEAD)" "$(git rev-parse --short HEAD)" +# printf "r%s.%s" "$(git rev-list --count HEAD)" "$(git rev-parse --short HEAD)" # Mercurial -# printf "r%s.%s" "$(hg identify -n)" "$(hg identify -i)" +# printf "r%s.%s" "$(hg identify -n)" "$(hg identify -i)" # Subversion -# printf "r%s" "$(svnversion | tr -d 'A-z')" +# printf "r%s" "$(svnversion | tr -d 'A-z')" - ( set -o pipefail + ( set -o pipefail git describe --long --tags 2>/dev/null | sed 's/\([^-]*-g\)/r\1/;s/-/./g' || printf "r%s.%s" "$(git rev-list --count HEAD)" "$(git rev-parse --short HEAD)" ) } prepare() { - cd "$srcdir/${pkgname%-git}" - mkdir Release + cd "$srcdir/${pkgname%-git}" + mkdir Release } build() { - cd "$srcdir/${pkgname%-git}/Release" - cmake ../ - make + cd "$srcdir/${pkgname%-git}/Release" + cmake ../ + make } check() { - cd "$srcdir/${pkgname%-git}/Release" - #make -k check + cd "$srcdir/${pkgname%-git}/Release" + #make -k check } package() { - cd "$srcdir/${pkgname%-git}/Release" - make DESTDIR="$pkgdir/" install + cd "$srcdir/${pkgname%-git}/Release" + make DESTDIR="$pkgdir/" install } diff --git a/CMakeLists.txt b/CMakeLists.txt index c07a05b..13f1927 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,4 +1,4 @@ -cmake_minimum_required (VERSION 2.8.4) +cmake_minimum_required (VERSION 3.0.0) include(CMakeDependentOption) # since MacOS 10.14 (XCode 10.0), default includes are no longer installed to /usr/include or /usr/local/include @@ -165,7 +165,7 @@ if(WIN32) Delete \\\"$SYSDIR\\\\dace.dll\\\" ") - unset(CPACK_PACKAGING_INSTALL_PREFIX) # replace by whatever it needs to be + unset(CPACK_PACKAGING_INSTALL_PREFIX) # replace by whatever it needs to be elseif(APPLE) set(CPACK_GENERATOR "TGZ;productbuild") set(CPACK_PACKAGE_VENDOR "DACEDevelopmentGroup") # no spaces allowed for packagebuild diff --git a/Doxyfile b/Doxyfile index 9099c5a..a5fb8b3 100644 --- a/Doxyfile +++ b/Doxyfile @@ -1,4 +1,4 @@ -# Doxyfile 1.8.16 +# Doxyfile 1.9.3 # This file describes the settings to be used by the documentation system # doxygen (www.doxygen.org) for a project. @@ -93,14 +93,6 @@ ALLOW_UNICODE_NAMES = NO OUTPUT_LANGUAGE = English -# The OUTPUT_TEXT_DIRECTION tag is used to specify the direction in which all -# documentation generated by doxygen is written. Doxygen will use this -# information to generate all generated output in the proper direction. -# Possible values are: None, LTR, RTL and Context. -# The default value is: None. - -OUTPUT_TEXT_DIRECTION = None - # If the BRIEF_MEMBER_DESC tag is set to YES, doxygen will include brief member # descriptions after the members that are listed in the file and class # documentation (similar to Javadoc). Set to NO to disable this. @@ -219,6 +211,14 @@ QT_AUTOBRIEF = NO MULTILINE_CPP_IS_BRIEF = NO +# By default Python docstrings are displayed as preformatted text and doxygen's +# special commands cannot be used. By setting PYTHON_DOCSTRING to NO the +# doxygen's special commands can be used and the contents of the docstring +# documentation blocks is shown as doxygen documentation. +# The default value is: YES. + +PYTHON_DOCSTRING = YES + # If the INHERIT_DOCS tag is set to YES then an undocumented member inherits the # documentation from any documented member that it re-implements. # The default value is: YES. @@ -242,25 +242,19 @@ TAB_SIZE = 4 # the documentation. An alias has the form: # name=value # For example adding -# "sideeffect=@par Side Effects:\n" +# "sideeffect=@par Side Effects:^^" # will allow you to put the command \sideeffect (or @sideeffect) in the # documentation, which will result in a user-defined paragraph with heading -# "Side Effects:". You can put \n's in the value part of an alias to insert -# newlines (in the resulting output). You can put ^^ in the value part of an -# alias to insert a newline as if a physical newline was in the original file. -# When you need a literal { or } or , in the value part of an alias you have to -# escape them by means of a backslash (\), this can lead to conflicts with the -# commands \{ and \} for these it is advised to use the version @{ and @} or use -# a double escape (\\{ and \\}) +# "Side Effects:". Note that you cannot put \n's in the value part of an alias +# to insert newlines (in the resulting output). You can put ^^ in the value part +# of an alias to insert a newline as if a physical newline was in the original +# file. When you need a literal { or } or , in the value part of an alias you +# have to escape them by means of a backslash (\), this can lead to conflicts +# with the commands \{ and \} for these it is advised to use the version @{ and +# @} or use a double escape (\\{ and \\}) ALIASES = -# This tag can be used to specify a number of word-keyword mappings (TCL only). -# A mapping has the form "name=value". For example adding "class=itcl::class" -# will allow you to use the command class in the itcl::class meaning. - -TCL_SUBST = - # Set the OPTIMIZE_OUTPUT_FOR_C tag to YES if your project consists of C sources # only. Doxygen will then generate output that is more tailored for C. For # instance, some of the names that are used will be different. The list of all @@ -301,19 +295,22 @@ OPTIMIZE_OUTPUT_SLICE = NO # parses. With this tag you can assign which parser to use for a given # extension. Doxygen has a built-in mapping, but you can override or extend it # using this tag. The format is ext=language, where ext is a file extension, and -# language is one of the parsers supported by doxygen: IDL, Java, Javascript, -# Csharp (C#), C, C++, D, PHP, md (Markdown), Objective-C, Python, Slice, -# Fortran (fixed format Fortran: FortranFixed, free formatted Fortran: +# language is one of the parsers supported by doxygen: IDL, Java, JavaScript, +# Csharp (C#), C, C++, Lex, D, PHP, md (Markdown), Objective-C, Python, Slice, +# VHDL, Fortran (fixed format Fortran: FortranFixed, free formatted Fortran: # FortranFree, unknown formatted Fortran: Fortran. In the later case the parser # tries to guess whether the code is fixed or free formatted code, this is the -# default for Fortran type files), VHDL, tcl. For instance to make doxygen treat -# .inc files as Fortran files (default is PHP), and .f files as C (default is -# Fortran), use: inc=Fortran f=C. +# default for Fortran type files). For instance to make doxygen treat .inc files +# as Fortran files (default is PHP), and .f files as C (default is Fortran), +# use: inc=Fortran f=C. # # Note: For files without extension you can use no_extension as a placeholder. # # Note that for custom extensions you also need to set FILE_PATTERNS otherwise -# the files are not read by doxygen. +# the files are not read by doxygen. When specifying no_extension you should add +# * to the FILE_PATTERNS. +# +# Note see also the list of default file extension mappings. EXTENSION_MAPPING = @@ -447,6 +444,19 @@ TYPEDEF_HIDES_STRUCT = NO LOOKUP_CACHE_SIZE = 0 +# The NUM_PROC_THREADS specifies the number threads doxygen is allowed to use +# during processing. When set to 0 doxygen will based this on the number of +# cores available in the system. You can set it explicitly to a value larger +# than 0 to get more control over the balance between CPU load and processing +# speed. At this moment only the input processing can be done using multiple +# threads. Since this is still an experimental feature the default is set to 1, +# which effectively disables parallel processing. Please report any issues you +# encounter. Generating dot graphs in parallel is controlled by the +# DOT_NUM_THREADS setting. +# Minimum value: 0, maximum value: 32, default value: 1. + +NUM_PROC_THREADS = 1 + #--------------------------------------------------------------------------- # Build related configuration options #--------------------------------------------------------------------------- @@ -510,6 +520,13 @@ EXTRACT_LOCAL_METHODS = NO EXTRACT_ANON_NSPACES = NO +# If this flag is set to YES, the name of an unnamed parameter in a declaration +# will be determined by the corresponding definition. By default unnamed +# parameters remain unnamed in the output. +# The default value is: YES. + +RESOLVE_UNNAMED_PARAMS = YES + # If the HIDE_UNDOC_MEMBERS tag is set to YES, doxygen will hide all # undocumented members inside documented classes or files. If set to NO these # members will be included in the various overviews, but no documentation @@ -527,8 +544,8 @@ HIDE_UNDOC_MEMBERS = NO HIDE_UNDOC_CLASSES = NO # If the HIDE_FRIEND_COMPOUNDS tag is set to YES, doxygen will hide all friend -# (class|struct|union) declarations. If set to NO, these declarations will be -# included in the documentation. +# declarations. If set to NO, these declarations will be included in the +# documentation. # The default value is: NO. HIDE_FRIEND_COMPOUNDS = NO @@ -547,11 +564,18 @@ HIDE_IN_BODY_DOCS = NO INTERNAL_DOCS = NO -# If the CASE_SENSE_NAMES tag is set to NO then doxygen will only generate file -# names in lower-case letters. If set to YES, upper-case letters are also -# allowed. This is useful if you have classes or files whose names only differ -# in case and if your file system supports case sensitive file names. Windows -# (including Cygwin) ands Mac users are advised to set this option to NO. +# With the correct setting of option CASE_SENSE_NAMES doxygen will better be +# able to match the capabilities of the underlying filesystem. In case the +# filesystem is case sensitive (i.e. it supports files in the same directory +# whose names only differ in casing), the option must be set to YES to properly +# deal with such files in case they appear in the input. For filesystems that +# are not case sensitive the option should be be set to NO to properly deal with +# output files written for symbols that only differ in casing, such as for two +# classes, one named CLASS and the other named Class, and to also support +# references to files without having to specify the exact matching casing. On +# Windows (including Cygwin) and MacOS, users should typically set this option +# to NO, whereas on Linux or other Unix flavors it should typically be set to +# YES. # The default value is: system dependent. CASE_SENSE_NAMES = YES @@ -570,6 +594,12 @@ HIDE_SCOPE_NAMES = NO HIDE_COMPOUND_REFERENCE= NO +# If the SHOW_HEADERFILE tag is set to YES then the documentation for a class +# will show which file needs to be included to use the class. +# The default value is: YES. + +SHOW_HEADERFILE = YES + # If the SHOW_INCLUDE_FILES tag is set to YES then doxygen will put a list of # the files that are included by a file in the documentation of that file. # The default value is: YES. @@ -727,7 +757,8 @@ FILE_VERSION_FILTER = # output files in an output format independent way. To create the layout file # that represents doxygen's defaults, run doxygen with the -l option. You can # optionally specify a file name after the option, if omitted DoxygenLayout.xml -# will be used as the name of the layout file. +# will be used as the name of the layout file. See also section "Changing the +# layout of pages" for information. # # Note that if you run doxygen from a directory containing a file called # DoxygenLayout.xml, doxygen will parse it automatically even if the LAYOUT_FILE @@ -773,24 +804,35 @@ WARNINGS = YES WARN_IF_UNDOCUMENTED = YES # If the WARN_IF_DOC_ERROR tag is set to YES, doxygen will generate warnings for -# potential errors in the documentation, such as not documenting some parameters -# in a documented function, or documenting parameters that don't exist or using -# markup commands wrongly. +# potential errors in the documentation, such as documenting some parameters in +# a documented function twice, or documenting parameters that don't exist or +# using markup commands wrongly. # The default value is: YES. WARN_IF_DOC_ERROR = YES +# If WARN_IF_INCOMPLETE_DOC is set to YES, doxygen will warn about incomplete +# function parameter documentation. If set to NO, doxygen will accept that some +# parameters have no documentation without warning. +# The default value is: YES. + +WARN_IF_INCOMPLETE_DOC = YES + # This WARN_NO_PARAMDOC option can be enabled to get warnings for functions that # are documented, but have no documentation for their parameters or return -# value. If set to NO, doxygen will only warn about wrong or incomplete -# parameter documentation, but not about the absence of documentation. If -# EXTRACT_ALL is set to YES then this flag will automatically be disabled. +# value. If set to NO, doxygen will only warn about wrong parameter +# documentation, but not about the absence of documentation. If EXTRACT_ALL is +# set to YES then this flag will automatically be disabled. See also +# WARN_IF_INCOMPLETE_DOC # The default value is: NO. WARN_NO_PARAMDOC = NO # If the WARN_AS_ERROR tag is set to YES then doxygen will immediately stop when -# a warning is encountered. +# a warning is encountered. If the WARN_AS_ERROR tag is set to FAIL_ON_WARNINGS +# then doxygen will continue running as if WARN_AS_ERROR tag is set to NO, but +# at the end of the doxygen process doxygen will return with a non-zero status. +# Possible values are: NO, YES and FAIL_ON_WARNINGS. # The default value is: NO. WARN_AS_ERROR = NO @@ -807,7 +849,10 @@ WARN_FORMAT = "$file:$line: $text" # The WARN_LOGFILE tag can be used to specify a file to which warning and error # messages should be written. If left blank the output is written to standard -# error (stderr). +# error (stderr). In case the file specified cannot be opened for writing the +# warning and error messages are written to standard error. When as file - is +# specified the warning and error messages are written to standard output +# (stdout). WARN_LOGFILE = @@ -827,8 +872,8 @@ INPUT = interfaces/cxx \ # This tag can be used to specify the character encoding of the source files # that doxygen parses. Internally doxygen uses the UTF-8 encoding. Doxygen uses # libiconv (or the iconv built into libc) for the transcoding. See the libiconv -# documentation (see: https://www.gnu.org/software/libiconv/) for the list of -# possible encodings. +# documentation (see: +# https://www.gnu.org/software/libiconv/) for the list of possible encodings. # The default value is: UTF-8. INPUT_ENCODING = UTF-8 @@ -841,13 +886,20 @@ INPUT_ENCODING = UTF-8 # need to set EXTENSION_MAPPING for the extension otherwise the files are not # read by doxygen. # +# Note the list of default checked file patterns might differ from the list of +# default file extension mappings. +# # If left blank the following patterns are tested:*.c, *.cc, *.cxx, *.cpp, # *.c++, *.java, *.ii, *.ixx, *.ipp, *.i++, *.inl, *.idl, *.ddl, *.odl, *.h, -# *.hh, *.hxx, *.hpp, *.h++, *.cs, *.d, *.php, *.php4, *.php5, *.phtml, *.inc, -# *.m, *.markdown, *.md, *.mm, *.dox, *.py, *.pyw, *.f90, *.f95, *.f03, *.f08, -# *.f, *.for, *.tcl, *.vhd, *.vhdl, *.ucf, *.qsf and *.ice. +# *.hh, *.hxx, *.hpp, *.h++, *.l, *.cs, *.d, *.php, *.php4, *.php5, *.phtml, +# *.inc, *.m, *.markdown, *.md, *.mm, *.dox (to be provided as doxygen C +# comment), *.py, *.pyw, *.f90, *.f95, *.f03, *.f08, *.f18, *.f, *.for, *.vhd, +# *.vhdl, *.ucf, *.qsf and *.ice. -FILE_PATTERNS = +FILE_PATTERNS = *.cpp \ + *.h \ + *.c \ + *.f # The RECURSIVE tag can be used to specify whether or not subdirectories should # be searched for input files as well. @@ -884,7 +936,7 @@ EXCLUDE_PATTERNS = # (namespaces, classes, functions, etc.) that should be excluded from the # output. The symbol name can be a fully qualified name, a word, or if the # wildcard * is used, a substring. Examples: ANamespace, AClass, -# AClass::ANamespace, ANamespace::*Test +# ANamespace::AClass, ANamespace::*Test # # Note that the wildcards are matched against the file with absolute path, so to # exclude all test directories use the pattern */test/* @@ -1070,13 +1122,6 @@ VERBATIM_HEADERS = YES ALPHABETICAL_INDEX = YES -# The COLS_IN_ALPHA_INDEX tag can be used to specify the number of columns in -# which the alphabetical index list will be split. -# Minimum value: 1, maximum value: 20, default value: 5. -# This tag requires that the tag ALPHABETICAL_INDEX is set to YES. - -COLS_IN_ALPHA_INDEX = 5 - # In case all classes in a project start with a common prefix, all classes will # be put under the same header in the alphabetical index. The IGNORE_PREFIX tag # can be used to specify a prefix (or a list of prefixes) that should be ignored @@ -1176,7 +1221,7 @@ HTML_EXTRA_FILES = # The HTML_COLORSTYLE_HUE tag controls the color of the HTML output. Doxygen # will adjust the colors in the style sheet and background images according to -# this color. Hue is specified as an angle on a colorwheel, see +# this color. Hue is specified as an angle on a color-wheel, see # https://en.wikipedia.org/wiki/Hue for more information. For instance the value # 0 represents red, 60 is yellow, 120 is green, 180 is cyan, 240 is blue, 300 # purple, and 360 is red again. @@ -1186,7 +1231,7 @@ HTML_EXTRA_FILES = HTML_COLORSTYLE_HUE = 220 # The HTML_COLORSTYLE_SAT tag controls the purity (or saturation) of the colors -# in the HTML output. For a value of 0 the output will use grayscales only. A +# in the HTML output. For a value of 0 the output will use gray-scales only. A # value of 255 will produce the most vivid colors. # Minimum value: 0, maximum value: 255, default value: 100. # This tag requires that the tag GENERATE_HTML is set to YES. @@ -1215,9 +1260,9 @@ HTML_TIMESTAMP = YES # If the HTML_DYNAMIC_MENUS tag is set to YES then the generated HTML # documentation will contain a main index with vertical navigation menus that -# are dynamically created via Javascript. If disabled, the navigation index will +# are dynamically created via JavaScript. If disabled, the navigation index will # consists of multiple levels of tabs that are statically embedded in every HTML -# page. Disable this option to support browsers that do not have Javascript, +# page. Disable this option to support browsers that do not have JavaScript, # like the Qt help browser. # The default value is: YES. # This tag requires that the tag GENERATE_HTML is set to YES. @@ -1247,10 +1292,11 @@ HTML_INDEX_NUM_ENTRIES = 100 # If the GENERATE_DOCSET tag is set to YES, additional index files will be # generated that can be used as input for Apple's Xcode 3 integrated development -# environment (see: https://developer.apple.com/xcode/), introduced with OSX -# 10.5 (Leopard). To create a documentation set, doxygen will generate a -# Makefile in the HTML output directory. Running make will produce the docset in -# that directory and running make install will install the docset in +# environment (see: +# https://developer.apple.com/xcode/), introduced with OSX 10.5 (Leopard). To +# create a documentation set, doxygen will generate a Makefile in the HTML +# output directory. Running make will produce the docset in that directory and +# running make install will install the docset in # ~/Library/Developer/Shared/Documentation/DocSets so that Xcode will find it at # startup. See https://developer.apple.com/library/archive/featuredarticles/Doxy # genXcode/_index.html for more information. @@ -1267,6 +1313,13 @@ GENERATE_DOCSET = YES DOCSET_FEEDNAME = DACE +# This tag determines the URL of the docset feed. A documentation feed provides +# an umbrella under which multiple documentation sets from a single provider +# (such as a company or product suite) can be grouped. +# This tag requires that the tag GENERATE_DOCSET is set to YES. + +DOCSET_FEEDURL = + # This tag specifies a string that should uniquely identify the documentation # set bundle. This should be a reverse domain-name style string, e.g. # com.mycompany.MyDocSet. Doxygen will append .docset to the name. @@ -1292,8 +1345,12 @@ DOCSET_PUBLISHER_NAME = Team # If the GENERATE_HTMLHELP tag is set to YES then doxygen generates three # additional HTML index files: index.hhp, index.hhc, and index.hhk. The # index.hhp is a project file that can be read by Microsoft's HTML Help Workshop -# (see: https://www.microsoft.com/en-us/download/details.aspx?id=21138) on -# Windows. +# on Windows. In the beginning of 2021 Microsoft took the original page, with +# a.o. the download links, offline the HTML help workshop was already many years +# in maintenance mode). You can download the HTML help workshop from the web +# archives at Installation executable (see: +# http://web.archive.org/web/20160201063255/http://download.microsoft.com/downlo +# ad/0/A/9/0A939EF6-E31C-430F-A3DF-DFAE7960D564/htmlhelp.exe). # # The HTML Help Workshop contains a compiler that can convert all HTML output # generated by doxygen into a single compiled HTML file (.chm). Compiled HTML @@ -1323,7 +1380,7 @@ CHM_FILE = DACE.chm HHC_LOCATION = # The GENERATE_CHI flag controls if a separate .chi index file is generated -# (YES) or that it should be included in the master .chm file (NO). +# (YES) or that it should be included in the main .chm file (NO). # The default value is: NO. # This tag requires that the tag GENERATE_HTMLHELP is set to YES. @@ -1368,7 +1425,8 @@ QCH_FILE = # The QHP_NAMESPACE tag specifies the namespace to use when generating Qt Help # Project output. For more information please see Qt Help Project / Namespace -# (see: https://doc.qt.io/archives/qt-4.8/qthelpproject.html#namespace). +# (see: +# https://doc.qt.io/archives/qt-4.8/qthelpproject.html#namespace). # The default value is: org.doxygen.Project. # This tag requires that the tag GENERATE_QHP is set to YES. @@ -1376,8 +1434,8 @@ QHP_NAMESPACE = org.doxygen.Project # The QHP_VIRTUAL_FOLDER tag specifies the namespace to use when generating Qt # Help Project output. For more information please see Qt Help Project / Virtual -# Folders (see: https://doc.qt.io/archives/qt-4.8/qthelpproject.html#virtual- -# folders). +# Folders (see: +# https://doc.qt.io/archives/qt-4.8/qthelpproject.html#virtual-folders). # The default value is: doc. # This tag requires that the tag GENERATE_QHP is set to YES. @@ -1385,16 +1443,16 @@ QHP_VIRTUAL_FOLDER = doc # If the QHP_CUST_FILTER_NAME tag is set, it specifies the name of a custom # filter to add. For more information please see Qt Help Project / Custom -# Filters (see: https://doc.qt.io/archives/qt-4.8/qthelpproject.html#custom- -# filters). +# Filters (see: +# https://doc.qt.io/archives/qt-4.8/qthelpproject.html#custom-filters). # This tag requires that the tag GENERATE_QHP is set to YES. QHP_CUST_FILTER_NAME = # The QHP_CUST_FILTER_ATTRS tag specifies the list of the attributes of the # custom filter to add. For more information please see Qt Help Project / Custom -# Filters (see: https://doc.qt.io/archives/qt-4.8/qthelpproject.html#custom- -# filters). +# Filters (see: +# https://doc.qt.io/archives/qt-4.8/qthelpproject.html#custom-filters). # This tag requires that the tag GENERATE_QHP is set to YES. QHP_CUST_FILTER_ATTRS = @@ -1406,9 +1464,9 @@ QHP_CUST_FILTER_ATTRS = QHP_SECT_FILTER_ATTRS = -# The QHG_LOCATION tag can be used to specify the location of Qt's -# qhelpgenerator. If non-empty doxygen will try to run qhelpgenerator on the -# generated .qhp file. +# The QHG_LOCATION tag can be used to specify the location (absolute path +# including file name) of Qt's qhelpgenerator. If non-empty doxygen will try to +# run qhelpgenerator on the generated .qhp file. # This tag requires that the tag GENERATE_QHP is set to YES. QHG_LOCATION = @@ -1451,16 +1509,28 @@ DISABLE_INDEX = NO # to work a browser that supports JavaScript, DHTML, CSS and frames is required # (i.e. any modern browser). Windows users are probably better off using the # HTML help feature. Via custom style sheets (see HTML_EXTRA_STYLESHEET) one can -# further fine-tune the look of the index. As an example, the default style -# sheet generated by doxygen has an example that shows how to put an image at -# the root of the tree instead of the PROJECT_NAME. Since the tree basically has -# the same information as the tab index, you could consider setting -# DISABLE_INDEX to YES when enabling this option. +# further fine tune the look of the index (see "Fine-tuning the output"). As an +# example, the default style sheet generated by doxygen has an example that +# shows how to put an image at the root of the tree instead of the PROJECT_NAME. +# Since the tree basically has the same information as the tab index, you could +# consider setting DISABLE_INDEX to YES when enabling this option. # The default value is: NO. # This tag requires that the tag GENERATE_HTML is set to YES. GENERATE_TREEVIEW = NO +# When both GENERATE_TREEVIEW and DISABLE_INDEX are set to YES, then the +# FULL_SIDEBAR option determines if the side bar is limited to only the treeview +# area (value NO) or if it should extend to the full height of the window (value +# YES). Setting this to YES gives a layout similar to +# https://docs.readthedocs.io with more room for contents, but less room for the +# project logo, title, and description. If either GENERATE_TREEVIEW or +# DISABLE_INDEX is set to NO, this option has no effect. +# The default value is: NO. +# This tag requires that the tag GENERATE_HTML is set to YES. + +FULL_SIDEBAR = NO + # The ENUM_VALUES_PER_LINE tag can be used to set the number of enum values that # doxygen will group on one line in the generated HTML documentation. # @@ -1485,6 +1555,24 @@ TREEVIEW_WIDTH = 250 EXT_LINKS_IN_WINDOW = NO +# If the OBFUSCATE_EMAILS tag is set to YES, doxygen will obfuscate email +# addresses. +# The default value is: YES. +# This tag requires that the tag GENERATE_HTML is set to YES. + +OBFUSCATE_EMAILS = YES + +# If the HTML_FORMULA_FORMAT option is set to svg, doxygen will use the pdf2svg +# tool (see https://github.com/dawbarton/pdf2svg) or inkscape (see +# https://inkscape.org) to generate formulas as SVG images instead of PNGs for +# the HTML output. These images will generally look nicer at scaled resolutions. +# Possible values are: png (the default) and svg (looks nicer but requires the +# pdf2svg or inkscape tool). +# The default value is: png. +# This tag requires that the tag GENERATE_HTML is set to YES. + +HTML_FORMULA_FORMAT = png + # Use this tag to change the font size of LaTeX formulas included as images in # the HTML documentation. When you change the font size after a successful # doxygen run you need to manually remove any form_*.png images from the HTML @@ -1505,8 +1593,14 @@ FORMULA_FONTSIZE = 10 FORMULA_TRANSPARENT = YES +# The FORMULA_MACROFILE can contain LaTeX \newcommand and \renewcommand commands +# to create new LaTeX commands to be used in formulas as building blocks. See +# the section "Including formulas" for details. + +FORMULA_MACROFILE = + # Enable the USE_MATHJAX option to render LaTeX formulas using MathJax (see -# https://www.mathjax.org) which uses client side Javascript for the rendering +# https://www.mathjax.org) which uses client side JavaScript for the rendering # instead of using pre-rendered bitmaps. Use this if you do not have LaTeX # installed or if you want to formulas look prettier in the HTML output. When # enabled you may also need to install MathJax separately and configure the path @@ -1516,11 +1610,29 @@ FORMULA_TRANSPARENT = YES USE_MATHJAX = NO +# With MATHJAX_VERSION it is possible to specify the MathJax version to be used. +# Note that the different versions of MathJax have different requirements with +# regards to the different settings, so it is possible that also other MathJax +# settings have to be changed when switching between the different MathJax +# versions. +# Possible values are: MathJax_2 and MathJax_3. +# The default value is: MathJax_2. +# This tag requires that the tag USE_MATHJAX is set to YES. + +MATHJAX_VERSION = MathJax_2 + # When MathJax is enabled you can set the default output format to be used for -# the MathJax output. See the MathJax site (see: -# http://docs.mathjax.org/en/latest/output.html) for more details. +# the MathJax output. For more details about the output format see MathJax +# version 2 (see: +# http://docs.mathjax.org/en/v2.7-latest/output.html) and MathJax version 3 +# (see: +# http://docs.mathjax.org/en/latest/web/components/output.html). # Possible values are: HTML-CSS (which is slower, but has the best -# compatibility), NativeMML (i.e. MathML) and SVG. +# compatibility. This is the name for Mathjax version 2, for MathJax version 3 +# this will be translated into chtml), NativeMML (i.e. MathML. Only supported +# for NathJax 2. For MathJax version 3 chtml will be used instead.), chtml (This +# is the name for Mathjax version 3, for MathJax version 2 this will be +# translated into HTML-CSS) and SVG. # The default value is: HTML-CSS. # This tag requires that the tag USE_MATHJAX is set to YES. @@ -1533,22 +1645,29 @@ MATHJAX_FORMAT = HTML-CSS # MATHJAX_RELPATH should be ../mathjax. The default value points to the MathJax # Content Delivery Network so you can quickly see the result without installing # MathJax. However, it is strongly recommended to install a local copy of -# MathJax from https://www.mathjax.org before deployment. -# The default value is: https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.5/. +# MathJax from https://www.mathjax.org before deployment. The default value is: +# - in case of MathJax version 2: https://cdn.jsdelivr.net/npm/mathjax@2 +# - in case of MathJax version 3: https://cdn.jsdelivr.net/npm/mathjax@3 # This tag requires that the tag USE_MATHJAX is set to YES. MATHJAX_RELPATH = http://cdn.mathjax.org/mathjax/latest # The MATHJAX_EXTENSIONS tag can be used to specify one or more MathJax # extension names that should be enabled during MathJax rendering. For example +# for MathJax version 2 (see +# https://docs.mathjax.org/en/v2.7-latest/tex.html#tex-and-latex-extensions): # MATHJAX_EXTENSIONS = TeX/AMSmath TeX/AMSsymbols +# For example for MathJax version 3 (see +# http://docs.mathjax.org/en/latest/input/tex/extensions/index.html): +# MATHJAX_EXTENSIONS = ams # This tag requires that the tag USE_MATHJAX is set to YES. MATHJAX_EXTENSIONS = # The MATHJAX_CODEFILE tag can be used to specify a file with javascript pieces # of code that will be used on startup of the MathJax code. See the MathJax site -# (see: http://docs.mathjax.org/en/latest/output.html) for more details. For an +# (see: +# http://docs.mathjax.org/en/v2.7-latest/output.html) for more details. For an # example see the documentation. # This tag requires that the tag USE_MATHJAX is set to YES. @@ -1576,7 +1695,7 @@ MATHJAX_CODEFILE = SEARCHENGINE = YES # When the SERVER_BASED_SEARCH tag is enabled the search engine will be -# implemented using a web server instead of a web client using Javascript. There +# implemented using a web server instead of a web client using JavaScript. There # are two flavors of web server based searching depending on the EXTERNAL_SEARCH # setting. When disabled, doxygen will generate a PHP script for searching and # an index file used by the script. When EXTERNAL_SEARCH is enabled the indexing @@ -1595,7 +1714,8 @@ SERVER_BASED_SEARCH = NO # # Doxygen ships with an example indexer (doxyindexer) and search engine # (doxysearch.cgi) which are based on the open source search engine library -# Xapian (see: https://xapian.org/). +# Xapian (see: +# https://xapian.org/). # # See the section "External Indexing and Searching" for details. # The default value is: NO. @@ -1608,8 +1728,9 @@ EXTERNAL_SEARCH = NO # # Doxygen ships with an example indexer (doxyindexer) and search engine # (doxysearch.cgi) which are based on the open source search engine library -# Xapian (see: https://xapian.org/). See the section "External Indexing and -# Searching" for details. +# Xapian (see: +# https://xapian.org/). See the section "External Indexing and Searching" for +# details. # This tag requires that the tag SEARCHENGINE is set to YES. SEARCHENGINE_URL = @@ -1718,29 +1839,31 @@ PAPER_TYPE = a4 EXTRA_PACKAGES = float -# The LATEX_HEADER tag can be used to specify a personal LaTeX header for the -# generated LaTeX document. The header should contain everything until the first -# chapter. If it is left blank doxygen will generate a standard header. See -# section "Doxygen usage" for information on how to let doxygen write the -# default header to a separate file. +# The LATEX_HEADER tag can be used to specify a user-defined LaTeX header for +# the generated LaTeX document. The header should contain everything until the +# first chapter. If it is left blank doxygen will generate a standard header. It +# is highly recommended to start with a default header using +# doxygen -w latex new_header.tex new_footer.tex new_stylesheet.sty +# and then modify the file new_header.tex. See also section "Doxygen usage" for +# information on how to generate the default header that doxygen normally uses. # -# Note: Only use a user-defined header if you know what you are doing! The -# following commands have a special meaning inside the header: $title, -# $datetime, $date, $doxygenversion, $projectname, $projectnumber, -# $projectbrief, $projectlogo. Doxygen will replace $title with the empty -# string, for the replacement values of the other commands the user is referred -# to HTML_HEADER. +# Note: Only use a user-defined header if you know what you are doing! +# Note: The header is subject to change so you typically have to regenerate the +# default header when upgrading to a newer version of doxygen. The following +# commands have a special meaning inside the header (and footer): For a +# description of the possible markers and block names see the documentation. # This tag requires that the tag GENERATE_LATEX is set to YES. LATEX_HEADER = -# The LATEX_FOOTER tag can be used to specify a personal LaTeX footer for the -# generated LaTeX document. The footer should contain everything after the last -# chapter. If it is left blank doxygen will generate a standard footer. See +# The LATEX_FOOTER tag can be used to specify a user-defined LaTeX footer for +# the generated LaTeX document. The footer should contain everything after the +# last chapter. If it is left blank doxygen will generate a standard footer. See # LATEX_HEADER for more information on how to generate a default footer and what -# special commands can be used inside the footer. -# -# Note: Only use a user-defined footer if you know what you are doing! +# special commands can be used inside the footer. See also section "Doxygen +# usage" for information on how to generate the default footer that doxygen +# normally uses. Note: Only use a user-defined footer if you know what you are +# doing! # This tag requires that the tag GENERATE_LATEX is set to YES. LATEX_FOOTER = @@ -1773,9 +1896,11 @@ LATEX_EXTRA_FILES = PDF_HYPERLINKS = YES -# If the USE_PDFLATEX tag is set to YES, doxygen will use pdflatex to generate -# the PDF file directly from the LaTeX files. Set this option to YES, to get a -# higher quality PDF documentation. +# If the USE_PDFLATEX tag is set to YES, doxygen will use the engine as +# specified with LATEX_CMD_NAME to generate the PDF file directly from the LaTeX +# files. Set this option to YES, to get a higher quality PDF documentation. +# +# See also section LATEX_CMD_NAME for selecting the engine. # The default value is: YES. # This tag requires that the tag GENERATE_LATEX is set to YES. @@ -1783,8 +1908,7 @@ USE_PDFLATEX = YES # If the LATEX_BATCHMODE tag is set to YES, doxygen will add the \batchmode # command to the generated LaTeX files. This will instruct LaTeX to keep running -# if errors occur, instead of asking the user for help. This option is also used -# when generating formulas in HTML. +# if errors occur, instead of asking the user for help. # The default value is: NO. # This tag requires that the tag GENERATE_LATEX is set to YES. @@ -1797,16 +1921,6 @@ LATEX_BATCHMODE = NO LATEX_HIDE_INDICES = YES -# If the LATEX_SOURCE_CODE tag is set to YES then doxygen will include source -# code with syntax highlighting in the LaTeX output. -# -# Note that which sources are shown also depends on other settings such as -# SOURCE_BROWSER. -# The default value is: NO. -# This tag requires that the tag GENERATE_LATEX is set to YES. - -LATEX_SOURCE_CODE = NO - # The LATEX_BIB_STYLE tag can be used to specify the style to use for the # bibliography, e.g. plainnat, or ieeetr. See # https://en.wikipedia.org/wiki/BibTeX and \cite for more info. @@ -1887,16 +2001,6 @@ RTF_STYLESHEET_FILE = RTF_EXTENSIONS_FILE = -# If the RTF_SOURCE_CODE tag is set to YES then doxygen will include source code -# with syntax highlighting in the RTF output. -# -# Note that which sources are shown also depends on other settings such as -# SOURCE_BROWSER. -# The default value is: NO. -# This tag requires that the tag GENERATE_RTF is set to YES. - -RTF_SOURCE_CODE = NO - #--------------------------------------------------------------------------- # Configuration options related to the man page output #--------------------------------------------------------------------------- @@ -1993,15 +2097,6 @@ GENERATE_DOCBOOK = NO DOCBOOK_OUTPUT = docbook -# If the DOCBOOK_PROGRAMLISTING tag is set to YES, doxygen will include the -# program listings (including syntax highlighting and cross-referencing -# information) to the DOCBOOK output. Note that enabling this will significantly -# increase the size of the DOCBOOK output. -# The default value is: NO. -# This tag requires that the tag GENERATE_DOCBOOK is set to YES. - -DOCBOOK_PROGRAMLISTING = NO - #--------------------------------------------------------------------------- # Configuration options for the AutoGen Definitions output #--------------------------------------------------------------------------- @@ -2180,15 +2275,6 @@ EXTERNAL_PAGES = YES # Configuration options related to the dot tool #--------------------------------------------------------------------------- -# If the CLASS_DIAGRAMS tag is set to YES, doxygen will generate a class diagram -# (in HTML and LaTeX) for classes with base or super classes. Setting the tag to -# NO turns the diagrams off. Note that this option also works with HAVE_DOT -# disabled, but it is recommended to install and use dot, since it yields more -# powerful graphs. -# The default value is: YES. - -CLASS_DIAGRAMS = YES - # You can include diagrams made with dia in doxygen documentation. Doxygen will # then run dia to produce the diagram and insert it in the documentation. The # DIA_PATH tag allows you to specify the directory where the dia binary resides. @@ -2245,11 +2331,14 @@ DOT_FONTSIZE = 10 DOT_FONTPATH = -# If the CLASS_GRAPH tag is set to YES then doxygen will generate a graph for -# each documented class showing the direct and indirect inheritance relations. -# Setting this tag to YES will force the CLASS_DIAGRAMS tag to NO. +# If the CLASS_GRAPH tag is set to YES (or GRAPH) then doxygen will generate a +# graph for each documented class showing the direct and indirect inheritance +# relations. In case HAVE_DOT is set as well dot will be used to draw the graph, +# otherwise the built-in generator will be used. If the CLASS_GRAPH tag is set +# to TEXT the direct and indirect inheritance relations will be shown as texts / +# links. +# Possible values are: NO, YES, TEXT and GRAPH. # The default value is: YES. -# This tag requires that the tag HAVE_DOT is set to YES. CLASS_GRAPH = YES @@ -2286,10 +2375,32 @@ UML_LOOK = NO # but if the number exceeds 15, the total amount of fields shown is limited to # 10. # Minimum value: 0, maximum value: 100, default value: 10. -# This tag requires that the tag HAVE_DOT is set to YES. +# This tag requires that the tag UML_LOOK is set to YES. UML_LIMIT_NUM_FIELDS = 10 +# If the DOT_UML_DETAILS tag is set to NO, doxygen will show attributes and +# methods without types and arguments in the UML graphs. If the DOT_UML_DETAILS +# tag is set to YES, doxygen will add type and arguments for attributes and +# methods in the UML graphs. If the DOT_UML_DETAILS tag is set to NONE, doxygen +# will not generate fields with class member information in the UML graphs. The +# class diagrams will look similar to the default class diagrams but using UML +# notation for the relationships. +# Possible values are: NO, YES and NONE. +# The default value is: NO. +# This tag requires that the tag UML_LOOK is set to YES. + +DOT_UML_DETAILS = NO + +# The DOT_WRAP_THRESHOLD tag can be used to set the maximum number of characters +# to display on a single line. If the actual line length exceeds this threshold +# significantly it will wrapped across multiple lines. Some heuristics are apply +# to avoid ugly line breaks. +# Minimum value: 0, maximum value: 1000, default value: 17. +# This tag requires that the tag HAVE_DOT is set to YES. + +DOT_WRAP_THRESHOLD = 17 + # If the TEMPLATE_RELATIONS tag is set to YES then the inheritance and # collaboration graphs will show the relations between templates and their # instances. @@ -2356,6 +2467,13 @@ GRAPHICAL_HIERARCHY = YES DIRECTORY_GRAPH = YES +# The DIR_GRAPH_MAX_DEPTH tag can be used to limit the maximum number of levels +# of child directories generated in directory dependency graphs by dot. +# Minimum value: 1, maximum value: 25, default value: 1. +# This tag requires that the tag DIRECTORY_GRAPH is set to YES. + +DIR_GRAPH_MAX_DEPTH = 1 + # The DOT_IMAGE_FORMAT tag can be used to set the image format of the images # generated by dot. For an explanation of the image formats see the section # output formats in the documentation of the dot tool (Graphviz (see: @@ -2409,10 +2527,10 @@ MSCFILE_DIRS = DIAFILE_DIRS = # When using plantuml, the PLANTUML_JAR_PATH tag should be used to specify the -# path where java can find the plantuml.jar file. If left blank, it is assumed -# PlantUML is not used or called during a preprocessing step. Doxygen will -# generate a warning when it encounters a \startuml command in this case and -# will not generate output for the diagram. +# path where java can find the plantuml.jar file or to the filename of jar file +# to be used. If left blank, it is assumed PlantUML is not used or called during +# a preprocessing step. Doxygen will generate a warning when it encounters a +# \startuml command in this case and will not generate output for the diagram. PLANTUML_JAR_PATH = @@ -2474,14 +2592,18 @@ DOT_MULTI_TARGETS = YES # If the GENERATE_LEGEND tag is set to YES doxygen will generate a legend page # explaining the meaning of the various boxes and arrows in the dot generated # graphs. +# Note: This tag requires that UML_LOOK isn't set, i.e. the doxygen internal +# graphical representation for inheritance and collaboration diagrams is used. # The default value is: YES. # This tag requires that the tag HAVE_DOT is set to YES. GENERATE_LEGEND = YES -# If the DOT_CLEANUP tag is set to YES, doxygen will remove the intermediate dot +# If the DOT_CLEANUP tag is set to YES, doxygen will remove the intermediate # files that are used to generate the various graphs. +# +# Note: This setting is not only used for dot files but also for msc temporary +# files. # The default value is: YES. -# This tag requires that the tag HAVE_DOT is set to YES. DOT_CLEANUP = YES diff --git a/Tutorials/Tutorial2/1Basics-Ex.cpp b/Tutorials/Tutorial2/1Basics-Ex.cpp index d5650ec..ec6a766 100644 --- a/Tutorials/Tutorial2/1Basics-Ex.cpp +++ b/Tutorials/Tutorial2/1Basics-Ex.cpp @@ -70,5 +70,5 @@ int main( void ) cout << "Exercise 1.2.2: Breaking bad\n" << ex1_2_2( 0.0, 0.0 ) << endl << ex1_2_2( x, y ) << endl; - } catch( exception &ex ) { cout << ex.what( ) << endl; } + } catch( exception const& ex ) { cout << ex.what( ) << endl; } } \ No newline at end of file diff --git a/core/contrib/dgamma.c b/core/contrib/dgamma.c index c98ea08..5ea4303 100644 --- a/core/contrib/dgamma.c +++ b/core/contrib/dgamma.c @@ -1,16 +1,16 @@ /* dgamma.f -- translated by f2c (version 20100827). You must link the resulting object file with libf2c: - on Microsoft Windows system, link with libf2c.lib; - on Linux or Unix systems, link with .../path/to/libf2c.a -lm - or, if you install libf2c.a in a standard place, with -lf2c -lm - -- in that order, at the end of the command line, as in - cc *.o -lf2c -lm - Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., - - http://www.netlib.org/f2c/libf2c.zip + on Microsoft Windows system, link with libf2c.lib; + on Linux or Unix systems, link with .../path/to/libf2c.a -lm + or, if you install libf2c.a in a standard place, with -lf2c -lm + -- in that order, at the end of the command line, as in + cc *.o -lf2c -lm + Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., + + http://www.netlib.org/f2c/libf2c.zip */ -/** \addtogroup DACE Contrib +/** \addtogroup DACEContrib Contrib * @{ */ @@ -27,19 +27,19 @@ doublereal dgamma_(const doublereal *x) static doublereal eps = 2.22e-16; static doublereal xinf = 1.79e308; static doublereal p[8] = { -1.71618513886549492533811, - 24.7656508055759199108314,-379.804256470945635097577, - 629.331155312818442661052,866.966202790413211295064, - -31451.2729688483675254357,-36144.4134186911729807069, - 66456.1438202405440627855 }; + 24.7656508055759199108314,-379.804256470945635097577, + 629.331155312818442661052,866.966202790413211295064, + -31451.2729688483675254357,-36144.4134186911729807069, + 66456.1438202405440627855 }; static doublereal q[8] = { -30.8402300119738975254353, - 315.350626979604161529144,-1015.15636749021914166146, - -3107.77167157231109440444,22538.1184209801510330112, - 4755.84627752788110767815,-134659.959864969306392456, - -115132.259675553483497211 }; + 315.350626979604161529144,-1015.15636749021914166146, + -3107.77167157231109440444,22538.1184209801510330112, + 4755.84627752788110767815,-134659.959864969306392456, + -115132.259675553483497211 }; static doublereal c__[7] = { -.001910444077728,8.4171387781295e-4, - -5.952379913043012e-4,7.93650793500350248e-4, - -.002777777777777681622553,.08333333333333333331554247, - .0057083835261 }; + -5.952379913043012e-4,7.93650793500350248e-4, + -.002777777777777681622553,.08333333333333333331554247, + .0057083835261 }; static doublereal half = .5; static doublereal twelve = 12.; static doublereal two = 2.; @@ -55,7 +55,7 @@ doublereal dgamma_(const doublereal *x) /* Builtin functions */ double d_int(doublereal *), sin(doublereal), log(doublereal), exp( - doublereal); + doublereal); /* Local variables */ integer i__, n; @@ -194,20 +194,20 @@ doublereal dgamma_(const doublereal *x) /* ---------------------------------------------------------------------- */ /* Argument is negative */ /* ---------------------------------------------------------------------- */ - y = -(*x); - y1 = d_int(&y); - res = y - y1; - if (res != zero) { - d__1 = y1 * half; - if (y1 != d_int(&d__1) * two) { - parity = TRUE_; - } - fact = -pi / sin(pi * res); - y += one; - } else { - res = xinf; - goto L900; - } + y = -(*x); + y1 = d_int(&y); + res = y - y1; + if (res != zero) { + d__1 = y1 * half; + if (y1 != d_int(&d__1) * two) { + parity = TRUE_; + } + fact = -pi / sin(pi * res); + y += one; + } else { + res = xinf; + goto L900; + } } /* ---------------------------------------------------------------------- */ /* Argument is positive */ @@ -216,82 +216,82 @@ doublereal dgamma_(const doublereal *x) /* ---------------------------------------------------------------------- */ /* Argument .LT. EPS */ /* ---------------------------------------------------------------------- */ - if (y >= xminin) { - res = one / y; - } else { - res = xinf; - goto L900; - } + if (y >= xminin) { + res = one / y; + } else { + res = xinf; + goto L900; + } } else if (y < twelve) { - y1 = y; - if (y < one) { + y1 = y; + if (y < one) { /* ---------------------------------------------------------------------- */ /* 0.0 .LT. argument .LT. 1.0 */ /* ---------------------------------------------------------------------- */ - z__ = y; - y += one; - } else { + z__ = y; + y += one; + } else { /* ---------------------------------------------------------------------- */ /* 1.0 .LT. argument .LT. 12.0, reduce argument if necessary */ /* ---------------------------------------------------------------------- */ - n = (integer) y - 1; - y -= (doublereal) n; - z__ = y - one; - } + n = (integer) y - 1; + y -= (doublereal) n; + z__ = y - one; + } /* ---------------------------------------------------------------------- */ /* Evaluate approximation for 1.0 .LT. argument .LT. 2.0 */ /* ---------------------------------------------------------------------- */ - xnum = zero; - xden = one; - for (i__ = 1; i__ <= 8; ++i__) { - xnum = (xnum + p[i__ - 1]) * z__; - xden = xden * z__ + q[i__ - 1]; + xnum = zero; + xden = one; + for (i__ = 1; i__ <= 8; ++i__) { + xnum = (xnum + p[i__ - 1]) * z__; + xden = xden * z__ + q[i__ - 1]; /* L260: */ - } - res = xnum / xden + one; - if (y1 < y) { + } + res = xnum / xden + one; + if (y1 < y) { /* ---------------------------------------------------------------------- */ /* Adjust result for case 0.0 .LT. argument .LT. 1.0 */ /* ---------------------------------------------------------------------- */ - res /= y1; - } else if (y1 > y) { + res /= y1; + } else if (y1 > y) { /* ---------------------------------------------------------------------- */ /* Adjust result for case 2.0 .LT. argument .LT. 12.0 */ /* ---------------------------------------------------------------------- */ - i__1 = n; - for (i__ = 1; i__ <= i__1; ++i__) { - res *= y; - y += one; + i__1 = n; + for (i__ = 1; i__ <= i__1; ++i__) { + res *= y; + y += one; /* L290: */ - } - } + } + } } else { /* ---------------------------------------------------------------------- */ /* Evaluate for argument .GE. 12.0, */ /* ---------------------------------------------------------------------- */ - if (y <= xbig) { - ysq = y * y; - sum = c__[6]; - for (i__ = 1; i__ <= 6; ++i__) { - sum = sum / ysq + c__[i__ - 1]; + if (y <= xbig) { + ysq = y * y; + sum = c__[6]; + for (i__ = 1; i__ <= 6; ++i__) { + sum = sum / ysq + c__[i__ - 1]; /* L350: */ - } - sum = sum / y - y + sqrtpi; - sum += (y - half) * log(y); - res = exp(sum); - } else { - res = xinf; - goto L900; - } + } + sum = sum / y - y + sqrtpi; + sum += (y - half) * log(y); + res = exp(sum); + } else { + res = xinf; + goto L900; + } } /* ---------------------------------------------------------------------- */ /* Final adjustments and return */ /* ---------------------------------------------------------------------- */ if (parity) { - res = -res; + res = -res; } if (fact != one) { - res = fact / res; + res = fact / res; } /* S900 GAMMA = RES */ L900: diff --git a/core/contrib/include/dacecontrib.h b/core/contrib/include/dacecontrib.h index 8504c6f..dbc9874 100644 --- a/core/contrib/include/dacecontrib.h +++ b/core/contrib/include/dacecontrib.h @@ -30,7 +30,7 @@ This file contains all contributed routines to the DACE core. It is never included publicly by DACE users or high level interfaces. */ -/** \addtogroup DACE Contrib +/** \addtogroup DACEContrib Contrib * @{ */ diff --git a/core/contrib/include/f2c.h b/core/contrib/include/f2c.h index f97272e..0c7b7f3 100644 --- a/core/contrib/include/f2c.h +++ b/core/contrib/include/f2c.h @@ -2,7 +2,7 @@ /** barf [ba:rf] 2. "He suggested using FORTRAN, and everybody barfed." - - From The Shogakukan DICTIONARY OF NEW ENGLISH (Second edition) */ + - From The Shogakukan DICTIONARY OF NEW ENGLISH (Second edition) */ #ifndef F2C_INCLUDE #define F2C_INCLUDE @@ -19,11 +19,11 @@ typedef long int logical; typedef short int shortlogical; typedef char logical1; typedef char integer1; -#ifdef INTEGER_STAR_8 /* Adjust for integer*8. */ -typedef long long longint; /* system-dependent */ -typedef unsigned long long ulongint; /* system-dependent */ -#define qbit_clear(a,b) ((a) & ~((ulongint)1 << (b))) -#define qbit_set(a,b) ((a) | ((ulongint)1 << (b))) +#ifdef INTEGER_STAR_8 /* Adjust for integer*8. */ +typedef long long longint; /* system-dependent */ +typedef unsigned long long ulongint; /* system-dependent */ +#define qbit_clear(a,b) ((a) & ~((ulongint)1 << (b))) +#define qbit_set(a,b) ((a) | ((ulongint)1 << (b))) #endif #define TRUE_ (1) @@ -49,109 +49,109 @@ typedef long int ftnint; /*external read, write*/ typedef struct -{ flag cierr; - ftnint ciunit; - flag ciend; - char *cifmt; - ftnint cirec; +{ flag cierr; + ftnint ciunit; + flag ciend; + char *cifmt; + ftnint cirec; } cilist; /*internal read, write*/ typedef struct -{ flag icierr; - char *iciunit; - flag iciend; - char *icifmt; - ftnint icirlen; - ftnint icirnum; +{ flag icierr; + char *iciunit; + flag iciend; + char *icifmt; + ftnint icirlen; + ftnint icirnum; } icilist; /*open*/ typedef struct -{ flag oerr; - ftnint ounit; - char *ofnm; - ftnlen ofnmlen; - char *osta; - char *oacc; - char *ofm; - ftnint orl; - char *oblnk; +{ flag oerr; + ftnint ounit; + char *ofnm; + ftnlen ofnmlen; + char *osta; + char *oacc; + char *ofm; + ftnint orl; + char *oblnk; } olist; /*close*/ typedef struct -{ flag cerr; - ftnint cunit; - char *csta; +{ flag cerr; + ftnint cunit; + char *csta; } cllist; /*rewind, backspace, endfile*/ typedef struct -{ flag aerr; - ftnint aunit; +{ flag aerr; + ftnint aunit; } alist; /* inquire */ typedef struct -{ flag inerr; - ftnint inunit; - char *infile; - ftnlen infilen; - ftnint *inex; /*parameters in standard's order*/ - ftnint *inopen; - ftnint *innum; - ftnint *innamed; - char *inname; - ftnlen innamlen; - char *inacc; - ftnlen inacclen; - char *inseq; - ftnlen inseqlen; - char *indir; - ftnlen indirlen; - char *infmt; - ftnlen infmtlen; - char *inform; - ftnint informlen; - char *inunf; - ftnlen inunflen; - ftnint *inrecl; - ftnint *innrec; - char *inblank; - ftnlen inblanklen; +{ flag inerr; + ftnint inunit; + char *infile; + ftnlen infilen; + ftnint *inex; /*parameters in standard's order*/ + ftnint *inopen; + ftnint *innum; + ftnint *innamed; + char *inname; + ftnlen innamlen; + char *inacc; + ftnlen inacclen; + char *inseq; + ftnlen inseqlen; + char *indir; + ftnlen indirlen; + char *infmt; + ftnlen infmtlen; + char *inform; + ftnint informlen; + char *inunf; + ftnlen inunflen; + ftnint *inrecl; + ftnint *innrec; + char *inblank; + ftnlen inblanklen; } inlist; #define VOID void -union Multitype { /* for multiple entry points */ - integer1 g; - shortint h; - integer i; - /* longint j; */ - real r; - doublereal d; - complex c; - doublecomplex z; - }; +union Multitype { /* for multiple entry points */ + integer1 g; + shortint h; + integer i; + /* longint j; */ + real r; + doublereal d; + complex c; + doublecomplex z; + }; typedef union Multitype Multitype; -/*typedef long int Long;*/ /* No longer used; formerly in Namelist */ +/*typedef long int Long;*/ /* No longer used; formerly in Namelist */ -struct Vardesc { /* for Namelist */ - char *name; - char *addr; - ftnlen *dims; - int type; - }; +struct Vardesc { /* for Namelist */ + char *name; + char *addr; + ftnlen *dims; + int type; + }; typedef struct Vardesc Vardesc; struct Namelist { - char *name; - Vardesc **vars; - int nvars; - }; + char *name; + Vardesc **vars; + int nvars; + }; typedef struct Namelist Namelist; #define abs(x) ((x) >= 0 ? (x) : -(x)) @@ -160,9 +160,9 @@ typedef struct Namelist Namelist; #define max(a,b) ((a) >= (b) ? (a) : (b)) #define dmin(a,b) (doublereal)min(a,b) #define dmax(a,b) (doublereal)max(a,b) -#define bit_test(a,b) ((a) >> (b) & 1) -#define bit_clear(a,b) ((a) & ~((uinteger)1 << (b))) -#define bit_set(a,b) ((a) | ((uinteger)1 << (b))) +#define bit_test(a,b) ((a) >> (b) & 1) +#define bit_clear(a,b) ((a) & ~((uinteger)1 << (b))) +#define bit_set(a,b) ((a) | ((uinteger)1 << (b))) /* procedure parameter types for -A and -C++ */ @@ -193,10 +193,10 @@ typedef /* Character */ VOID (*H_fp)(); typedef /* Subroutine */ int (*S_fp)(); #endif /* E_fp is for real functions when -R is not specified */ -typedef VOID C_f; /* complex function */ -typedef VOID H_f; /* character function */ -typedef VOID Z_f; /* double complex function */ -typedef doublereal E_f; /* real function with -R not specified */ +typedef VOID C_f; /* complex function */ +typedef VOID H_f; /* character function */ +typedef VOID Z_f; /* double complex function */ +typedef doublereal E_f; /* real function with -R not specified */ /* undef any lower-case symbols that your C compiler predefines, e.g.: */ @@ -381,5 +381,5 @@ extern void z_exp(doublecomplex *, doublecomplex *); extern void z_log(doublecomplex *, doublecomplex *); extern void z_sin(doublecomplex *, doublecomplex *); extern void z_sqrt(doublecomplex *, doublecomplex *); - } + } #endif diff --git a/core/contrib/libf2c.c b/core/contrib/libf2c.c index 59f005c..d66ac2b 100644 --- a/core/contrib/libf2c.c +++ b/core/contrib/libf2c.c @@ -3,10 +3,10 @@ Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., - http://www.netlib.org/f2c/libf2c.zip + http://www.netlib.org/f2c/libf2c.zip */ -/** \addtogroup DACE Contrib +/** \addtogroup DACEContrib Contrib * @{ */ @@ -40,22 +40,22 @@ x = *ap; n = *bp; if(n != 0) - { - if(n < 0) - { - n = -n; - x = 1/x; - } - for(u = n; ; ) - { - if(u & 01) - pow *= x; - if(u >>= 1) - x *= x; - else - break; - } - } + { + if(n < 0) + { + n = -n; + x = 1/x; + } + for(u = n; ; ) + { + if(u & 01) + pow *= x; + if(u >>= 1) + x *= x; + else + break; + } + } return(pow); } #ifdef __cplusplus diff --git a/core/contrib/netlib/zeta.c b/core/contrib/netlib/zeta.c index 87886a3..1ae0e67 100644 --- a/core/contrib/netlib/zeta.c +++ b/core/contrib/netlib/zeta.c @@ -1,6 +1,6 @@ -/* zeta.c +/* zeta.c * - * Riemann zeta function of two arguments + * Riemann zeta function of two arguments * * * @@ -18,7 +18,7 @@ * * inf. * - -x - * zeta(x,q) = > (k+q) + * zeta(x,q) = > (k+q) * - * k=0 * @@ -26,11 +26,11 @@ * The Euler-Maclaurin summation formula is used to obtain * the expansion * - * n + * n * - -x - * zeta(x,q) = > (k+q) - * - - * k=1 + * zeta(x,q) = > (k+q) + * - + * k=1 * * 1-x inf. B x(x+1)...(x+2j) * (n+q) 1 - 2j @@ -98,26 +98,26 @@ int i; double a, b, k, s, t, w; if( x == 1.0 ) - goto retinf; + goto retinf; if( x < 1.0 ) - { + { domerr: - mtherr( "zeta", DOMAIN ); - return(0.0); - } + mtherr( "zeta", DOMAIN ); + return(0.0); + } if( q <= 0.0 ) - { - if(q == floor(q)) - { - mtherr( "zeta", SING ); + { + if(q == floor(q)) + { + mtherr( "zeta", SING ); retinf: - return( MAXNUM ); - } - if( x != floor(x) ) - goto domerr; /* because q^-x not defined */ - } + return( MAXNUM ); + } + if( x != floor(x) ) + goto domerr; /* because q^-x not defined */ + } /* Euler-Maclaurin summation formula */ /* @@ -134,14 +134,14 @@ a = q; i = 0; b = 0.0; while( (i < 9) || (a <= 9.0) ) - { - i += 1; - a += 1.0; - b = pow( a, -x ); - s += b; - if( fabs(b/s) < MACHEP ) - goto done; - } + { + i += 1; + a += 1.0; + b = pow( a, -x ); + s += b; + if( fabs(b/s) < MACHEP ) + goto done; + } w = a; s += b*w/(x-1.0); @@ -149,19 +149,19 @@ s -= 0.5 * b; a = 1.0; k = 0.0; for( i=0; i<12; i++ ) - { - a *= x + k; - b /= w; - t = a*b/A[i]; - s = s + t; - t = fabs(t/s); - if( t < MACHEP ) - goto done; - k += 1.0; - a *= x + k; - b /= w; - k += 1.0; - } + { + a *= x + k; + b /= w; + t = a*b/A[i]; + s = s + t; + t = fabs(t/s); + if( t < MACHEP ) + goto done; + k += 1.0; + a *= x + k; + b /= w; + k += 1.0; + } done: return(s); } @@ -175,11 +175,11 @@ return(s); s = pow( q, -x ); a = q; do - { - a += 2.0; - b = pow( a, -x ); - s += b; - } + { + a += 2.0; + b = pow( a, -x ); + s += b; + } while( b/s > MACHEP ); b = pow( 2.0, -x ); diff --git a/core/contrib/psi.c b/core/contrib/psi.c index d4a9e94..2d848e9 100644 --- a/core/contrib/psi.c +++ b/core/contrib/psi.c @@ -1,15 +1,21 @@ /* netlib/psi.f -- translated by f2c (version 20100827). You must link the resulting object file with libf2c: - on Microsoft Windows system, link with libf2c.lib; - on Linux or Unix systems, link with .../path/to/libf2c.a -lm - or, if you install libf2c.a in a standard place, with -lf2c -lm - -- in that order, at the end of the command line, as in - cc *.o -lf2c -lm - Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., - - http://www.netlib.org/f2c/libf2c.zip + on Microsoft Windows system, link with libf2c.lib; + on Linux or Unix systems, link with .../path/to/libf2c.a -lm + or, if you install libf2c.a in a standard place, with -lf2c -lm + -- in that order, at the end of the command line, as in + cc *.o -lf2c -lm + Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., + + http://www.netlib.org/f2c/libf2c.zip */ +/** \addtogroup DACEContrib Contrib + * @{ + */ + +/// @cond + #include "f2c.h" doublereal psi_(const doublereal *xx) @@ -24,20 +30,20 @@ doublereal psi_(const doublereal *xx) static doublereal x01d = 128.; static doublereal x02 = 6.9464496836234126266e-4; static doublereal p1[9] = { .004510468124576293416,5.4932855833000385356, - 376.46693175929276856,7952.5490849151998065,71451.59581895193321, - 306559.76301987365674,636069.97788964458797,580413.12783537569993, - 165856.95029761022321 }; + 376.46693175929276856,7952.5490849151998065,71451.59581895193321, + 306559.76301987365674,636069.97788964458797,580413.12783537569993, + 165856.95029761022321 }; static doublereal q1[8] = { 96.141654774222358525,2628.771579058119333, - 29862.49702225027792,162065.66091533671639,434878.80712768329037, - 542563.84537269993733,242421.85002017985252, - 6.4155223783576225996e-8 }; + 29862.49702225027792,162065.66091533671639,434878.80712768329037, + 542563.84537269993733,242421.85002017985252, + 6.4155223783576225996e-8 }; static doublereal p2[7] = { -2.7103228277757834192,-15.166271776896121383, - -19.784554148719218667,-8.8100958828312219821, - -1.4479614616899842986,-.073689600332394549911, - -6.5135387732718171306e-21 }; + -19.784554148719218667,-8.8100958828312219821, + -1.4479614616899842986,-.073689600332394549911, + -6.5135387732718171306e-21 }; static doublereal q2[6] = { 44.992760373789365846,202.40955312679931159, - 247.36979003315290057,107.42543875702278326,17.463965060678569906, - .88427520398873480342 }; + 247.36979003315290057,107.42543875702278326,17.463965060678569906, + .88427520398873480342 }; static doublereal fourth = .25; static doublereal half = .5; static doublereal one = 1.; @@ -209,25 +215,25 @@ doublereal psi_(const doublereal *xx) /* Check for valid arguments, then branch to appropriate algorithm */ /* ---------------------------------------------------------------------- */ if (-x >= xmax1 || w < xmin1) { - goto L410; + goto L410; } else if (x >= half) { - goto L200; + goto L200; /* ---------------------------------------------------------------------- */ /* X < 0.5, use reflection formula: psi(1-x) = psi(x) + pi * cot(pi*x) */ /* Use 1/X for PI*COTAN(PI*X) when XMIN1 < |X| <= XSMALL. */ /* ---------------------------------------------------------------------- */ } else if (w <= xsmall) { - aug = -one / x; - goto L150; + aug = -one / x; + goto L150; } /* ---------------------------------------------------------------------- */ /* Argument reduction for cot */ /* ---------------------------------------------------------------------- */ /* L100: */ if (x < zero) { - sgn = piov4; + sgn = piov4; } else { - sgn = -piov4; + sgn = -piov4; } w -= d_int(&w); nq = (integer) (w * four); @@ -239,11 +245,11 @@ doublereal psi_(const doublereal *xx) /* ---------------------------------------------------------------------- */ n = nq / 2; if (n + n != nq) { - w = one - w; + w = one - w; } z__ = piov4 * w; if (n % 2 != 0) { - sgn = -sgn; + sgn = -sgn; } /* ---------------------------------------------------------------------- */ /* determine the final value for -pi * cotan(pi*x) */ @@ -253,18 +259,18 @@ doublereal psi_(const doublereal *xx) /* ---------------------------------------------------------------------- */ /* Check for singularity */ /* ---------------------------------------------------------------------- */ - if (z__ == zero) { - goto L410; - } - aug = sgn * (four / tan(z__)); + if (z__ == zero) { + goto L410; + } + aug = sgn * (four / tan(z__)); } else { - aug = sgn * (four * tan(z__)); + aug = sgn * (four * tan(z__)); } L150: x = one - x; L200: if (x > three) { - goto L300; + goto L300; } /* ---------------------------------------------------------------------- */ /* 0.5 <= X <= 3.0 */ @@ -272,8 +278,8 @@ doublereal psi_(const doublereal *xx) den = x; upper = p1[0] * x; for (i__ = 1; i__ <= 7; ++i__) { - den = (den + q1[i__ - 1]) * x; - upper = (upper + p1[i__]) * x; + den = (den + q1[i__ - 1]) * x; + upper = (upper + p1[i__]) * x; /* L210: */ } den = (upper + p1[8]) / (den + q1[7]); @@ -285,15 +291,15 @@ doublereal psi_(const doublereal *xx) /* ---------------------------------------------------------------------- */ L300: if (x < xlarge) { - w = one / (x * x); - den = w; - upper = p2[0] * w; - for (i__ = 1; i__ <= 5; ++i__) { - den = (den + q2[i__ - 1]) * w; - upper = (upper + p2[i__]) * w; + w = one / (x * x); + den = w; + upper = p2[0] * w; + for (i__ = 1; i__ <= 5; ++i__) { + den = (den + q2[i__ - 1]) * w; + upper = (upper + p2[i__]) * w; /* L310: */ - } - aug = (upper + p2[6]) / (den + q2[5]) - half / x + aug; + } + aug = (upper + p2[6]) / (den + q2[5]) - half / x + aug; } ret_val = aug + log(x); goto L500; @@ -303,10 +309,13 @@ doublereal psi_(const doublereal *xx) L410: ret_val = xinf; if (x > zero) { - ret_val = -xinf; + ret_val = -xinf; } L500: return ret_val; /* ---------- Last card of PSI ---------- */ } /* psi_ */ + +/// @endcond +/** @}*/ diff --git a/core/contrib/ribesl.c b/core/contrib/ribesl.c index 745683f..5e3637c 100644 --- a/core/contrib/ribesl.c +++ b/core/contrib/ribesl.c @@ -1,16 +1,16 @@ /* ribesl.f -- translated by f2c (version 20100827). You must link the resulting object file with libf2c: - on Microsoft Windows system, link with libf2c.lib; - on Linux or Unix systems, link with .../path/to/libf2c.a -lm - or, if you install libf2c.a in a standard place, with -lf2c -lm - -- in that order, at the end of the command line, as in - cc *.o -lf2c -lm - Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., - - http://www.netlib.org/f2c/libf2c.zip + on Microsoft Windows system, link with libf2c.lib; + on Linux or Unix systems, link with .../path/to/libf2c.a -lm + or, if you install libf2c.a in a standard place, with -lf2c -lm + -- in that order, at the end of the command line, as in + cc *.o -lf2c -lm + Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., + + http://www.netlib.org/f2c/libf2c.zip */ -/** \addtogroup DACE Contrib +/** \addtogroup DACEContrib Contrib * @{ */ @@ -19,7 +19,7 @@ #include "f2c.h" /* Subroutine */ int ribesl_(const doublereal *x, const doublereal *alpha, const integer *nb, - const integer *ize, doublereal *b, integer *ncalc) + const integer *ize, doublereal *b, integer *ncalc) { /* Initialized data */ @@ -42,7 +42,7 @@ /* Builtin functions */ double sqrt(doublereal), pow_di(const doublereal *, const integer *), pow_dd( - const doublereal *, const doublereal *), exp(doublereal); + const doublereal *, const doublereal *), exp(doublereal); /* Local variables */ integer k, l, n; @@ -51,7 +51,7 @@ doublereal pold; integer nbmx; doublereal test, empal, halfx, tempa, tempb, tempc, psave, plast, tover, - emp2al; + emp2al; extern doublereal dgamma_(doublereal *); doublereal psavel; integer nstart; @@ -250,303 +250,303 @@ /* Check for X, NB, OR IZE out of range. */ /* ------------------------------------------------------------------- */ if (*nb > 0 && *x >= zero && *alpha >= zero && *alpha < one && ((*ize == 1 - && *x <= exparg) || (*ize == 2 && *x <= xlarge))) { + && *x <= exparg) || (*ize == 2 && *x <= xlarge))) { /* ------------------------------------------------------------------- */ /* Use 2-term ascending series for small X */ /* ------------------------------------------------------------------- */ - *ncalc = *nb; - magx = (integer) (*x); - if (*x >= rtnsig) { + *ncalc = *nb; + magx = (integer) (*x); + if (*x >= rtnsig) { /* ------------------------------------------------------------------- */ /* Initialize the forward sweep, the P-sequence of Olver */ /* ------------------------------------------------------------------- */ - nbmx = *nb - magx; - n = magx + 1; - i__1 = n + n; - en = (doublereal) i__1 + (*alpha + *alpha); - plast = one; - p = en / *x; + nbmx = *nb - magx; + n = magx + 1; + i__1 = n + n; + en = (doublereal) i__1 + (*alpha + *alpha); + plast = one; + p = en / *x; /* ------------------------------------------------------------------- */ /* Calculate general significance test */ /* ------------------------------------------------------------------- */ - test = ensig + ensig; - if (magx << 1 > nsig * 5) { - test = sqrt(test * p); - } else { - test /= pow_di(&const__, &magx); - } - if (nbmx >= 3) { + test = ensig + ensig; + if (magx << 1 > nsig * 5) { + test = sqrt(test * p); + } else { + test /= pow_di(&const__, &magx); + } + if (nbmx >= 3) { /* ------------------------------------------------------------------- */ /* Calculate P-sequence until N = NB-1. Check for possible overflow. */ /* ------------------------------------------------------------------- */ - tover = enten / ensig; - nstart = magx + 2; - nend = *nb - 1; - i__1 = nend; - for (k = nstart; k <= i__1; ++k) { - n = k; - en += two; - pold = plast; - plast = p; - p = en * plast / *x + pold; - if (p > tover) { + tover = enten / ensig; + nstart = magx + 2; + nend = *nb - 1; + i__1 = nend; + for (k = nstart; k <= i__1; ++k) { + n = k; + en += two; + pold = plast; + plast = p; + p = en * plast / *x + pold; + if (p > tover) { /* ------------------------------------------------------------------- */ /* To avoid overflow, divide P-sequence by TOVER. Calculate */ /* P-sequence until ABS(P) .GT. 1. */ /* ------------------------------------------------------------------- */ - tover = enten; - p /= tover; - plast /= tover; - psave = p; - psavel = plast; - nstart = n + 1; + tover = enten; + p /= tover; + plast /= tover; + psave = p; + psavel = plast; + nstart = n + 1; L60: - ++n; - en += two; - pold = plast; - plast = p; - p = en * plast / *x + pold; - if (p <= one) { - goto L60; - } - tempb = en / *x; + ++n; + en += two; + pold = plast; + plast = p; + p = en * plast / *x + pold; + if (p <= one) { + goto L60; + } + tempb = en / *x; /* ------------------------------------------------------------------- */ /* Calculate backward test, and find NCALC, the highest N */ /* such that the test is passed. */ /* ------------------------------------------------------------------- */ - test = pold * plast / ensig; - test *= half - half / (tempb * tempb); - p = plast * tover; - --n; - en -= two; - nend = min(*nb,n); - i__2 = nend; - for (l = nstart; l <= i__2; ++l) { - *ncalc = l; - pold = psavel; - psavel = psave; - psave = en * psavel / *x + pold; - if (psave * psavel > test) { - goto L90; - } + test = pold * plast / ensig; + test *= half - half / (tempb * tempb); + p = plast * tover; + --n; + en -= two; + nend = min(*nb,n); + i__2 = nend; + for (l = nstart; l <= i__2; ++l) { + *ncalc = l; + pold = psavel; + psavel = psave; + psave = en * psavel / *x + pold; + if (psave * psavel > test) { + goto L90; + } /* L80: */ - } - *ncalc = nend + 1; + } + *ncalc = nend + 1; L90: - --(*ncalc); - goto L120; - } + --(*ncalc); + goto L120; + } /* L100: */ - } - n = nend; - i__1 = n + n; - en = (doublereal) i__1 + (*alpha + *alpha); + } + n = nend; + i__1 = n + n; + en = (doublereal) i__1 + (*alpha + *alpha); /* ------------------------------------------------------------------- */ /* Calculate special significance test for NBMX .GT. 2. */ /* ------------------------------------------------------------------- */ /* Computing MAX */ - d__1 = test, d__2 = sqrt(plast * ensig) * sqrt(p + p); - test = max(d__1,d__2); - } + d__1 = test, d__2 = sqrt(plast * ensig) * sqrt(p + p); + test = max(d__1,d__2); + } /* ------------------------------------------------------------------- */ /* Calculate P-sequence until significance test passed. */ /* ------------------------------------------------------------------- */ L110: - ++n; - en += two; - pold = plast; - plast = p; - p = en * plast / *x + pold; - if (p < test) { - goto L110; - } + ++n; + en += two; + pold = plast; + plast = p; + p = en * plast / *x + pold; + if (p < test) { + goto L110; + } /* ------------------------------------------------------------------- */ /* Initialize the backward recursion and the normalization sum. */ /* ------------------------------------------------------------------- */ L120: - ++n; - en += two; - tempb = zero; - tempa = one / p; - em = (doublereal) n - one; - empal = em + *alpha; - emp2al = em - one + (*alpha + *alpha); - sum = tempa * empal * emp2al / em; - nend = n - *nb; - if (nend < 0) { + ++n; + en += two; + tempb = zero; + tempa = one / p; + em = (doublereal) n - one; + empal = em + *alpha; + emp2al = em - one + (*alpha + *alpha); + sum = tempa * empal * emp2al / em; + nend = n - *nb; + if (nend < 0) { /* ------------------------------------------------------------------- */ /* N .LT. NB, so store B(N) and set higher orders to zero. */ /* ------------------------------------------------------------------- */ - b[n] = tempa; - nend = -nend; - i__1 = nend; - for (l = 1; l <= i__1; ++l) { + b[n] = tempa; + nend = -nend; + i__1 = nend; + for (l = 1; l <= i__1; ++l) { /* L130: */ - b[n + l] = zero; - } - } else { - if (nend > 0) { + b[n + l] = zero; + } + } else { + if (nend > 0) { /* ------------------------------------------------------------------- */ /* Recur backward via difference equation, calculating (but */ /* not storing) B(N), until N = NB. */ /* ------------------------------------------------------------------- */ - i__1 = nend; - for (l = 1; l <= i__1; ++l) { - --n; - en -= two; - tempc = tempb; - tempb = tempa; - tempa = en * tempb / *x + tempc; - em -= one; - emp2al -= one; - if (n == 1) { - goto L150; - } - if (n == 2) { - emp2al = one; - } - empal -= one; - sum = (sum + tempa * empal) * emp2al / em; + i__1 = nend; + for (l = 1; l <= i__1; ++l) { + --n; + en -= two; + tempc = tempb; + tempb = tempa; + tempa = en * tempb / *x + tempc; + em -= one; + emp2al -= one; + if (n == 1) { + goto L150; + } + if (n == 2) { + emp2al = one; + } + empal -= one; + sum = (sum + tempa * empal) * emp2al / em; /* L140: */ - } - } + } + } /* ------------------------------------------------------------------- */ /* Store B(NB) */ /* ------------------------------------------------------------------- */ L150: - b[n] = tempa; - if (*nb <= 1) { - sum = sum + sum + tempa; - goto L230; - } + b[n] = tempa; + if (*nb <= 1) { + sum = sum + sum + tempa; + goto L230; + } /* ------------------------------------------------------------------- */ /* Calculate and Store B(NB-1) */ /* ------------------------------------------------------------------- */ - --n; - en -= two; - b[n] = en * tempa / *x + tempb; - if (n == 1) { - goto L220; - } - em -= one; - emp2al -= one; - if (n == 2) { - emp2al = one; - } - empal -= one; - sum = (sum + b[n] * empal) * emp2al / em; - } - nend = n - 2; - if (nend > 0) { + --n; + en -= two; + b[n] = en * tempa / *x + tempb; + if (n == 1) { + goto L220; + } + em -= one; + emp2al -= one; + if (n == 2) { + emp2al = one; + } + empal -= one; + sum = (sum + b[n] * empal) * emp2al / em; + } + nend = n - 2; + if (nend > 0) { /* ------------------------------------------------------------------- */ /* Calculate via difference equation and store B(N), until N = 2. */ /* ------------------------------------------------------------------- */ - i__1 = nend; - for (l = 1; l <= i__1; ++l) { - --n; - en -= two; - b[n] = en * b[n + 1] / *x + b[n + 2]; - em -= one; - emp2al -= one; - if (n == 2) { - emp2al = one; - } - empal -= one; - sum = (sum + b[n] * empal) * emp2al / em; + i__1 = nend; + for (l = 1; l <= i__1; ++l) { + --n; + en -= two; + b[n] = en * b[n + 1] / *x + b[n + 2]; + em -= one; + emp2al -= one; + if (n == 2) { + emp2al = one; + } + empal -= one; + sum = (sum + b[n] * empal) * emp2al / em; /* L200: */ - } - } + } + } /* ------------------------------------------------------------------- */ /* Calculate B(1) */ /* ------------------------------------------------------------------- */ - b[1] = two * empal * b[2] / *x + b[3]; + b[1] = two * empal * b[2] / *x + b[3]; L220: - sum = sum + sum + b[1]; + sum = sum + sum + b[1]; /* ------------------------------------------------------------------- */ /* Normalize. Divide all B(N) by sum. */ /* ------------------------------------------------------------------- */ L230: - if (*alpha != zero) { - d__1 = one + *alpha; - d__2 = *x * half; - d__3 = -(*alpha); - sum = sum * dgamma_(&d__1) * pow_dd(&d__2, &d__3); - } - if (*ize == 1) { - sum *= exp(-(*x)); - } - tempa = enmten; - if (sum > one) { - tempa *= sum; - } - i__1 = *nb; - for (n = 1; n <= i__1; ++n) { - if (b[n] < tempa) { - b[n] = zero; - } - b[n] /= sum; + if (*alpha != zero) { + d__1 = one + *alpha; + d__2 = *x * half; + d__3 = -(*alpha); + sum = sum * dgamma_(&d__1) * pow_dd(&d__2, &d__3); + } + if (*ize == 1) { + sum *= exp(-(*x)); + } + tempa = enmten; + if (sum > one) { + tempa *= sum; + } + i__1 = *nb; + for (n = 1; n <= i__1; ++n) { + if (b[n] < tempa) { + b[n] = zero; + } + b[n] /= sum; /* L260: */ - } - return 0; + } + return 0; /* ------------------------------------------------------------------- */ /* Two-term ascending series for small X. */ /* ------------------------------------------------------------------- */ - } else { - tempa = one; - empal = one + *alpha; - halfx = zero; - if (*x > enmten) { - halfx = half * *x; - } - if (*alpha != zero) { - tempa = pow_dd(&halfx, alpha) / dgamma_(&empal); - } - if (*ize == 2) { - tempa *= exp(-(*x)); - } - tempb = zero; - if (*x + one > one) { - tempb = halfx * halfx; - } - b[1] = tempa + tempa * tempb / empal; - if (*x != zero && b[1] == zero) { - *ncalc = 0; - } - if (*nb > 1) { - if (*x == zero) { - i__1 = *nb; - for (n = 2; n <= i__1; ++n) { - b[n] = zero; + } else { + tempa = one; + empal = one + *alpha; + halfx = zero; + if (*x > enmten) { + halfx = half * *x; + } + if (*alpha != zero) { + tempa = pow_dd(&halfx, alpha) / dgamma_(&empal); + } + if (*ize == 2) { + tempa *= exp(-(*x)); + } + tempb = zero; + if (*x + one > one) { + tempb = halfx * halfx; + } + b[1] = tempa + tempa * tempb / empal; + if (*x != zero && b[1] == zero) { + *ncalc = 0; + } + if (*nb > 1) { + if (*x == zero) { + i__1 = *nb; + for (n = 2; n <= i__1; ++n) { + b[n] = zero; /* L310: */ - } - } else { + } + } else { /* ------------------------------------------------------------------- */ /* Calculate higher-order functions. */ /* ------------------------------------------------------------------- */ - tempc = halfx; - tover = (enmten + enmten) / *x; - if (tempb != zero) { - tover = enmten / tempb; - } - i__1 = *nb; - for (n = 2; n <= i__1; ++n) { - tempa /= empal; - empal += one; - tempa *= tempc; - if (tempa <= tover * empal) { - tempa = zero; - } - b[n] = tempa + tempa * tempb / empal; - if (b[n] == zero && *ncalc > n) { - *ncalc = n - 1; - } + tempc = halfx; + tover = (enmten + enmten) / *x; + if (tempb != zero) { + tover = enmten / tempb; + } + i__1 = *nb; + for (n = 2; n <= i__1; ++n) { + tempa /= empal; + empal += one; + tempa *= tempc; + if (tempa <= tover * empal) { + tempa = zero; + } + b[n] = tempa + tempa * tempb / empal; + if (b[n] == zero && *ncalc > n) { + *ncalc = n - 1; + } /* L340: */ - } - } - } - } + } + } + } + } } else { - *ncalc = min(*nb,0) - 1; + *ncalc = min(*nb,0) - 1; } return 0; /* ---------- Last line of RIBESL ---------- */ diff --git a/core/contrib/rjbesl.c b/core/contrib/rjbesl.c index 0acbc0a..d2f3cc4 100644 --- a/core/contrib/rjbesl.c +++ b/core/contrib/rjbesl.c @@ -1,16 +1,16 @@ /* rjbesl.f -- translated by f2c (version 20100827). You must link the resulting object file with libf2c: - on Microsoft Windows system, link with libf2c.lib; - on Linux or Unix systems, link with .../path/to/libf2c.a -lm - or, if you install libf2c.a in a standard place, with -lf2c -lm - -- in that order, at the end of the command line, as in - cc *.o -lf2c -lm - Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., - - http://www.netlib.org/f2c/libf2c.zip + on Microsoft Windows system, link with libf2c.lib; + on Linux or Unix systems, link with .../path/to/libf2c.a -lm + or, if you install libf2c.a in a standard place, with -lf2c -lm + -- in that order, at the end of the command line, as in + cc *.o -lf2c -lm + Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., + + http://www.netlib.org/f2c/libf2c.zip */ -/** \addtogroup DACE Contrib +/** \addtogroup DACEContrib Contrib * @{ */ @@ -19,7 +19,7 @@ #include "f2c.h" /* Subroutine */ int rjbesl_(const doublereal *x, const doublereal *alpha, const integer *nb, - doublereal *b, integer *ncalc) + doublereal *b, integer *ncalc) { /* Initialized data */ @@ -34,11 +34,11 @@ static doublereal enmten = 8.9e-308; static doublereal xlarge = 1e4; static doublereal fact[25] = { 1.,1.,2.,6.,24.,120.,720.,5040.,40320., - 362880.,3628800.,39916800.,479001600.,6227020800.,87178291200., - 1.307674368e12,2.0922789888e13,3.55687428096e14,6.402373705728e15, - 1.21645100408832e17,2.43290200817664e18,5.109094217170944e19, - 1.12400072777760768e21,2.585201673888497664e22, - 6.2044840173323943936e23 }; + 362880.,3628800.,39916800.,479001600.,6227020800.,87178291200., + 1.307674368e12,2.0922789888e13,3.55687428096e14,6.402373705728e15, + 1.21645100408832e17,2.43290200817664e18,5.109094217170944e19, + 1.12400072777760768e21,2.585201673888497664e22, + 6.2044840173323943936e23 }; static doublereal twopi1 = 6.28125; static doublereal twopi2 = .001935307179586476925286767; static doublereal zero = 0.; @@ -54,7 +54,7 @@ /* Builtin functions */ double pow_dd(const doublereal *, const doublereal *), sqrt(doublereal), d_int( - const doublereal *), sin(doublereal), cos(doublereal); + const doublereal *), sin(doublereal), cos(doublereal); /* Local variables */ integer i__, j, k, l, m, n; @@ -65,7 +65,7 @@ doublereal pold; integer nbmx; doublereal vcos, test, vsin, alpem, halfx, tempa, tempb, tempc, psave, - plast, tover, alp2em; + plast, tover, alp2em; extern doublereal dgamma_(const doublereal *); doublereal psavel; integer nstart; @@ -250,385 +250,385 @@ /* --------------------------------------------------------------------- */ magx = (integer) (*x); if (*nb > 0 && *x >= zero && *x <= xlarge && *alpha >= zero && *alpha < - one) { + one) { /* --------------------------------------------------------------------- */ /* Initialize result array to zero. */ /* --------------------------------------------------------------------- */ - *ncalc = *nb; - i__1 = *nb; - for (i__ = 1; i__ <= i__1; ++i__) { - b[i__] = zero; + *ncalc = *nb; + i__1 = *nb; + for (i__ = 1; i__ <= i__1; ++i__) { + b[i__] = zero; /* L20: */ - } + } /* --------------------------------------------------------------------- */ /* Branch to use 2-term ascending series for small X and asymptotic */ /* form for large X when NB is not too large. */ /* --------------------------------------------------------------------- */ - if (*x < rtnsig) { + if (*x < rtnsig) { /* --------------------------------------------------------------------- */ /* Two-term ascending series for small X. */ /* --------------------------------------------------------------------- */ - tempa = one; - alpem = one + *alpha; - halfx = zero; - if (*x > enmten) { - halfx = half * *x; - } - if (*alpha != zero) { - tempa = pow_dd(&halfx, alpha) / (*alpha * dgamma_(alpha)); - } - tempb = zero; - if (*x + one > one) { - tempb = -halfx * halfx; - } - b[1] = tempa + tempa * tempb / alpem; - if (*x != zero && b[1] == zero) { - *ncalc = 0; - } - if (*nb != 1) { - if (*x <= zero) { - i__1 = *nb; - for (n = 2; n <= i__1; ++n) { - b[n] = zero; + tempa = one; + alpem = one + *alpha; + halfx = zero; + if (*x > enmten) { + halfx = half * *x; + } + if (*alpha != zero) { + tempa = pow_dd(&halfx, alpha) / (*alpha * dgamma_(alpha)); + } + tempb = zero; + if (*x + one > one) { + tempb = -halfx * halfx; + } + b[1] = tempa + tempa * tempb / alpem; + if (*x != zero && b[1] == zero) { + *ncalc = 0; + } + if (*nb != 1) { + if (*x <= zero) { + i__1 = *nb; + for (n = 2; n <= i__1; ++n) { + b[n] = zero; /* L30: */ - } - } else { + } + } else { /* --------------------------------------------------------------------- */ /* Calculate higher order functions. */ /* --------------------------------------------------------------------- */ - tempc = halfx; - tover = (enmten + enmten) / *x; - if (tempb != zero) { - tover = enmten / tempb; - } - i__1 = *nb; - for (n = 2; n <= i__1; ++n) { - tempa /= alpem; - alpem += one; - tempa *= tempc; - if (tempa <= tover * alpem) { - tempa = zero; - } - b[n] = tempa + tempa * tempb / alpem; - if (b[n] == zero && *ncalc > n) { - *ncalc = n - 1; - } + tempc = halfx; + tover = (enmten + enmten) / *x; + if (tempb != zero) { + tover = enmten / tempb; + } + i__1 = *nb; + for (n = 2; n <= i__1; ++n) { + tempa /= alpem; + alpem += one; + tempa *= tempc; + if (tempa <= tover * alpem) { + tempa = zero; + } + b[n] = tempa + tempa * tempb / alpem; + if (b[n] == zero && *ncalc > n) { + *ncalc = n - 1; + } /* L50: */ - } - } - } - } else if (*x > twofiv && *nb <= magx + 1) { + } + } + } + } else if (*x > twofiv && *nb <= magx + 1) { /* --------------------------------------------------------------------- */ /* Asymptotic series for X .GT. 21.0. */ /* --------------------------------------------------------------------- */ - xc = sqrt(pi2 / *x); + xc = sqrt(pi2 / *x); /* Computing 2nd power */ - d__1 = eighth / *x; - xin = d__1 * d__1; - m = 11; - if (*x >= three5) { - m = 8; - } - if (*x >= one30) { - m = 4; - } - xm = four * (doublereal) m; + d__1 = eighth / *x; + xin = d__1 * d__1; + m = 11; + if (*x >= three5) { + m = 8; + } + if (*x >= one30) { + m = 4; + } + xm = four * (doublereal) m; /* --------------------------------------------------------------------- */ /* Argument reduction for SIN and COS routines. */ /* --------------------------------------------------------------------- */ - d__1 = *x / (twopi1 + twopi2) + half; - t = d_int(&d__1); - z__ = *x - t * twopi1 - t * twopi2 - (*alpha + half) / pi2; - vsin = sin(z__); - vcos = cos(z__); - gnu = *alpha + *alpha; - for (i__ = 1; i__ <= 2; ++i__) { - s = (xm - one - gnu) * (xm - one + gnu) * xin * half; - t = (gnu - (xm - three)) * (gnu + (xm - three)); - capp = s * t / fact[m * 2]; - t1 = (gnu - (xm + one)) * (gnu + (xm + one)); - capq = s * t1 / fact[(m << 1) + 1]; - xk = xm; - k = m + m; - t1 = t; - i__1 = m; - for (j = 2; j <= i__1; ++j) { - xk -= four; - s = (xk - one - gnu) * (xk - one + gnu); - t = (gnu - (xk - three)) * (gnu + (xk - three)); - capp = (capp + one / fact[k - 2]) * s * t * xin; - capq = (capq + one / fact[k - 1]) * s * t1 * xin; - k += -2; - t1 = t; + d__1 = *x / (twopi1 + twopi2) + half; + t = d_int(&d__1); + z__ = *x - t * twopi1 - t * twopi2 - (*alpha + half) / pi2; + vsin = sin(z__); + vcos = cos(z__); + gnu = *alpha + *alpha; + for (i__ = 1; i__ <= 2; ++i__) { + s = (xm - one - gnu) * (xm - one + gnu) * xin * half; + t = (gnu - (xm - three)) * (gnu + (xm - three)); + capp = s * t / fact[m * 2]; + t1 = (gnu - (xm + one)) * (gnu + (xm + one)); + capq = s * t1 / fact[(m << 1) + 1]; + xk = xm; + k = m + m; + t1 = t; + i__1 = m; + for (j = 2; j <= i__1; ++j) { + xk -= four; + s = (xk - one - gnu) * (xk - one + gnu); + t = (gnu - (xk - three)) * (gnu + (xk - three)); + capp = (capp + one / fact[k - 2]) * s * t * xin; + capq = (capq + one / fact[k - 1]) * s * t1 * xin; + k += -2; + t1 = t; /* L70: */ - } - capp += one; - capq = (capq + one) * (gnu * gnu - one) * (eighth / *x); - b[i__] = xc * (capp * vcos - capq * vsin); - if (*nb == 1) { - goto L300; - } - t = vsin; - vsin = -vcos; - vcos = t; - gnu += two; + } + capp += one; + capq = (capq + one) * (gnu * gnu - one) * (eighth / *x); + b[i__] = xc * (capp * vcos - capq * vsin); + if (*nb == 1) { + goto L300; + } + t = vsin; + vsin = -vcos; + vcos = t; + gnu += two; /* L80: */ - } + } /* --------------------------------------------------------------------- */ /* If NB .GT. 2, compute J(X,ORDER+I) I = 2, NB-1 */ /* --------------------------------------------------------------------- */ - if (*nb > 2) { - gnu = *alpha + *alpha + two; - i__1 = *nb; - for (j = 3; j <= i__1; ++j) { - b[j] = gnu * b[j - 1] / *x - b[j - 2]; - gnu += two; + if (*nb > 2) { + gnu = *alpha + *alpha + two; + i__1 = *nb; + for (j = 3; j <= i__1; ++j) { + b[j] = gnu * b[j - 1] / *x - b[j - 2]; + gnu += two; /* L90: */ - } - } + } + } /* --------------------------------------------------------------------- */ /* Use recurrence to generate results. First initialize the */ /* calculation of P*S. */ /* --------------------------------------------------------------------- */ - } else { - nbmx = *nb - magx; - n = magx + 1; - i__1 = n + n; - en = (doublereal) i__1 + (*alpha + *alpha); - plast = one; - p = en / *x; + } else { + nbmx = *nb - magx; + n = magx + 1; + i__1 = n + n; + en = (doublereal) i__1 + (*alpha + *alpha); + plast = one; + p = en / *x; /* --------------------------------------------------------------------- */ /* Calculate general significance test. */ /* --------------------------------------------------------------------- */ - test = ensig + ensig; - if (nbmx >= 3) { + test = ensig + ensig; + if (nbmx >= 3) { /* --------------------------------------------------------------------- */ /* Calculate P*S until N = NB-1. Check for possible overflow. */ /* --------------------------------------------------------------------- */ - tover = enten / ensig; - nstart = magx + 2; - nend = *nb - 1; - i__1 = nstart + nstart; - en = (doublereal) i__1 - two + (*alpha + *alpha); - i__1 = nend; - for (k = nstart; k <= i__1; ++k) { - n = k; - en += two; - pold = plast; - plast = p; - p = en * plast / *x - pold; - if (p > tover) { + tover = enten / ensig; + nstart = magx + 2; + nend = *nb - 1; + i__1 = nstart + nstart; + en = (doublereal) i__1 - two + (*alpha + *alpha); + i__1 = nend; + for (k = nstart; k <= i__1; ++k) { + n = k; + en += two; + pold = plast; + plast = p; + p = en * plast / *x - pold; + if (p > tover) { /* --------------------------------------------------------------------- */ /* To avoid overflow, divide P*S by TOVER. Calculate P*S until */ /* ABS(P) .GT. 1. */ /* --------------------------------------------------------------------- */ - tover = enten; - p /= tover; - plast /= tover; - psave = p; - psavel = plast; - nstart = n + 1; + tover = enten; + p /= tover; + plast /= tover; + psave = p; + psavel = plast; + nstart = n + 1; L100: - ++n; - en += two; - pold = plast; - plast = p; - p = en * plast / *x - pold; - if (p <= one) { - goto L100; - } - tempb = en / *x; + ++n; + en += two; + pold = plast; + plast = p; + p = en * plast / *x - pold; + if (p <= one) { + goto L100; + } + tempb = en / *x; /* --------------------------------------------------------------------- */ /* Calculate backward test and find NCALC, the highest N such that */ /* the test is passed. */ /* --------------------------------------------------------------------- */ - test = pold * plast * (half - half / (tempb * tempb)); - test /= ensig; - p = plast * tover; - --n; - en -= two; - nend = min(*nb,n); - i__2 = nend; - for (l = nstart; l <= i__2; ++l) { - pold = psavel; - psavel = psave; - psave = en * psavel / *x - pold; - if (psave * psavel > test) { - *ncalc = l - 1; - goto L190; - } + test = pold * plast * (half - half / (tempb * tempb)); + test /= ensig; + p = plast * tover; + --n; + en -= two; + nend = min(*nb,n); + i__2 = nend; + for (l = nstart; l <= i__2; ++l) { + pold = psavel; + psavel = psave; + psave = en * psavel / *x - pold; + if (psave * psavel > test) { + *ncalc = l - 1; + goto L190; + } /* L110: */ - } - *ncalc = nend; - goto L190; - } + } + *ncalc = nend; + goto L190; + } /* L130: */ - } - n = nend; - i__1 = n + n; - en = (doublereal) i__1 + (*alpha + *alpha); + } + n = nend; + i__1 = n + n; + en = (doublereal) i__1 + (*alpha + *alpha); /* --------------------------------------------------------------------- */ /* Calculate special significance test for NBMX .GT. 2. */ /* --------------------------------------------------------------------- */ /* Computing MAX */ - d__1 = test, d__2 = sqrt(plast * ensig) * sqrt(p + p); - test = max(d__1,d__2); - } + d__1 = test, d__2 = sqrt(plast * ensig) * sqrt(p + p); + test = max(d__1,d__2); + } /* --------------------------------------------------------------------- */ /* Calculate P*S until significance test passes. */ /* --------------------------------------------------------------------- */ L140: - ++n; - en += two; - pold = plast; - plast = p; - p = en * plast / *x - pold; - if (p < test) { - goto L140; - } + ++n; + en += two; + pold = plast; + plast = p; + p = en * plast / *x - pold; + if (p < test) { + goto L140; + } /* --------------------------------------------------------------------- */ /* Initialize the backward recursion and the normalization sum. */ /* --------------------------------------------------------------------- */ L190: - ++n; - en += two; - tempb = zero; - tempa = one / p; - m = (n << 1) - (n / 2 << 2); - sum = zero; - i__1 = n / 2; - em = (doublereal) i__1; - alpem = em - one + *alpha; - alp2em = em + em + *alpha; - if (m != 0) { - sum = tempa * alpem * alp2em / em; - } - nend = n - *nb; - if (nend > 0) { + ++n; + en += two; + tempb = zero; + tempa = one / p; + m = (n << 1) - (n / 2 << 2); + sum = zero; + i__1 = n / 2; + em = (doublereal) i__1; + alpem = em - one + *alpha; + alp2em = em + em + *alpha; + if (m != 0) { + sum = tempa * alpem * alp2em / em; + } + nend = n - *nb; + if (nend > 0) { /* --------------------------------------------------------------------- */ /* Recur backward via difference equation, calculating (but not */ /* storing) B(N), until N = NB. */ /* --------------------------------------------------------------------- */ - i__1 = nend; - for (l = 1; l <= i__1; ++l) { - --n; - en -= two; - tempc = tempb; - tempb = tempa; - tempa = en * tempb / *x - tempc; - m = 2 - m; - if (m != 0) { - em -= one; - alp2em = em + em + *alpha; - if (n == 1) { - goto L210; - } - alpem = em - one + *alpha; - if (alpem == zero) { - alpem = one; - } - sum = (sum + tempa * alp2em) * alpem / em; - } + i__1 = nend; + for (l = 1; l <= i__1; ++l) { + --n; + en -= two; + tempc = tempb; + tempb = tempa; + tempa = en * tempb / *x - tempc; + m = 2 - m; + if (m != 0) { + em -= one; + alp2em = em + em + *alpha; + if (n == 1) { + goto L210; + } + alpem = em - one + *alpha; + if (alpem == zero) { + alpem = one; + } + sum = (sum + tempa * alp2em) * alpem / em; + } /* L200: */ - } - } + } + } /* --------------------------------------------------------------------- */ /* Store B(NB). */ /* --------------------------------------------------------------------- */ L210: - b[n] = tempa; - if (nend >= 0) { - if (*nb <= 1) { - alp2em = *alpha; - if (*alpha + one == one) { - alp2em = one; - } - sum += b[1] * alp2em; - goto L250; - } else { + b[n] = tempa; + if (nend >= 0) { + if (*nb <= 1) { + alp2em = *alpha; + if (*alpha + one == one) { + alp2em = one; + } + sum += b[1] * alp2em; + goto L250; + } else { /* --------------------------------------------------------------------- */ /* Calculate and store B(NB-1). */ /* --------------------------------------------------------------------- */ - --n; - en -= two; - b[n] = en * tempa / *x - tempb; - if (n == 1) { - goto L240; - } - m = 2 - m; - if (m != 0) { - em -= one; - alp2em = em + em + *alpha; - alpem = em - one + *alpha; - if (alpem == zero) { - alpem = one; - } - sum = (sum + b[n] * alp2em) * alpem / em; - } - } - } - nend = n - 2; - if (nend != 0) { + --n; + en -= two; + b[n] = en * tempa / *x - tempb; + if (n == 1) { + goto L240; + } + m = 2 - m; + if (m != 0) { + em -= one; + alp2em = em + em + *alpha; + alpem = em - one + *alpha; + if (alpem == zero) { + alpem = one; + } + sum = (sum + b[n] * alp2em) * alpem / em; + } + } + } + nend = n - 2; + if (nend != 0) { /* --------------------------------------------------------------------- */ /* Calculate via difference equation and store B(N), until N = 2. */ /* --------------------------------------------------------------------- */ - i__1 = nend; - for (l = 1; l <= i__1; ++l) { - --n; - en -= two; - b[n] = en * b[n + 1] / *x - b[n + 2]; - m = 2 - m; - if (m != 0) { - em -= one; - alp2em = em + em + *alpha; - alpem = em - one + *alpha; - if (alpem == zero) { - alpem = one; - } - sum = (sum + b[n] * alp2em) * alpem / em; - } + i__1 = nend; + for (l = 1; l <= i__1; ++l) { + --n; + en -= two; + b[n] = en * b[n + 1] / *x - b[n + 2]; + m = 2 - m; + if (m != 0) { + em -= one; + alp2em = em + em + *alpha; + alpem = em - one + *alpha; + if (alpem == zero) { + alpem = one; + } + sum = (sum + b[n] * alp2em) * alpem / em; + } /* L230: */ - } - } + } + } /* --------------------------------------------------------------------- */ /* Calculate B(1). */ /* --------------------------------------------------------------------- */ - b[1] = two * (*alpha + one) * b[2] / *x - b[3]; + b[1] = two * (*alpha + one) * b[2] / *x - b[3]; L240: - em -= one; - alp2em = em + em + *alpha; - if (alp2em == zero) { - alp2em = one; - } - sum += b[1] * alp2em; + em -= one; + alp2em = em + em + *alpha; + if (alp2em == zero) { + alp2em = one; + } + sum += b[1] * alp2em; /* --------------------------------------------------------------------- */ /* Normalize. Divide all B(N) by sum. */ /* --------------------------------------------------------------------- */ L250: - if (*alpha + one != one) { - d__1 = *x * half; - d__2 = -(*alpha); - sum = sum * dgamma_(alpha) * pow_dd(&d__1, &d__2); - } - tempa = enmten; - if (sum > one) { - tempa *= sum; - } - i__1 = *nb; - for (n = 1; n <= i__1; ++n) { - if ((d__1 = b[n], abs(d__1)) < tempa) { - b[n] = zero; - } - b[n] /= sum; + if (*alpha + one != one) { + d__1 = *x * half; + d__2 = -(*alpha); + sum = sum * dgamma_(alpha) * pow_dd(&d__1, &d__2); + } + tempa = enmten; + if (sum > one) { + tempa *= sum; + } + i__1 = *nb; + for (n = 1; n <= i__1; ++n) { + if ((d__1 = b[n], abs(d__1)) < tempa) { + b[n] = zero; + } + b[n] /= sum; /* L260: */ - } - } + } + } /* --------------------------------------------------------------------- */ /* Error return -- X, NB, or ALPHA is out of range. */ /* --------------------------------------------------------------------- */ } else { - b[1] = zero; - *ncalc = min(*nb,0) - 1; + b[1] = zero; + *ncalc = min(*nb,0) - 1; } /* --------------------------------------------------------------------- */ /* Exit */ diff --git a/core/contrib/rkbesl.c b/core/contrib/rkbesl.c index b58cfd2..7138045 100644 --- a/core/contrib/rkbesl.c +++ b/core/contrib/rkbesl.c @@ -1,16 +1,16 @@ /* rkbesl.f -- translated by f2c (version 20100827). You must link the resulting object file with libf2c: - on Microsoft Windows system, link with libf2c.lib; - on Linux or Unix systems, link with .../path/to/libf2c.a -lm - or, if you install libf2c.a in a standard place, with -lf2c -lm - -- in that order, at the end of the command line, as in - cc *.o -lf2c -lm - Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., - - http://www.netlib.org/f2c/libf2c.zip + on Microsoft Windows system, link with libf2c.lib; + on Linux or Unix systems, link with .../path/to/libf2c.a -lm + or, if you install libf2c.a in a standard place, with -lf2c -lm + -- in that order, at the end of the command line, as in + cc *.o -lf2c -lm + Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., + + http://www.netlib.org/f2c/libf2c.zip */ -/** \addtogroup DACE Contrib +/** \addtogroup DACEContrib Contrib * @{ */ @@ -19,7 +19,7 @@ #include "f2c.h" /* Subroutine */ int rkbesl_(const doublereal *x, const doublereal *alpha, const integer *nb, - const integer *ize, doublereal *bk, integer *ncalc) + const integer *ize, doublereal *bk, integer *ncalc) { /* Initialized data */ @@ -29,25 +29,25 @@ static doublereal xmin = 2.23e-308; static doublereal xmax = 705.342; static doublereal p[8] = { .805629875690432845,20.4045500205365151, - 157.705605106676174,536.671116469207504,900.382759291288778, - 730.923886650660393,229.299301509425145,.822467033424113231 }; + 157.705605106676174,536.671116469207504,900.382759291288778, + 730.923886650660393,229.299301509425145,.822467033424113231 }; static doublereal q[7] = { 29.4601986247850434,277.577868510221208, - 1206.70325591027438,2762.91444159791519,3443.74050506564618, - 2210.63190113378647,572.267338359892221 }; + 1206.70325591027438,2762.91444159791519,3443.74050506564618, + 2210.63190113378647,572.267338359892221 }; static doublereal r__[5] = { -.48672575865218401848,13.079485869097804016, - -101.96490580880537526,347.65409106507813131, - 3.495898124521934782e-4 }; + -101.96490580880537526,347.65409106507813131, + 3.495898124521934782e-4 }; static doublereal s[4] = { -25.579105509976461286,212.57260432226544008, - -610.69018684944109624,422.69668805777760407 }; + -610.69018684944109624,422.69668805777760407 }; static doublereal t[6] = { 1.6125990452916363814e-10, - 2.5051878502858255354e-8,2.7557319615147964774e-6, - 1.9841269840928373686e-4,.0083333333333334751799, - .16666666666666666446 }; + 2.5051878502858255354e-8,2.7557319615147964774e-6, + 1.9841269840928373686e-4,.0083333333333334751799, + .16666666666666666446 }; static doublereal estm[6] = { 52.0583,5.7607,2.7782,14.4303,185.3004, - 9.3715 }; + 9.3715 }; static doublereal one = 1.; static doublereal estf[7] = { 41.8341,7.1075,6.4306,42.511,1.35633, - 84.5096,20. }; + 84.5096,20. }; static doublereal two = 2.; static doublereal zero = 0.; static doublereal four = 4.; @@ -62,7 +62,7 @@ /* Builtin functions */ double log(doublereal), exp(doublereal), sinh(doublereal), sqrt( - doublereal), d_int(const doublereal *); + doublereal), d_int(const doublereal *); /* Local variables */ doublereal c__; @@ -267,342 +267,342 @@ enu = *alpha; *ncalc = min(*nb,0) - 2; if (*nb > 0 && (enu >= zero && enu < one) && (*ize >= 1 && *ize <= 2) && ( - *ize != 1 || ex <= xmax) && ex > zero) { - k = 0; - if (enu < sqxmin) { - enu = zero; - } - if (enu > half) { - k = 1; - enu -= one; - } - twonu = enu + enu; - iend = *nb + k - 1; - c__ = enu * enu; - d3 = -c__; - if (ex <= one) { + *ize != 1 || ex <= xmax) && ex > zero) { + k = 0; + if (enu < sqxmin) { + enu = zero; + } + if (enu > half) { + k = 1; + enu -= one; + } + twonu = enu + enu; + iend = *nb + k - 1; + c__ = enu * enu; + d3 = -c__; + if (ex <= one) { /* --------------------------------------------------------------------- */ /* Calculation of P0 = GAMMA(1+ALPHA) * (2/X)**ALPHA */ /* Q0 = GAMMA(1-ALPHA) * (X/2)**ALPHA */ /* --------------------------------------------------------------------- */ - d1 = zero; - d2 = p[0]; - t1 = one; - t2 = q[0]; - for (i__ = 2; i__ <= 7; i__ += 2) { - d1 = c__ * d1 + p[i__ - 1]; - d2 = c__ * d2 + p[i__]; - t1 = c__ * t1 + q[i__ - 1]; - t2 = c__ * t2 + q[i__]; + d1 = zero; + d2 = p[0]; + t1 = one; + t2 = q[0]; + for (i__ = 2; i__ <= 7; i__ += 2) { + d1 = c__ * d1 + p[i__ - 1]; + d2 = c__ * d2 + p[i__]; + t1 = c__ * t1 + q[i__ - 1]; + t2 = c__ * t2 + q[i__]; /* L10: */ - } - d1 = enu * d1; - t1 = enu * t1; - f1 = log(ex); - f0 = a + enu * (p[7] - enu * (d1 + d2) / (t1 + t2)) - f1; - q0 = exp(-enu * (a - enu * (p[7] + enu * (d1 - d2) / (t1 - t2)) - - f1)); - f1 = enu * f0; - p0 = exp(f1); + } + d1 = enu * d1; + t1 = enu * t1; + f1 = log(ex); + f0 = a + enu * (p[7] - enu * (d1 + d2) / (t1 + t2)) - f1; + q0 = exp(-enu * (a - enu * (p[7] + enu * (d1 - d2) / (t1 - t2)) - + f1)); + f1 = enu * f0; + p0 = exp(f1); /* --------------------------------------------------------------------- */ /* Calculation of F0 = */ /* --------------------------------------------------------------------- */ - d1 = r__[4]; - t1 = one; - for (i__ = 1; i__ <= 4; ++i__) { - d1 = c__ * d1 + r__[i__ - 1]; - t1 = c__ * t1 + s[i__ - 1]; + d1 = r__[4]; + t1 = one; + for (i__ = 1; i__ <= 4; ++i__) { + d1 = c__ * d1 + r__[i__ - 1]; + t1 = c__ * t1 + s[i__ - 1]; /* L20: */ - } - if (abs(f1) <= half) { - f1 *= f1; - d2 = zero; - for (i__ = 1; i__ <= 6; ++i__) { - d2 = f1 * d2 + t[i__ - 1]; + } + if (abs(f1) <= half) { + f1 *= f1; + d2 = zero; + for (i__ = 1; i__ <= 6; ++i__) { + d2 = f1 * d2 + t[i__ - 1]; /* L30: */ - } - d2 = f0 + f0 * f1 * d2; - } else { - d2 = sinh(f1) / enu; - } - f0 = d2 - enu * d1 / (t1 * p0); - if (ex <= tinyx) { + } + d2 = f0 + f0 * f1 * d2; + } else { + d2 = sinh(f1) / enu; + } + f0 = d2 - enu * d1 / (t1 * p0); + if (ex <= tinyx) { /* -------------------------------------------------------------------- */ /* X.LE.1.0E-10 */ /* Calculation of K(ALPHA,X) and X*K(ALPHA+1,X)/K(ALPHA,X) */ /* -------------------------------------------------------------------- */ - bk[1] = f0 + ex * f0; - if (*ize == 1) { - bk[1] -= ex * bk[1]; - } - ratio = p0 / f0; - c__ = ex * xinf; - if (k != 0) { + bk[1] = f0 + ex * f0; + if (*ize == 1) { + bk[1] -= ex * bk[1]; + } + ratio = p0 / f0; + c__ = ex * xinf; + if (k != 0) { /* -------------------------------------------------------------------- */ /* Calculation of K(ALPHA,X) and X*K(ALPHA+1,X)/K(ALPHA,X), */ /* ALPHA .GE. 1/2 */ /* -------------------------------------------------------------------- */ - *ncalc = -1; - if (bk[1] >= c__ / ratio) { - goto L500; - } - bk[1] = ratio * bk[1] / ex; - twonu += two; - ratio = twonu; - } - *ncalc = 1; - if (*nb == 1) { - goto L500; - } + *ncalc = -1; + if (bk[1] >= c__ / ratio) { + goto L500; + } + bk[1] = ratio * bk[1] / ex; + twonu += two; + ratio = twonu; + } + *ncalc = 1; + if (*nb == 1) { + goto L500; + } /* -------------------------------------------------------------------- */ /* Calculate K(ALPHA+L,X)/K(ALPHA+L-1,X), L = 1, 2, ... , NB-1 */ /* -------------------------------------------------------------------- */ - *ncalc = -1; - i__1 = *nb; - for (i__ = 2; i__ <= i__1; ++i__) { - if (ratio >= c__) { - goto L500; - } - bk[i__] = ratio / ex; - twonu += two; - ratio = twonu; + *ncalc = -1; + i__1 = *nb; + for (i__ = 2; i__ <= i__1; ++i__) { + if (ratio >= c__) { + goto L500; + } + bk[i__] = ratio / ex; + twonu += two; + ratio = twonu; /* L80: */ - } - *ncalc = 1; - goto L420; - } else { + } + *ncalc = 1; + goto L420; + } else { /* -------------------------------------------------------------------- */ /* 1.0E-10 .LT. X .LE. 1.0 */ /* -------------------------------------------------------------------- */ - c__ = one; - x2by4 = ex * ex / four; - p0 = half * p0; - q0 = half * q0; - d1 = -one; - d2 = zero; - bk1 = zero; - bk2 = zero; - f1 = f0; - f2 = p0; + c__ = one; + x2by4 = ex * ex / four; + p0 = half * p0; + q0 = half * q0; + d1 = -one; + d2 = zero; + bk1 = zero; + bk2 = zero; + f1 = f0; + f2 = p0; L100: - d1 += two; - d2 += one; - d3 = d1 + d3; - c__ = x2by4 * c__ / d2; - f0 = (d2 * f0 + p0 + q0) / d3; - p0 /= d2 - enu; - q0 /= d2 + enu; - t1 = c__ * f0; - t2 = c__ * (p0 - d2 * f0); - bk1 += t1; - bk2 += t2; - if ((d__1 = t1 / (f1 + bk1), abs(d__1)) > eps || (d__2 = t2 / - (f2 + bk2), abs(d__2)) > eps) { - goto L100; - } - bk1 = f1 + bk1; - bk2 = two * (f2 + bk2) / ex; - if (*ize == 2) { - d1 = exp(ex); - bk1 *= d1; - bk2 *= d1; - } - wminf = estf[0] * ex + estf[1]; - } - } else if (eps * ex > one) { + d1 += two; + d2 += one; + d3 = d1 + d3; + c__ = x2by4 * c__ / d2; + f0 = (d2 * f0 + p0 + q0) / d3; + p0 /= d2 - enu; + q0 /= d2 + enu; + t1 = c__ * f0; + t2 = c__ * (p0 - d2 * f0); + bk1 += t1; + bk2 += t2; + if ((d__1 = t1 / (f1 + bk1), abs(d__1)) > eps || (d__2 = t2 / + (f2 + bk2), abs(d__2)) > eps) { + goto L100; + } + bk1 = f1 + bk1; + bk2 = two * (f2 + bk2) / ex; + if (*ize == 2) { + d1 = exp(ex); + bk1 *= d1; + bk2 *= d1; + } + wminf = estf[0] * ex + estf[1]; + } + } else if (eps * ex > one) { /* -------------------------------------------------------------------- */ /* X .GT. ONE/EPS */ /* -------------------------------------------------------------------- */ - *ncalc = *nb; - bk1 = one / (d__ * sqrt(ex)); - i__1 = *nb; - for (i__ = 1; i__ <= i__1; ++i__) { - bk[i__] = bk1; + *ncalc = *nb; + bk1 = one / (d__ * sqrt(ex)); + i__1 = *nb; + for (i__ = 1; i__ <= i__1; ++i__) { + bk[i__] = bk1; /* L110: */ - } - goto L500; - } else { + } + goto L500; + } else { /* -------------------------------------------------------------------- */ /* X .GT. 1.0 */ /* -------------------------------------------------------------------- */ - twox = ex + ex; - blpha = zero; - ratio = zero; - if (ex <= four) { + twox = ex + ex; + blpha = zero; + ratio = zero; + if (ex <= four) { /* -------------------------------------------------------------------- */ /* Calculation of K(ALPHA+1,X)/K(ALPHA,X), 1.0 .LE. X .LE. 4.0 */ /* -------------------------------------------------------------------- */ - d__1 = estm[0] / ex + estm[1]; - d2 = d_int(&d__1); - m = (integer) d2; - d1 = d2 + d2; - d2 -= half; - d2 *= d2; - i__1 = m; - for (i__ = 2; i__ <= i__1; ++i__) { - d1 -= two; - d2 -= d1; - ratio = (d3 + d2) / (twox + d1 - ratio); + d__1 = estm[0] / ex + estm[1]; + d2 = d_int(&d__1); + m = (integer) d2; + d1 = d2 + d2; + d2 -= half; + d2 *= d2; + i__1 = m; + for (i__ = 2; i__ <= i__1; ++i__) { + d1 -= two; + d2 -= d1; + ratio = (d3 + d2) / (twox + d1 - ratio); /* L120: */ - } + } /* -------------------------------------------------------------------- */ /* Calculation of I(|ALPHA|,X) and I(|ALPHA|+1,X) by backward */ /* recurrence and K(ALPHA,X) from the wronskian */ /* -------------------------------------------------------------------- */ - d__1 = estm[2] * ex + estm[3]; - d2 = d_int(&d__1); - m = (integer) d2; - c__ = abs(enu); - d3 = c__ + c__; - d1 = d3 - one; - f1 = xmin; - f0 = (two * (c__ + d2) / ex + half * ex / (c__ + d2 + one)) * - xmin; - i__1 = m; - for (i__ = 3; i__ <= i__1; ++i__) { - d2 -= one; - f2 = (d3 + d2 + d2) * f0; - blpha = (one + d1 / d2) * (f2 + blpha); - f2 = f2 / ex + f1; - f1 = f0; - f0 = f2; + d__1 = estm[2] * ex + estm[3]; + d2 = d_int(&d__1); + m = (integer) d2; + c__ = abs(enu); + d3 = c__ + c__; + d1 = d3 - one; + f1 = xmin; + f0 = (two * (c__ + d2) / ex + half * ex / (c__ + d2 + one)) * + xmin; + i__1 = m; + for (i__ = 3; i__ <= i__1; ++i__) { + d2 -= one; + f2 = (d3 + d2 + d2) * f0; + blpha = (one + d1 / d2) * (f2 + blpha); + f2 = f2 / ex + f1; + f1 = f0; + f0 = f2; /* L130: */ - } - f1 = (d3 + two) * f0 / ex + f1; - d1 = zero; - t1 = one; - for (i__ = 1; i__ <= 7; ++i__) { - d1 = c__ * d1 + p[i__ - 1]; - t1 = c__ * t1 + q[i__ - 1]; + } + f1 = (d3 + two) * f0 / ex + f1; + d1 = zero; + t1 = one; + for (i__ = 1; i__ <= 7; ++i__) { + d1 = c__ * d1 + p[i__ - 1]; + t1 = c__ * t1 + q[i__ - 1]; /* L140: */ - } - p0 = exp(c__ * (a + c__ * (p[7] - c__ * d1 / t1) - log(ex))) / - ex; - f2 = (c__ + half - ratio) * f1 / ex; - bk1 = p0 + (d3 * f0 - f2 + f0 + blpha) / (f2 + f1 + f0) * p0; - if (*ize == 1) { - bk1 *= exp(-ex); - } - wminf = estf[2] * ex + estf[3]; - } else { + } + p0 = exp(c__ * (a + c__ * (p[7] - c__ * d1 / t1) - log(ex))) / + ex; + f2 = (c__ + half - ratio) * f1 / ex; + bk1 = p0 + (d3 * f0 - f2 + f0 + blpha) / (f2 + f1 + f0) * p0; + if (*ize == 1) { + bk1 *= exp(-ex); + } + wminf = estf[2] * ex + estf[3]; + } else { /* -------------------------------------------------------------------- */ /* Calculation of K(ALPHA,X) and K(ALPHA+1,X)/K(ALPHA,X), by backward */ /* recurrence, for X .GT. 4.0 */ /* -------------------------------------------------------------------- */ - d__1 = estm[4] / ex + estm[5]; - dm = d_int(&d__1); - m = (integer) dm; - d2 = dm - half; - d2 *= d2; - d1 = dm + dm; - i__1 = m; - for (i__ = 2; i__ <= i__1; ++i__) { - dm -= one; - d1 -= two; - d2 -= d1; - ratio = (d3 + d2) / (twox + d1 - ratio); - blpha = (ratio + ratio * blpha) / dm; + d__1 = estm[4] / ex + estm[5]; + dm = d_int(&d__1); + m = (integer) dm; + d2 = dm - half; + d2 *= d2; + d1 = dm + dm; + i__1 = m; + for (i__ = 2; i__ <= i__1; ++i__) { + dm -= one; + d1 -= two; + d2 -= d1; + ratio = (d3 + d2) / (twox + d1 - ratio); + blpha = (ratio + ratio * blpha) / dm; /* L160: */ - } - bk1 = one / ((d__ + d__ * blpha) * sqrt(ex)); - if (*ize == 1) { - bk1 *= exp(-ex); - } - wminf = estf[4] * (ex - (d__1 = ex - estf[6], abs(d__1))) + - estf[5]; - } + } + bk1 = one / ((d__ + d__ * blpha) * sqrt(ex)); + if (*ize == 1) { + bk1 *= exp(-ex); + } + wminf = estf[4] * (ex - (d__1 = ex - estf[6], abs(d__1))) + + estf[5]; + } /* -------------------------------------------------------------------- */ /* Calculation of K(ALPHA+1,X) from K(ALPHA,X) and */ /* K(ALPHA+1,X)/K(ALPHA,X) */ /* -------------------------------------------------------------------- */ - bk2 = bk1 + bk1 * (enu + half - ratio) / ex; - } + bk2 = bk1 + bk1 * (enu + half - ratio) / ex; + } /* -------------------------------------------------------------------- */ /* Calculation of 'NCALC', K(ALPHA+I,X), I = 0, 1, ... , NCALC-1, */ /* K(ALPHA+I,X)/K(ALPHA+I-1,X), I = NCALC, NCALC+1, ... , NB-1 */ /* -------------------------------------------------------------------- */ - *ncalc = *nb; - bk[1] = bk1; - if (iend == 0) { - goto L500; - } - j = 2 - k; - if (j > 0) { - bk[j] = bk2; - } - if (iend == 1) { - goto L500; - } + *ncalc = *nb; + bk[1] = bk1; + if (iend == 0) { + goto L500; + } + j = 2 - k; + if (j > 0) { + bk[j] = bk2; + } + if (iend == 1) { + goto L500; + } /* Computing MIN */ - i__1 = (integer) (wminf - enu); - m = min(i__1,iend); - i__1 = m; - for (i__ = 2; i__ <= i__1; ++i__) { - t1 = bk1; - bk1 = bk2; - twonu += two; - if (ex < one) { - if (bk1 >= xinf / twonu * ex) { - goto L195; - } - goto L187; - } else { - if (bk1 / ex >= xinf / twonu) { - goto L195; - } - } + i__1 = (integer) (wminf - enu); + m = min(i__1,iend); + i__1 = m; + for (i__ = 2; i__ <= i__1; ++i__) { + t1 = bk1; + bk1 = bk2; + twonu += two; + if (ex < one) { + if (bk1 >= xinf / twonu * ex) { + goto L195; + } + goto L187; + } else { + if (bk1 / ex >= xinf / twonu) { + goto L195; + } + } L187: - bk2 = twonu / ex * bk1 + t1; - itemp = i__; - ++j; - if (j > 0) { - bk[j] = bk2; - } + bk2 = twonu / ex * bk1 + t1; + itemp = i__; + ++j; + if (j > 0) { + bk[j] = bk2; + } /* L190: */ - } + } L195: - m = itemp; - if (m == iend) { - goto L500; - } - ratio = bk2 / bk1; - mplus1 = m + 1; - *ncalc = -1; - i__1 = iend; - for (i__ = mplus1; i__ <= i__1; ++i__) { - twonu += two; - ratio = twonu / ex + one / ratio; - ++j; - if (j > 1) { - bk[j] = ratio; - } else { - if (bk2 >= xinf / ratio) { - goto L500; - } - bk2 = ratio * bk2; - } + m = itemp; + if (m == iend) { + goto L500; + } + ratio = bk2 / bk1; + mplus1 = m + 1; + *ncalc = -1; + i__1 = iend; + for (i__ = mplus1; i__ <= i__1; ++i__) { + twonu += two; + ratio = twonu / ex + one / ratio; + ++j; + if (j > 1) { + bk[j] = ratio; + } else { + if (bk2 >= xinf / ratio) { + goto L500; + } + bk2 = ratio * bk2; + } /* L410: */ - } + } /* Computing MAX */ - i__1 = mplus1 - k; - *ncalc = max(i__1,1); - if (*ncalc == 1) { - bk[1] = bk2; - } - if (*nb == 1) { - goto L500; - } + i__1 = mplus1 - k; + *ncalc = max(i__1,1); + if (*ncalc == 1) { + bk[1] = bk2; + } + if (*nb == 1) { + goto L500; + } L420: - j = *ncalc + 1; - i__1 = *nb; - for (i__ = j; i__ <= i__1; ++i__) { - if (bk[*ncalc] >= xinf / bk[i__]) { - goto L500; - } - bk[i__] = bk[*ncalc] * bk[i__]; - *ncalc = i__; + j = *ncalc + 1; + i__1 = *nb; + for (i__ = j; i__ <= i__1; ++i__) { + if (bk[*ncalc] >= xinf / bk[i__]) { + goto L500; + } + bk[i__] = bk[*ncalc] * bk[i__]; + *ncalc = i__; /* L430: */ - } + } } L500: return 0; diff --git a/core/contrib/rybesl.c b/core/contrib/rybesl.c index 6e37451..ad4d119 100644 --- a/core/contrib/rybesl.c +++ b/core/contrib/rybesl.c @@ -1,16 +1,16 @@ /* rybesl.f -- translated by f2c (version 20100827). You must link the resulting object file with libf2c: - on Microsoft Windows system, link with libf2c.lib; - on Linux or Unix systems, link with .../path/to/libf2c.a -lm - or, if you install libf2c.a in a standard place, with -lf2c -lm - -- in that order, at the end of the command line, as in - cc *.o -lf2c -lm - Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., - - http://www.netlib.org/f2c/libf2c.zip + on Microsoft Windows system, link with libf2c.lib; + on Linux or Unix systems, link with .../path/to/libf2c.a -lm + or, if you install libf2c.a in a standard place, with -lf2c -lm + -- in that order, at the end of the command line, as in + cc *.o -lf2c -lm + Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., + + http://www.netlib.org/f2c/libf2c.zip */ -/** \addtogroup DACE Contrib +/** \addtogroup DACEContrib Contrib * @{ */ @@ -19,7 +19,7 @@ #include "f2c.h" /* Subroutine */ int rybesl_(const doublereal *x, const doublereal *alpha, const integer *nb, - doublereal *by, integer *ncalc) + doublereal *by, integer *ncalc) { /* Initialized data */ @@ -37,16 +37,16 @@ static doublereal half = .5; static doublereal xlarge = 1e8; static doublereal ch[21] = { -6.7735241822398840964e-24, - -6.1455180116049879894e-23,2.9017595056104745456e-21, - 1.3639417919073099464e-19,2.3826220476859635824e-18, - -9.0642907957550702534e-18,-1.4943667065169001769e-15, - -3.3919078305362211264e-14,-1.7023776642512729175e-13, - 9.1609750938768647911e-12,2.4230957900482704055e-10, - 1.7451364971382984243e-9,-3.3126119768180852711e-8, - -8.6592079961391259661e-7,-4.9717367041957398581e-6, - 7.6309597585908126618e-5,.0012719271366545622927, - .0017063050710955562222,-.07685284084478667369, - -.28387654227602353814,.92187029365045265648 }; + -6.1455180116049879894e-23,2.9017595056104745456e-21, + 1.3639417919073099464e-19,2.3826220476859635824e-18, + -9.0642907957550702534e-18,-1.4943667065169001769e-15, + -3.3919078305362211264e-14,-1.7023776642512729175e-13, + 9.1609750938768647911e-12,2.4230957900482704055e-10, + 1.7451364971382984243e-9,-3.3126119768180852711e-8, + -8.6592079961391259661e-7,-4.9717367041957398581e-6, + 7.6309597585908126618e-5,.0012719271366545622927, + .0017063050710955562222,-.07685284084478667369, + -.28387654227602353814,.92187029365045265648 }; static doublereal one = 1.; static doublereal two = 2.; static doublereal three = 3.; @@ -61,7 +61,7 @@ /* Builtin functions */ double d_int(const doublereal *), sqrt(doublereal), sin(doublereal), cos( - doublereal), log(doublereal), pow_dd(const doublereal *, const doublereal *); + doublereal), log(doublereal), pow_dd(const doublereal *, const doublereal *); /* Local variables */ doublereal b, c__, d__, e, f, g, h__; @@ -69,8 +69,8 @@ doublereal p, q, r__, s, d1, d2, q0, x2; integer na; doublereal pa, qa, en, ya, ex, pa1, qa1, en1, ya1, den, odd, aye, div, - dmu, xna, enu, alfa, ddiv, even, term, gamma, cosmu, sinmu, - twobyx; + dmu, xna, enu, alfa, ddiv, even, term, gamma, cosmu, sinmu, + twobyx; /* ---------------------------------------------------------------------- */ @@ -249,266 +249,266 @@ ex = *x; enu = *alpha; if (*nb > 0 && *x >= xmin && ex < xlarge && enu >= zero && enu < one) { - d__1 = enu + half; - xna = d_int(&d__1); - na = (integer) xna; - if (na == 1) { - enu -= xna; - } - if (enu == -half) { - p = sq2bpi / sqrt(ex); - ya = p * sin(ex); - ya1 = -p * cos(ex); - } else if (ex < three) { + d__1 = enu + half; + xna = d_int(&d__1); + na = (integer) xna; + if (na == 1) { + enu -= xna; + } + if (enu == -half) { + p = sq2bpi / sqrt(ex); + ya = p * sin(ex); + ya1 = -p * cos(ex); + } else if (ex < three) { /* ---------------------------------------------------------------------- */ /* Use Temme's scheme for small X */ /* ---------------------------------------------------------------------- */ - b = ex * half; - d__ = -log(b); - f = enu * d__; - d__1 = -enu; - e = pow_dd(&b, &d__1); - if (abs(enu) < del) { - c__ = onbpi; - } else { - c__ = enu / sin(enu * pi); - } + b = ex * half; + d__ = -log(b); + f = enu * d__; + d__1 = -enu; + e = pow_dd(&b, &d__1); + if (abs(enu) < del) { + c__ = onbpi; + } else { + c__ = enu / sin(enu * pi); + } /* ---------------------------------------------------------------------- */ /* Computation of sinh(f)/f */ /* ---------------------------------------------------------------------- */ - if (abs(f) < one) { - x2 = f * f; - en = ten9; - s = one; - for (i__ = 1; i__ <= 9; ++i__) { - s = s * x2 / en / (en - one) + one; - en -= two; + if (abs(f) < one) { + x2 = f * f; + en = ten9; + s = one; + for (i__ = 1; i__ <= 9; ++i__) { + s = s * x2 / en / (en - one) + one; + en -= two; /* L80: */ - } - } else { - s = (e - one / e) * half / f; - } + } + } else { + s = (e - one / e) * half / f; + } /* ---------------------------------------------------------------------- */ /* Computation of 1/gamma(1-a) using Chebyshev polynomials */ /* ---------------------------------------------------------------------- */ - x2 = enu * enu * eight; - aye = ch[0]; - even = zero; - alfa = ch[1]; - odd = zero; - for (i__ = 3; i__ <= 19; i__ += 2) { - even = -(aye + aye + even); - aye = -even * x2 - aye + ch[i__ - 1]; - odd = -(alfa + alfa + odd); - alfa = -odd * x2 - alfa + ch[i__]; + x2 = enu * enu * eight; + aye = ch[0]; + even = zero; + alfa = ch[1]; + odd = zero; + for (i__ = 3; i__ <= 19; i__ += 2) { + even = -(aye + aye + even); + aye = -even * x2 - aye + ch[i__ - 1]; + odd = -(alfa + alfa + odd); + alfa = -odd * x2 - alfa + ch[i__]; /* L40: */ - } - even = (even * half + aye) * x2 - aye + ch[20]; - odd = (odd + alfa) * two; - gamma = odd * enu + even; + } + even = (even * half + aye) * x2 - aye + ch[20]; + odd = (odd + alfa) * two; + gamma = odd * enu + even; /* ---------------------------------------------------------------------- */ /* End of computation of 1/gamma(1-a) */ /* ---------------------------------------------------------------------- */ - g = e * gamma; - e = (e + one / e) * half; - f = two * c__ * (odd * e + even * s * d__); - e = enu * enu; - p = g * c__; - q = onbpi / g; - c__ = enu * piby2; - if (abs(c__) < del) { - r__ = one; - } else { - r__ = sin(c__) / c__; - } - r__ = pi * c__ * r__ * r__; - c__ = one; - d__ = -b * b; - h__ = zero; - ya = f + r__ * q; - ya1 = p; - en = zero; + g = e * gamma; + e = (e + one / e) * half; + f = two * c__ * (odd * e + even * s * d__); + e = enu * enu; + p = g * c__; + q = onbpi / g; + c__ = enu * piby2; + if (abs(c__) < del) { + r__ = one; + } else { + r__ = sin(c__) / c__; + } + r__ = pi * c__ * r__ * r__; + c__ = one; + d__ = -b * b; + h__ = zero; + ya = f + r__ * q; + ya1 = p; + en = zero; L100: - en += one; - if ((d__1 = g / (one + abs(ya)), abs(d__1)) + (d__2 = h__ / (one - + abs(ya1)), abs(d__2)) > eps) { - f = (f * en + p + q) / (en * en - e); - c__ = c__ * d__ / en; - p /= en - enu; - q /= en + enu; - g = c__ * (f + r__ * q); - h__ = c__ * p - en * g; - ya += g; - ya1 += h__; - goto L100; - } - ya = -ya; - ya1 = -ya1 / b; - } else if (ex < thresh) { + en += one; + if ((d__1 = g / (one + abs(ya)), abs(d__1)) + (d__2 = h__ / (one + + abs(ya1)), abs(d__2)) > eps) { + f = (f * en + p + q) / (en * en - e); + c__ = c__ * d__ / en; + p /= en - enu; + q /= en + enu; + g = c__ * (f + r__ * q); + h__ = c__ * p - en * g; + ya += g; + ya1 += h__; + goto L100; + } + ya = -ya; + ya1 = -ya1 / b; + } else if (ex < thresh) { /* ---------------------------------------------------------------------- */ /* Use Temme's scheme for moderate X */ /* ---------------------------------------------------------------------- */ - c__ = (half - enu) * (half + enu); - b = ex + ex; - e = ex * onbpi * cos(enu * pi) / eps; - e *= e; - p = one; - q = -ex; - r__ = one + ex * ex; - s = r__; - en = two; + c__ = (half - enu) * (half + enu); + b = ex + ex; + e = ex * onbpi * cos(enu * pi) / eps; + e *= e; + p = one; + q = -ex; + r__ = one + ex * ex; + s = r__; + en = two; L200: - if (r__ * en * en < e) { - en1 = en + one; - d__ = (en - one + c__ / en) / s; - p = (en + en - p * d__) / en1; - q = (-b + q * d__) / en1; - s = p * p + q * q; - r__ *= s; - en = en1; - goto L200; - } - f = p / s; - p = f; - g = -q / s; - q = g; + if (r__ * en * en < e) { + en1 = en + one; + d__ = (en - one + c__ / en) / s; + p = (en + en - p * d__) / en1; + q = (-b + q * d__) / en1; + s = p * p + q * q; + r__ *= s; + en = en1; + goto L200; + } + f = p / s; + p = f; + g = -q / s; + q = g; L220: - en -= one; - if (en > zero) { - r__ = en1 * (two - p) - two; - s = b + en1 * q; - d__ = (en - one + c__ / en) / (r__ * r__ + s * s); - p = d__ * r__; - q = d__ * s; - e = f + one; - f = p * e - g * q; - g = q * e + p * g; - en1 = en; - goto L220; - } - f = one + f; - d__ = f * f + g * g; - pa = f / d__; - qa = -g / d__; - d__ = enu + half - p; - q += ex; - pa1 = (pa * q - qa * d__) / ex; - qa1 = (qa * q + pa * d__) / ex; - b = ex - piby2 * (enu + half); - c__ = cos(b); - s = sin(b); - d__ = sq2bpi / sqrt(ex); - ya = d__ * (pa * s + qa * c__); - ya1 = d__ * (qa1 * s - pa1 * c__); - } else { + en -= one; + if (en > zero) { + r__ = en1 * (two - p) - two; + s = b + en1 * q; + d__ = (en - one + c__ / en) / (r__ * r__ + s * s); + p = d__ * r__; + q = d__ * s; + e = f + one; + f = p * e - g * q; + g = q * e + p * g; + en1 = en; + goto L220; + } + f = one + f; + d__ = f * f + g * g; + pa = f / d__; + qa = -g / d__; + d__ = enu + half - p; + q += ex; + pa1 = (pa * q - qa * d__) / ex; + qa1 = (qa * q + pa * d__) / ex; + b = ex - piby2 * (enu + half); + c__ = cos(b); + s = sin(b); + d__ = sq2bpi / sqrt(ex); + ya = d__ * (pa * s + qa * c__); + ya1 = d__ * (qa1 * s - pa1 * c__); + } else { /* ---------------------------------------------------------------------- */ /* Use Campbell's asymptotic scheme. */ /* ---------------------------------------------------------------------- */ - na = 0; - d__1 = ex / fivpi; - d1 = d_int(&d__1); - i__ = (integer) d1; - dmu = ex - one5 * d1 - d1 * pim5 - (*alpha + half) * piby2; - if (i__ - (i__ / 2 << 1) == 0) { - cosmu = cos(dmu); - sinmu = sin(dmu); - } else { - cosmu = -cos(dmu); - sinmu = -sin(dmu); - } - ddiv = eight * ex; - dmu = *alpha; - den = sqrt(ex); - for (k = 1; k <= 2; ++k) { - p = cosmu; - cosmu = sinmu; - sinmu = -p; - d1 = (two * dmu - one) * (two * dmu + one); - d2 = zero; - div = ddiv; - p = zero; - q = zero; - q0 = d1 / div; - term = q0; - for (i__ = 2; i__ <= 20; ++i__) { - d2 += eight; - d1 -= d2; - div += ddiv; - term = -term * d1 / div; - p += term; - d2 += eight; - d1 -= d2; - div += ddiv; - term = term * d1 / div; - q += term; - if (abs(term) <= eps) { - goto L320; - } + na = 0; + d__1 = ex / fivpi; + d1 = d_int(&d__1); + i__ = (integer) d1; + dmu = ex - one5 * d1 - d1 * pim5 - (*alpha + half) * piby2; + if (i__ - (i__ / 2 << 1) == 0) { + cosmu = cos(dmu); + sinmu = sin(dmu); + } else { + cosmu = -cos(dmu); + sinmu = -sin(dmu); + } + ddiv = eight * ex; + dmu = *alpha; + den = sqrt(ex); + for (k = 1; k <= 2; ++k) { + p = cosmu; + cosmu = sinmu; + sinmu = -p; + d1 = (two * dmu - one) * (two * dmu + one); + d2 = zero; + div = ddiv; + p = zero; + q = zero; + q0 = d1 / div; + term = q0; + for (i__ = 2; i__ <= 20; ++i__) { + d2 += eight; + d1 -= d2; + div += ddiv; + term = -term * d1 / div; + p += term; + d2 += eight; + d1 -= d2; + div += ddiv; + term = term * d1 / div; + q += term; + if (abs(term) <= eps) { + goto L320; + } /* L310: */ - } + } L320: - p += one; - q += q0; - if (k == 1) { - ya = sq2bpi * (p * cosmu - q * sinmu) / den; - } else { - ya1 = sq2bpi * (p * cosmu - q * sinmu) / den; - } - dmu += one; + p += one; + q += q0; + if (k == 1) { + ya = sq2bpi * (p * cosmu - q * sinmu) / den; + } else { + ya1 = sq2bpi * (p * cosmu - q * sinmu) / den; + } + dmu += one; /* L350: */ - } - } - if (na == 1) { - h__ = two * (enu + one) / ex; - if (h__ > one) { - if (abs(ya1) > xinf / h__) { - h__ = zero; - ya = zero; - } - } - h__ = h__ * ya1 - ya; - ya = ya1; - ya1 = h__; - } + } + } + if (na == 1) { + h__ = two * (enu + one) / ex; + if (h__ > one) { + if (abs(ya1) > xinf / h__) { + h__ = zero; + ya = zero; + } + } + h__ = h__ * ya1 - ya; + ya = ya1; + ya1 = h__; + } /* ---------------------------------------------------------------------- */ /* Now have first one or two Y's */ /* ---------------------------------------------------------------------- */ - by[1] = ya; - by[2] = ya1; - if (ya1 == zero) { - *ncalc = 1; - } else { - aye = one + *alpha; - twobyx = two / ex; - *ncalc = 2; - i__1 = *nb; - for (i__ = 3; i__ <= i__1; ++i__) { - if (twobyx < one) { - if ((d__1 = by[i__ - 1], abs(d__1)) * twobyx >= xinf / - aye) { - goto L450; - } - } else { - if ((d__1 = by[i__ - 1], abs(d__1)) >= xinf / aye / - twobyx) { - goto L450; - } - } - by[i__] = twobyx * aye * by[i__ - 1] - by[i__ - 2]; - aye += one; - ++(*ncalc); + by[1] = ya; + by[2] = ya1; + if (ya1 == zero) { + *ncalc = 1; + } else { + aye = one + *alpha; + twobyx = two / ex; + *ncalc = 2; + i__1 = *nb; + for (i__ = 3; i__ <= i__1; ++i__) { + if (twobyx < one) { + if ((d__1 = by[i__ - 1], abs(d__1)) * twobyx >= xinf / + aye) { + goto L450; + } + } else { + if ((d__1 = by[i__ - 1], abs(d__1)) >= xinf / aye / + twobyx) { + goto L450; + } + } + by[i__] = twobyx * aye * by[i__ - 1] - by[i__ - 2]; + aye += one; + ++(*ncalc); /* L400: */ - } - } + } + } L450: - i__1 = *nb; - for (i__ = *ncalc + 1; i__ <= i__1; ++i__) { - by[i__] = zero; + i__1 = *nb; + for (i__ = *ncalc + 1; i__ <= i__1; ++i__) { + by[i__] = zero; /* L460: */ - } + } } else { - by[1] = zero; - *ncalc = min(*nb,0) - 1; + by[1] = zero; + *ncalc = min(*nb,0) - 1; } /* L900: */ return 0; diff --git a/core/contrib/zeta.c b/core/contrib/zeta.c index 0400f18..e7f2609 100644 --- a/core/contrib/zeta.c +++ b/core/contrib/zeta.c @@ -1,6 +1,6 @@ -/* zeta.c +/* zeta.c * - * Riemann zeta function of two arguments + * Riemann zeta function of two arguments * * * @@ -53,7 +53,7 @@ * Series, and Products, p. 1073; Academic Press, 1980. * */ - + /* Cephes Math Library Release 2.8: June, 2000 Copyright 1984, 1987, 2000 by Stephen L. Moshier @@ -61,6 +61,12 @@ Copyright 1984, 1987, 2000 by Stephen L. Moshier Modified interface by Alexander Wittig, November 2019. */ +/** \addtogroup DACEContrib Contrib + * @{ + */ + +/// @cond + #include /* Expansion coefficients @@ -88,68 +94,71 @@ static const double MACHEP = 1.11022302462515654042E-16; double zeta_(const double x, const double q, unsigned int *err) { - double a, b, k, s, t, w; + double a, b, k, s, t, w; - if( err ) *err = 0; + if( err ) *err = 0; - if( x == 1.0 ) - goto retinf; + if( x == 1.0 ) + goto retinf; - if( x < 1.0 ) - { + if( x < 1.0 ) + { domerr: - if( err ) *err = 1; - return(NAN); - } + if( err ) *err = 1; + return(NAN); + } - if( q <= 0.0 ) - { - if(q == floor(q)) - { + if( q <= 0.0 ) + { + if(q == floor(q)) + { retinf: - if( err ) *err = 2; - return(INFINITY); - } - if( x != floor(x) ) - goto domerr; /* because q^-x not defined */ - } + if( err ) *err = 2; + return(INFINITY); + } + if( x != floor(x) ) + goto domerr; /* because q^-x not defined */ + } - /* Euler-Maclaurin summation formula */ - /* Permit negative q but continue sum until n+q > +9 . - * This case should be handled by a reflection formula. - * If q<0 and x is an integer, there is a relation to - * the polygamma function. - */ - s = pow( q, -x ); - a = q; - b = 0.0; - for( unsigned int i = 0; (i < 9) || (a <= 9.0); i++ ) - { - a += 1.0; - b = pow( a, -x ); - s += b; - if( fabs(b/s) < MACHEP ) - goto done; - } + /* Euler-Maclaurin summation formula */ + /* Permit negative q but continue sum until n+q > +9 . + * This case should be handled by a reflection formula. + * If q<0 and x is an integer, there is a relation to + * the polygamma function. + */ + s = pow( q, -x ); + a = q; + b = 0.0; + for( unsigned int i = 0; (i < 9) || (a <= 9.0); i++ ) + { + a += 1.0; + b = pow( a, -x ); + s += b; + if( fabs(b/s) < MACHEP ) + goto done; + } - w = a; - s += b*w/(x-1.0); - s -= 0.5 * b; - a = 1.0; - k = 0.0; - for( unsigned int i=0; i<12; i++ ) - { - a *= x + k; - b /= w; - t = a*b/A[i]; - s = s + t; - if( fabs(t/s) < MACHEP ) - goto done; - k += 1.0; - a *= x + k; - b /= w; - k += 1.0; - } + w = a; + s += b*w/(x-1.0); + s -= 0.5 * b; + a = 1.0; + k = 0.0; + for( unsigned int i=0; i<12; i++ ) + { + a *= x + k; + b /= w; + t = a*b/A[i]; + s = s + t; + if( fabs(t/s) < MACHEP ) + goto done; + k += 1.0; + a *= x + k; + b /= w; + k += 1.0; + } done: - return(s); + return(s); } + +/// @endcond +/** @}*/ diff --git a/core/daceaux.c b/core/daceaux.c index b58092e..7b27b01 100644 --- a/core/daceaux.c +++ b/core/daceaux.c @@ -26,7 +26,7 @@ * Author: Politecnico di Milano */ -/** \addtogroup DACE Core +/** \addtogroup DACE Core * @{ */ @@ -214,7 +214,7 @@ unsigned int daceDecodeExponents(unsigned int ic, const unsigned int no, const u order monomial). Then call it repeatedly with the previous state of p[] and it will put the next monomial in p[] as well as return the order of that monomial. - \note When there are no more monomials left, the funcion will cycle back to + \note When there are no more monomials left, the funcion will cycle back to the constant monomial of order zero. \note This function handles the cases nv=0 and no=0 correctly (i.e. returning 0) */ @@ -250,7 +250,7 @@ unsigned int daceNextMonomial(unsigned int p[], const unsigned int no, const uns order monomial). Then call it repeatedly with the previous state of p[] and it will put the next monomial in p[] as well as return the order of that monomial. - \note When there are no more monomials left, the funcion will cycle back to + \note When there are no more monomials left, the funcion will cycle back to the constant monomial of order zero. \note This function handles the cases nv=0 and no=0 correctly (i.e. returning 0) */ @@ -334,7 +334,7 @@ void daceDecode(const unsigned int jc, unsigned int jj[]) \param[in] cc C array of nmmax monomials \param[in] inc Pointer to DA object to pack the monomials into */ -void dacePack(double cc[], DACEDA *inc) +void dacePack(double cc[restrict], DACEDA *restrict inc) { monomial *ipoc; unsigned int ilmc, illc; @@ -345,34 +345,27 @@ void dacePack(double cc[], DACEDA *inc) if(LIKELY(DACECom.lfi == 0)) #endif { + // provide loop without error checking for big enough target DA if(LIKELY(ilmc >= DACECom.nmmax)) { for(unsigned int i = 0; i < DACECom.nmmax; i++) { -//#define OPTIMIZE // not clear if this actually helps in any way -#ifndef OPTIMIZE - if(fabs(cc[i]) >= DACECom_t.eps) // && (DACECom.ieo[i] <= DACECom_t.nocut) is avoided for performance + if(LIKELY(!(fabs(cc[i]) <= DACECom_t.eps))) // && (DACECom.ieo[i] <= DACECom_t.nocut) is avoided for performance { ic->ii = i; ic->cc = cc[i]; ic++; } cc[i] = 0.0; -#else - ic->ii = i; - ic->cc = cc[i]; - ic += (fabs(cc[i]) >= DACECom_t.eps); - cc[i] = 0.0; -#endif } } else { for(unsigned int i = 0; i < DACECom.nmmax; i++) { - if(fabs(cc[i]) >= DACECom_t.eps && DACECom.ieo[i] <= DACECom_t.nocut) + if(LIKELY(!(fabs(cc[i]) <= DACECom_t.eps) && DACECom.ieo[i] <= DACECom_t.nocut)) // here we remove also cut orders to save as much space as possible { - if(ic >= ipoc+ilmc) + if(UNLIKELY(ic >= ipoc+ilmc)) { daceSetError(__func__, DACE_ERROR, 21); for(unsigned int j = i; j < DACECom.nmmax; j++) cc[j] = 0.0; @@ -391,9 +384,9 @@ void dacePack(double cc[], DACEDA *inc) { for(unsigned int i = 0; i < DACECom.nmmax; i++) { - if((DACECom.ifi[i] != 0) && (fabs(cc[i]) >= DACECom_t.eps) && (DACECom.ieo[i] <= DACECom_t.nocut)) + if(LIKELY((DACECom.ifi[i] != 0) && !(fabs(cc[i]) <= DACECom_t.eps) && (DACECom.ieo[i] <= DACECom_t.nocut))) { - if(ic >= ipoc+ilmc) + if(UNLIKELY(ic >= ipoc+ilmc)) { daceSetError(__func__, DACE_ERROR, 21); for(unsigned int j = 0; j < DACECom.nmmax; j++) cc[j] = 0.0; diff --git a/core/dacebasic.c b/core/dacebasic.c index 71ebbb9..25bc75a 100644 --- a/core/dacebasic.c +++ b/core/dacebasic.c @@ -26,7 +26,7 @@ * Author: Politecnico di Milano */ -/** \addtogroup DACE Core +/** \addtogroup DACE Core * @{ */ @@ -65,7 +65,7 @@ void daceCreateVariable(DACEDA *ina, const unsigned int i, const double ckon) return; } - if(fabs(ckon) < DACECom_t.eps) + if(fabs(ckon) <= DACECom_t.eps) { return; } @@ -117,15 +117,15 @@ void daceCreateMonomial(DACEDA *ina, const unsigned int jj[], const double ckon) return; } - if(fabs(ckon) >= DACECom_t.eps) + if(fabs(ckon) <= DACECom_t.eps) { - ipoa->ii = daceEncode(jj); - ipoa->cc = ckon; - daceSetLength(ina, 1); + daceSetLength(ina, 0); } else { - daceSetLength(ina, 0); + ipoa->ii = daceEncode(jj); + ipoa->cc = ckon; + daceSetLength(ina, 1); } } @@ -323,10 +323,10 @@ void daceSetCoefficient(DACEDA *ina, const unsigned int jj[], const double cjj) daceSetCoefficient0(ina, daceEncode(jj), cjj); } -/*! Extract coefficient of a monomial in a DA object. - \param[in] ina Pointer to DA object to extract monomial coefficient from +/*! Set coefficient of a monomial in a DA object. + \param[in] ina Pointer to DA object to set monomial in \param[in] ic DA coding integer of the monomial to set - \return The coefficient of the given monomial in the DA object + \param[in] cjj Value of the corresponding coefficient */ void daceSetCoefficient0(DACEDA *ina, const unsigned int ic, const double cjj) { @@ -390,7 +390,7 @@ void daceSetCoefficient0(DACEDA *ina, const unsigned int ic, const double cjj) if(insert) { // insert a new monomial before the i-th monomial - if(fabs(cjj) < DACECom_t.eps) return; + if(fabs(cjj) <= DACECom_t.eps) return; if(illa+1 > ilma) { @@ -410,17 +410,17 @@ void daceSetCoefficient0(DACEDA *ina, const unsigned int ic, const double cjj) else { // replace the i-th monomial - if(fabs(cjj) < DACECom_t.eps) + if(!(fabs(cjj) <= DACECom_t.eps)) + { + i->cc = cjj; + } + else { //memmove(i, i+1, (ipoa+illa - i-1)*sizeof(monomial)); for(monomial *ii = i; ii < ipoa+(illa-1); ii++) *ii = *(ii+1); daceSetLength(ina, illa-1); } - else - { - i->cc = cjj; - } } } @@ -506,7 +506,7 @@ void daceCopyFiltering(const DACEDA *ina, DACEDA *inb) { for(monomial *ia = ipoa; ia < ipoa+illa; ia++) { - if(fabs(ia->cc) < DACECom_t.eps || DACECom.ieo[ia->ii] > DACECom_t.nocut) + if(fabs(ia->cc) <= DACECom_t.eps || DACECom.ieo[ia->ii] > DACECom_t.nocut) continue; *ib = *ia; ib++; @@ -516,7 +516,7 @@ void daceCopyFiltering(const DACEDA *ina, DACEDA *inb) { for(monomial *ia = ipoa; ia < ipoa+illa; ia++) { - if(fabs(ia->cc) < DACECom_t.eps || DACECom.ieo[ia->ii] > DACECom_t.nocut) + if(fabs(ia->cc) <= DACECom_t.eps || DACECom.ieo[ia->ii] > DACECom_t.nocut) continue; if(ib >= ipob+ilmb) { @@ -535,7 +535,7 @@ void daceCopyFiltering(const DACEDA *ina, DACEDA *inb) and less than or equal imax. \param[in] ina Pointer to DA object to trim \param[in] imin Minimum order to keep - \param[in] imin Maximum order to keep + \param[in] imax Maximum order to keep \param[in] inc Pointer to DA object to store the truncated result in */ void daceTrim(const DACEDA *ina, const unsigned int imin, const unsigned int imax, DACEDA *inc) @@ -568,8 +568,8 @@ void daceTrim(const DACEDA *ina, const unsigned int imin, const unsigned int ima } /*! Copy monomials from a DA object ina to DA object inb if the same monomial - is non-zero in DA vector inc, while filtering out terms below the current - threshold + is non-zero in DA object inc, while filtering out terms below the current + cutoff \param[in] ina Pointer to DA object to filter \param[in] inb Pointer to DA object to store the filtered result in \param[in] inc Pointer to DA object providing the filter template @@ -595,7 +595,7 @@ void daceFilter(const DACEDA *ina, DACEDA *inb, const DACEDA *inc) ic++; if(ic >= ipoc+illc) break; - if(ia->ii < ic->ii || fabs(ia->cc) < DACECom_t.eps || DACECom.ieo[ia->ii] > DACECom_t.nocut) + if(ia->ii < ic->ii || fabs(ia->cc) <= DACECom_t.eps || DACECom.ieo[ia->ii] > DACECom_t.nocut) continue; *ib = *ia; @@ -611,7 +611,7 @@ void daceFilter(const DACEDA *ina, DACEDA *inb, const DACEDA *inc) ic++; if(ic >= ipoc+illc) break; - if(ia->ii < ic->ii || fabs(ia->cc) < DACECom_t.eps || DACECom.ieo[ia->ii] > DACECom_t.nocut) + if(ia->ii < ic->ii || fabs(ia->cc) <= DACECom_t.eps || DACECom.ieo[ia->ii] > DACECom_t.nocut) continue; if(ib >= ipob+ilmb) diff --git a/core/dacecompat.c b/core/dacecompat.c index 5b9f9f1..a45453f 100644 --- a/core/dacecompat.c +++ b/core/dacecompat.c @@ -26,7 +26,7 @@ * Author: Politecnico di Milano */ -/** \addtogroup DACE Core +/** \addtogroup DACE Core * @{ */ @@ -41,7 +41,7 @@ /*! Return the number of non-zero monomials in a DA object. \param[in] ina Pointer to DA object to get length of - \return Number of non-zero monomials + \param[out] size Number of non-zero monomials \deprecated Has been renamed to daceGetLength() \sa daceGetLength() */ diff --git a/core/daceerror.c b/core/daceerror.c index a7f4f05..d07e041 100644 --- a/core/daceerror.c +++ b/core/daceerror.c @@ -26,7 +26,7 @@ * Author: Politecnico di Milano */ -/** \addtogroup DACE Core +/** \addtogroup DACE Core * @{ */ @@ -106,30 +106,30 @@ void daceClearError() *DACEDbg.msg = '\0'; } -/*! This subroutine serves as an error handler for errors within the dace. +/*! This subroutine serves as an error handler for errors within the dace. It is intended mostly for development and debugging. More descriptive error messages should be displayed by the user interface. - + The error codes are defined as XYY with X indicating the severity and YY corresponding to the actual error code - + Severity Levels X 1 = Info: Informative, no action required - + 3 = Warning: Serious, possibly incorrect use of DACE routines - + 6 = Error: Recoverable, result may not be correct or assumptions have been made - + 9 = Error: Unrecoverable, new call to DACEINI is required to reinitialize DACE, DACE objects are no longer valid - + 10 = Critical: Crash in the DACE, just printing as much as possible and dying - + Currently used error codes XYY are defined in daceerror.h - + \param[in] c an error string representing the error \param[in] ix is the error severity code \param[in] iyy is the error code diff --git a/core/daceeval.c b/core/daceeval.c index 3f19491..831db19 100644 --- a/core/daceeval.c +++ b/core/daceeval.c @@ -26,7 +26,7 @@ * Author: Politecnico di Milano */ -/** \addtogroup DACE Core +/** \addtogroup DACE Core * @{ */ @@ -160,7 +160,7 @@ void daceReplaceVariable(const DACEDA *ina, const unsigned int from, const unsig if(from < 1 || from > DACECom.nvmax || to < 1 || to > DACECom.nvmax) { - daceSetError(__func__, DACE_ERROR, 24); + daceSetError(__func__, DACE_ERROR, 24); daceCreateConstant(inc, 0.0); return; } @@ -216,7 +216,7 @@ void daceScaleVariable(const DACEDA *ina, const unsigned int nvar, const double if(nvar < 1 || nvar > DACECom.nvmax) { - daceSetError(__func__, DACE_ERROR, 24); + daceSetError(__func__, DACE_ERROR, 24); daceCreateConstant(inc, 0.0); return; } @@ -278,7 +278,7 @@ void daceTranslateVariable(const DACEDA *ina, const unsigned int nvar, const dou if(nvar < 1 || nvar > DACECom.nvmax) { - daceSetError(__func__, DACE_ERROR, 24); + daceSetError(__func__, DACE_ERROR, 24); daceCreateConstant(inc, 0.0); return; } @@ -294,7 +294,7 @@ void daceTranslateVariable(const DACEDA *ina, const unsigned int nvar, const dou double *cc = dacecalloc(DACECom.nmmax, sizeof(double)); double *powa = dacecalloc(DACECom.nomax+1, sizeof(double)); double *powc = dacecalloc(DACECom.nomax+1, sizeof(double)); - double *binomial = dacecalloc((DACECom.nomax+1)*(DACECom.nomax+1), sizeof(double)); + double *binomial = dacecalloc(((size_t)DACECom.nomax+1)*((size_t)DACECom.nomax+1), sizeof(double)); #endif // precompute powers of a and c @@ -475,7 +475,7 @@ void daceEvalTree(const DACEDA *das[], const unsigned int count, double ac[], un exit(1); } #endif - + #if DACE_MEMORY_MODEL != DACE_MEMORY_STATIC dacefree(nc); dacefree(p); diff --git a/core/daceinit.c b/core/daceinit.c index 3b59d03..3ecb712 100644 --- a/core/daceinit.c +++ b/core/daceinit.c @@ -26,7 +26,7 @@ * Author: Politecnico di Milano */ -/** \addtogroup DACE Core +/** \addtogroup DACE Core * @{ */ @@ -46,8 +46,12 @@ /*! Set up the ordering and addressing arrays in the common data structure and initialize DA memory. \note MUST BE CALLED BEFORE ANY OTHER DA ROUTINE CAN BE USED. + \note Also initializes the truncation order to the maximum computation order + and disables the DA epsilon cutoff by setting it to 0.0. \param[in] no order of the Taylor polynomials; \param[in] nv number of variables considered. + \sa daceTruncationOrder + \sa daceSetEpsilon */ void daceInitialize(unsigned int no, unsigned int nv) { @@ -59,7 +63,7 @@ void daceInitialize(unsigned int no, unsigned int nv) daceSetError(__func__, DACE_INFO, 67); no = 1; } - + if(nv < 1) { daceSetError(__func__, DACE_INFO, 68); nv = 1; @@ -72,7 +76,7 @@ void daceInitialize(unsigned int no, unsigned int nv) DACECom.epsmac = DACECom.epsmac/2.0; } DACECom.epsmac = DACECom.epsmac*2.0; - + // Reset memory, purging all previous DA objects, if any. daceFreeMemory(); @@ -146,7 +150,7 @@ void daceInitialize(unsigned int no, unsigned int nv) i++; // increase total monomial count } while((no2 = daceNextOrderedMonomial(p2, no-no1, nv2)) > 0); } while((no1 = daceNextOrderedMonomial(p1, no, nv1)) > 0); - + // free memory #if DACE_MEMORY_MODEL != DACE_MEMORY_STATIC dacefree(p1); @@ -185,7 +189,13 @@ void daceInitialize(unsigned int no, unsigned int nv) \note The main thread must call daceInitialize once before spawning new threads. All spawned threads must then call daceInitializeThread to initialize the thread before performing any other operations. + \note daceInitialize MUST NOT be called again by any thread while other threads + are active. + \note Also initializes the truncation order to the maximum computation order + and disables the DA epsilon cutoff by setting it to 0.0. \sa daceInitialize + \sa daceTruncationOrder + \sa daceSetEpsilon */ void daceInitializeThread() { @@ -209,11 +219,13 @@ void daceCleanupThread() } /*! Set up thread local data structures without resetting error. - */ + Also initializes the truncation order to the maximum computation order + and disables the DA epsilon cutoff by setting it to 0.0. +*/ void daceInitializeThread0() { - // Set cutoff to (arbitrarily) be 2^10 smaller than machine epsilon - DACECom_t.eps = DACECom.epsmac/pown(2.0, 10); + // Set initial cutoff to zero, disabling the optimization + DACECom_t.eps = 0.0; // Set truncation order to full order DACECom_t.nocut = DACECom.nomax; @@ -244,9 +256,15 @@ void daceGetVersion(int *imaj, int *imin, int *ipat) } /*! Set cutoff value to eps and return the previous value. - \param[in] eps New cutoff value below which coefficients may be flushed to + \param[in] eps New cutoff value at or below which coefficients can be flushed to zero for efficiency purposes. \return The previous value of the cutoff + \note This feature can have severe unintended consequences if used incorrectly! + Flushing occurs for any intermediate result also within the DACE, and can result + in wrong solutions whenever DA coefficients become very small relative to epsilon. + For example, division by a large DA divisor can cause the (internally calculated) + multiplicative inverse to be entirely flushed to zero, resulting in a zero DA + quotient independently of the size of the dividend. \sa daceGetEpsilon */ double daceSetEpsilon(const double eps) @@ -266,7 +284,7 @@ double daceGetEpsilon() } /*! Get machine epsilon value. - \return The experimentally computed machine epsilon + \return The experimentally determined machine epsilon */ double daceGetMachineEpsilon() { diff --git a/core/daceio.c b/core/daceio.c index a5b838c..8d23a82 100644 --- a/core/daceio.c +++ b/core/daceio.c @@ -26,7 +26,7 @@ * Author: Politecnico di Milano */ -/** \addtogroup DACE Core +/** \addtogroup DACE Core * @{ */ @@ -461,7 +461,7 @@ unsigned int daceBlobSize(const void *blob) /*! Import a DA object in a binary format. \param[in] blob Pointer to memory where the data is stored - \param[in] ina The DA object to import into + \param[in] inc The DA object to import into \note This routine will silently truncate orders above the currently set maximum computation order as well as any extra variables present. */ @@ -471,7 +471,7 @@ void daceImportBlob(const void *blob, DACEDA *inc) if(data->magic != DACE_BINARY_MAGIC) { - daceSetError(__func__, DACE_ERROR, 31); + daceSetError(__func__, DACE_ERROR, 31); daceCreateConstant(inc, 0.0); return; } @@ -485,7 +485,7 @@ void daceImportBlob(const void *blob, DACEDA *inc) #if DACE_MEMORY_MODEL == DACE_MEMORY_STATIC if(nv > DACE_STATIC_NVMAX) { - daceSetError(__func__, DACE_ERROR, 23); + daceSetError(__func__, DACE_ERROR, 23); daceCreateConstant(inc, 0.0); return; } diff --git a/core/dacemath.c b/core/dacemath.c index 179d77b..a56c89f 100755 --- a/core/dacemath.c +++ b/core/dacemath.c @@ -126,7 +126,7 @@ void daceMultiply(const DACEDA *ina, const DACEDA *inb, DACEDA *inc) ipbeg[i] = emb + daceCountMonomials(i - 1, DACECom.nvmax); } #else - static DACE_THREAD_LOCAL double *cc = NULL; + static DACE_THREAD_LOCAL double *restrict cc = NULL; static DACE_THREAD_LOCAL extended_monomial *emb = NULL; static DACE_THREAD_LOCAL extended_monomial **ipbeg = NULL; static DACE_THREAD_LOCAL extended_monomial **ipend = NULL; @@ -159,7 +159,7 @@ void daceMultiply(const DACEDA *ina, const DACEDA *inb, DACEDA *inc) daceVariableInformation(inb, &ipob, &ilmb, &illb); // sort so that ina is the short DA vector - if(illa>illb) + if(illa > illb) { unsigned int t1; t1 = illb; illb = illa; illa = t1; @@ -171,10 +171,10 @@ void daceMultiply(const DACEDA *ina, const DACEDA *inb, DACEDA *inc) for(unsigned int i = 0; i <= DACECom_t.nocut; i++) ipend[i] = ipbeg[i]; // sort vector b by order - for(monomial *ib = ipob; ib < ipob+illb; ib++) + for(const monomial *ib = ipob; ib < ipob+illb; ib++) { const unsigned int noib = DACECom.ieo[ib->ii]; - if(noib > DACECom_t.nocut) continue; + if(UNLIKELY(noib > DACECom_t.nocut)) continue; ipend[noib]->i1 = DACECom.ie1[ib->ii]; ipend[noib]->i2 = DACECom.ie2[ib->ii]; ipend[noib]->cc = ib->cc; @@ -182,7 +182,7 @@ void daceMultiply(const DACEDA *ina, const DACEDA *inb, DACEDA *inc) } // perform actual multiplication - for(monomial *ia = ipoa; ia < ipoa+illa; ia++) + for(const monomial *ia = ipoa; ia < ipoa+illa; ia++) { const unsigned int i1ia = DACECom.ie1[ia->ii]; const unsigned int i2ia = DACECom.ie2[ia->ii]; @@ -191,7 +191,7 @@ void daceMultiply(const DACEDA *ina, const DACEDA *inb, DACEDA *inc) //#pragma omp parallel for for(int noib = DACECom_t.nocut-DACECom.ieo[ia->ii]; noib >= 0; noib--) { - for(extended_monomial *ib = ipbeg[noib]; ib < ipend[noib]; ib++) + for(const extended_monomial *ib = ipbeg[noib]; ib < ipend[noib]; ib++) { const unsigned int ic = DACECom.ia1[i1ia+ib->i1] + DACECom.ia2[i2ia+ib->i2]; cc[ic] += ccia*ib->cc; @@ -259,7 +259,7 @@ void daceDivide(const DACEDA *ina, const DACEDA *inb, DACEDA *inc) /*! Square a DA object. \param[in] ina Pointer to the DA object to square - \param[out] inc Pointer to the DA object to store the result in + \param[out] inb Pointer to the DA object to store the result in \note This routine is aliasing safe, i.e. inc can be the same as ina. */ void daceSquare(const DACEDA *ina, DACEDA *inb) @@ -328,12 +328,12 @@ void daceMultiplyDouble(const DACEDA *ina, const double ckon, DACEDA *inb) continue; const double c = ia->cc*ckon; - if(fabs(c) < DACECom_t.eps) - continue; - - ib->cc = c; - ib->ii = ia->ii; - ib++; + if(!(fabs(c) <= DACECom_t.eps)) + { + ib->cc = c; + ib->ii = ia->ii; + ib++; + } } } else @@ -345,17 +345,17 @@ void daceMultiplyDouble(const DACEDA *ina, const double ckon, DACEDA *inb) continue; const double c = ia->cc*ckon; - if(fabs(c) < DACECom_t.eps) - continue; - - if(ib >= ibmax) + if(!(fabs(c) <= DACECom_t.eps)) { - daceSetError(__func__, DACE_ERROR, 21); - break; + if(ib >= ibmax) + { + daceSetError(__func__, DACE_ERROR, 21); + break; + } + ib->cc = c; + ib->ii = ia->ii; + ib++; } - ib->cc = c; - ib->ii = ia->ii; - ib++; } } @@ -568,7 +568,7 @@ void daceDifferentiate(const unsigned int idif, const DACEDA *ina, DACEDA *inc) } /*! Integral of DA object with respect to a given independent variable. - \param[in] idif Number of the independent variable with respect to which the + \param[in] iint Number of the independent variable with respect to which the integral is taken \param[in] ina Pointer to the DA object to operate on \param[out] inc Pointer to the DA object to store the result in @@ -608,16 +608,17 @@ void daceIntegrate(const unsigned int iint, const DACEDA *ina, DACEDA *inc) const unsigned int ic2 = DACECom.ie2[i->ii]; const unsigned int ipow = (ic2/idiv)%ibase; const double ccc = i->cc/(ipow+1); - if(fabs(ccc) < DACECom_t.eps) - continue; - if(ic >= icmax) + if(!(fabs(ccc) <= DACECom_t.eps)) { - daceSetError(__func__, DACE_ERROR, 21); - break; + if(ic >= icmax) + { + daceSetError(__func__, DACE_ERROR, 21); + break; + } + ic->ii = DACECom.ia1[ic1] + DACECom.ia2[ic2+idiv]; + ic->cc = ccc; + ic = ic+1; } - ic->ii = DACECom.ia1[ic1] + DACECom.ia2[ic2+idiv]; - ic->cc = ccc; - ic = ic+1; } } else @@ -630,16 +631,17 @@ void daceIntegrate(const unsigned int iint, const DACEDA *ina, DACEDA *inc) const unsigned int ic2 = DACECom.ie2[i->ii]; const unsigned int ipow = (ic1/idiv)%ibase; const double ccc = i->cc/(ipow+1); - if(fabs(ccc) < DACECom_t.eps) - continue; - if(ic >= icmax) + if(!(fabs(ccc) <= DACECom_t.eps)) { - daceSetError(__func__, DACE_ERROR, 21); - break; + if(ic >= icmax) + { + daceSetError(__func__, DACE_ERROR, 21); + break; + } + ic->ii = DACECom.ia1[ic1+idiv] + DACECom.ia2[ic2]; + ic->cc = ccc; + ic = ic+1; } - ic->ii = DACECom.ia1[ic1+idiv] + DACECom.ia2[ic2]; - ic->cc = ccc; - ic = ic+1; } } @@ -732,7 +734,7 @@ void dacePowerDouble(const DACEDA *ina, const double p, DACEDA *inc) /*! Raise a DA object to the p-th integer power. \param[in] ina Pointer to the DA object to operate on - \param[in] p Power to which to raise the DA object + \param[in] np Power to which to raise the DA object \param[out] inc Pointer to the DA object to store the result in \note This routine is aliasing safe, i.e. inc can be the same as ina. */ @@ -1230,7 +1232,7 @@ void daceArcTangent(const DACEDA *ina, DACEDA *inc) /*! Arctangent of ina/inb with proper sign in [-pi, pi]. This function follows the C standard atan2(y,x) function syntax. \param[in] ina Pointer to the first DA object to operate on - \param[in] ina Pointer to the second DA object to operate on + \param[in] inb Pointer to the second DA object to operate on \param[out] inc Pointer to the DA object to store the result in \note This routine is aliasing safe, i.e. inc can be the same as ina. */ @@ -1341,11 +1343,25 @@ void daceHyperbolicCosine(const DACEDA *ina, DACEDA *inc) void daceHyperbolicTangent(const DACEDA *ina, DACEDA *inc) { DACEDA itemp; + const double a0 = daceGetConstant(ina); daceAllocateDA(&itemp, 0); - daceHyperbolicSine(ina, &itemp); - daceHyperbolicCosine(ina, inc); - daceDivide(&itemp, inc, inc); + if(a0 > 0.0) + { + daceMultiplyDouble(ina, -2.0, &itemp); + daceExponential(&itemp, &itemp); + daceAddDouble(&itemp, 1.0, inc); + daceDoubleSubtract(&itemp, 1.0, &itemp); + daceDivide(&itemp, inc, inc); + } + else + { + daceMultiplyDouble(ina, 2.0, &itemp); + daceExponential(&itemp, &itemp); + daceAddDouble(&itemp, 1.0, inc); + daceAddDouble(&itemp, -1.0, &itemp); + daceDivide(&itemp, inc, inc); + } daceFreeDA(&itemp); } @@ -1859,7 +1875,8 @@ void daceEvaluateScaledModifiedBesselFunction(const DACEDA *ina, const double bz } /*! Compute the partial Logarithmic Gamma function of a DA object (without constant part). - \param[in] ina Pointer to the DA object to operate on (constant part != 0, -1, -2, ...) + \param[in] ina Pointer to the DA object to operate on + \param[in] a0 Constant part \param[out] inc Pointer to the DA object to store the result in \note This routine is aliasing safe, i.e. inc can be the same as ina. \note No argument checking is performed to ensure values are within allowable range. @@ -2043,7 +2060,7 @@ void daceWeightedSum(const DACEDA *ina, const double afac, const DACEDA *inb, co if(DACECom.ieo[ja] <= DACECom_t.nocut) { const double ccc = ia->cc*afac + ib->cc*bfac; - if(fabs(ccc) >= DACECom_t.eps) + if(!(fabs(ccc) <= DACECom_t.eps)) { if(ic >= icmax) { @@ -2067,7 +2084,7 @@ void daceWeightedSum(const DACEDA *ina, const double afac, const DACEDA *inb, co if(DACECom.ieo[ja] <= DACECom_t.nocut) { const double ccc = ia->cc*afac; - if(fabs(ccc) >= DACECom_t.eps) + if(!(fabs(ccc) <= DACECom_t.eps)) { if(ic >= icmax) { @@ -2090,7 +2107,7 @@ void daceWeightedSum(const DACEDA *ina, const double afac, const DACEDA *inb, co if(DACECom.ieo[jb] <= DACECom_t.nocut) { const double ccc = ib->cc*bfac; - if(fabs(ccc) >= DACECom_t.eps) + if(!(fabs(ccc) <= DACECom_t.eps)) { if(ic >= icmax) { @@ -2131,7 +2148,7 @@ void daceWeightedSum(const DACEDA *ina, const double afac, const DACEDA *inb, co if(DACECom.ieo[is->ii] <= DACECom_t.nocut) { const double ccc = is->cc*fac; - if(fabs(ccc) >= DACECom_t.eps) + if(!(fabs(ccc) <= DACECom_t.eps)) { if(ic >= icmax) { diff --git a/core/dacememory.c b/core/dacememory.c index 2e4d48e..2b6dcc6 100644 --- a/core/dacememory.c +++ b/core/dacememory.c @@ -26,7 +26,7 @@ * Author: Politecnico di Milano */ -/** \addtogroup DACE Core +/** \addtogroup DACE Core * @{ */ @@ -71,8 +71,8 @@ typedef struct dvariable { // Global memory array typedef struct dmem { - monomial *mem; - variable *var; + monomial *mem; + variable *var; unsigned int nda, mda, nst, lda, lst; } dacemem; @@ -111,7 +111,7 @@ void daceReallocateMemory(const unsigned int nvar, const unsigned int nmem) daceSetError(__func__, DACE_PANIC, 2); exit(1); } - + // static memory blocks static monomial dace_static_mem[DACE_STATIC_MEM_SIZE]; static variable dace_static_var[DACE_STATIC_VAR_SIZE]; @@ -290,7 +290,7 @@ void daceMemoryDump() /*! Extract internal information about a DA object. \param[in] inc Pointer to the DA object to extract information from \param[out] ipoc Pointer to an array of monomials allocated for this variable - \param[out] ilmc Pointer where to store the maximum number of monomials allocated in this DA object + \param[out] ilmc Pointer where to store the maximum number of monomials allocated in this DA object \param[out] illc Pointer where to store the currently used length of this DA object */ void daceVariableInformation(const DACEDA *inc, monomial **ipoc, unsigned int *ilmc, unsigned int *illc) @@ -414,7 +414,7 @@ void daceMemoryDump() /*! Extract internal information about a DA object. \param[in] inc Pointer to the DA object to extract information from \param[out] ipoc Pointer to an array of monomials allocated for this variable - \param[out] ilmc Pointer where to store the maximum number of monomials allocated in this DA object + \param[out] ilmc Pointer where to store the maximum number of monomials allocated in this DA object \param[out] illc Pointer where to store the currently used length of this DA object */ void daceVariableInformation(const DACEDA *inc, monomial **ipoc, unsigned int *ilmc, unsigned int *illc) @@ -438,7 +438,7 @@ void daceSetLength(DACEDA *inc, const size_t len) { if(UNLIKELY(inc->max < len)) { - daceSetError(__func__, DACE_PANIC, 7); + daceSetError(__func__, DACE_PANIC, 7); exit(1); // catastrophic error because we may have written past the end of the variable and contaminated memory there } diff --git a/core/dacenorm.c b/core/dacenorm.c index 75d1877..64f5d16 100644 --- a/core/dacenorm.c +++ b/core/dacenorm.c @@ -26,7 +26,7 @@ * Author: Politecnico di Milano */ -/** \addtogroup DACE Core +/** \addtogroup DACE Core * @{ */ @@ -108,7 +108,7 @@ void daceOrderedNorm(const DACEDA *ina, const unsigned int ivar, const unsigned if(ivar > DACECom.nvmax) { - daceSetError(__func__, DACE_ERROR, 24); + daceSetError(__func__, DACE_ERROR, 24); return; } @@ -197,6 +197,7 @@ void daceOrderedNorm(const DACEDA *ina, const unsigned int ivar, const unsigned 0 = max norm 1 = sum norm >1 = corresponding vector norm + \param[in] nc Maximum order to estimate \param[out] c C array of length nc+1 containing the grouped estimates \param[out] err C array of length min(nc, nomax)+1 containing the residuals of the exponential fit at each order. If NULL is passed in, no residuals @@ -228,7 +229,7 @@ void daceEstimate(const DACEDA *ina, const unsigned int ivar, const unsigned int double xtx[2][2] = {{0.0}, {0.0}}; for(unsigned int i = 1; i <= DACECom.nomax; i++) { - if(onorm[i] > DACECom_t.eps) + if(!(onorm[i] <= DACECom_t.eps)) { xtx[0][0] += i*i; xtx[0][1] -= i; diff --git a/core/include/dace/daceaux.h b/core/include/dace/daceaux.h index 2c4a7c4..9ee1198 100644 --- a/core/include/dace/daceaux.h +++ b/core/include/dace/daceaux.h @@ -81,14 +81,14 @@ typedef struct ddbg { // Basic memory structure of a monomial typedef struct dmonomial { double cc; - unsigned int ii; + unsigned int ii; } monomial; // Basic memory structure of an extended monomial (these are written and read from binary files, so enforce no padding!) #pragma pack(push,1) typedef struct dextendedmonomial { - unsigned int i1, i2; - double cc; + unsigned int i1, i2; + double cc; } extended_monomial; #pragma pack(pop) @@ -143,7 +143,7 @@ unsigned int daceNextOrderedMonomial(unsigned int p[], const unsigned int no, co // internal routines void daceInitializeThread0(); void daceSetError(const char *c, const unsigned int ix, const unsigned int iyy); -void dacePack(double cc[], DACEDA *inc); +void dacePack(double cc[restrict], DACEDA *restrict inc); void daceMultiplicativeInverse0(const DACEDA *ina, DACEDA *inc, const double a0); int BesselWrapper(const double x, const int n0, const int n1, const int type, double *bz); int ModifiedBesselWrapper(const double x, const int n0, const int n1, const int type, double *bz); diff --git a/core/include/dace/dacecompat.h b/core/include/dace/dacecompat.h index 5166f39..f7d98e4 100644 --- a/core/include/dace/dacecompat.h +++ b/core/include/dace/dacecompat.h @@ -31,7 +31,7 @@ Where necessary, small shims with the old interface are provided in dacecompat.c. */ -/** \addtogroup DACE Core +/** \addtogroup DACE Core * @{ */ #ifndef DINAMICA_DACECOMPAT_H_ @@ -52,7 +52,7 @@ #define dacegetnot() daceGetTruncationOrder() #define dacesetnot(fnot) daceSetTruncationOrder(fnot) -//#define dacegeterr() daceGetError() // do not redefine this routine, it causes problems with the inlining of a function of the same name we do in C++ interface! +//#define dacegeterr() daceGetError() // do not redefine this routine, it causes problems with the inlining of a function of the same name we do in C++ interface! #define dacegetxerr() daceGetErrorX() #define dacegetyyerr() daceGetErrorYY() #define daceclrerr() daceClearError() diff --git a/core/include/dace/dacecore.h b/core/include/dace/dacecore.h index 07618fc..d3fbc0c 100644 --- a/core/include/dace/dacecore.h +++ b/core/include/dace/dacecore.h @@ -30,7 +30,7 @@ User interface header for DACE core library. Includes all relevant headers with public interfaces to the DACE core. */ -/** \addtogroup DACE Core +/** \addtogroup DACE Core * @{ */ diff --git a/core/include/dace/dacecore_s.h b/core/include/dace/dacecore_s.h index 85dd06a..aa438d1 100644 --- a/core/include/dace/dacecore_s.h +++ b/core/include/dace/dacecore_s.h @@ -31,7 +31,7 @@ Includes all relevant headers with public interfaces to the DACE core with the correct API decorations for Windows static linking. */ -/** \addtogroup DACE Core +/** \addtogroup DACE Core * @{ */ diff --git a/core/include/dace/daceerror.h b/core/include/dace/daceerror.h index 85e935e..e227d44 100644 --- a/core/include/dace/daceerror.h +++ b/core/include/dace/daceerror.h @@ -29,7 +29,7 @@ #ifndef DINAMICA_DACEERROR_H #define DINAMICA_DACEERROR_H -/** \addtogroup DACE Core +/** \addtogroup DACE Core * @{ */ typedef struct { @@ -100,7 +100,7 @@ DACE_API const errstrings DACEerr[] = { { 57, ""}, { 58, ""}, { 59, ""}, - { 60, ""}, + { 60, ""}, { 161, "Free or invalid variable"}, { 162, "Truncation order too high"}, { 163, "Inacurate estimate"}, diff --git a/interfaces/cxx/DA.cpp b/interfaces/cxx/DA.cpp index cd69548..7765539 100644 --- a/interfaces/cxx/DA.cpp +++ b/interfaces/cxx/DA.cpp @@ -66,7 +66,7 @@ void DA::init(const unsigned int ord, const unsigned int nvar) { */ try { checkVersion(); - } catch(DACEException ex) { + } catch(DACEException const& ex) { std::cerr << ex << std::endl; std::terminate(); } daceInitialize(ord,nvar); @@ -85,8 +85,7 @@ bool DA::isInitialized() { void DA::version(int &maj, int &min, int &patch) { /*! Return the major and minor version number of the DACE - along with the "COSY flag" indicating if the core library - is COSY or DACE. + along with the patch level of the library. \param[out] maj major DACE version number; \param[out] min minor DACE version number; \param[out] patch patch level of DACE version. @@ -366,16 +365,16 @@ double DA::getCoefficient(const std::vector &jj) const { double coeff; const unsigned int nvar = daceGetMaxVariables(); - if(jj.size() >= nvar) - { - coeff = daceGetCoefficient(m_index, jj.data()); - } - else - { - std::vector temp(jj); - temp.resize(nvar, 0); - coeff = daceGetCoefficient(m_index, temp.data()); - } + if(jj.size() >= nvar) + { + coeff = daceGetCoefficient(m_index, jj.data()); + } + else + { + std::vector temp(jj); + temp.resize(nvar, 0); + coeff = daceGetCoefficient(m_index, temp.data()); + } if(daceGetError()) DACEException(); @@ -391,16 +390,16 @@ void DA::setCoefficient(const std::vector &jj, const double coeff) // check arguments const unsigned int nvar = daceGetMaxVariables(); - if(jj.size() >= nvar) - { - daceSetCoefficient(m_index, jj.data(), coeff); - } - else - { - std::vector temp(jj); - temp.resize(nvar, 0); - daceSetCoefficient(m_index, temp.data(), coeff); - } + if(jj.size() >= nvar) + { + daceSetCoefficient(m_index, jj.data(), coeff); + } + else + { + std::vector temp(jj); + temp.resize(nvar, 0); + daceSetCoefficient(m_index, temp.data(), coeff); + } if(daceGetError()) DACEException(); } @@ -421,17 +420,17 @@ Monomial DA::getMonomial(const unsigned int npos) const{ \sa DA::getMonomials */ Monomial m; - getMonomial(npos, m); + getMonomial(npos, m); return m; } void DA::getMonomial(const unsigned int npos, Monomial &m) const { /*! Return the Monomial corresponding to the non-zero coefficient at the given - position in the DA object (monomials use one based indexing!). + position in the DA object (monomials use one based indexing!). \param[in] npos position within the DA object. The ordering of the Monomials - within a DA object is implementation dependent and does not correspond - to the order in which Monomials are listed in the ASCII output routines. + within a DA object is implementation dependent and does not correspond + to the order in which Monomials are listed in the ASCII output routines. \param[out] m the monomial object in which to store the corresponding monomial. \throw DACE::DACEException \sa Monomial @@ -831,17 +830,17 @@ DA operator/(const double c, const DA &da){ DA DA::multiplyMonomials(const DA &da) const { /*! Multiply the DA vector with another DA vector da monomial by monomial. - This is the equivalent of coefficient-wise multiplication (like in DA addition). + This is the equivalent of coefficient-wise multiplication (like in DA addition). \param[in] da DA vector to multiply with coefficient-wise \return A new DA object containing the result of the operation. \throw DACE::DACEException \sa DA::evalMonomial */ - DA temp; - daceMultiplyMonomials(m_index, da.m_index, temp.m_index); - if (daceGetError()) DACEException(); + DA temp; + daceMultiplyMonomials(m_index, da.m_index, temp.m_index); + if (daceGetError()) DACEException(); - return temp; + return temp; } DA DA::divide(const unsigned int var, const unsigned int p) const{ @@ -1655,14 +1654,14 @@ double DA::evalMonomials(const DA &values) const { /*! Evaluates the DA vector using the coefficients in argument values as the values for each monomial. This is equivalent to a monomial-wise dot product of two DA vectors. \param[in] values DA vector containing the values of each monomial - \return The resutl of the evaluation. + \return The result of the evaluation. \throw DACE::DACEException \sa DA::multiplyMonomial */ - const double res = daceEvalMonomials(m_index, values.m_index); - if (daceGetError()) DACEException(); + const double res = daceEvalMonomials(m_index, values.m_index); + if (daceGetError()) DACEException(); - return res; + return res; } DA DA::replaceVariable(const unsigned int from, const unsigned int to, const double val) const{ @@ -1769,16 +1768,16 @@ std::istream& operator>>(std::istream &in, DA &da){ \throw DACE::DACEException \note When using binary IO operations, make sure the stream is opened in ios_base::binary mode! Some C++ libraries are known to mangle the input otherwise which will break the ability to read binary DA objects. - Setting the binary flag for all IO (also text based) does not affect the output and is recommended. + Setting the binary flag for all IO (also text based) does not affect the output and is recommended. \sa DA::fromString */ storedDA sda(in); if(sda.isValid()) - { + { da = sda; // automatically cast sDA to DA return in; - } + } else { const std::string endstr = "------------------------------------------------"; // from daceio.c @@ -1898,9 +1897,9 @@ DA DA::fromString(const std::vector &str){ DA DA::read(std::istream &is){ /*! Read a binary representation of a DA from is. \throw DACE::DACEException - \note When using binary IO operations, make sure the stream is opened in ios_base::binary mode! - Some C++ libraries are known to mangle the input otherwise which will break the ability to read binary DA objects. - Setting the binary flag for all IO (also text based) does not affect the output and is recommended. + \note When using binary IO operations, make sure the stream is opened in ios_base::binary mode! + Some C++ libraries are known to mangle the input otherwise which will break the ability to read binary DA objects. + Setting the binary flag for all IO (also text based) does not affect the output and is recommended. */ storedDA sda(is); return sda; // automatically cast to DA @@ -1948,8 +1947,9 @@ DA divide(const DA &da, const unsigned int var, const unsigned int p){ /*! Divide by independent variable var raised to power p. The result is copied in a new DA object. \param[in] da DA object. - \param[in] i variable with respect to which the derivative is calculated. - \return A new DA object containing the result of the derivation. + \param[in] var variable number to divide by. + \param[in] p power of variable var to divide by. + \return A new DA object containing the result of the division. \throw DACE::DACEException \sa DA::divide */ @@ -2148,9 +2148,10 @@ DA icrt(const DA &da){ return da.icrt();} DA hypot(const DA &da1, const DA &da2){ -/*! Compute the hypotenuse (sqrt(a*a + b*b)) of two DA objects. +/*! Compute the hypotenuse (sqrt(da1*da1 + da2*da2)) of two DA objects. The result is copied in a new DA object. - \param[in] da a given DA object. + \param[in] da1 first DA object. + \param[in] da2 second DA object. \return A new DA object containing the result of the operation. \throw DACE::DACEException \sa DA::hypot @@ -2194,7 +2195,6 @@ DA log10(const DA &da){ /*! Compute the 10 based logarithm of a DA object. The result is copied in a new DA object. \param[in] da a given DA object. - \param[in] b base with respect to which the logarithm is computed (base 10 set as default base). \return A new DA object containing the result of the operation. \throw DACE::DACEException \sa DA::log10 @@ -2205,7 +2205,6 @@ DA log2(const DA &da){ /*! Compute the 2 based logarithm of a DA object. The result is copied in a new DA object. \param[in] da a given DA object. - \param[in] b base with respect to which the logarithm is computed (base 10 set as default base). \return A new DA object containing the result of the operation. \throw DACE::DACEException \sa DA::log2 @@ -2488,6 +2487,7 @@ DA LogGammaFunction(const DA &da){ DA PsiFunction(const unsigned int n, const DA &da){ /*! Compute the n-th Psi function of a DA object. The result is copied in a new DA object. + \param[in] n order of the Psi function to compute. \param[in] da a given DA object. \return A new DA object containing the result of the operation. \throw DACE::DACEException @@ -2687,11 +2687,11 @@ const unsigned int storedDA::headerSize = daceBlobSize(NULL); // create new storedDA from an existing DA storedDA::storedDA(const DA &da){ - unsigned int len; + unsigned int len; - daceExportBlob(da.m_index, NULL, len); - resize(len); - daceExportBlob(da.m_index, data(), len); + daceExportBlob(da.m_index, NULL, len); + resize(len); + daceExportBlob(da.m_index, data(), len); if(daceGetError()) DACEException(); } @@ -2703,11 +2703,11 @@ storedDA::storedDA(const std::vector &data) : std::vector(data){ storedDA::storedDA(std::istream &is) : std::vector(storedDA::headerSize){ // read blob header is.read(data(), headerSize); - if(is.gcount() != headerSize) - { + if(is.gcount() != headerSize) + { resize((size_t)is.gcount()); return; - } + } // check validity of blob header and read the rest const unsigned int len = daceBlobSize(data()); @@ -2716,15 +2716,15 @@ storedDA::storedDA(std::istream &is) : std::vector(storedDA::headerSize){ return; } else if(len > headerSize) - { + { resize(len); - is.read(data()+headerSize, len-headerSize); - if(is.gcount() != (len-headerSize)) - { + is.read(data()+headerSize, len-headerSize); + if(is.gcount() != (len-headerSize)) + { resize((size_t)(headerSize+is.gcount())); - return; - } - } + return; + } + } } // return if this storedDA data appears to be valid diff --git a/interfaces/cxx/DACEException.cpp b/interfaces/cxx/DACEException.cpp index b1e3ef1..ddd2464 100644 --- a/interfaces/cxx/DACEException.cpp +++ b/interfaces/cxx/DACEException.cpp @@ -161,14 +161,14 @@ namespace DACE{ const int id = m_x*100+m_yy; int i = 0; std::stringstream s; - + if(m_x > 10) { for(i=length-1; (i>0)&&(DACEerr[i].ID != id); i--); - + s << DACEerr[i].msg << " (ID: " << DACEerr[i].ID << ")" ; - + } else { @@ -230,7 +230,7 @@ namespace DACE{ /******************************************************************************** * Friend functions *********************************************************************************/ - std::ostream& operator<< (std::ostream &out, const DACEException &ex){ + std::ostream& operator<< (std::ostream &out, const DACEException &ex){ /*! Overload of std::operator<< in iostream. \param[in] out standard output stream. \param[in] ex Exception to be printed in the stream diff --git a/interfaces/cxx/compiledDA.cpp b/interfaces/cxx/compiledDA.cpp index d2400e2..fae366b 100644 --- a/interfaces/cxx/compiledDA.cpp +++ b/interfaces/cxx/compiledDA.cpp @@ -153,7 +153,7 @@ template<> void compiledDA::eval(const std::vector &args, std::vector &r daceAllocateDA(tmp,0); // prepare temporary powers daceCreateConstant(xm[0],1.0); - + // constant part for(unsigned int i=0; i template V AlgebraicVector::eval(const V &args) const template<> template AlgebraicVector AlgebraicVector::eval(const std::initializer_list l) const{ /*! Evaluate a vector of polynomials with an braced initializer list of type U and return an AlgebraicVector of type U with the results. - \param[in] args Braced initializer list containing the arguments. + \param[in] l Braced initializer list containing the arguments. \return A new AlgebraicVector of type U containing the results of the evaluation. \note C++ is not able to derive the type of elements of an initializer list automatically. That means eval() must be called explicitly as e.g. eval({1.0, 2.0, 3.0}) when @@ -937,7 +937,7 @@ template std::istream& operator>>(std::istream &in, AlgebraicVector< // fill the AlgebraicVector for (size_t i = 0; in.good() && (i < vec_size); i++) { in >> obj[i]; - + // check the next character if (in.peek() == '\n') // the previous operator>> does not consume the \n character when an AlgebraicVector (with T != DA) is considered in.ignore(); // ignore the next character @@ -1246,7 +1246,7 @@ template AlgebraicVector eval(const AlgebraicVector &obj, con /*! Evaluate an AlgebraicVector with an braced initializer list of type T and return an AlgebraicVector of type T with the results. \param[in] obj An AlgebraicVector. - \param[in] args Braced initializer list containing the arguments. + \param[in] l Braced initializer list containing the arguments. \return A new AlgebraicVector of type T containing the results of the evaluation. \note C++ is not able to derive the type of elements of an initializer list automatically. That means eval() must be called explicitly as e.g. eval(x, {1.0, 2.0, 3.0}) when diff --git a/interfaces/cxx/include/dace/DA.h b/interfaces/cxx/include/dace/DA.h index 4cf8af1..ddb206e 100644 --- a/interfaces/cxx/include/dace/DA.h +++ b/interfaces/cxx/include/dace/DA.h @@ -56,14 +56,14 @@ template class AlgebraicVector; /*! Basic DA class representing a single polynomial. */ class DACE_API DA { - friend class compiledDA; - friend class storedDA; + friend class compiledDA; + friend class storedDA; private: static bool initialized; //!< Indicates if DA::init() was called static std::stack TOstack; //!< Truncation order stack - // XXX: Mauro, MSVC is bitching around (also for other templated classes, e.g. std::string): Warning C4251 'DACE::DA::TOstack': class 'std::stack>>' needs to have dll - interface to be used by clients of class 'DACE::DA' - DACEDA m_index; //!< Index to the DA vector + // XXX: Mauro, MSVC is bitching around (also for other templated classes, e.g. std::string): Warning C4251 'DACE::DA::TOstack': class 'std::stack>>' needs to have dll - interface to be used by clients of class 'DACE::DA' + DACEDA m_index; //!< Index to the DA vector public: /******************************************************************************** @@ -71,7 +71,7 @@ class DACE_API DA *********************************************************************************/ static void init(const unsigned int ord, const unsigned int nvar); //!< DACE initialization static bool isInitialized(); //!< Get DACE initialization status - static void version(int &maj, int &min, int &cosy_flag); //!< DACE core routines version + static void version(int &maj, int &min, int &patch); //!< DACE core routines version static void checkVersion(); //!< Check DACE C++ interface and core routines version for compatibility static double setEps(const double eps); //!< Set truncation epsilon static double getEps(); //!< Get truncation epsilon @@ -103,10 +103,10 @@ class DACE_API DA double cons() const; //!< Get constant part of a DA AlgebraicVector linear() const; //!< Get linear part of a DA AlgebraicVector gradient() const; //!< Gradient vector with respect to all independent DA variables - double getCoefficient(const std::vector &jj) const; //!< Get specific coefficient + double getCoefficient(const std::vector &jj) const; //!< Get specific coefficient void setCoefficient(const std::vector &jj, const double coeff); //!< Set specific coefficient - Monomial getMonomial(const unsigned int npos) const; //!< Get the Monomial at given position - void getMonomial(const unsigned int npos, Monomial &m) const; //!< Extract the Monomial at given position + Monomial getMonomial(const unsigned int npos) const; //!< Get the Monomial at given position + void getMonomial(const unsigned int npos, Monomial &m) const; //!< Extract the Monomial at given position std::vector getMonomials() const; //!< Get std::vector of all non-zero Monomials /******************************************************************************** @@ -134,27 +134,27 @@ class DACE_API DA *********************************************************************************/ DA operator-() const; //!< Negate this DA - friend DA DACE_API operator+(const DA &da1, const DA &da2); //!< Addition between two DAs - friend DA DACE_API operator+(const DA &da, const double c); //!< Addition between a DA and a constant - friend DA DACE_API operator+(const double c, const DA &da); //!< Addition between a constant and a DA + friend DA DACE_API operator+(const DA &da1, const DA &da2); //!< Addition between two DAs + friend DA DACE_API operator+(const DA &da, const double c); //!< Addition between a DA and a constant + friend DA DACE_API operator+(const double c, const DA &da); //!< Addition between a constant and a DA - friend DA DACE_API operator-(const DA &da1, const DA &da2); //!< Subtraction between two DAs - friend DA DACE_API operator-(const DA &da, const double c); //!< Subtraction between a DA and a constant - friend DA DACE_API operator-(const double c, const DA &da); //!< Subtraction between a constant and a DA + friend DA DACE_API operator-(const DA &da1, const DA &da2); //!< Subtraction between two DAs + friend DA DACE_API operator-(const DA &da, const double c); //!< Subtraction between a DA and a constant + friend DA DACE_API operator-(const double c, const DA &da); //!< Subtraction between a constant and a DA - friend DA DACE_API operator*(const DA &da1, const DA &da2); //!< Multiplication between two DAs - friend DA DACE_API operator*(const DA &da, const double c); //!< Multiplication between a DA and a constant - friend DA DACE_API operator*(const double c, const DA &da); //!< Multiplication between a constant and a DA + friend DA DACE_API operator*(const DA &da1, const DA &da2); //!< Multiplication between two DAs + friend DA DACE_API operator*(const DA &da, const double c); //!< Multiplication between a DA and a constant + friend DA DACE_API operator*(const double c, const DA &da); //!< Multiplication between a constant and a DA - friend DA DACE_API operator/(const DA &da1, const DA &da2); //!< Division between two DAs - friend DA DACE_API operator/(const DA &da, const double c); //!< Division between a DA and a constant - friend DA DACE_API operator/(const double c, const DA &da); //!< Division between a constant and a DA + friend DA DACE_API operator/(const DA &da1, const DA &da2); //!< Division between two DAs + friend DA DACE_API operator/(const DA &da, const double c); //!< Division between a DA and a constant + friend DA DACE_API operator/(const double c, const DA &da); //!< Division between a constant and a DA /******************************************************************************** * Math routines *********************************************************************************/ - DA multiplyMonomials(const DA &da) const; //!< Multiply the DA with the argument da monomial by monomial (i.e. coefficient-wise) - DA divide(const unsigned int var, const unsigned int p = 1) const; //!< Divide by an independent variable to some power + DA multiplyMonomials(const DA &da) const; //!< Multiply the DA with the argument da monomial by monomial (i.e. coefficient-wise) + DA divide(const unsigned int var, const unsigned int p = 1) const; //!< Divide by an independent variable to some power DA deriv(const unsigned int i) const; //!< Derivative with respect to given variable DA deriv(const std::vector ind) const; //!< Derivative with respect to given variables DA integ(const unsigned int i) const; //!< Integral with respect to given variable @@ -225,8 +225,8 @@ class DACE_API DA template T evalScalar(const T &arg) const; //!< Evaluation with a single arithmetic type T argument (not efficient for repeated evaluation!) compiledDA compile() const; //!< Compile current DA for efficient repeated evaluation DA plug(const unsigned int var, const double val = 0.0) const; //!< Partial evaluation to replace given independent DA variable by value val - double evalMonomials(const DA &da) const; //!< evaluate the DA providing the value of every monomial in da - DA replaceVariable(const unsigned int from = 0, const unsigned int to = 0, const double val = 1.0) const; + double evalMonomials(const DA &values) const; //!< evaluate the DA providing the value of every monomial in da + DA replaceVariable(const unsigned int from = 0, const unsigned int to = 0, const double val = 1.0) const; //!< Replace variable number from by val times variable number to DA scaleVariable(const unsigned int var = 0, const double val = 1.0) const; //!< Scale variable number var by val @@ -238,8 +238,8 @@ class DACE_API DA *********************************************************************************/ std::string toString() const; //!< Convert to string representation void write(std::ostream &os) const; //!< Write binary representation of DA to stream - friend DACE_API std::ostream& operator<< (std::ostream &out, const DA &da); //!< Output to C++ stream in text form - friend DACE_API std::istream& operator>> (std::istream &in, DA &da); //!< Input from C++ stream in text form + friend DACE_API std::ostream& operator<< (std::ostream &out, const DA &da); //!< Output to C++ stream in text form + friend DACE_API std::istream& operator>> (std::istream &in, DA &da); //!< Input from C++ stream in text form /******************************************************************************** * Static factory routines @@ -344,16 +344,16 @@ class DACE_API storedDA : std::vector static const unsigned int headerSize; public: - storedDA(const DA &da); //!< Constructor from a DA. - storedDA(const std::vector &data); //!< Constructor from binary data. - storedDA(std::istream &is); //!< Constructor from stream. + storedDA(const DA &da); //!< Constructor from a DA. + storedDA(const std::vector &data); //!< Constructor from binary data. + storedDA(std::istream &is); //!< Constructor from stream. bool isValid() const; //!< Is the data a valid DACE blob operator DA() const; //!< Cast to DA operator std::string() const; //!< Cast to std::string - friend DACE_API std::ostream& operator<<(std::ostream &out, const storedDA &sda); //!< Output to C++ stream in binary form + friend DACE_API std::ostream& operator<<(std::ostream &out, const storedDA &sda); //!< Output to C++ stream in binary form }; } diff --git a/interfaces/cxx/include/dace/DACEException.h b/interfaces/cxx/include/dace/DACEException.h index 6f9b5b1..8040ae6 100644 --- a/interfaces/cxx/include/dace/DACEException.h +++ b/interfaces/cxx/include/dace/DACEException.h @@ -42,7 +42,7 @@ class DACE_API DACEException : public std::exception private: int m_x; //!< Severity code int m_yy; //!< Error code - std::string msg; //!< Error message + std::string msg; //!< Error message static int severity; //!< Default severity code static bool warning; //!< Default warning status void execute() const; //!< Execute the exception @@ -57,7 +57,7 @@ class DACE_API DACEException : public std::exception static void setSeverity(const int n); //!< Select the desired severity code static void setWarning(const bool w); //!< Select the warning status - friend DACE_API std::ostream& operator<< (std::ostream &out, const DACEException &ex); //!< Overload output stream operator + friend DACE_API std::ostream& operator<< (std::ostream &out, const DACEException &ex); //!< Overload output stream operator }; }