Replace GROMACS-specific math functions with C++11 versions
authorErik Lindahl <erik@kth.se>
Sat, 10 Oct 2015 11:17:47 +0000 (13:17 +0200)
committerMark Abraham <mark.j.abraham@gmail.com>
Tue, 13 Oct 2015 13:01:32 +0000 (15:01 +0200)
Change-Id: I8e65fd6154273cbb8a18f5ba0b4eaa6c82729971

34 files changed:
CMakeLists.txt
cmake/TestFloatFormat.cpp [deleted file]
cmake/gmxManageBlueGene.cmake
cmake/gmxTestFloatFormat.cmake [deleted file]
cmake/gmxTestIsfinite.cmake [deleted file]
src/config.h.cmakein
src/gromacs/correlationfunctions/expfit.cpp
src/gromacs/ewald/long-range-correction.cpp
src/gromacs/ewald/pme-load-balancing.cpp
src/gromacs/ewald/pme-solve.cpp
src/gromacs/fileio/matio.cpp
src/gromacs/fileio/pdbio.cpp
src/gromacs/gmxana/gmx_anaeig.cpp
src/gromacs/gmxana/gmx_analyze.cpp
src/gromacs/gmxana/gmx_dipoles.cpp
src/gromacs/gmxana/gmx_dos.cpp
src/gromacs/gmxana/gmx_genion.cpp
src/gromacs/gmxana/gmx_hydorder.cpp
src/gromacs/gmxana/gmx_order.cpp
src/gromacs/gmxana/gmx_pme_error.cpp
src/gromacs/gmxana/gmx_trjcat.cpp
src/gromacs/gmxana/gmx_wham.cpp
src/gromacs/gmxpreprocess/pdb2top.cpp
src/gromacs/gmxpreprocess/readir.cpp
src/gromacs/math/calculate-ewald-splitting-coefficient.cpp
src/gromacs/math/utilities.cpp
src/gromacs/math/utilities.h
src/gromacs/mdlib/calc_verletbuf.cpp
src/gromacs/mdlib/clincs.cpp
src/gromacs/mdlib/forcerec.cpp
src/gromacs/mdlib/nbnxn_kernels/nbnxn_kernel_gpu_ref.cpp
src/gromacs/simd/tests/simd_math.cpp
src/gromacs/tables/forcetable.cpp
src/gromacs/topology/atomprop.cpp

index 31de403c44dc32e21f4dd8b71a4307c403daf0f8..fb81bd193ad606effa28dedeb2ea63c5c63b0b2c 100644 (file)
@@ -594,11 +594,6 @@ if(GMX_USE_RDTSCP)
     set(HAVE_RDTSCP 1)
 endif()
 
-include(gmxTestFloatFormat)
-gmx_test_float_format(GMX_FLOAT_FORMAT_IEEE754
-                      GMX_IEEE754_BIG_ENDIAN_BYTE_ORDER
-                      GMX_IEEE754_BIG_ENDIAN_WORD_ORDER)
-
 include(gmxTestLargeFiles)
 gmx_test_large_files(GMX_LARGEFILES)
 
@@ -608,11 +603,6 @@ gmx_test_sigusr1(HAVE_SIGUSR1)
 include(gmxTestPipes)
 gmx_test_pipes(HAVE_PIPES)
 
-include(gmxTestIsfinite)
-gmx_test_isfinite(HAVE_ISFINITE)
-gmx_test__isfinite(HAVE__ISFINITE)
-gmx_test__finite(HAVE__FINITE)
-
 include(gmxTestCXX11)
 gmx_test_cxx11(GMX_CXX11_SUPPORTED GMX_CXX11_FLAGS)
 if(NOT GMX_CXX11_SUPPORTED)
diff --git a/cmake/TestFloatFormat.cpp b/cmake/TestFloatFormat.cpp
deleted file mode 100644 (file)
index 5c1e224..0000000
+++ /dev/null
@@ -1,68 +0,0 @@
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * Copyright (c) 2014, by the GROMACS development team, led by
- * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
- * and including many others, as listed in the AUTHORS file in the
- * top-level source directory and at http://www.gromacs.org.
- *
- * GROMACS is free software; you can redistribute it and/or
- * modify it under the terms of the GNU Lesser General Public License
- * as published by the Free Software Foundation; either version 2.1
- * of the License, or (at your option) any later version.
- *
- * GROMACS is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
- * Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public
- * License along with GROMACS; if not, see
- * http://www.gnu.org/licenses, or write to the Free Software Foundation,
- * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
- *
- * If you want to redistribute modifications to GROMACS, please
- * consider that scientific software is very special. Version
- * control is crucial - bugs must be traceable. We will be happy to
- * consider code for inclusion in the official distribution, but
- * derived work must not be called official GROMACS. Details are found
- * in the README & COPYING files - if they are missing, get the
- * official version at http://www.gromacs.org.
- *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the research papers on the package. Check out http://www.gromacs.org.
- */
-#include <stdio.h>
-
-volatile const double abc [10] = {
-    /* Zero-terminated strings encoded as floating-point numbers */
-    /* "GROMACSX" in ascii    */
-    (double)  3.80279098314984902657e+35, (double) 0.0,
-    /* "GROMACSX" in ebcdic   */
-    (double) -1.37384666579378297437e+38, (double) 0.0,
-    /* "D__float" (vax)       */
-    (double)  3.53802595280598432000e+18, (double) 0.0,
-    /* "IBMHEXFP" s390/ascii  */
-    (double)  1.77977764695171661377e+10, (double) 0.0,
-    /* "IBMHEXFP" s390/ebcdic */
-    (double) -5.22995989424860458374e+10, (double) 0.0,
-};
-
-/* Check that a double is 8 bytes - compilation dies if it isnt */
-extern char xyz [sizeof(double) == 8 ? 1 : -1];
-
-int
-main()
-{
-    int         i;
-    double      d;
-
-    /* Make sure some compilers do not optimize away the entire structure
-     * with floating-point data by using it to produce a return value.
-     */
-    for (i = 0, d = 0; i < 10; i++)
-    {
-        d += abc[i];
-    }
-    return (d == 12345.0);
-}
index 300519789b685aa9ada7237ec8f347e52eb95f99..aa9a97ffa577a31ad080633f7a5e09893b901c0b 100644 (file)
@@ -1,7 +1,7 @@
 #
 # This file is part of the GROMACS molecular simulation package.
 #
-# Copyright (c) 2012,2013,2014, by the GROMACS development team, led by
+# Copyright (c) 2012,2013,2014,2015, by the GROMACS development team, led by
 # Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
 # and including many others, as listed in the AUTHORS file in the
 # top-level source directory and at http://www.gromacs.org.
@@ -67,7 +67,3 @@ set(GMX_MPI ON CACHE BOOL "MPI is required on BlueGene" FORCE)
 # warnings about harmless things in src/gromacs/utility/cstringutil.h.
 set(HAVE_PWD_H OFF)
 
-# The automatic testing for endianness does not work for the BlueGene cross-compiler
-set(GMX_FLOAT_FORMAT_IEEE754 1 CACHE INTERNAL "" FORCE)
-set(GMX_IEEE754_BIG_ENDIAN_BYTE_ORDER 1 CACHE INTERNAL "BlueGene has big-endian floating-point byte order (by default)" FORCE)
-set(GMX_IEEE754_BIG_ENDIAN_WORD_ORDER 1 CACHE INTERNAL "BlueGene has big-endian floating-point word order (by default)" FORCE)
diff --git a/cmake/gmxTestFloatFormat.cmake b/cmake/gmxTestFloatFormat.cmake
deleted file mode 100644 (file)
index 91f6d45..0000000
+++ /dev/null
@@ -1,113 +0,0 @@
-#
-# This file is part of the GROMACS molecular simulation package.
-#
-# Copyright (c) 2009,2014, by the GROMACS development team, led by
-# Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
-# and including many others, as listed in the AUTHORS file in the
-# top-level source directory and at http://www.gromacs.org.
-#
-# GROMACS is free software; you can redistribute it and/or
-# modify it under the terms of the GNU Lesser General Public License
-# as published by the Free Software Foundation; either version 2.1
-# of the License, or (at your option) any later version.
-#
-# GROMACS is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-# Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public
-# License along with GROMACS; if not, see
-# http://www.gnu.org/licenses, or write to the Free Software Foundation,
-# Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
-#
-# If you want to redistribute modifications to GROMACS, please
-# consider that scientific software is very special. Version
-# control is crucial - bugs must be traceable. We will be happy to
-# consider code for inclusion in the official distribution, but
-# derived work must not be called official GROMACS. Details are found
-# in the README & COPYING files - if they are missing, get the
-# official version at http://www.gromacs.org.
-#
-# To help us fund GROMACS development, we humbly ask that you cite
-# the research papers on the package. Check out http://www.gromacs.org.
-
-# - Define macro to determine floating point format properties
-#
-#  GMX_TEST_FLOAT_FORMAT(FP_IEEE754 FP_BIG_ENDIAN_BYTE FP_BIG_ENDIAN_WORD)
-#
-#  The thee variables are set to true when:
-#  FP_IEEE754          Floating-point numbers are stored in IEEE754 format
-#  FP_BIG_ENDIAN_BYTE  The order of bytes in each FP word (4 bytes) is big endian
-#  FP_BIG_ENDIAN_WORD  The order of FP words in double precision dwords (2 words) is big endian
-#
-#  On *most* platforms the two last tests will be the same as the integer endian,
-#  big e.g. ARM processors have different byte/word order for floating-point storage,
-#  so we figured it is a good idea to test both before relying on the format.
-
-MACRO(GMX_TEST_FLOAT_FORMAT FP_IEEE754 FP_BIG_ENDIAN_BYTE FP_BIG_ENDIAN_WORD)
-    IF(NOT DEFINED HAVE_${FP_IEEE754})
-        MESSAGE(STATUS "Checking floating point format")
-
-        TRY_COMPILE(HAVE_${FP_IEEE754} "${CMAKE_BINARY_DIR}"    
-                    "${CMAKE_SOURCE_DIR}/cmake/TestFloatFormat.cpp"
-                    COPY_FILE "${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/TestFloatFormat.bin")  
-
-        if(HAVE_${FP_IEEE754})
-
-            # dont match first/last letter because of string rounding errors :-)
-            FILE(STRINGS "${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/TestFloatFormat.bin"
-                 GMX_IEEE754_BB_BW LIMIT_COUNT 1 REGEX "ROMACS")
-            FILE(STRINGS "${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/TestFloatFormat.bin"
-                 GMX_IEEE754_BB_LW LIMIT_COUNT 1 REGEX "CSXGRO")
-            FILE(STRINGS "${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/TestFloatFormat.bin"
-                 GMX_IEEE754_LB_BW LIMIT_COUNT 1 REGEX "ORGXSC") 
-            FILE(STRINGS "${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/TestFloatFormat.bin"
-                 GMX_IEEE754_LB_LW REGEX "SCAMOR")
-           FILE(STRINGS "${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/TestFloatFormat.bin"
-                 GMX_IEEE754 REGEX "GROMACS|CSXGRO|ORGXSC|SCAMOR")
-
-            # OS X Universal binaries will contain both strings, set it to the host
-            IF(GMX_IEEE754_BB_BW  AND  GMX_IEEE754_LB_LW)
-                IF(CMAKE_SYSTEM_PROCESSOR MATCHES powerpc)
-                    SET(GMX_IEEE754_BB_BW TRUE)
-                    SET(GMX_IEEE754_LB_LW FALSE)
-                ELSE()
-                    SET(GMX_IEEE754_BB_BW FALSE)
-                    SET(GMX_IEEE754_LB_LW TRUE)
-                ENDIF()
-                MESSAGE(STATUS "GMX_TEST_IEEE754_FORMAT found different results, consider setting CMAKE_OSX_ARCHITECTURES or CMAKE_TRY_COMPILE_OSX_ARCHITECTURES to one or no architecture !")
-            ENDIF()
-
-            IF(GMX_IEEE754)
-                SET(${FP_IEEE754} 1 CACHE INTERNAL "Result of test for IEEE754 FP format" FORCE)
-            ENDIF()
-
-            IF(GMX_IEEE754_BB_BW  OR  GMX_IEEE754_BB_LW)
-                SET(${FP_BIG_ENDIAN_BYTE} 1 CACHE INTERNAL "Result of test for big endian FP byte order" FORCE)
-            ENDIF()
-
-            IF(GMX_IEEE754_BB_BW  OR  GMX_IEEE754_LB_BW)
-                SET(${FP_BIG_ENDIAN_WORD} 1 CACHE INTERNAL "Result of test for big endian FP word order" FORCE)
-            ENDIF()
-
-        endif()
-
-        # just some informational output for the user
-        if(GMX_IEEE754_BB_BW)
-            MESSAGE(STATUS "Checking floating point format - IEEE754 (BE byte, BE word)")
-        elseif(GMX_IEEE754_BB_LW)
-            MESSAGE(STATUS "Checking floating point format - IEEE754 (BE byte, LE word)")
-        elseif(GMX_IEEE754_LB_BW)
-            MESSAGE(STATUS "Checking floating point format - IEEE754 (LE byte, BE word)")
-        elseif(GMX_IEEE754_LB_LW)
-            MESSAGE(STATUS "Checking floating point format - IEEE754 (LE byte, LE word)")
-        else()
-            MESSAGE(STATUS "Checking floating point format - unknown")
-            MESSAGE(WARNING "Cannot detect your floating-point format. It is extremely unlikely to be anything else than IEEE754, but if we do not know the endian we need to rely on your OS providing the math functions erfd() and erfcd() rather than using our built-in ones.")
-        endif()
-    ENDIF()
-ENDMACRO(GMX_TEST_FLOAT_FORMAT FP_IEEE754 FP_BIG_ENDIAN_BYTE FP_BIG_ENDIAN_WORD)
-
-
-
diff --git a/cmake/gmxTestIsfinite.cmake b/cmake/gmxTestIsfinite.cmake
deleted file mode 100644 (file)
index fbccdba..0000000
+++ /dev/null
@@ -1,145 +0,0 @@
-#
-# This file is part of the GROMACS molecular simulation package.
-#
-# Copyright (c) 2012,2014,2015, by the GROMACS development team, led by
-# Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
-# and including many others, as listed in the AUTHORS file in the
-# top-level source directory and at http://www.gromacs.org.
-#
-# GROMACS is free software; you can redistribute it and/or
-# modify it under the terms of the GNU Lesser General Public License
-# as published by the Free Software Foundation; either version 2.1
-# of the License, or (at your option) any later version.
-#
-# GROMACS is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-# Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public
-# License along with GROMACS; if not, see
-# http://www.gnu.org/licenses, or write to the Free Software Foundation,
-# Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
-#
-# If you want to redistribute modifications to GROMACS, please
-# consider that scientific software is very special. Version
-# control is crucial - bugs must be traceable. We will be happy to
-# consider code for inclusion in the official distribution, but
-# derived work must not be called official GROMACS. Details are found
-# in the README & COPYING files - if they are missing, get the
-# official version at http://www.gromacs.org.
-#
-# To help us fund GROMACS development, we humbly ask that you cite
-# the research papers on the package. Check out http://www.gromacs.org.
-
-# - Define macro to check if isfinite or _isfinite exists
-#
-#  gmx_test_isfinite(VARIABLE)
-#
-#  VARIABLE will be set to true if isfinite exists
-#
-#  gmx_test__isfinite(VARIABLE)
-#
-#  VARIABLE will be set to true if _isfinite exists
-#
-#  gmx_test__finite(VARIABLE) - disabled since it doesn't seem to work the way the MSVC docs suggest
-#
-#  VARIABLE will be set to true if _finite exists
-#
-
-MACRO(gmx_test_isfinite VARIABLE)
-
-  if(NOT DEFINED isfinite_compile_ok)
-    MESSAGE(STATUS "Checking for isfinite")
-
-    set(CMAKE_REQUIRED_INCLUDES "math.h")
-    if (HAVE_LIBM)
-        set(CMAKE_REQUIRED_LIBRARIES "m")
-    endif()
-    check_c_source_compiles(
-      "#include <math.h>
-int main(void) {
-  float f;
-  isfinite(f);
-  return 0;
-}" isfinite_compile_ok)
-
-    if(isfinite_compile_ok)
-        MESSAGE(STATUS "Checking for isfinite - yes")
-    else()
-        MESSAGE(STATUS "Checking for isfinite - no")
-    endif()
-    set(isfinite_compile_ok "${isfinite_compile_ok}" CACHE INTERNAL "Result of isfinite check")
-    set(CMAKE_REQUIRED_INCLUDES)
-    set(CMAKE_REQUIRED_LIBRARIES)
-  endif()
-
-  if(isfinite_compile_ok)
-    set(${VARIABLE} ${isfinite_compile_ok}
-                "Result of test for isfinite")
-  endif()
-
-ENDMACRO(gmx_test_isfinite VARIABLE)
-
-MACRO(gmx_test__isfinite VARIABLE)
-
-  if(NOT DEFINED _isfinite_compile_ok)
-    MESSAGE(STATUS "Checking for _isfinite")
-
-    set(CMAKE_REQUIRED_INCLUDES "math.h")
-    if (HAVE_LIBM)
-        set(CMAKE_REQUIRED_LIBRARIES "m")
-    endif()
-    check_c_source_compiles(
-      "#include <math.h>
-int main(void) {
-  float f;
-  _isfinite(f);
-}" _isfinite_compile_ok)
-
-    if(_isfinite_compile_ok)
-        MESSAGE(STATUS "Checking for _isfinite - yes")
-    else()
-        MESSAGE(STATUS "Checking for _isfinite - no")
-    endif()
-    set(_isfinite_compile_ok "${_isfinite_compile_ok}" CACHE INTERNAL "Result of _isfinite check")
-    set(CMAKE_REQUIRED_INCLUDES)
-    set(CMAKE_REQUIRED_LIBRARIES)
-  endif()
-
-  if(_isfinite_compile_ok)
-    set(${VARIABLE} ${_isfinite_compile_ok}
-                "Result of test for _isfinite")
-  endif()
-
-ENDMACRO(gmx_test__isfinite VARIABLE)
-
-# Necessary for MSVC
-MACRO(gmx_test__finite VARIABLE)
-
-  if(NOT DEFINED _finite_compile_ok)
-    MESSAGE(STATUS "Checking for _finite")
-
-    set(CMAKE_REQUIRED_INCLUDES "float.h")
-    check_c_source_compiles(
-      "#include <float.h>
-int main(void) {
-  float f;
-  _finite(f);
-}" _finite_compile_ok)
-
-    if(_finite_compile_ok)
-        MESSAGE(STATUS "Checking for _finite - yes")
-    else()
-        MESSAGE(STATUS "Checking for _finite - no")
-    endif()
-    set(_finite_compile_ok "${_finite_compile_ok}" CACHE INTERNAL "Result of _finite check")
-    set(CMAKE_REQUIRED_INCLUDES)
-  endif()
-
-  if(_finite_compile_ok)
-    set(${VARIABLE} ${_finite_compile_ok}
-                "Result of test for _finite")
-  endif()
-
-ENDMACRO(gmx_test__finite VARIABLE)
index e09647aedba441edda7df5cf636c703cf9caadb3..b6ef7769993aa00febd7c26477109c9fe1d46d8b 100644 (file)
 /* TODO: For now, disable Doxygen warnings from here */
 /*! \cond */
 
-/* IEEE754 floating-point format. Memory layout is defined by macros
- * GMX_IEEE754_BIG_ENDIAN_BYTE_ORDER and GMX_IEEE754_BIG_ENDIAN_WORD_ORDER. 
- */
-#cmakedefine01 GMX_FLOAT_FORMAT_IEEE754
-
 /* Work around broken calloc() */
 #cmakedefine GMX_BROKEN_CALLOC
 
 /* Define to 1 if you have the all the affinity functions in sched.h */
 #cmakedefine HAVE_SCHED_AFFINITY
 
-/* Bytes in IEEE fp word are in big-endian order if set, little-endian if not.
-   Only relevant when FLOAT_FORMAT_IEEE754 is defined. */
-#cmakedefine01 GMX_IEEE754_BIG_ENDIAN_BYTE_ORDER
-
-/* The two words in a double precision variable are in b ig-endian order if
-   set, little-endian if not. Do NOT assume this is the same as the byte
-   order! Only relevant when FLOAT_FORMAT_IEEE754 is defined. */
-#cmakedefine01 GMX_IEEE754_BIG_ENDIAN_WORD_ORDER
-
 /* Define if SIGUSR1 is present */
 #cmakedefine01 HAVE_SIGUSR1
 
index 379bacc25b78d2aaeca0e8142a7c942a585f13d2..ef743abbb5d3a69fe180193069a8a9a4ade8c44f 100644 (file)
 
 #include "expfit.h"
 
-#include <math.h>
 #include <string.h>
 
+#include <cmath>
+
 #include <algorithm>
 
 #include "external/lmfit/lmcurve.h"
@@ -150,7 +151,7 @@ static double lmc_erffit (double x, const double *a)
     if (a[3] != 0)
     {
         erfarg = (x-a[2])/(a[3]*a[3]);
-        myerf  = gmx_erfd(erfarg);
+        myerf  = std::erf(erfarg);
     }
     else
     {
@@ -201,7 +202,7 @@ static double safe_expm1(double x)
     }
     else
     {
-        return gmx_expm1(x);
+        return std::expm1(x);
     }
 }
 
index 408f742650ef34bc9a6618e09b8a32964629351f..b7e4f6e94d818358e27eace4c90864ee3dbf1b2e 100644 (file)
@@ -38,7 +38,7 @@
 
 #include "long-range-correction.h"
 
-#include <math.h>
+#include <cmath>
 
 #include "gromacs/legacyheaders/names.h"
 #include "gromacs/legacyheaders/types/commrec.h"
@@ -226,7 +226,7 @@ void ewald_LRcorrection(int start, int end,
 
                                     dr       = 1.0/rinv;
                                     ewcdr    = ewc_q*dr;
-                                    vc       = qqA*gmx_erf(ewcdr)*rinv;
+                                    vc       = qqA*std::erf(ewcdr)*rinv;
                                     Vexcl_q += vc;
 #ifdef GMX_DOUBLE
                                     /* Relative accuracy at R_ERF_R_INACC of 3e-10 */
@@ -388,7 +388,7 @@ void ewald_LRcorrection(int start, int end,
                                     real dr;
 
                                     dr           = 1.0/rinv;
-                                    v            = gmx_erf(ewc_q*dr)*rinv;
+                                    v            = std::erf(ewc_q*dr)*rinv;
                                     vc           = qqL*v;
                                     Vexcl_q     += vc;
                                     /* fscal is the scalar force pre-multiplied by rinv,
index 42fe06195eeca255986c8053edf6aeb4ae694bc6..73902b5ae0ecfcbdd6697eea559c0d8f69cdb702 100644 (file)
@@ -800,7 +800,7 @@ pme_load_balance(pme_load_balancing_t      *pme_lb,
     /* TODO: centralize the code that sets the potentials shifts */
     if (ic->coulomb_modifier == eintmodPOTSHIFT)
     {
-        ic->sh_ewald = gmx_erfc(ic->ewaldcoeff_q*ic->rcoulomb);
+        ic->sh_ewald = std::erfc(ic->ewaldcoeff_q*ic->rcoulomb);
     }
     if (EVDW_PME(ic->vdwtype))
     {
index 60936bb666696e7c3bced726e99df1c7664c3bd8..af9975f9034b7a311564063aae508ecbb28fd001 100644 (file)
@@ -39,7 +39,7 @@
 
 #include "pme-solve.h"
 
-#include <math.h>
+#include <cmath>
 
 #include "gromacs/fft/parallel_3dfft.h"
 #include "gromacs/math/vec.h"
@@ -279,7 +279,7 @@ gmx_inline static void calc_exponentials_lj(int start, int end, real *r, real *t
     for (kx = start; kx < end; kx++)
     {
         mk       = tmp2[kx];
-        tmp2[kx] = sqrt(M_PI)*mk*gmx_erfc(mk);
+        tmp2[kx] = sqrt(M_PI)*mk*std::erfc(mk);
     }
 }
 #endif
index 8b1e7a976a1ce6e37246be813e5c2055fbae4254..964bd942c4762626285fa9db1e9a03e8be5a4acc 100644 (file)
@@ -39,6 +39,7 @@
 #include "matio.h"
 
 #include <cctype>
+#include <cmath>
 #include <cstdio>
 #include <cstring>
 
@@ -54,8 +55,6 @@
 #include "gromacs/utility/programcontext.h"
 #include "gromacs/utility/smalloc.h"
 
-#define round(a) (int)(a+0.5)
-
 static const char mapper[] = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789!@#$%^&*()-_=+{}|;:',<.>/?";
 #define NMAP (long int)strlen(mapper)
 
@@ -895,7 +894,7 @@ static void write_xpm_data(FILE *out, int n_x, int n_y, real **mat,
         fprintf(out, "\"");
         for (i = 0; (i < n_x); i++)
         {
-            c = gmx_nint((mat[i][j]-lo)*invlevel);
+            c = std::round((mat[i][j]-lo)*invlevel);
             if (c < 0)
             {
                 c = 0;
@@ -945,11 +944,11 @@ static void write_xpm_data3(FILE *out, int n_x, int n_y, real **mat,
         {
             if (mat[i][j] >= mid)
             {
-                c = nmid+gmx_nint((mat[i][j]-mid)*invlev_hi);
+                c = nmid+std::round((mat[i][j]-mid)*invlev_hi);
             }
             else if (mat[i][j] >= lo)
             {
-                c = gmx_nint((mat[i][j]-lo)*invlev_lo);
+                c = std::round((mat[i][j]-lo)*invlev_lo);
             }
             else
             {
index f4c5c23d6e8dfe6ed76e398340bda3c664a8196b..f742a9b0b8f45641786664a33c59cbfd1c9c0ea5 100644 (file)
@@ -528,14 +528,14 @@ void get_pdb_atomnumber(t_atoms *atoms, gmx_atomprop_t aps)
             anm_copy[2] = nc;
             if (gmx_atomprop_query(aps, epropElement, "???", anm_copy, &eval))
             {
-                atomnumber = gmx_nint(eval);
+                atomnumber = std::round(eval);
             }
             else
             {
                 anm_copy[1] = nc;
                 if (gmx_atomprop_query(aps, epropElement, "???", anm_copy, &eval))
                 {
-                    atomnumber = gmx_nint(eval);
+                    atomnumber = std::round(eval);
                 }
             }
         }
@@ -550,7 +550,7 @@ void get_pdb_atomnumber(t_atoms *atoms, gmx_atomprop_t aps)
             anm_copy[1] = nc;
             if (gmx_atomprop_query(aps, epropElement, "???", anm_copy, &eval))
             {
-                atomnumber = gmx_nint(eval);
+                atomnumber = std::round(eval);
             }
         }
         atoms->atom[i].atomnumber = atomnumber;
index 92d5f55e5fb0665d7497aafa302d40e2e9b060ef..0a1a1e99f4296ffce4ae45ea752805f714af0bf6 100644 (file)
@@ -81,7 +81,7 @@ static void calc_entropy_qh(FILE *fp, int n, real eigval[], real temp, int nskip
             lambda = eigval[i]*AMU;
             w      = std::sqrt(BOLTZMANN*temp/lambda)/NANO;
             hwkT   = (hbar*w)/(BOLTZMANN*temp);
-            dS     = (hwkT/gmx_expm1(hwkT) - gmx_log1p(-std::exp(-hwkT)));
+            dS     = (hwkT/std::expm1(hwkT) - std::log1p(-std::exp(-hwkT)));
             S     += dS;
             if (debug)
             {
index 6704ca6e93b3cb832f9c5e5d4b71fec11a4f31b1..d61c913a5a06507ee2734d762eb15c681bf67dcf 100644 (file)
@@ -89,7 +89,7 @@ static void power_fit(int n, int nset, real **val, real *t)
         fprintf(stdout, "First time is not larger than 0, using index number as time for power fit\n");
         for (i = 0; i < n; i++)
         {
-            x[i] = gmx_log1p(static_cast<real>(i));
+            x[i] = std::log1p(static_cast<real>(i));
         }
     }
 
index b4313049bc304bdf75997079ed1f0c7903248de3..a96cb4be770c0ea987705ae2faeac6e1eb917ce9 100644 (file)
@@ -134,7 +134,7 @@ static void done_gkrbin(t_gkrbin **gb)
 
 static void add2gkr(t_gkrbin *gb, real r, real cosa, real phi)
 {
-    int  cy, index = gmx_nint(r/gb->spacing);
+    int  cy, index = std::round(r/gb->spacing);
     real alpha;
 
     if (index < gb->nelem)
@@ -253,7 +253,7 @@ void do_gkr(t_gkrbin *gb, int ncos, int *ngrp, int *molindex[],
                     cosa = cos_angle(mu[gi], mu[gj]);
                     phi  = 0;
                 }
-                if (debug || gmx_isnan(cosa))
+                if (debug || std::isnan(cosa))
                 {
                     fprintf(debug ? debug : stderr,
                             "mu[%d] = %5.2f %5.2f %5.2f |mi| = %5.2f, mu[%d] = %5.2f %5.2f %5.2f |mj| = %5.2f rr = %5.2f cosa = %5.2f\n",
index cef2d16969de6badef8a4d3a50097bb9be41c263..c980bd57f1cd0e2c38db7fae7cae65be671cc6ff 100644 (file)
@@ -215,7 +215,7 @@ static real wSsolid(real nu, real beta)
     }
     else
     {
-        return bhn/gmx_expm1(bhn) - gmx_log1p(-std::exp(-bhn));
+        return bhn/std::expm1(bhn) - std::log1p(-std::exp(-bhn));
     }
 }
 
@@ -243,7 +243,7 @@ static real wEsolid(real nu, real beta)
     }
     else
     {
-        return bhn/2 + bhn/gmx_expm1(bhn)-1;
+        return bhn/2 + bhn/std::expm1(bhn)-1;
     }
 }
 
index 81cf326297f6aa9d7de0b8a692cd35532df2dac6..63a1197fc6e56f216a7023899c15e8c05612ad1c 100644 (file)
@@ -37,6 +37,7 @@
 #include "gmxpre.h"
 
 #include <cctype>
+#include <cmath>
 #include <cstdlib>
 #include <cstring>
 
@@ -425,14 +426,14 @@ int gmx_genion(int argc, char *argv[])
     {
         qtot += atoms.atom[i].q;
     }
-    iqtot = gmx_nint(qtot);
+    iqtot = std::round(qtot);
 
 
     if (conc > 0)
     {
         /* Compute number of ions to be added */
         vol   = det(box);
-        nsalt = gmx_nint(conc*vol*AVOGADRO/1e24);
+        nsalt = std::round(conc*vol*AVOGADRO/1e24);
         p_num = abs(nsalt*n_q);
         n_num = abs(nsalt*p_q);
     }
index 351d81b44589e2cb9a8b751867711756bed4018c..a23c1c817ac58a63e93ef70e907b4658bcca7fed 100644 (file)
@@ -232,9 +232,9 @@ static void find_tetra_order_grid(t_topology top, int ePBC,
         *skmean += skmol[i];
 
         /* Compute sliced stuff in x y z*/
-        slindex_x = gmx_nint((1+x[i][XX]/box[XX][XX])*nslicex) % nslicex;
-        slindex_y = gmx_nint((1+x[i][YY]/box[YY][YY])*nslicey) % nslicey;
-        slindex_z = gmx_nint((1+x[i][ZZ]/box[ZZ][ZZ])*nslicez) % nslicez;
+        slindex_x = static_cast<int>(std::round((1+x[i][XX]/box[XX][XX])*nslicex)) % nslicex;
+        slindex_y = static_cast<int>(std::round((1+x[i][YY]/box[YY][YY])*nslicey)) % nslicey;
+        slindex_z = static_cast<int>(std::round((1+x[i][ZZ]/box[ZZ][ZZ])*nslicez)) % nslicez;
         sggrid[slindex_x][slindex_y][slindex_z] += sgmol[i];
         skgrid[slindex_x][slindex_y][slindex_z] += skmol[i];
         (sl_count[slindex_x][slindex_y][slindex_z])++;
index 3b6f8e5016e48dea84f8a8b68cd2c93f28e4d3bf..476ba4538d8e5a6295a25cf4f1e3d7787db622d6 100644 (file)
@@ -220,7 +220,7 @@ static void find_nearest_neighbours(int ePBC,
         *skmean += skmol[i];
 
         /* Compute sliced stuff */
-        sl_index           = gmx_nint((1+x[i][slice_dim]/box[slice_dim][slice_dim])*nslice) % nslice;
+        sl_index           = static_cast<int>(std::round((1+x[i][slice_dim]/box[slice_dim][slice_dim])*nslice)) % nslice;
         sgslice[sl_index] += sgmol[i];
         skslice[sl_index] += skmol[i];
         sl_count[sl_index]++;
index ef2303d1de02e635853c9244e0ad07e505c5e4bd..3b76c432ed3526f493b2e0b474f86e1b29d5a5c8 100644 (file)
@@ -949,7 +949,7 @@ static void estimate_PME_error(t_inputinfo *info, t_state *state,
     if (MASTER(cr))
     {
         calc_q2all(mtop, &(info->q2all), &(info->q2allnr));
-        info->ewald_rtol[0] = gmx_erfc(info->rcoulomb[0]*info->ewald_beta[0]);
+        info->ewald_rtol[0] = std::erfc(info->rcoulomb[0]*info->ewald_beta[0]);
         /* Write some info to log file */
         fprintf(fp_out, "Box volume              : %g nm^3\n", info->volume);
         fprintf(fp_out, "Number of charged atoms : %d (total atoms %d)\n", ncharges, info->natoms);
@@ -1055,7 +1055,7 @@ static void estimate_PME_error(t_inputinfo *info, t_state *state,
             }
         }
 
-        info->ewald_rtol[0] = gmx_erfc(info->rcoulomb[0]*info->ewald_beta[0]);
+        info->ewald_rtol[0] = std::erfc(info->rcoulomb[0]*info->ewald_beta[0]);
 
         if (MASTER(cr))
         {
index be4dcb24813ec0a12fa16895b9a6f58534a018c9..85961695b11e2252999652f093f54e816f043b1d 100644 (file)
@@ -352,7 +352,7 @@ static void do_demux(int nset, char *fnms[], char *fnms_out[], int nval,
         fp_out[i] = open_trx(fnms_out[i], "w");
     }
     k = 0;
-    if (gmx_nint(time[k] - t) != 0)
+    if (std::round(time[k] - t) != 0)
     {
         gmx_fatal(FARGS, "First time in demuxing table does not match trajectories");
     }
@@ -372,7 +372,7 @@ static void do_demux(int nset, char *fnms[], char *fnms_out[], int nval,
         }
         for (i = 0; (i < nset); i++)
         {
-            j = gmx_nint(value[i][k]);
+            j = std::round(value[i][k]);
             range_check(j, 0, nset);
             if (bSet[j])
             {
@@ -541,7 +541,7 @@ int gmx_trjcat(int argc, char *argv[])
                 fprintf(debug, "%10g", t[i]);
                 for (j = 0; (j < nset); j++)
                 {
-                    fprintf(debug, "  %3d", gmx_nint(val[j][i]));
+                    fprintf(debug, "  %3d", static_cast<int>(std::round(val[j][i])));
                 }
                 fprintf(debug, "\n");
             }
index 29de95b70176f1bba11ac424741f7499bed3d6d5..e2ef5455bb2ca6976cb96d962417c4fdd894c882 100644 (file)
@@ -1433,7 +1433,7 @@ void create_synthetic_histo(t_UmbrellaWindow *synthWindow, t_UmbrellaWindow *thi
             /* get flat distribution in [0,1] using cumulative distribution function of Gauusian
                Note: CDF(Gaussian) = 0.5*{1+erf[x/sqrt(2)]}
              */
-            r = 0.5*(1+gmx_erf(x*invsqrt2));
+            r = 0.5*(1+std::erf(x*invsqrt2));
             searchCumulative(thisWindow->cum[pullid], nbins+1, r, &r_index);
             synthWindow->Histo[0][r_index] += 1.;
         }
index be585755f93f683bb4b9110f34c8a9cc2a97803a..05592f9ed70c2d37f77c7948e4a800e00d54b517 100644 (file)
@@ -122,7 +122,7 @@ gmx_bool is_int(double x)
     {
         x = -x;
     }
-    ix = gmx_nint(x);
+    ix = std::round(x);
 
     return (fabs(x-ix) < tol);
 }
index 4745447dac35cef6e86815f8748fdc922cc3fc13..21999e6014194680618848165733f4bdd430905b 100644 (file)
@@ -163,7 +163,7 @@ static void GetSimTemps(int ntemps, t_simtemp *simtemp, double *temperature_lamb
         }
         else if (simtemp->eSimTempScale == esimtempEXPONENTIAL)
         {
-            simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*(gmx_expm1(temperature_lambdas[i])/gmx_expm1(1.0));
+            simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*(std::expm1(temperature_lambdas[i])/std::expm1(1.0));
         }
         else
         {
index f5b4b762691d78975b8f05de170701e75a16f6bf..d653d9ffea553f1bb394a74795ab952b9e9046fe 100644 (file)
@@ -3,7 +3,7 @@
  *
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
@@ -53,7 +53,7 @@ real calc_ewaldcoeff_q(real rc, real rtol)
         i++;
         beta *= 2;
     }
-    while (gmx_erfc(beta*rc) > rtol);
+    while (std::erfc(beta*rc) > rtol);
 
     /* Do a binary search with tolerance 2^-60 */
     n    = i+60;
@@ -62,7 +62,7 @@ real calc_ewaldcoeff_q(real rc, real rtol)
     for (i = 0; i < n; i++)
     {
         beta = (low+high)/2;
-        if (gmx_erfc(beta*rc) > rtol)
+        if (std::erfc(beta*rc) > rtol)
         {
             low = beta;
         }
index ae263828e57cd564eed5704b1c2442a9494510bc..aab8a36078403c773849f9cfc0acde0a7ce4d48e 100644 (file)
 
 #include <cmath>
 #ifdef GMX_NATIVE_WINDOWS
+// for _BitScanReverse()
 #include <intrin.h>
 #endif
 
-#ifdef HAVE__FINITE
-#include <float.h>
-#endif
 #ifndef GMX_NATIVE_WINDOWS
+// for fp exception control stuff
 #include <fenv.h>
 #endif
 
-int gmx_nint(real a)
-{
-    const real half = .5;
-    int        result;
-
-    result = (a < 0.) ? ((int)(a - half)) : ((int)(a + half));
-    return result;
-}
-
-real cuberoot(real x)
-{
-    if (x < 0)
-    {
-        return (-std::pow(-x, 1.0/3.0));
-    }
-    else
-    {
-        return (std::pow(x, 1.0/3.0));
-    }
-}
-
-real sign(real x, real y)
-{
-    if (y < 0)
-    {
-        return -fabs(x);
-    }
-    else
-    {
-        return +fabs(x);
-    }
-}
-
-/* Double and single precision erf() and erfc() from
- * the Sun Freely Distributable Math Library FDLIBM.
- * See http://www.netlib.org/fdlibm
- * Specific file used: s_erf.c, version 1.3 95/01/18
- */
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunSoft, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
-
-static const double
-    tiny        = 1e-300,
-    half        =  5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
-    one         =  1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
-    two         =  2.00000000000000000000e+00, /* 0x40000000, 0x00000000 */
-/* c = (float)0.84506291151 */
-    erx =  8.45062911510467529297e-01,         /* 0x3FEB0AC1, 0x60000000 */
-/*
- * Coefficients for approximation to  erf on [0,0.84375]
- */
-    efx  =  1.28379167095512586316e-01,        /* 0x3FC06EBA, 0x8214DB69 */
-    efx8 =  1.02703333676410069053e+00,        /* 0x3FF06EBA, 0x8214DB69 */
-    pp0  =  1.28379167095512558561e-01,        /* 0x3FC06EBA, 0x8214DB68 */
-    pp1  = -3.25042107247001499370e-01,        /* 0xBFD4CD7D, 0x691CB913 */
-    pp2  = -2.84817495755985104766e-02,        /* 0xBF9D2A51, 0xDBD7194F */
-    pp3  = -5.77027029648944159157e-03,        /* 0xBF77A291, 0x236668E4 */
-    pp4  = -2.37630166566501626084e-05,        /* 0xBEF8EAD6, 0x120016AC */
-    qq1  =  3.97917223959155352819e-01,        /* 0x3FD97779, 0xCDDADC09 */
-    qq2  =  6.50222499887672944485e-02,        /* 0x3FB0A54C, 0x5536CEBA */
-    qq3  =  5.08130628187576562776e-03,        /* 0x3F74D022, 0xC4D36B0F */
-    qq4  =  1.32494738004321644526e-04,        /* 0x3F215DC9, 0x221C1A10 */
-    qq5  = -3.96022827877536812320e-06,        /* 0xBED09C43, 0x42A26120 */
-/*
- * Coefficients for approximation to  erf  in [0.84375,1.25]
- */
-    pa0  = -2.36211856075265944077e-03,        /* 0xBF6359B8, 0xBEF77538 */
-    pa1  =  4.14856118683748331666e-01,        /* 0x3FDA8D00, 0xAD92B34D */
-    pa2  = -3.72207876035701323847e-01,        /* 0xBFD7D240, 0xFBB8C3F1 */
-    pa3  =  3.18346619901161753674e-01,        /* 0x3FD45FCA, 0x805120E4 */
-    pa4  = -1.10894694282396677476e-01,        /* 0xBFBC6398, 0x3D3E28EC */
-    pa5  =  3.54783043256182359371e-02,        /* 0x3FA22A36, 0x599795EB */
-    pa6  = -2.16637559486879084300e-03,        /* 0xBF61BF38, 0x0A96073F */
-    qa1  =  1.06420880400844228286e-01,        /* 0x3FBB3E66, 0x18EEE323 */
-    qa2  =  5.40397917702171048937e-01,        /* 0x3FE14AF0, 0x92EB6F33 */
-    qa3  =  7.18286544141962662868e-02,        /* 0x3FB2635C, 0xD99FE9A7 */
-    qa4  =  1.26171219808761642112e-01,        /* 0x3FC02660, 0xE763351F */
-    qa5  =  1.36370839120290507362e-02,        /* 0x3F8BEDC2, 0x6B51DD1C */
-    qa6  =  1.19844998467991074170e-02,        /* 0x3F888B54, 0x5735151D */
-/*
- * Coefficients for approximation to  erfc in [1.25,1/0.35]
- */
-    ra0  = -9.86494403484714822705e-03,        /* 0xBF843412, 0x600D6435 */
-    ra1  = -6.93858572707181764372e-01,        /* 0xBFE63416, 0xE4BA7360 */
-    ra2  = -1.05586262253232909814e+01,        /* 0xC0251E04, 0x41B0E726 */
-    ra3  = -6.23753324503260060396e+01,        /* 0xC04F300A, 0xE4CBA38D */
-    ra4  = -1.62396669462573470355e+02,        /* 0xC0644CB1, 0x84282266 */
-    ra5  = -1.84605092906711035994e+02,        /* 0xC067135C, 0xEBCCABB2 */
-    ra6  = -8.12874355063065934246e+01,        /* 0xC0545265, 0x57E4D2F2 */
-    ra7  = -9.81432934416914548592e+00,        /* 0xC023A0EF, 0xC69AC25C */
-    sa1  =  1.96512716674392571292e+01,        /* 0x4033A6B9, 0xBD707687 */
-    sa2  =  1.37657754143519042600e+02,        /* 0x4061350C, 0x526AE721 */
-    sa3  =  4.34565877475229228821e+02,        /* 0x407B290D, 0xD58A1A71 */
-    sa4  =  6.45387271733267880336e+02,        /* 0x40842B19, 0x21EC2868 */
-    sa5  =  4.29008140027567833386e+02,        /* 0x407AD021, 0x57700314 */
-    sa6  =  1.08635005541779435134e+02,        /* 0x405B28A3, 0xEE48AE2C */
-    sa7  =  6.57024977031928170135e+00,        /* 0x401A47EF, 0x8E484A93 */
-    sa8  = -6.04244152148580987438e-02,        /* 0xBFAEEFF2, 0xEE749A62 */
-/*
- * Coefficients for approximation to  erfc in [1/.35,28]
- */
-    rb0  = -9.86494292470009928597e-03,        /* 0xBF843412, 0x39E86F4A */
-    rb1  = -7.99283237680523006574e-01,        /* 0xBFE993BA, 0x70C285DE */
-    rb2  = -1.77579549177547519889e+01,        /* 0xC031C209, 0x555F995A */
-    rb3  = -1.60636384855821916062e+02,        /* 0xC064145D, 0x43C5ED98 */
-    rb4  = -6.37566443368389627722e+02,        /* 0xC083EC88, 0x1375F228 */
-    rb5  = -1.02509513161107724954e+03,        /* 0xC0900461, 0x6A2E5992 */
-    rb6  = -4.83519191608651397019e+02,        /* 0xC07E384E, 0x9BDC383F */
-    sb1  =  3.03380607434824582924e+01,        /* 0x403E568B, 0x261D5190 */
-    sb2  =  3.25792512996573918826e+02,        /* 0x40745CAE, 0x221B9F0A */
-    sb3  =  1.53672958608443695994e+03,        /* 0x409802EB, 0x189D5118 */
-    sb4  =  3.19985821950859553908e+03,        /* 0x40A8FFB7, 0x688C246A */
-    sb5  =  2.55305040643316442583e+03,        /* 0x40A3F219, 0xCEDF3BE6 */
-    sb6  =  4.74528541206955367215e+02,        /* 0x407DA874, 0xE79FE763 */
-    sb7  = -2.24409524465858183362e+01;        /* 0xC03670E2, 0x42712D62 */
-
-double gmx_erfd(double x)
-{
-#if GMX_FLOAT_FORMAT_IEEE754
-    gmx_int32_t hx, ix, i;
-    double      R, S, P, Q, s, y, z, r;
-
-    union
-    {
-        double d;
-        int    i[2];
-    }
-    conv;
-
-    conv.d = x;
-
-#if GMX_IEEE754_BIG_ENDIAN_WORD_ORDER
-    hx = conv.i[0];
-#else
-    hx = conv.i[1];
-#endif
-
-    ix = hx&0x7fffffff;
-    if (ix >= 0x7ff00000)
-    {
-        /* erf(nan)=nan */
-        i = ((gmx_uint32_t)hx>>31)<<1;
-        return (double)(1-i)+one/x; /* erf(+-inf)=+-1 */
-    }
-
-    if (ix < 0x3feb0000)
-    {
-        /* |x|<0.84375 */
-        if (ix < 0x3e300000)
-        {
-            /* |x|<2**-28 */
-            if (ix < 0x00800000)
-            {
-                return 0.125*(8.0*x+efx8*x);  /*avoid underflow */
-            }
-            return x + efx*x;
-        }
-        z = x*x;
-        r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
-        s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
-        y = r/s;
-        return x + x*y;
-    }
-    if (ix < 0x3ff40000)
-    {
-        /* 0.84375 <= |x| < 1.25 */
-        s = fabs(x)-one;
-        P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
-        Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
-        if (hx >= 0)
-        {
-            return erx + P/Q;
-        }
-        else
-        {
-            return -erx - P/Q;
-        }
-    }
-    if (ix >= 0x40180000)
-    {
-        /* inf>|x|>=6 */
-        if (hx >= 0)
-        {
-            return one-tiny;
-        }
-        else
-        {
-            return tiny-one;
-        }
-    }
-    x = fabs(x);
-    s = one/(x*x);
-    if (ix < 0x4006DB6E)
-    {
-        /* |x| < 1/0.35 */
-        R = ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(ra5+s*(ra6+s*ra7))))));
-        S = one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(sa5+s*(sa6+s*(sa7+s*sa8)))))));
-    }
-    else
-    {
-        /* |x| >= 1/0.35 */
-        R = rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(rb5+s*rb6)))));
-        S = one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(sb5+s*(sb6+s*sb7))))));
-    }
-
-    conv.d = x;
-
-#if GMX_IEEE754_BIG_ENDIAN_WORD_ORDER
-    conv.i[1] = 0;
-#else
-    conv.i[0] = 0;
-#endif
-
-    z = conv.d;
-
-    r  =  exp(-z*z-0.5625)*exp((z-x)*(z+x)+R/S);
-    if (hx >= 0)
-    {
-        return one-r/x;
-    }
-    else
-    {
-        return r/x-one;
-    }
-#else
-    /* No IEEE754 information. We need to trust that the OS provides erf(). */
-    return erf(x);
-#endif
-}
-
-
-double gmx_erfcd(double x)
-{
-#if GMX_FLOAT_FORMAT_IEEE754
-    gmx_int32_t hx, ix;
-    double      R, S, P, Q, s, y, z, r;
-
-    union
-    {
-        double d;
-        int    i[2];
-    }
-    conv;
-
-    conv.d = x;
-
-#if GMX_IEEE754_BIG_ENDIAN_WORD_ORDER
-    hx = conv.i[0];
-#else
-    hx = conv.i[1];
-#endif
-
-    ix = hx&0x7fffffff;
-    if (ix >= 0x7ff00000)
-    {
-        /* erfc(nan)=nan */
-        /* erfc(+-inf)=0,2 */
-        return (double)(((gmx_uint32_t)hx>>31)<<1)+one/x;
-    }
-
-    if (ix < 0x3feb0000)
-    {
-        /* |x|<0.84375 */
-        if (ix < 0x3c700000)     /* |x|<2**-56 */
-        {
-            return one-x;
-        }
-        z = x*x;
-        r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
-        s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
-        y = r/s;
-        if (hx < 0x3fd00000)
-        {
-            /* x<1/4 */
-            return one-(x+x*y);
-        }
-        else
-        {
-            r  = x*y;
-            r += (x-half);
-            return half - r;
-        }
-    }
-
-    if (ix < 0x3ff40000)
-    {
-        /* 0.84375 <= |x| < 1.25 */
-        s = fabs(x)-one;
-        P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
-        Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
-        if (hx >= 0)
-        {
-            z  = one-erx; return z - P/Q;
-        }
-        else
-        {
-            z = erx+P/Q; return one+z;
-        }
-    }
-    if (ix < 0x403c0000)
-    {
-        /* |x|<28 */
-        x = fabs(x);
-        s = one/(x*x);
-        if (ix < 0x4006DB6D)
-        {
-            /* |x| < 1/.35 ~ 2.857143*/
-            R = ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(ra5+s*(ra6+s*ra7))))));
-            S = one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(sa5+s*(sa6+s*(sa7+s*sa8)))))));
-        }
-        else
-        {
-            /* |x| >= 1/.35 ~ 2.857143 */
-            if (hx < 0 && ix >= 0x40180000)
-            {
-                return two-tiny; /* x < -6 */
-            }
-            R = rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(rb5+s*rb6)))));
-            S = one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(sb5+s*(sb6+s*sb7))))));
-        }
-
-        conv.d = x;
-
-#if GMX_IEEE754_BIG_ENDIAN_WORD_ORDER
-        conv.i[1] = 0;
-#else
-        conv.i[0] = 0;
-#endif
-
-        z = conv.d;
-
-        r  =  exp(-z*z-0.5625)*exp((z-x)*(z+x)+R/S);
-
-        if (hx > 0)
-        {
-            return r/x;
-        }
-        else
-        {
-            return two-r/x;
-        }
-    }
-    else
-    {
-        if (hx > 0)
-        {
-            return tiny*tiny;
-        }
-        else
-        {
-            return two-tiny;
-        }
-    }
-#else
-    /* No IEEE754 information. We need to trust that the OS provides erfc(). */
-    return erfc(x);
-#endif
-}
-
-
-static const float
-    tinyf =  1e-30,
-    halff =  5.0000000000e-01, /* 0x3F000000 */
-    onef  =  1.0000000000e+00, /* 0x3F800000 */
-    twof  =  2.0000000000e+00, /* 0x40000000 */
-/* c = (subfloat)0.84506291151 */
-    erxf =  8.4506291151e-01,  /* 0x3f58560b */
-/*
- * Coefficients for approximation to  erf on [0,0.84375]
- */
-    efxf  =  1.2837916613e-01, /* 0x3e0375d4 */
-    efx8f =  1.0270333290e+00, /* 0x3f8375d4 */
-    pp0f  =  1.2837916613e-01, /* 0x3e0375d4 */
-    pp1f  = -3.2504209876e-01, /* 0xbea66beb */
-    pp2f  = -2.8481749818e-02, /* 0xbce9528f */
-    pp3f  = -5.7702702470e-03, /* 0xbbbd1489 */
-    pp4f  = -2.3763017452e-05, /* 0xb7c756b1 */
-    qq1f  =  3.9791721106e-01, /* 0x3ecbbbce */
-    qq2f  =  6.5022252500e-02, /* 0x3d852a63 */
-    qq3f  =  5.0813062117e-03, /* 0x3ba68116 */
-    qq4f  =  1.3249473704e-04, /* 0x390aee49 */
-    qq5f  = -3.9602282413e-06, /* 0xb684e21a */
-/*
- * Coefficients for approximation to  erf  in [0.84375,1.25]
- */
-    pa0f = -2.3621185683e-03,  /* 0xbb1acdc6 */
-    pa1f =  4.1485610604e-01,  /* 0x3ed46805 */
-    pa2f = -3.7220788002e-01,  /* 0xbebe9208 */
-    pa3f =  3.1834661961e-01,  /* 0x3ea2fe54 */
-    pa4f = -1.1089469492e-01,  /* 0xbde31cc2 */
-    pa5f =  3.5478305072e-02,  /* 0x3d1151b3 */
-    pa6f = -2.1663755178e-03,  /* 0xbb0df9c0 */
-    qa1f =  1.0642088205e-01,  /* 0x3dd9f331 */
-    qa2f =  5.4039794207e-01,  /* 0x3f0a5785 */
-    qa3f =  7.1828655899e-02,  /* 0x3d931ae7 */
-    qa4f =  1.2617121637e-01,  /* 0x3e013307 */
-    qa5f =  1.3637083583e-02,  /* 0x3c5f6e13 */
-    qa6f =  1.1984500103e-02,  /* 0x3c445aa3 */
-/*
- * Coefficients for approximation to  erfc in [1.25,1/0.35]
- */
-    ra0f = -9.8649440333e-03,  /* 0xbc21a093 */
-    ra1f = -6.9385856390e-01,  /* 0xbf31a0b7 */
-    ra2f = -1.0558626175e+01,  /* 0xc128f022 */
-    ra3f = -6.2375331879e+01,  /* 0xc2798057 */
-    ra4f = -1.6239666748e+02,  /* 0xc322658c */
-    ra5f = -1.8460508728e+02,  /* 0xc3389ae7 */
-    ra6f = -8.1287437439e+01,  /* 0xc2a2932b */
-    ra7f = -9.8143291473e+00,  /* 0xc11d077e */
-    sa1f =  1.9651271820e+01,  /* 0x419d35ce */
-    sa2f =  1.3765776062e+02,  /* 0x4309a863 */
-    sa3f =  4.3456588745e+02,  /* 0x43d9486f */
-    sa4f =  6.4538726807e+02,  /* 0x442158c9 */
-    sa5f =  4.2900814819e+02,  /* 0x43d6810b */
-    sa6f =  1.0863500214e+02,  /* 0x42d9451f */
-    sa7f =  6.5702495575e+00,  /* 0x40d23f7c */
-    sa8f = -6.0424413532e-02,  /* 0xbd777f97 */
-/*
- * Coefficients for approximation to  erfc in [1/.35,28]
- */
-    rb0f = -9.8649431020e-03,  /* 0xbc21a092 */
-    rb1f = -7.9928326607e-01,  /* 0xbf4c9dd4 */
-    rb2f = -1.7757955551e+01,  /* 0xc18e104b */
-    rb3f = -1.6063638306e+02,  /* 0xc320a2ea */
-    rb4f = -6.3756646729e+02,  /* 0xc41f6441 */
-    rb5f = -1.0250950928e+03,  /* 0xc480230b */
-    rb6f = -4.8351919556e+02,  /* 0xc3f1c275 */
-    sb1f =  3.0338060379e+01,  /* 0x41f2b459 */
-    sb2f =  3.2579251099e+02,  /* 0x43a2e571 */
-    sb3f =  1.5367296143e+03,  /* 0x44c01759 */
-    sb4f =  3.1998581543e+03,  /* 0x4547fdbb */
-    sb5f =  2.5530502930e+03,  /* 0x451f90ce */
-    sb6f =  4.7452853394e+02,  /* 0x43ed43a7 */
-    sb7f = -2.2440952301e+01;  /* 0xc1b38712 */
-
-
-float gmx_erff(float x)
-{
-    gmx_int32_t hx, ix, i;
-    float       R, S, P, Q, s, y, z, r;
-
-    union
-    {
-        float  f;
-        int    i;
-    }
-    conv;
-
-    conv.f = x;
-    hx     = conv.i;
-
-    ix = hx&0x7fffffff;
-    if (ix >= 0x7f800000)
-    {
-        /* erf(nan)=nan */
-        i = ((gmx_uint32_t)hx>>31)<<1;
-        return (float)(1-i)+onef/x; /* erf(+-inf)=+-1 */
-    }
-
-    if (ix < 0x3f580000)
-    {
-        /* |x|<0.84375 */
-        if (ix < 0x31800000)
-        {
-            /* |x|<2**-28 */
-            if (ix < 0x04000000)
-            {
-                return (float)0.125*((float)8.0*x+efx8f*x);             /*avoid underflow */
-            }
-            return x + efxf*x;
-        }
-        z = x*x;
-        r = pp0f+z*(pp1f+z*(pp2f+z*(pp3f+z*pp4f)));
-        s = onef+z*(qq1f+z*(qq2f+z*(qq3f+z*(qq4f+z*qq5f))));
-        y = r/s;
-        return x + x*y;
-    }
-    if (ix < 0x3fa00000)
-    {
-        /* 0.84375 <= |x| < 1.25 */
-        s = fabs(x)-onef;
-        P = pa0f+s*(pa1f+s*(pa2f+s*(pa3f+s*(pa4f+s*(pa5f+s*pa6f)))));
-        Q = onef+s*(qa1f+s*(qa2f+s*(qa3f+s*(qa4f+s*(qa5f+s*qa6f)))));
-        if (hx >= 0)
-        {
-            return erxf + P/Q;
-        }
-        else
-        {
-            return -erxf - P/Q;
-        }
-    }
-    if (ix >= 0x40c00000)
-    {
-        /* inf>|x|>=6 */
-        if (hx >= 0)
-        {
-            return onef-tinyf;
-        }
-        else
-        {
-            return tinyf-onef;
-        }
-    }
-    x = fabs(x);
-    s = onef/(x*x);
-    if (ix < 0x4036DB6E)
-    {
-        /* |x| < 1/0.35 */
-        R = ra0f+s*(ra1f+s*(ra2f+s*(ra3f+s*(ra4f+s*(ra5f+s*(ra6f+s*ra7f))))));
-        S = onef+s*(sa1f+s*(sa2f+s*(sa3f+s*(sa4f+s*(sa5f+s*(sa6f+s*(sa7f+s*sa8f)))))));
-    }
-    else
-    {
-        /* |x| >= 1/0.35 */
-        R = rb0f+s*(rb1f+s*(rb2f+s*(rb3f+s*(rb4f+s*(rb5f+s*rb6f)))));
-        S = onef+s*(sb1f+s*(sb2f+s*(sb3f+s*(sb4f+s*(sb5f+s*(sb6f+s*sb7f))))));
-    }
-
-    conv.f = x;
-    conv.i = conv.i & 0xfffff000;
-    z      = conv.f;
-
-    r  =  exp(-z*z-(float)0.5625)*exp((z-x)*(z+x)+R/S);
-    if (hx >= 0)
-    {
-        return onef-r/x;
-    }
-    else
-    {
-        return r/x-onef;
-    }
-}
-
-float gmx_erfcf(float x)
-{
-    gmx_int32_t hx, ix;
-    float       R, S, P, Q, s, y, z, r;
-
-    union
-    {
-        float  f;
-        int    i;
-    }
-    conv;
-
-    conv.f = x;
-    hx     = conv.i;
-
-    ix = hx&0x7fffffff;
-    if (ix >= 0x7f800000)
-    {
-        /* erfc(nan)=nan */
-        /* erfc(+-inf)=0,2 */
-        return (float)(((gmx_uint32_t)hx>>31)<<1)+onef/x;
-    }
-
-    if (ix < 0x3f580000)
-    {
-        /* |x|<0.84375 */
-        if (ix < 0x23800000)
-        {
-            return onef-x;  /* |x|<2**-56 */
-        }
-        z = x*x;
-        r = pp0f+z*(pp1f+z*(pp2f+z*(pp3f+z*pp4f)));
-        s = onef+z*(qq1f+z*(qq2f+z*(qq3f+z*(qq4f+z*qq5f))));
-        y = r/s;
-        if (hx < 0x3e800000)
-        {
-            /* x<1/4 */
-            return onef-(x+x*y);
-        }
-        else
-        {
-            r  = x*y;
-            r += (x-halff);
-            return halff - r;
-        }
-    }
-    if (ix < 0x3fa00000)
-    {
-        /* 0.84375 <= |x| < 1.25 */
-        s = fabs(x)-onef;
-        P = pa0f+s*(pa1f+s*(pa2f+s*(pa3f+s*(pa4f+s*(pa5f+s*pa6f)))));
-        Q = onef+s*(qa1f+s*(qa2f+s*(qa3f+s*(qa4f+s*(qa5f+s*qa6f)))));
-        if (hx >= 0)
-        {
-            z  = onef-erxf; return z - P/Q;
-        }
-        else
-        {
-            z = erxf+P/Q; return onef+z;
-        }
-    }
-    if (ix < 0x41e00000)
-    {
-        /* |x|<28 */
-        x = fabs(x);
-        s = onef/(x*x);
-        if (ix < 0x4036DB6D)
-        {
-            /* |x| < 1/.35 ~ 2.857143*/
-            R = ra0f+s*(ra1f+s*(ra2f+s*(ra3f+s*(ra4f+s*(ra5f+s*(ra6f+s*ra7f))))));
-            S = onef+s*(sa1f+s*(sa2f+s*(sa3f+s*(sa4f+s*(sa5f+s*(sa6f+s*(sa7f+s*sa8f)))))));
-        }
-        else
-        {
-            /* |x| >= 1/.35 ~ 2.857143 */
-            if (hx < 0 && ix >= 0x40c00000)
-            {
-                return twof-tinyf;                     /* x < -6 */
-            }
-            R = rb0f+s*(rb1f+s*(rb2f+s*(rb3f+s*(rb4f+s*(rb5f+s*rb6f)))));
-            S = onef+s*(sb1f+s*(sb2f+s*(sb3f+s*(sb4f+s*(sb5f+s*(sb6f+s*sb7f))))));
-        }
-
-        conv.f = x;
-        conv.i = conv.i & 0xfffff000;
-        z      = conv.f;
-
-        r  =  exp(-z*z-(float)0.5625)*exp((z-x)*(z+x)+R/S);
-        if (hx > 0)
-        {
-            return r/x;
-        }
-        else
-        {
-            return twof-r/x;
-        }
-    }
-    else
-    {
-        if (hx > 0)
-        {
-            return tinyf*tinyf;
-        }
-        else
-        {
-            return twof-tinyf;
-        }
-    }
-}
-
-
-gmx_bool gmx_isfinite(real gmx_unused x)
-{
-    gmx_bool returnval;
-
-#ifdef HAVE__FINITE
-    returnval = _finite(x);
-#elif defined HAVE_ISFINITE
-    returnval = std::isfinite(x);
-#elif defined HAVE__ISFINITE
-    returnval = _isfinite(x);
-#else
-    /* If no suitable function was found, assume the value is
-     * finite. */
-    returnval = TRUE;
-#endif
-    return returnval;
-}
-
-gmx_bool gmx_isnan(real x)
-{
-    return x != x;
-}
-
 int
 gmx_within_tol(double   f1,
                double   f2,
index 467430d9eaaea76ff34057cb4c2823c08a0c0bc3..2f65fef687d9c7b35ea1534e7ff77cac161945dd 100644 (file)
@@ -82,33 +82,6 @@ extern "C" {
 #define M_2_SQRTPI  1.128379167095513
 #endif
 
-int     gmx_nint(real a);
-real    sign(real x, real y);
-
-real    cuberoot (real a);
-double  gmx_erfd(double x);
-double  gmx_erfcd(double x);
-float   gmx_erff(float x);
-float   gmx_erfcf(float x);
-#ifdef GMX_DOUBLE
-#define gmx_erf(x)   gmx_erfd(x)
-#define gmx_erfc(x)  gmx_erfcd(x)
-#else
-#define gmx_erf(x)   gmx_erff(x)
-#define gmx_erfc(x)  gmx_erfcf(x)
-#endif
-
-#if defined(_MSC_VER) && _MSC_VER < 1800
-#define gmx_expm1(x) (exp(x)-1)
-#define gmx_log1p(x) log(1+x)
-#else
-#define gmx_expm1 expm1
-#define gmx_log1p log1p
-#endif
-
-gmx_bool gmx_isfinite(real x);
-gmx_bool gmx_isnan(real x);
-
 /*! \brief Check if two numbers are within a tolerance
  *
  *  This routine checks if the relative difference between two numbers is
index fcb0713a687cc6c5dc06a876a9ae9daef75e9677..6a1362f2ea838da179f1b749e939e3835b459673 100644 (file)
 #include "calc_verletbuf.h"
 
 #include <assert.h>
-#include <math.h>
 #include <stdlib.h>
 
+#include <cmath>
+
 #include <algorithm>
 
 #include "gromacs/legacyheaders/typedefs.h"
@@ -564,7 +565,7 @@ static void approx_2dof(real s2, real x, real *shift, real *scale)
     real ex, er;
 
     ex = exp(-x*x/(2*s2));
-    er = gmx_erfc(x/sqrt(2*s2));
+    er = std::erfc(x/sqrt(2*s2));
 
     *shift = -x + sqrt(2*s2/M_PI)*ex/er;
     *scale = 0.5*M_PI*exp(ex*ex/(M_PI*er*er))*er;
@@ -678,7 +679,7 @@ static real ener_drift(const verletbuf_atomtype_t *att, int natt,
                  * atom density still needs to be put in.
                  */
                 c_exp  = exp(-rsh*rsh/(2*s2))/sqrt(2*M_PI);
-                c_erfc = 0.5*gmx_erfc(rsh/(sqrt(2*s2)));
+                c_erfc = 0.5*std::erfc(rsh/(sqrt(2*s2)));
             }
             s      = sqrt(s2);
             rsh2   = rsh*rsh;
@@ -988,8 +989,8 @@ void calc_verlet_buffer_size(const gmx_mtop_t *mtop, real boxvol,
         b      = calc_ewaldcoeff_q(ir->rcoulomb, ir->ewald_rtol);
         rc     = ir->rcoulomb;
         br     = b*rc;
-        md1_el = elfac*(b*exp(-br*br)*M_2_SQRTPI/rc + gmx_erfc(br)/(rc*rc));
-        d2_el  = elfac/(rc*rc)*(2*b*(1 + br*br)*exp(-br*br)*M_2_SQRTPI + 2*gmx_erfc(br)/rc);
+        md1_el = elfac*(b*exp(-br*br)*M_2_SQRTPI/rc + std::erfc(br)/(rc*rc));
+        d2_el  = elfac/(rc*rc)*(2*b*(1 + br*br)*exp(-br*br)*M_2_SQRTPI + 2*std::erfc(br)/rc);
     }
     else
     {
index 5a57aa8d44bafd519b862fe6a99c584c22bfdcb5..5a6c4649fb0f027a8d7f922bbec85781778e92e9 100644 (file)
 #include "config.h"
 
 #include <assert.h>
-#include <math.h>
 #include <stdlib.h>
 
+#include <cmath>
+
 #include <algorithm>
 
 #include "gromacs/domdec/domdec.h"
@@ -2436,7 +2437,7 @@ static void lincs_warning(FILE *fplog,
             {
                 fprintf(fplog, "%s", buf);
             }
-            if (!gmx_isfinite(d1))
+            if (!std::isfinite(d1))
             {
                 gmx_fatal(FARGS, "Bond length not finite.");
             }
index 114dc25bea697fd8efa5f1f60d8db41ab9105345..7b8c2baf52863ce9d38eee3a6797c4742cf69015 100644 (file)
 #include "config.h"
 
 #include <assert.h>
-#include <math.h>
 #include <stdlib.h>
 #include <string.h>
 
+#include <cmath>
+
 #include <algorithm>
 
 #include "gromacs/domdec/domdec.h"
@@ -2052,7 +2053,7 @@ init_interaction_const(FILE                       *fp,
 
     if (fr->coulomb_modifier == eintmodPOTSHIFT)
     {
-        ic->sh_ewald = gmx_erfc(ic->ewaldcoeff_q*ic->rcoulomb);
+        ic->sh_ewald = std::erfc(ic->ewaldcoeff_q*ic->rcoulomb);
     }
     else
     {
index bab0003464886153fd708fbee47345dd7124976f..362977f06fab9963235974cd514081848a0eecf3 100644 (file)
@@ -284,7 +284,7 @@ nbnxn_kernel_gpu_ref(const nbnxn_pairlist_t     *nbl,
 
                                     if (bEner)
                                     {
-                                        vcoul = qq*((int_bit - gmx_erf(iconst->ewaldcoeff_q*r))*rinv - int_bit*iconst->sh_ewald);
+                                        vcoul = qq*((int_bit - std::erf(iconst->ewaldcoeff_q*r))*rinv - int_bit*iconst->sh_ewald);
                                     }
                                 }
 
index 5b3a44b66ebe8fe6a9d5cd0be763f27d38173d5a..e8f40936e01c7e139a5112f8fed8b4df56a11ddd 100644 (file)
 
 #include "config.h"
 
+#include <cmath>
+
 #include <vector>
 
-#include "gromacs/math/utilities.h"
 #include "gromacs/options/basicoptions.h"
 #include "gromacs/simd/simd.h"
 
@@ -272,18 +273,10 @@ TEST_F(SimdMathTest, gmxSimdLogR)
     GMX_EXPECT_SIMD_FUNC_NEAR(ref_log, gmx_simd_log_r);
 }
 
-// MSVC does not support exp2(), so we have no reference to test against
-#ifndef _MSC_VER
-/*! \brief Function wrapper for exp2(x), with argument/return in default Gromacs precision */
-real ref_exp2(real x)
-{
-    return exp2(x);
-}
-
 TEST_F(SimdMathTest, gmxSimdExp2R)
 {
     setRange(-100, 100);
-    GMX_EXPECT_SIMD_FUNC_NEAR(ref_exp2, gmx_simd_exp2_r);
+    GMX_EXPECT_SIMD_FUNC_NEAR(std::exp2, gmx_simd_exp2_r);
 
     // We do not care about the SIMD implementation getting denormal values right,
     // but they must be clamped to zero rather than producing garbage.
@@ -291,7 +284,7 @@ TEST_F(SimdMathTest, gmxSimdExp2R)
     setAbsTol(GMX_REAL_EPS);
 
     // First two values will have denormal results in single, third value in double too.
-    GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(ref_exp2(-150.0), ref_exp2(-300.0), ref_exp2(-1050.0)),
+    GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(std::exp2(-150.0), std::exp2(-300.0), std::exp2(-1050.0)),
                               gmx_simd_exp2_r(setSimdRealFrom3R(-150.0, -300.0, -1050.0)));
 
     // Reset absolute tolerance to enforce ULP checking
@@ -299,10 +292,9 @@ TEST_F(SimdMathTest, gmxSimdExp2R)
 
     // Make sure that underflowing values are set to zero.
     // First two values underflow in single, third value in double too.
-    GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(ref_exp2(-200.0), ref_exp2(-600.0), ref_exp2(-1500.0)),
+    GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(std::exp2(-200.0), std::exp2(-600.0), std::exp2(-1500.0)),
                               gmx_simd_exp2_r(setSimdRealFrom3R(-200.0, -600.0, -1500.0)));
 }
-#endif
 
 /*! \brief Function wrapper for exp(x), with argument/return in default Gromacs precision */
 real ref_exp(real x)
@@ -334,12 +326,12 @@ TEST_F(SimdMathTest, gmxSimdExpR)
 
 /*! \brief Function wrapper for erf(x), with argument/return in default Gromacs precision.
  *
- * \note The single-precision gmx_erff() in gmxlib is slightly lower precision
- * than the SIMD flavor, so we use double for reference.
+ * \note Single-precision erf() in some libraries can be slightly lower precision
+ * than the SIMD flavor, so we use a cast to force double precision for reference.
  */
 real ref_erf(real x)
 {
-    return gmx_erfd(x);
+    return std::erf(static_cast<double>(x));
 }
 
 TEST_F(SimdMathTest, gmxSimdErfR)
@@ -351,12 +343,12 @@ TEST_F(SimdMathTest, gmxSimdErfR)
 
 /*! \brief Function wrapper for erfc(x), with argument/return in default Gromacs precision.
  *
- * \note The single-precision gmx_erfcf() in gmxlib is slightly lower precision
- * than the SIMD flavor, so we use double for reference.
+ * \note Single-precision erfc() in some libraries can be slightly lower precision
+ * than the SIMD flavor, so we use a cast to force double precision for reference.
  */
 real ref_erfc(real x)
 {
-    return gmx_erfcd(x);
+    return std::erfc(static_cast<double>(x));
 }
 
 TEST_F(SimdMathTest, gmxSimdErfcR)
@@ -488,7 +480,7 @@ real ref_pmecorrF(real x)
     if (x != 0)
     {
         real y = sqrt(x);
-        return 2*exp(-x)/(sqrt(M_PI)*x) - gmx_erfd(y)/(x*y);
+        return 2*exp(-x)/(sqrt(M_PI)*x) - std::erf(static_cast<double>(y))/(x*y);
     }
     else
     {
@@ -515,7 +507,7 @@ TEST_F(SimdMathTest, gmxSimdPmecorrForceR)
 real ref_pmecorrV(real x)
 {
     real y = sqrt(x);
-    return gmx_erfd(y)/y;
+    return std::erf(static_cast<double>(y))/y;
 }
 
 // The PME corrections will be added to ~1/r, so absolute tolerance of EPS is fine.
@@ -608,17 +600,14 @@ TEST_F(SimdMathTest, gmxSimdLogSingleaccuracyR)
     GMX_EXPECT_SIMD_FUNC_NEAR(ref_log, gmx_simd_log_singleaccuracy_r);
 }
 
-// MSVC does not support exp2(), so we have no reference to test against
-#ifndef _MSC_VER
 TEST_F(SimdMathTest, gmxSimdExp2SingleaccuracyR)
 {
     /* Increase the allowed error by the difference between the actual precision and single */
     setUlpTol(ulpTol_ * (1LL << (std::numeric_limits<real>::digits-std::numeric_limits<float>::digits)));
 
     setRange(-100, 100);
-    GMX_EXPECT_SIMD_FUNC_NEAR(ref_exp2, gmx_simd_exp2_singleaccuracy_r);
+    GMX_EXPECT_SIMD_FUNC_NEAR(std::exp2, gmx_simd_exp2_singleaccuracy_r);
 }
-#endif
 
 TEST_F(SimdMathTest, gmxSimdExpSingleaccuracyR)
 {
index 018403e45abadfc31a04c3085efdd1f6aa402f00..1b18e4a2084508b32a1f9f764bae90c9a2cb1e7c 100644 (file)
@@ -136,7 +136,7 @@ double v_q_ewald_lr(double beta, double r)
     }
     else
     {
-        return gmx_erfd(beta*r)/r;
+        return std::erf(beta*r)/r;
     }
 }
 
@@ -353,7 +353,7 @@ real ewald_spline3_table_scale(const interaction_const_t *ic)
         real   sc_q;
 
         /* Energy tolerance: 0.1 times the cut-off jump */
-        etol  = 0.1*gmx_erfc(ic->ewaldcoeff_q*ic->rcoulomb);
+        etol  = 0.1*std::erfc(ic->ewaldcoeff_q*ic->rcoulomb);
 
         sc_q  = spline3_table_scale(erf_x_d3, ic->ewaldcoeff_q, etol);
 
@@ -925,11 +925,11 @@ static void fill_table(t_tabledata *td, int tp, const t_forcerec *fr,
                 break;
             case etabEwald:
             case etabEwaldSwitch:
-                Vcut  = gmx_erfc(ewc*rc)/rc;
+                Vcut  = std::erfc(ewc*rc)/rc;
                 break;
             case etabEwaldUser:
                 /* Only calculate minus the reciprocal space contribution */
-                Vcut  = -gmx_erf(ewc*rc)/rc;
+                Vcut  = -std::erf(ewc*rc)/rc;
                 break;
             case etabRF:
             case etabRF_ZERO:
@@ -1068,14 +1068,14 @@ static void fill_table(t_tabledata *td, int tp, const t_forcerec *fr,
                 break;
             case etabEwald:
             case etabEwaldSwitch:
-                Vtab  = gmx_erfc(ewc*r)/r;
-                Ftab  = gmx_erfc(ewc*r)/r2+exp(-(ewc*ewc*r2))*ewc*M_2_SQRTPI/r;
+                Vtab  = std::erfc(ewc*r)/r;
+                Ftab  = std::erfc(ewc*r)/r2+exp(-(ewc*ewc*r2))*ewc*M_2_SQRTPI/r;
                 break;
             case etabEwaldUser:
             case etabEwaldUserSwitch:
                 /* Only calculate the negative of the reciprocal space contribution */
-                Vtab  = -gmx_erf(ewc*r)/r;
-                Ftab  = -gmx_erf(ewc*r)/r2+exp(-(ewc*ewc*r2))*ewc*M_2_SQRTPI/r;
+                Vtab  = -std::erf(ewc*r)/r;
+                Ftab  = -std::erf(ewc*r)/r2+exp(-(ewc*ewc*r2))*ewc*M_2_SQRTPI/r;
                 break;
             case etabLJ6Ewald:
                 Vtab  = -r6*exp(-ewclj*ewclj*r2)*(1 + ewclj*ewclj*r2 + pow4(ewclj)*r2*r2/2);
index 0c237a122006200f8b8841650c0323769941b9bf..418a07c6f1ea53e45dac099198d13e22316e32a7 100644 (file)
@@ -3,7 +3,7 @@
  *
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
@@ -42,6 +42,8 @@
 #include <stdio.h>
 #include <string.h>
 
+#include <cmath>
+
 #include "gromacs/fileio/strdb.h"
 #include "gromacs/legacyheaders/copyrite.h"
 #include "gromacs/math/utilities.h"
@@ -406,7 +408,7 @@ char *gmx_atomprop_element(gmx_atomprop_t aps, int atomnumber)
     set_prop(aps, epropElement);
     for (i = 0; (i < ap->prop[epropElement].nprop); i++)
     {
-        if (gmx_nint(ap->prop[epropElement].value[i]) == atomnumber)
+        if (std::round(ap->prop[epropElement].value[i]) == atomnumber)
         {
             return ap->prop[epropElement].atomnm[i];
         }
@@ -424,7 +426,7 @@ int gmx_atomprop_atomnumber(gmx_atomprop_t aps, const char *elem)
     {
         if (gmx_strcasecmp(ap->prop[epropElement].atomnm[i], elem) == 0)
         {
-            return gmx_nint(ap->prop[epropElement].value[i]);
+            return std::round(ap->prop[epropElement].value[i]);
         }
     }
     return -1;