Optionally use CORE-MATH correctly rounded elementary functions

CMakeLists.txt

@@ -74,6 +74,9 @@ option(USE_SIGNALS
 option(USE_GMP
        "Use GNU MP Bignum library (LGPL)"
        ON)
+option(USE_CRMATH
+       "Use correctly rounded CORE_MATH elementary functions"
+       OFF)
 option(USE_TCMALLOC
        "Use Google tcmalloc instead of default malloc"
        ${DEFAULT_USE_TCMALLOC})

cmake/Config.cmake

@@ -123,6 +123,36 @@ set(CMAKE_EXTRA_INCLUDE_FILES ${CMAKE_EXTRA_INCLUDE_FILES} math.h wchar.h)
 #  set(CMAKE_REQUIRED_FLAGS ${CMAKE_REQUIRED_FLAGS} -Wno-builtin-declaration-mismatch)
 #endif()
 
+if(USE_CRMATH) 
+  if(NOT TARGET crmath)  # already built?
+    # Repo existance check
+    set(CRMATH_REPO "https://github.com/indigobio/core-math-binary64.git")
+    set(CRMATH_TAG  "main")
+    execute_process(
+      COMMAND git ls-remote --exit-code --heads --tags "${CRMATH_REPO}" "${CRMATH_TAG}"
+      RESULT_VARIABLE git_check
+      OUTPUT_QUIET
+      ERROR_QUIET
+    )
+    if(git_check EQUAL 0)
+      include(FetchContent)
+      FetchContent_Declare(
+        crmath_repo
+        GIT_REPOSITORY "${CRMATH_REPO}"
+        GIT_TAG        "${CRMATH_TAG}"
+      )
+      FetchContent_Populate(crmath_repo)
+      message(STATUS "Looking for " ${CRMATH_REPO} " - found")
+      # define source directory for building libcrmath
+      set(CRMATH_REPO_SOURCE_DIR "${CMAKE_BINARY_DIR}/_deps/crmath_repo-src")
+      set(HAVE_CRMATH 1)
+    else()
+      message(STATUS "Looking for " ${CRMATH_REPO} " - not found")
+      set(HAVE_CRMATH 0)
+    endif()  # git_check EQUAL 0
+  endif()  # NOT TARGET crmath
+endif() # USE_CRMATH
+
 ################
 # Types
 SET(CMAKE_TRY_COMPILE_TARGET_TYPE_SAVE ${CMAKE_TRY_COMPILE_TARGET_TYPE})

src/CMakeLists.txt

@@ -492,6 +492,68 @@ if(APPLE)
   set(LIBSWIPL_LIBRARIES ${LIBSWIPL_LIBRARIES} "-framework CoreFoundation")
 endif()
 
+if(HAVE_CRMATH)
+  if(NOT TARGET crmath)  # already built?
+    # CRMATH_REPO_SOURCE_DIR defined in cmake/Config.cmake
+	# Adds required subset from repo, unused functions in comments
+    add_library(crmath STATIC 
+  ${CRMATH_REPO_SOURCE_DIR}/src/acos/acos.c
+  ${CRMATH_REPO_SOURCE_DIR}/src/acosh/acosh.c
+  ${CRMATH_REPO_SOURCE_DIR}/src/asin/asin.c
+  ${CRMATH_REPO_SOURCE_DIR}/src/asinh/asinh.c
+#  src/asinpi/asinpi.c
+  ${CRMATH_REPO_SOURCE_DIR}/src/atan/atan.c
+  ${CRMATH_REPO_SOURCE_DIR}/src/atan2/atan2.c
+  ${CRMATH_REPO_SOURCE_DIR}/src/atan2/tint.h
+#  src/atan2pi/atan2pi.c
+#  src/atan2pi/tint.h
+  ${CRMATH_REPO_SOURCE_DIR}/src/atanh/atanh.c
+#  src/atanpi/atanpi.c
+#  src/cbrt/cbrt.c
+  ${CRMATH_REPO_SOURCE_DIR}/src/cos/cos.c
+  ${CRMATH_REPO_SOURCE_DIR}/src/cosh/cosh.c
+#  src/cospi/cospi.c
+  ${CRMATH_REPO_SOURCE_DIR}/src/erf/erf.c
+  ${CRMATH_REPO_SOURCE_DIR}/src/erfc/erfc.c
+  ${CRMATH_REPO_SOURCE_DIR}/src/exp/exp.c
+#  src/exp10/exp10.c
+#  src/exp10m1/exp10m1.c
+#  src/exp2/exp2.c
+#  src/exp2m1/exp2m1.c
+#  src/expm1/expm1.c
+#  src/hypot/hypot.c
+  ${CRMATH_REPO_SOURCE_DIR}/src/log/dint.h
+  ${CRMATH_REPO_SOURCE_DIR}/src/log/log.c
+  ${CRMATH_REPO_SOURCE_DIR}/src/log10/dint.h
+  ${CRMATH_REPO_SOURCE_DIR}/src/log10/log10.c
+#  src/log10p1/dint.h
+#  src/log10p1/log10p1.c
+#  src/log1p/dint.h
+#  src/log1p/log1p.c
+#  src/log2/log2.c
+#  src/log2p1/dint_log2p1.h
+#  src/log2p1/log2p1.c
+  ${CRMATH_REPO_SOURCE_DIR}/src/pow/dint.h
+  ${CRMATH_REPO_SOURCE_DIR}/src/pow/pow.c
+  ${CRMATH_REPO_SOURCE_DIR}/src/pow/pow.h
+  ${CRMATH_REPO_SOURCE_DIR}/src/pow/qint.h
+#  src/rsqrt/rsqrt.c
+  ${CRMATH_REPO_SOURCE_DIR}/src/sin/sin.c
+#  src/sincos/sincos.c
+  ${CRMATH_REPO_SOURCE_DIR}/src/sinh/sinh.c
+#  src/sinpi/sinpi.c
+  ${CRMATH_REPO_SOURCE_DIR}/src/tan/tan.c
+  ${CRMATH_REPO_SOURCE_DIR}/src/tanh/tanh.c
+#  src/tanpi/tanpi.c
+#  src/tgamma/tgamma.c
+    )
+    # not sure what compiler flags to use for various platforms/compilers so
+    # just using base set from Bob Burgers repo 
+    target_compile_options(crmath PRIVATE -O3 -ffp-contract=on -fno-math-errno)
+    set(LIBSWIPL_LIBRARIES ${LIBSWIPL_LIBRARIES} crmath)
+  endif()  # NOT TARGET crmath
+endif()  # HAVE_CRMATH
+
 # build swipl
 add_executable(swipl${PGO_SUFFIX} ${SWIPL_SRC})
 target_link_libraries(swipl${PGO_SUFFIX} ${SWIPL_LIBRARIES} libswipl${PGO_SUFFIX})

src/config.h.cmake

@@ -289,3 +289,6 @@
 #cmakedefine SWIPL_TMP_DIR "@SWIPL_TMP_DIR@"
 
 #cmakedefine MSYS2 @MSYS2@
+
+#cmakedefine HAVE_CRMATH @HAVE_CRMATH@
+

src/pl-arith.c

@@ -1923,6 +1923,10 @@ FE_TONEAREST mode is in effect which is not a proven fact.
 This set of correctly rounded functions (cr_exp, cr_log, cr_sin, cr_cos,
 etc.) are used to redefine the standard functions as used in the various
 'ar_...' arithmetic functions.
+
+There are two possible implementations of the cr_xxx elementary functions.
+If HAVE_CRMATH is true, they are part of the core library as defined in
+pl-crmath.h. Otherwise they are implemented here as described above. 
 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 
 
@@ -1942,6 +1946,9 @@ static double cr_##func(double in) \
   return result; \
 }
 
+#if HAVE_CRMATH
+#include <pl-crmath.h>
+#else
 // #define unary function with abs(result) =< 1.0
 #define CR_FUNC_1(func) \
 static double cr_##func(double in) \
@@ -1987,31 +1994,22 @@ cr_exp(double in)
   if ( roundMode != FE_TONEAREST ) fesetround(roundMode);
   return result;
 }
-#define exp(x) cr_exp(x)
 
 CR_FUNC(log)
-#define log(x) cr_log(x)
 
 CR_FUNC(log10)
-#define log10(x) cr_log10(x)
 
 CR_FUNC_1(sin)
-#define sin(x) cr_sin(x)
 
 CR_FUNC_1(cos)
-#define cos(x) cr_cos(x)
 
 CR_FUNC(tan)
-#define tan(x) cr_tan(x)
 
 CR_FUNC(asin)
-#define asin(x) cr_asin(x)
 
 CR_FUNC(acos)
-#define acos(x) cr_acos(x)
 
 CR_FUNC(atan)
-#define atan(x) cr_atan(x)
 
 // special case for binary atan2 function with standard rounding, finite result
 static double
@@ -2028,10 +2026,8 @@ cr_atan2(double y, double x)
   if ( roundMode != FE_TONEAREST ) fesetround(roundMode);
   return result;
 }
-#define atan2(y,x) cr_atan2(y,x)
 
 CR_FUNC(sinh)
-#define sinh(x) cr_sinh(x)
 
 // special case for unary cosh with lower limit of 1.0
 static double
@@ -2055,10 +2051,8 @@ cr_cosh(double in)
   if ( roundMode != FE_TONEAREST ) fesetround(roundMode);
   return result;
 }
-#define cosh(x) cr_cosh(x)
 
 CR_FUNC_1(tanh)
-#define tanh(x) cr_tanh(x)
 
 // special case for binary function pow with lower limit of 0.0
 static double
@@ -2084,7 +2078,6 @@ cr_pow(double base, double exp)
   if ( roundMode != FE_TONEAREST ) fesetround(roundMode);
   return result;
 }
-#define pow(b,e) cr_pow(b,e)
 
 // special case for unary erf with bounded by +/- 1.0
 static double
@@ -2111,19 +2104,34 @@ cr_erf(double in)
   if ( roundMode != FE_TONEAREST ) fesetround(roundMode);
   return result;
 }
-#define erf(x) cr_erf(x)
 
 static double cr_erfc(double in)
 { return 1.0 - cr_erf(in);
 }
-#define erfc(x) cr_erfc(x)
+#endif /*HAVE_CRMATH*/
 
-CR_FUNC(lgamma)
+CR_FUNC(lgamma)  // for both cases (not in libcrmath, Windows issue?)
+
+#define exp(x) cr_exp(x)
+#define log(x) cr_log(x)
+#define log10(x) cr_log10(x)
+#define sin(x) cr_sin(x)
+#define cos(x) cr_cos(x)
+#define tan(x) cr_tan(x)
+#define asin(x) cr_asin(x)
+#define acos(x) cr_acos(x)
+#define atan(x) cr_atan(x)
+#define atan2(y,x) cr_atan2(y,x)
+#define sinh(x) cr_sinh(x)
+#define cosh(x) cr_cosh(x)
+#define tanh(x) cr_tanh(x)
+#define pow(b,e) cr_pow(b,e)
+#define erf(x) cr_erf(x)
+#define erfc(x) cr_erfc(x)
 #define lgamma(x) cr_lgamma(x)
 
 #endif /*O_ROUND_UP_DOWN*/
 
-
 /* Unary functions requiring double argument */
 
 #define UNAIRY_FLOAT_FUNCTION(name, op) \

src/pl-crmath.h

@@ -0,0 +1,55 @@
+/*
+ *	Required subset from crmath repo. (License??)
+ *  Files not in subset are commented out.
+ */
+
+#pragma once
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+double cr_acos(double x);
+double cr_acosh(double x);
+//double cr_acospi(double x);
+double cr_asin(double x);
+double cr_asinh(double x);
+//double cr_asinpi(double x);
+double cr_atan(double x);
+double cr_atan2(double y0, double x0);
+//double cr_atan2pi(double y, double x);
+double cr_atanh(double x);
+//double cr_atanpi(double x);
+//double cr_cbrt(double x);
+double cr_cos(double x);
+double cr_cosh(double x);
+//double cr_cospi(double x);
+double cr_erf(double x);
+double cr_erfc(double x);
+double cr_exp(double x);
+//double cr_exp10(double x);
+//double cr_exp10m1(double x);
+//double cr_exp2(double x);
+//double cr_exp2m1(double x);
+//double cr_expm1(double x);
+//double cr_hypot(double x, double y);
+double cr_log(double x);
+double cr_log10(double x);
+//double cr_log10p1(double x);
+//double cr_log1p(double x);
+//double cr_log2(double x);
+//double cr_log2p1(double x);
+double cr_pow(double x, double y);
+//double cr_rsqrt(double x);
+double cr_sin(double x);
+//double cr_sincos(double x, double* s, double* c);
+double cr_sinh(double x);
+//double cr_sinpi(double x);
+double cr_tan(double x);
+double cr_tanh(double x);
+//double cr_tanpi(double x);
+//double cr_tgamma(double x);
+
+#ifdef __cplusplus
+}
+#endif