New cppduals from https://gitlab.com/tesch1/cppduals

downloaded may 1st, 2025
2025-05-01 21:09:19 +02:00 · 2025-05-01 21:09:19 +02:00 · 64284d0fda
parent 3680a8c24d
commit 64284d0fda
22 changed files with 2056 additions and 1082 deletions
--- a/src/include/cppduals/.appveyor.yml
+++ b/src/include/cppduals/.appveyor.yml
@ -3,9 +3,8 @@ clone_folder: c:\projects\cppduals
 clone_depth: 3
 image:
 #- Visual Studio 2013
 #- Visual Studio 2015
 - Visual Studio 2017
 #- Visual Studio 2022
 configuration:
 - Release
@ -14,12 +13,10 @@ configuration:
 # Do not build feature branch with open Pull Requests
 skip_branch_with_pr: true
 # skip unsupported combinations
 init:
 - echo %APPVEYOR_BUILD_WORKER_IMAGE%
 - if "%APPVEYOR_BUILD_WORKER_IMAGE%"=="Visual Studio 2013" ( set generator="Visual Studio 12 2013" )
 - if "%APPVEYOR_BUILD_WORKER_IMAGE%"=="Visual Studio 2015" ( set generator="Visual Studio 14 2015" )
 - if "%APPVEYOR_BUILD_WORKER_IMAGE%"=="Visual Studio 2017" ( set generator="Visual Studio 15 2017" )
 - if "%APPVEYOR_BUILD_WORKER_IMAGE%"=="Visual Studio 2022" ( set generator="Visual Studio 17 2022" )
 - echo %generator%
 before_build:
--- a/src/include/cppduals/.gitlab-ci.yml
+++ b/src/include/cppduals/.gitlab-ci.yml
@ -1,5 +1,4 @@
-#image: ubuntu:19.04
+image: fedora:41
 image: fedora:30
 variables:
  GIT_DEPTH: 3
@ -9,79 +8,105 @@ stages:
 - build
 - test
 - cover
- publish
+- pages
-
+- upload
-before_script:
+- release
 #- apt-get update --yes
 #- apt-get install --yes cmake g++ git doxygen lcov graphviz
 - dnf install -y gcc-c++ make cmake git doxygen lcov graphviz
 build:
  stage: build
 #  variables:
 #    CC: clang
 #    CXX: clang++
  script:
  - dnf install -y gcc-c++ make cmake git doxygen lcov graphviz
  - echo $CXX
-  - cmake -Bbuild -H. -DCPPDUALS_TESTING=ON
+  - gcc -march=native -dM -E - </dev/null | grep -i sse
  - cmake -Bbuild -H. -DCPPDUALS_TESTING=ON -DCMAKE_BUILD_TYPE=RelWithDebInfo
  - cmake --build build
  - cmake -Bbuild-latest -H. -DCPPDUALS_TESTING=ON -DCPPDUALS_EIGEN_LATEST=ON
  - cmake --build build-latest
  artifacts:
    expire_in: 1 week
    paths:
-    - build
+    - build/
-    - build-latest
+  only:
  - merge_requests
  - master
 test:
  stage: test
  script:
  - dnf install -y gcc-c++ make cmake git doxygen lcov graphviz
  - cmake --build build --target test
  - cmake --build build-latest --target test
  dependencies:
  - build
  only:
  - merge_requests
  - master
 cover:
  stage: cover
  script:
-  - cmake -Bbuild-cov -H. -DCODE_COVERAGE=ON -DCPPDUALS_TESTING=ON
+  - dnf install -y gcc-c++ make cmake git doxygen lcov graphviz
  - cmake -Bbuild-cov -H. -DCPPDUALS_TESTING=ON -DCPPDUALS_CODE_COVERAGE=ON -DCMAKE_BUILD_TYPE=RelWithDebInfo
  - cmake --build build-cov --target cov
  - cmake --build build-cov --target cov-html
-  coverage: '/Total:|\w*\d+\.\d+/'
+  coverage: "/Total:|\\w*\\d+\\.\\d+/"
  artifacts:
-    expire_in: 1 day
+    expire_in: 1 week
    paths:
-    - build-cov
+    - build-cov/
  only:
  - merge_requests
 pages:
  script:
-  - cmake -Bbuild -H. -DCODE_COVERAGE=ON -DCPPDUALS_TESTING=ON
+  - dnf install -y gcc-c++ make cmake git doxygen lcov graphviz
-  - cmake --build build --target cppduals_docs
+  - cmake -Bbuild-cov -H. -DCODE_COVERAGE=ON -DCPPDUALS_TESTING=ON
-  - cmake --build build --target cov-html
+  - cmake --build build-cov --target cov-html
-  - mv build/docs public/
+  - cmake --build build-cov --target cppduals_docs
-  - mv build/coverage public/
+  - mv build-cov/coverage public/
-  coverage: '/Total:|\w*\d+\.\d+/'
+  - mv build-cov/docs public/
  coverage: "/Total:|\\w*\\d+\\.\\d+/"
  artifacts:
    paths:
-    - public
+    - public/
  only:
  - master
-publish:
+variables:
-  stage: publish
+  # Package version can only contain numbers (0-9), and dots (.).
-  dependencies:
+  # Must be in the format of X.Y.Z, i.e. should match /\A\d+\.\d+\.\d+\z/ regular expresion.
-  - build
+  # See https://docs.gitlab.com/ee/user/packages/generic_packages/#publish-a-package-file
-  environment:
+  PACKAGE_VERSION: $CI_COMMIT_TAG
-    name: publish
+  PACKAGE_REGISTRY_URL: "${CI_API_V4_URL}/projects/${CI_PROJECT_ID}/packages/generic/release/${PACKAGE_VERSION}"
-  only:
+  HEADER_ONLY_PACKAGE: "cppduals-h-${CI_COMMIT_TAG#v}.tgz"
-  - /^v\d+\.\d+\.\d+$/
+
-  except:
+upload:
-  - branches
+  stage: upload
-  before_script:
+  image: curlimages/curl:latest
-  - dnf install -y python3-requests
+  rules:
    - if: $CI_COMMIT_TAG
  script:
-#  - ln -s cppduals-h-${CI_BUILD_TAG#v} .
+    - |
-#  - tar czvhf cppduals-h-${CI_BUILD_TAG#v}.tgz cppduals-h-${CI_BUILD_TAG#v}/duals cppduals-h-${CI_BUILD_TAG#v}/CMakeLists.txt
+      tar czvf ${HEADER_ONLY_PACKAGE} duals/
-  - tar czvf cppduals-h-${CI_BUILD_TAG#v}.tgz duals CMakeLists.txt
+      curl --header "JOB-TOKEN: ${CI_JOB_TOKEN}" \
-  - ./doc/gitlab-release --message "Release ${CI_BUILD_TAG}" cppduals-h-${CI_BUILD_TAG#v}.tgz
+        --upload-file ${HEADER_ONLY_PACKAGE} \
                      ${PACKAGE_REGISTRY_URL}/${HEADER_ONLY_PACKAGE}
 release:
  # Caution, as of 2021-02-02 these assets links require a login, see:
  # https://gitlab.com/gitlab-org/gitlab/-/issues/299384
  stage: release
  image: registry.gitlab.com/gitlab-org/release-cli:latest
  rules:
    - if: $CI_COMMIT_TAG
  script:
    - |
      release-cli create \
        --name "Release $CI_COMMIT_TAG" \
        --tag-name $CI_COMMIT_TAG \
        --assets-link "{\"name\":\"${HEADER_ONLY_PACKAGE}\",\"url\":\"${PACKAGE_REGISTRY_URL}/${HEADER_ONLY_PACKAGE}\"}"
 sast:
  variables:
    SAST_DEFAULT_ANALYZERS: flawfinder
  stage: test
 include:
 - template: Security/SAST.gitlab-ci.yml
--- a/src/include/cppduals/CMakeLists.txt
+++ b/src/include/cppduals/CMakeLists.txt
@ -1,14 +1,17 @@
 #
 # CMake for cppduals
 #
-cmake_minimum_required (VERSION 3.1)
+cmake_minimum_required (VERSION 3.14)
 project (cppduals
-  VERSION 0.3.1
+  VERSION 0.6.2
  LANGUAGES C CXX
  )
 include (GNUInstallDirs)
-set (CMAKE_CXX_STANDARD 11 CACHE STRING "Which C++ standard to test against.")
+if (NOT CMAKE_CXX_STANDARD)
  set (CMAKE_CXX_STANDARD 17 CACHE STRING "Which C++ standard to test against.")
 endif()
 message (STATUS "CXX_STANDARD: ${CMAKE_CXX_STANDARD}")
 set (CMAKE_CXX_STANDARD_REQUIRED ON)
 set (CMAKE_DISABLE_IN_SOURCE_BUILD ON)
 if (CMAKE_SOURCE_DIR STREQUAL CMAKE_CURRENT_SOURCE_DIR)
@ -38,11 +41,11 @@ if (CPPDUALS_STANDALONE AND
  message (STATUS "No install prefix specified; using '${CMAKE_INSTALL_PREFIX}'")
 endif ()
-set_property (CACHE CMAKE_CXX_STANDARD PROPERTY STRINGS 11 14 17 20)
+#set_property (CACHE CMAKE_CXX_STANDARD PROPERTY STRINGS 11 14 17 20)
 set_property (CACHE CMAKE_CXX_STANDARD PROPERTY STRINGS 20)
 option (CPPDUALS_TESTING "Enable testing" OFF)
 option (CPPDUALS_BENCHMARK "Enable benchmarking" OFF)
 option (CPPDUALS_EIGEN_LATEST "Eigen latest" OFF)
 option (CPPDUALS_USE_LIBCXX "When testing use flags for libc++" OFF)
 set (EIGEN3_INCLUDE_DIRS "" CACHE PATH "Path to Eigen includes" )
@ -91,7 +94,7 @@ add_library (cppduals::duals ALIAS cppduals_duals)
 #
 if (CPPDUALS_TESTING)
-  cmake_minimum_required (VERSION 3.10) # need gtest_discover_tests
+  cmake_minimum_required (VERSION 3.14) # need gtest_discover_tests
  file (MAKE_DIRECTORY "${PROJECT_BINARY_DIR}/thirdparty")
  # generator name
@ -109,7 +112,6 @@ if (CPPDUALS_TESTING)
                             "-DCMAKE_BUILD_TYPE=${CMAKE_BUILD_TYPE}"
                             "-DCMAKE_MAKE_PROGRAM=${CMAKE_MAKE_PROGRAM}"
                             "-DCPPDUALS_BENCHMARK=${CPPDUALS_BENCHMARK}"
                             "-DCPPDUALS_EIGEN_LATEST=${CPPDUALS_EIGEN_LATEST}"
                             "-DCPPDUALS_USE_LIBCXX=${CPPDUALS_USE_LIBCXX}"
                             "${PROJECT_SOURCE_DIR}/thirdparty"
    RESULT_VARIABLE DEPS_CONFIG_RESULT
@ -144,18 +146,27 @@ endif ()
 # Code Coverage Configuration
 #
 add_library (cppduals_coverage_config INTERFACE)
-option (CODE_COVERAGE "Enable coverage reporting" OFF)
+option (CPPDUALS_CODE_COVERAGE "Enable coverage reporting" OFF)
-if (CODE_COVERAGE AND CMAKE_CXX_COMPILER_ID MATCHES "GNU|Clang")
+if (CPPDUALS_CODE_COVERAGE AND NOT CPPDUALS_TESTING)
  message(FATAL_ERROR "CPPDUALS_CODE_COVERAGE requires CPPDUALS_TESTING to be enabled")
 endif()
 if (CPPDUALS_CODE_COVERAGE)
  if (CMAKE_CXX_COMPILER_ID MATCHES "GNU|Clang")
    message(STATUS "Enabling coverage reporting")
    # Add required flags (GCC & LLVM/Clang)
    target_compile_options (cppduals_coverage_config INTERFACE
      -O0        # no optimization
      -g         # generate debug info
      --coverage # sets all required flags
      -fprofile-arcs 
      -ftest-coverage
    )
-  if (CMAKE_VERSION VERSION_GREATER_EQUAL 3.13)
+    target_link_options (cppduals_coverage_config INTERFACE 
-    target_link_options (cppduals_coverage_config INTERFACE --coverage)
+      --coverage
-  else ()
+    )
-    target_link_libraries (cppduals_coverage_config INTERFACE --coverage)
+  else()
    # error out if coverage is enabled but compiler is not supported
    message(FATAL_ERROR "Coverage reporting is disabled for unknown compilers.")
  endif ()
 endif ()
@ -221,3 +232,8 @@ if (ETAGS)
    COMMAND ${ETAGS} --language=c++ --append `find ${PROJECT_BINARY_DIR}/thirdparty/eigenX/src/eigenX -type f`
    )
 endif (ETAGS)
 #
 # Packaging
 #
 include (CPack)
--- a/src/include/cppduals/README.md
+++ b/src/include/cppduals/README.md
@ -96,25 +96,6 @@ to specify library dependencies:
  target_link_libraries (your_target PRIVATE cppduals::duals)
 ```
 Older versions of CMake can achieve a similar result using the ``ExternalProject``
 family of commands and modifying the global preprocessor search path:
 ```cmake
  include(ExternalProject)
  # Have CMake download the library headers only
  set (CPPDUALS_TAG v0.4.1)
  set (CPPDUALS_MD5 7efe49496b8d0e3d3ffbcd3c68f542f3)
  ExternalProject_Add (cppduals
    URL https://gitlab.com/tesch1/cppduals/-/archive/${CPPDUALS_TAG}/cppduals-${CPPDUALS_TAG}.tar.bz2
    URL_HASH MD5=${CPPDUALS_MD5}
    CONFIGURE_COMMAND "" BUILD_COMMAND "" INSTALL_COMMAND "" )
  # Make include directory globally visible
  ExternalProject_Get_Property (cppduals source_dir)
  include_directories (${source_dir}/)
 ```
 Alternatively, `cppduals` supports installation and discovery via the
 `find_package` utility. First, download and install the library to a
 location of your choosing:
@ -299,6 +280,50 @@ thus licensed under [MPL-2](http://www.mozilla.org/MPL/2.0/FAQ.html) .
 ChangeLog
 =========
 v0.6.0
 ======
 - target at least c++17.
 - tested with eigen 3.3.7 & 3.3.8.
 - rearrange how we're (illegally using/abusing) std:: by moving many implementation details
  into duals::detail tricks - makes everything work with LIBCPP on osx now (as of Sequoia 15.2).
 - upgrade google test & benchmark libraries.
 - update to more modern cmake, at least 3.14 required now.
 - dont try to use Eigen's .exp() with duals, not supported.
 Known Issues
 - fmt library support is very out of date - should update to support std::format.
 - coverage is not working.
 - SSE/AVX not working - open to fixes but dont have a machine to test on at the moment.
 v0.5.4
 ======
 - upgrade google test library
 v0.5.3
 ======
 - fix some problem with pow()
 v0.5.2
 ======
 - change optional libfmt print support to fmt 7.1.3 (from 6.x)
 - change default build standard to c++14.
 v0.5.1
 ======
 - packaging cleanup
 v0.5.0
 ======
 - added ARM NEON support.  tested on apple M1 - note apple's BLAS gemm
  is ~2x faster than Eigen-generated matrix mul :(
 - fixed atan2 and pow
 v0.4.1
 ======
@ -351,4 +376,4 @@ Todo
 - Add multi-variate differentiation capability.
 - Non-x86_64 (CUDA/AltiVec/HIP/NEON/...) vectorization.
 - Higher-order derivatives.
-
+- finish NEON (is why clang failing on osx?  or just make it faster!)
--- a/src/include/cppduals/duals/Makefile.am
+++ b/src/include/cppduals/duals/Makefile.am
@ -1,6 +0,0 @@
 ## Process this file with automake to produce Makefile.in
 noinst_HEADERS = \
 	dual
 MAINTAINERCLEANFILES = Makefile.in
--- a/src/include/cppduals/duals/arch/SSE/Dual.h
+++ b/src/include/cppduals/duals/arch/SSE/Dual.h
@ -10,6 +10,7 @@
 #ifndef EIGEN_DUAL_SSE_H
 #define EIGEN_DUAL_SSE_H
 #include <pmmintrin.h> // SSE3
 namespace Eigen {
@ -150,17 +151,9 @@ template<> EIGEN_STRONG_INLINE void prefetch<duals::dual<float> >(const duals::d
 template<> EIGEN_STRONG_INLINE duals::dual<float>  pfirst<Packet2df>(const Packet2df& a)
 {
  #if EIGEN_GNUC_AT_MOST(4,3)
  // Workaround gcc 4.2 ICE - this is not performance wise ideal, but who cares...
  // This workaround also fix invalid code generation with gcc 4.3
  EIGEN_ALIGN16 duals::dual<float> res[2];
  _mm_store_ps((float*)res, a.v);
  return res[0];
  #else
  duals::dual<float> res;
  _mm_storel_pi((__m64*)&res, a.v);
  return res;
  #endif
 }
 template<> EIGEN_STRONG_INLINE Packet2df preverse(const Packet2df& a)
--- a/src/include/cppduals/duals/dual
+++ b/src/include/cppduals/duals/dual
--- a/src/include/cppduals/duals/dual_eigen
+++ b/src/include/cppduals/duals/dual_eigen
@ -141,7 +141,6 @@ namespace Eigen {
 template<typename T>
 struct NumTraits<duals::dual<T> > : GenericNumTraits<T>
 {
  typedef typename NumTraits<T>::Real ReallyReal;
  typedef duals::dual<T> Real;
  typedef duals::dual<T> Literal;
  typedef duals::dual<T> Nested;
@ -159,13 +158,13 @@ struct NumTraits<duals::dual<T> > : GenericNumTraits<T>
  EIGEN_DEVICE_FUNC
  static inline Real epsilon()        { return Real(NumTraits<T>::epsilon()); }
  EIGEN_DEVICE_FUNC
-  static inline ReallyReal dummy_precision()  { return NumTraits<T>::dummy_precision(); }
+  static inline Real dummy_precision()  { return NumTraits<T>::dummy_precision(); }
  EIGEN_DEVICE_FUNC
-  static inline ReallyReal highest()  { return NumTraits<T>::highest(); }
+  static inline Real highest()  { return NumTraits<T>::highest(); }
  EIGEN_DEVICE_FUNC
-  static inline ReallyReal lowest()   { return NumTraits<T>::lowest(); }
+  static inline Real lowest()   { return NumTraits<T>::lowest(); }
  EIGEN_DEVICE_FUNC
-  static inline ReallyReal digits10() { return NumTraits<T>::digits10(); }
+  static inline int digits10() { return NumTraits<T>::digits10(); }
 };
 #if !defined(CPPDUALS_NO_EIGEN_PROMOTION)
@ -218,6 +217,20 @@ using duals::dconj;
 namespace internal {
 #if 0
 // For MatrixExponential.h to treat duals::dual<T> as a known type and compile it
 template<typename T> struct is_exp_known_type;
 template<typename MatrixType, typename T> struct matrix_exp_computeUV;
 template<typename T> struct is_exp_known_type<duals::dual<T>> : is_exp_known_type<T> {};
 template <typename MatrixType, typename T>
 struct matrix_exp_computeUV<MatrixType, duals::dual<T> > : matrix_exp_computeUV<MatrixType, T>
 {
  typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
 };
 #endif
 #if 1
 // this is used by the packet math for SSE to copy raw duals around.
 template<typename T>
 struct real_impl<duals::dual<T> >
 {
@ -251,6 +264,43 @@ struct real_ref_retval<duals::dual<Scalar>>
  typedef typename NumTraits<Scalar>::Real & type;
 };
 #else ////
 template<typename T>
 struct real_impl<duals::dual<T> >
 {
  typedef duals::dual<T> RealScalar;
  EIGEN_DEVICE_FUNC
  static inline RealScalar run(const duals::dual<T>& x)
  {
    return x;
  }
 };
 template<typename Scalar>
 struct real_ref_impl<duals::dual<Scalar>>
 {
  typedef typename NumTraits<duals::dual<Scalar> >::Real RealScalar;
  EIGEN_DEVICE_FUNC
  static inline RealScalar & run(duals::dual<Scalar> & x)
  {
    return reinterpret_cast<RealScalar*>(&x)[0];
  }
  EIGEN_DEVICE_FUNC
  static inline const RealScalar & run(const duals::dual<Scalar> & x)
  {
    return reinterpret_cast<const RealScalar *>(&x)[0];
  }
 };
 template<typename Scalar>
 struct real_ref_retval<duals::dual<Scalar>>
 {
  typedef typename NumTraits<duals::dual<Scalar>>::Real & type;
 };
 #endif /// 0
 template<typename T> struct scalar_random_op<duals::dual<T>>
 {
  EIGEN_EMPTY_STRUCT_CTOR(scalar_random_op)
@ -384,7 +434,8 @@ template<typename RealScalar,bool Conj> struct dconj_helper<RealScalar, duals::d
 #elif defined(EIGEN_VECTORIZE_ALTIVEC) || defined(EIGEN_VECTORIZE_VSX)
 //  #include "duals/arch/AltiVec/Dual.h" // TODO
 #elif defined EIGEN_VECTORIZE_NEON
-//#include "duals/arch/NEON/Dual.h" // TODO
+  #include "duals/arch/NEON/Dual.h"
  #include "duals/arch/NEON/ComplexDual.h"
 #elif defined EIGEN_VECTORIZE_ZVECTOR
 //  #include "duals/arch/ZVector/Dual.h" // TODO
 #endif
--- a/src/include/cppduals/tests/CMakeLists.txt
+++ b/src/include/cppduals/tests/CMakeLists.txt
@ -10,47 +10,53 @@
 # Public License v. 2.0. If a copy of the MPL was not distributed
 # with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
-# gtest_discover_tests requires 3.10
+cmake_minimum_required (VERSION 3.14)
 cmake_minimum_required (VERSION 3.10)
 # Configure google-test as a downloadable library.
 include (GoogleTest)
 if (WIN32)
  add_definitions (-D_USE_MATH_DEFINES)
 endif ()
 include_directories ("${CMAKE_SOURCE_DIR}")
 include_directories ("${DEPS_ROOT}/include")
 include_directories ("${EIGEN3_INCLUDE_DIRS}")
 #include_directories ("${MPFR_INCLUDES}")
 include_directories ("${EXPOKIT_INCLUDE_DIR}")
 set (IOFORMAT "IOFormat(FullPrecision, DontAlignCols, \", \", \"\\$<SEMICOLON>\\n\", \"\", \"\", \"[\", \"]\")")
 add_definitions (-DEIGEN_DEFAULT_IO_FORMAT=${IOFORMAT})
 #add_definitions (-DEIGEN_DEFAULT_IO_FORMAT=EIGEN_IO_FORMAT_OCTAVE)
 #
 # Correctness & Coverage
 #
-if (NOT MSVC)
+message ("OSX_ARCHITECTURES: ${CMAKE_SYSTEM_PROCESSOR}")
-  set (OPT_FLAGS "-O2")
+set (OSX_ARCHITECTURES "arm64;x86_64")
-  set (OPT_FLAGS "-O3;-msse3")
+if (CMAKE_SYSTEM_PROCESSOR MATCHES "(x86)|(X86)|(amd64)|(AMD64)")
-  set (OPT_FLAGS "-O3;-mavx2;-mfma")
+  set (X86 TRUE)
  set (OPT_FLAGS "-O3;-march=native")
 else ()
-  set (OPT_FLAGS "/arch:IA32")
+  set (X86 FALSE)
-  set (OPT_FLAGS "/arch:SSE")
+endif ()
-  set (OPT_FLAGS "/arch:SSE2")
+
-  set (OPT_FLAGS "/arch:AVX2")
+# Set optimization flags only if coverage is not enabled
 if (NOT CPPDUALS_CODE_COVERAGE)
  if (X86)
    if (NOT MSVC)
      set (OPT_FLAGS "-O3;-march=native")
    else ()
      set (OPT_FLAGS "/arch:AVX2")
    endif ()
  else ()
    set (OPT_FLAGS "-O3")
  endif ()
 else()
  if (X86)
    if (NOT MSVC)
      set (OPT_FLAGS "-mavx2")
    else ()
      set (OPT_FLAGS "/arch:AVX2")
    endif ()
  else ()
    set (OPT_FLAGS "")
  endif ()
 endif ()
 #set (OPT_FLAGS "${OPT_FLAGS};-fsanitize=address;-fno-omit-frame-pointer")
 set (OPT_FLAGS "${OPT_FLAGS};-DCPPDUALS_VECTORIZE_CDUAL")
 #set (OPT_FLAGS "${OPT_FLAGS};-DCPPDUALS_DONT_VECTORIZE")
 #set (OPT_FLAGS "${OPT_FLAGS};-DEIGEN_DONT_VECTORIZE")
 set (ALL_TESTS
-  test_dual test_funcs test_eigen test_packets
+  test_dual test_cdual test_funcs test_eigen test_packets
  test_vectorize test_solve test_expm test_1 test_fmt
  example
  )
@ -70,33 +76,30 @@ foreach (TEST_ ${ALL_TESTS})
    add_executable (${TEST} ${TEST_}.cpp)
    set (PHASE_FLAGS ${OPT_FLAGS} -DPHASE_${PHASE})
-    string (REPLACE ";" ", " L2 "${CMAKE_CXX_FLAGS}")
+    string (REPLACE ";" ", " L2 "${OPT_FLAGS};${CMAKE_CXX_FLAGS}")
    target_compile_options (${TEST} PUBLIC ${PHASE_FLAGS})
    target_compile_definitions (${TEST} PRIVATE "OPT_FLAGS=${L2}")
-    target_link_libraries (${TEST} gtest_main cppduals_coverage_config)
+    
-    #target_link_libraries (${TEST} -lasan)
+    # Link with coverage config first to ensure flags are used
-    add_dependencies (${TEST} eigenX expokitX)
+    target_link_libraries (${TEST} 
      cppduals_coverage_config  # Link coverage first
      gtest_main 
      eigen 
      expokit
    )
    gtest_discover_tests (${TEST} TEST_LIST ${TEST}_targets)
    set_tests_properties (${${TEST}_targets} PROPERTIES TIMEOUT 10)
    # -ftest-coverage
  endforeach (PHASE)
 endforeach (TEST_)
 # special for fmt
 target_compile_features (test_fmt_1 PUBLIC cxx_std_14)
 target_link_libraries (test_fmt_1 fmt::fmt)
 if (CPPDUALS_BENCHMARK)
  #
  # Benchmarks
  #
  message ("searching: ${DEPS_ROOT} for google benchmark libs")
  find_library (BENCHMARK_LIBRARY benchmark PATHS ${DEPS_ROOT}/lib)
  find_library (BENCHMARKM_LIBRARY benchmark_main PATHS ${DEPS_ROOT}/lib)
  #find_library (PTHREAD_LIBRARY pthread)
  message ("BENCHMARK_LIBRARY: ${BENCHMARK_LIBRARY}")
  include_directories ("${BENCHMARK_INC_DIR}")
  if (Boost_FOUND AND NO)
    add_definitions (-DHAVE_BOOST=1)
    include_directories ("${Boost_INCLUDE_DIRS}")
@ -110,11 +113,14 @@ if (CPPDUALS_BENCHMARK)
      set (BLA_STATIC OFF)
    endif (NOT BLA_STATIC)
    set (BLA_VENDOR OpenBLAS)
-  endif (NOT APPLE AND NOT BLA_VENDOR)
+  endif ()
  find_package (BLAS REQUIRED)
  #find_package (LAPACK REQUIRED)
  add_definitions (-DHAVE_BLAS)
  #add_definitions (-DEIGEN_USE_BLAS)
  if (APPLE)
    add_definitions (-DACCELERATE_NEW_LAPACK -DACCELERATE_LAPACK_ILP64)
  endif ()
  # find lapacke.h cblas.h
  set (CBLAS_HINTS ${BLAS_DIR} ${LAPACK_DIR} /usr /usr/local /opt /opt/local)
@ -124,18 +130,21 @@ if (CPPDUALS_BENCHMARK)
    /opt
    /opt/local
    /usr/local/opt
-    /System/Library/Frameworks)
+    /System/Library/Frameworks
    ${BLAS_LIBRARIES}
    )
  # Finds the include directories for lapacke.h
  find_path (LAPACKE_INCLUDE_DIRS
-    NAMES lapacke.h
+    REQUIRED
    NAMES lapacke.h clapack.h
    HINTS ${CBLAS_HINTS}
    PATH_SUFFIXES
    include inc include/x86_64 include/x64
-    openblas/include
+    openblas/include openblas
-#    Accelerate.framework/Versions/Current/Frameworks/vecLib.framework/Versions/Current/Headers
+    Frameworks/vecLib.framework/Headers
    PATHS ${CBLAS_PATHS}
-    DOC "LAPACK(E) include header lapacke.h")
+    DOC "LAPACK(E) include header lapacke.h/clapack.h")
  mark_as_advanced (LAPACKE_INCLUDE_DIRS)
  if (LAPACKE_INCLUDE_DIRS)
    include_directories (${LAPACKE_INCLUDE_DIRS})
@ -145,12 +154,13 @@ if (CPPDUALS_BENCHMARK)
  # Finds the include directories for cblas*.h
  find_path (CBLAS_INCLUDE_DIRS
    REQUIRED
    NAMES cblas.h cblas_openblas.h cblas-openblas.h
    HINTS ${CBLAS_HINTS}
    PATH_SUFFIXES
    include inc include/x86_64 include/x64
-    openblas/include
+    openblas/include openblas
-#    Accelerate.framework/Versions/Current/Frameworks/vecLib.framework/Versions/Current/Headers
+    Frameworks/vecLib.framework/Headers
    PATHS ${CBLAS_PATHS}
    DOC "BLAS include header cblas.h")
  mark_as_advanced (CBLAS_INCLUDE_DIRS)
@ -161,80 +171,168 @@ if (CPPDUALS_BENCHMARK)
      break()
    endif (EXISTS "${CBLAS_INCLUDE_DIRS}/${cblas}")
  endforeach (cblas)
  message ("Found BLAS    : ${BLAS_LIBRARIES}")
  message ("Found cBLAS   : ${CBLAS_INCLUDE_DIRS}")
  message ("Found lapacke : ${LAPACKE_INCLUDE_DIRS}")
-  set (OPT_FLAGS "")
+  set (BMK_FLAGS "")
  if (NOT MSVC)
-    #set (OPT_FLAGS "-O3;-mavx")
+    #set (BMK_FLAGS "-O3;-mavx")
-    #set (OPT_FLAGS "-O3;-march=native;-fopenmp")
+    #set (BMK_FLAGS "-O3;-march=native;-fopenmp")
-    set (OPT_FLAGS "-O3;-msse3;-fopenmp")
+    if (X86)
-    set (OPT_FLAGS "-O3")
+      set (BMK_FLAGS "-msse3;-fopenmp")
-    set (OPT_FLAGS "-O3;-msse3")
+      set (BMK_FLAGS "-msse3")
-    set (OPT_FLAGS "-O3;-march=native;-funroll-loops")
+      set (BMK_FLAGS "-march=native;-funroll-loops")
-    set (OPT_FLAGS "-O3;-msse3;-mavx2;-mfma")
+      set (BMK_FLAGS "-msse3;-mavx2;-mfma")
    set (OPT_FLAGS "-O3;-march=native")
    #set (OPT_FLAGS "${OPT_FLAGS};-save-temps;-fverbose-asm")
    else ()
-    set (OPT_FLAGS "/arch:IA32")
+      set (BMK_FLAGS "-march=native")
-    set (OPT_FLAGS "/arch:SSE")
+    endif ()
-    set (OPT_FLAGS "/arch:SSE2")
+    #set (BMK_FLAGS "-O3;${BMK_FLAGS}")
-    set (OPT_FLAGS "/arch:AVX2")
+    #set (BMK_FLAGS "${BMK_FLAGS};-save-temps;-fverbose-asm")
  else ()
    set (BMK_FLAGS "/arch:IA32")
    set (BMK_FLAGS "/arch:SSE")
    set (BMK_FLAGS "/arch:SSE2")
    set (BMK_FLAGS "/arch:AVX2")
  endif ()
-  #set (OPT_FLAGS "${OPT_FLAGS};-DEIGEN_DONT_VECTORIZE")
+  #set (BMK_FLAGS "${BMK_FLAGS};-DEIGEN_DONT_VECTORIZE")
-  #set (OPT_FLAGS "${OPT_FLAGS};-DCPPDUALS_DONT_VECTORIZE")
+  #set (BMK_FLAGS "${BMK_FLAGS};-DCPPDUALS_DONT_VECTORIZE")
-  #set (OPT_FLAGS "${OPT_FLAGS};-DCPPDUALS_DONT_VECTORIZE_CDUAL")
+  #set (BMK_FLAGS "${BMK_FLAGS};-DCPPDUALS_DONT_VECTORIZE_CDUAL")
-  foreach (BENCH  bench_dual bench_eigen bench_gemm bench_example bench_fmt)
+  foreach (VECTORIZE YES NO)
-    add_executable (${BENCH} ${BENCH}.cpp)
+    foreach (BENCH  bench_dual bench_eigen bench_exp bench_gemm bench_example bench_fmt)
-    target_compile_options (${BENCH} PUBLIC ${OPT_FLAGS})
+      if (NOT VECTORIZE)
        set (BENCHE ${BENCH})
      else ()
        set (BENCHE ${BENCH}_novec)
      endif ()
      add_executable (${BENCHE} ${BENCH}.cpp)
      target_compile_options (${BENCHE} PUBLIC ${BMK_FLAGS})
      if (NOT VECTORIZE)
        target_compile_options (${BENCHE} PUBLIC "-DEIGEN_DONT_VECTORIZE")
      endif()
      #set_target_properties (${BENCH} PROPERTIES LINK_FLAGS -fopenmp)
-    #target_link_options (${BENCH} PUBLIC ${OPT_FLAGS})
+      #target_link_options (${BENCH} PUBLIC ${BMK_FLAGS})
-    string (REPLACE ";" ", " L2 "${OPT_FLAGS} ${CMAKE_CXX_FLAGS}")
+      string (REPLACE ";" ", " L2 "${BMK_FLAGS} ${CMAKE_CXX_FLAGS}")
-    target_compile_definitions (${BENCH} PRIVATE "OPT_FLAGS=${L2}")
+      target_compile_definitions (${BENCHE} PRIVATE "BMK_FLAGS=${L2}")
-    add_dependencies (${BENCH} benchmarkX eigenX expokitX)
+      target_link_libraries (${BENCHE} benchmark::benchmark ${BLAS_LIBRARIES} eigen expokit)
-    target_link_libraries (${BENCH} ${BENCHMARK_LIBRARY} -lpthread ${BLAS_LIBRARIES})
+    endforeach ()
-  endforeach (BENCH)
+  endforeach ()
  target_link_libraries (bench_fmt fmt::fmt)
  target_link_libraries (bench_fmt_novec fmt::fmt)
 endif (CPPDUALS_BENCHMARK)
 add_executable (sandbox sandbox.cpp)
-add_dependencies (sandbox eigenX expokitX ) # mpfrX mprealX
+#target_compile_options (sandbox PUBLIC ${BMK_FLAGS})
 #target_compile_options (sandbox PUBLIC ${OPT_FLAGS})
 target_compile_options (sandbox PUBLIC -DCPPDUALS_VECTORIZE_CDUAL)
 if (MSVC)
  if (X86)
    target_compile_options (sandbox PUBLIC /arch:AVX2)
  endif ()
 else ()
  if (X86)
    target_compile_options (sandbox PUBLIC -O1 -msse3 -mavx2 -mfma)
  else ()
    target_compile_options (sandbox PUBLIC -O1 )  # -mfpu=neon
  endif ()
 endif ()
 set_target_properties (sandbox PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR})
-target_link_libraries (sandbox PUBLIC cppduals_coverage_config)
+target_link_libraries (sandbox PUBLIC cppduals_coverage_config eigen expokit)
 #
 # Generate coverage reports
 #
-if (CODE_COVERAGE)
+if (CPPDUALS_CODE_COVERAGE)
-  add_custom_target (cov
+  # Find required tools
  get_filename_component(CMAKE_CXX_COMPILER_DIR "${CMAKE_CXX_COMPILER}" DIRECTORY)
  set(GENHTML_PATH ${CMAKE_CXX_COMPILER_DIR}/genhtml)
  if(NOT EXISTS ${GENHTML_PATH})
    find_program(GENHTML_PATH genhtml REQUIRED)
  endif()
  if (CMAKE_CXX_COMPILER_ID MATCHES "Clang")
    # Clang
    # often in the same directory as clang - search there too
    get_filename_component(CMAKE_CXX_COMPILER_DIR "${CMAKE_CXX_COMPILER}" DIRECTORY)
    set(LLVM_COV_PATH ${CMAKE_CXX_COMPILER_DIR}/llvm-cov)
    set(LLVM_PROFDATA_PATH ${CMAKE_CXX_COMPILER_DIR}/llvm-profdata)
    if(NOT EXISTS ${LLVM_COV_PATH})
      find_program(LLVM_COV_PATH llvm-cov REQUIRED)
    endif()
    if(NOT EXISTS ${LLVM_PROFDATA_PATH})
      find_program(LLVM_PROFDATA_PATH llvm-profdata REQUIRED)
    endif()
    add_custom_target(cov
      DEPENDS ${ALL_TEST_BINS}
-    COMMAND ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR} --target test
+      COMMAND ${CMAKE_COMMAND} -E remove_directory coverage/profraw
-    COMMAND $<TARGET_FILE:sandbox>
+      COMMAND ${CMAKE_COMMAND} -E make_directory coverage/profraw
-    COMMAND lcov --capture --directory . --output-file coverage.info
+      COMMAND ${CMAKE_COMMAND} --build . --target test "LLVM_PROFILE_FILE=coverage/profraw/%p.profraw"
-    COMMAND lcov --remove coverage.info '/usr/*' --output-file coverage.info
+      COMMAND ${LLVM_PROFDATA_PATH} merge -sparse coverage/profraw/*.profraw -o coverage/coverage.profdata
-    COMMAND lcov --remove coverage.info '*/thirdparty/*' --output-file coverage.info
+      COMMAND ${LLVM_COV_PATH} show
-    COMMAND lcov --remove coverage.info '*/googletest/*' --output-file coverage.info
+        -instr-profile=coverage/coverage.profdata
-    COMMAND lcov --list coverage.info
+        -format=html -output-dir=coverage/html
        -show-line-counts-or-regions
        -show-instantiation-summary
        ${ALL_TEST_BINS}
      COMMAND ${LLVM_COV_PATH} report
        -instr-profile=coverage/coverage.profdata
        ${ALL_TEST_BINS}
      WORKING_DIRECTORY ${CMAKE_BINARY_DIR}
    )
-  add_custom_target (cov-html
+  elseif(CMAKE_CXX_COMPILER_ID MATCHES "GNU")
    # GCC
    set(LCOV_PATH ${CMAKE_CXX_COMPILER_DIR}/lcov)
    if(NOT EXISTS ${LCOV_PATH})
      find_program(LCOV_PATH lcov REQUIRED)
    endif()
    # Extract GCC version number and use matching gcov version
    execute_process(
      COMMAND ${CMAKE_CXX_COMPILER} -dumpversion
      OUTPUT_VARIABLE GCC_VERSION
      OUTPUT_STRIP_TRAILING_WHITESPACE
    )
    set(GCOV_PATH ${CMAKE_CXX_COMPILER_DIR}/gcov-${GCC_VERSION})
    if(NOT EXISTS ${GCOV_PATH})
      find_program(GCOV_PATH gcov-${GCC_VERSION} REQUIRED)
    endif()
    add_custom_target(cov
      DEPENDS ${ALL_TEST_BINS}
      COMMAND ${CMAKE_COMMAND} -E make_directory coverage
      COMMAND ${LCOV_PATH} --directory . --zerocounters
      # Add initial capture before running tests
      COMMAND ${LCOV_PATH} --directory . --capture --initial --gcov-tool ${GCOV_PATH} --output-file coverage/coverage.base --ignore-errors inconsistent
      # Run tests and capture coverage data
      COMMAND ctest --output-on-failure
      COMMAND ${LCOV_PATH} --directory . --capture --gcov-tool ${GCOV_PATH} --output-file coverage/coverage.info --ignore-errors inconsistent
      # Combine baseline and test coverage data
      COMMAND ${LCOV_PATH} --add-tracefile coverage/coverage.base --add-tracefile coverage/coverage.info --gcov-tool ${GCOV_PATH} --output-file coverage/coverage.total
      # Only look at the coverage of the tests and duals library
      COMMAND ${LCOV_PATH} --extract coverage/coverage.total '*/tests/*' '*/duals/*' --gcov-tool ${GCOV_PATH} --output-file coverage/coverage.info
      # Remove unwanted paths
      #COMMAND ${LCOV_PATH} --remove coverage/coverage.total '/usr/*' --gcov-tool ${GCOV_PATH} --output-file coverage/coverage.info
      #COMMAND ${LCOV_PATH} --remove coverage/coverage.info '*/thirdparty/*' --gcov-tool ${GCOV_PATH} --output-file coverage/coverage.info
      #COMMAND ${LCOV_PATH} --remove coverage/coverage.info '*/googletest/*' --gcov-tool ${GCOV_PATH} --output-file coverage/coverage.info
      # Add --ignore-errors empty to prevent failure on empty coverage data
      COMMAND ${LCOV_PATH} --list coverage/coverage.info --ignore-errors empty --gcov-tool ${GCOV_PATH}
      WORKING_DIRECTORY ${CMAKE_BINARY_DIR}
    )
  else()
    message(FATAL_ERROR "No coverage tool found for ${CMAKE_CXX_COMPILER_ID}")
  endif()
  add_custom_target(cov-html
    DEPENDS cov
-    COMMAND rm -rf ../coverage
+    COMMAND ${CMAKE_COMMAND} -E remove_directory coverage/html
-    COMMAND genhtml --ignore-errors source coverage.info --legend --title "make cov"
+    COMMAND ${GENHTML_PATH} --ignore-errors source coverage/coverage.info --legend --title "cppduals coverage"
-    --output-directory=../coverage
+      --output-directory=coverage/html
-    COMMAND echo "output in coverage/index.html"
+    COMMAND ${CMAKE_COMMAND} -E echo "Coverage report generated at coverage/html/index.html"
    WORKING_DIRECTORY ${CMAKE_BINARY_DIR}
  )
 endif ()
--- a/src/include/cppduals/tests/bench_eigen.cpp
+++ b/src/include/cppduals/tests/bench_eigen.cpp
@ -11,13 +11,13 @@
 // (c)2019 Michael Tesch. tesch1@gmail.com
 //
 #include <duals/dual_eigen>
 #include <iostream>
 #include <fstream>
 #include <complex>
 #include <memory>
 #include "type_name.hpp"
 #include <duals/dual_eigen>
 #include <Eigen/Core>
 #include <Eigen/Dense>
@ -31,67 +31,6 @@ using namespace duals;
 template< class T > struct type_identity { typedef T type; };
 namespace Eigen {
 namespace internal {
 template<typename T> struct is_exp_known_type;
 template<typename T> struct is_exp_known_type<std::complex<T>> : is_exp_known_type<T> {};
 #if 0
 template <typename RealScalar> struct MatrixExponentialScalingOp;
 template <typename RealScalar>
 struct MatrixExponentialScalingOp<duals::dual<RealScalar>>
 {
  MatrixExponentialScalingOp(int squarings) : m_squarings(squarings) { }
  inline const duals::dual<RealScalar> operator() (const duals::dual<RealScalar> & x) const
  {
    using std::ldexp;
    return ldexp(x, -m_squarings);
  }
  typedef std::complex<duals::dual<RealScalar>> ComplexScalar;
  inline const ComplexScalar operator() (const ComplexScalar& x) const
  {
    using std::ldexp;
    return ComplexScalar(ldexp(x.real(), -m_squarings), ldexp(x.imag(), -m_squarings));
  }
  private:
    int m_squarings;
 };
 #endif
 }}
 #include <unsupported/Eigen/MatrixFunctions>
 namespace Eigen {
 namespace internal {
 template <typename MatrixType, typename T>
 struct matrix_exp_computeUV<MatrixType, duals::dual<T> >
 {
  typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
  template <typename ArgType>
  static void run(const ArgType& arg, MatrixType& U, MatrixType& V, int& squarings)
  {
    using std::frexp;
    using std::pow;
    const RealScalar l1norm = arg.cwiseAbs().colwise().sum().maxCoeff();
    squarings = 0;
    if (l1norm < 1.495585217958292e-002) {
      matrix_exp_pade3(arg, U, V);
    } else if (l1norm < 2.539398330063230e-001) {
      matrix_exp_pade5(arg, U, V);
    } else if (l1norm < 9.504178996162932e-001) {
      matrix_exp_pade7(arg, U, V);
    } else if (l1norm < 2.097847961257068e+000) {
      matrix_exp_pade9(arg, U, V);
    } else {
      const RealScalar maxnorm = 5.371920351148152;
      frexp(l1norm / maxnorm, &squarings);
      if (squarings < 0) squarings = 0;
      MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<RealScalar>(squarings));
      matrix_exp_pade13(A, U, V);
    }
  }
 };
 }}
 /* encode the type into an integer for benchmark output */
 template<typename Tp> struct type_num { /* should fail */ };
 template<> struct type_num<float>     { static constexpr int id = 1; };
@ -271,87 +210,6 @@ template <class Rt, class U> void B_MatVec(benchmark::State& state) {
  state.SetComplexityN(state.range(0));
 }
 template <class Rt> void B_Expm(benchmark::State& state)
 {
  int N = state.range(0);
  //Rt S(1);
  MatrixX<Rt> A = MatrixX<Rt>::Random(N, N);
  MatrixX<Rt> B = MatrixX<Rt>::Zero(N, N);
  //A = S * A / A.norm();
  for (auto _ : state) {
    B = A.exp();
    benchmark::ClobberMemory();
  }
  state.counters["type"] = type_num<Rt>::id;
  state.SetComplexityN(state.range(0));
 }
 template <class Rt> void B_ExpPadm(benchmark::State& state)
 {
  int N = state.range(0);
  //Rt S(1);
  MatrixX<Rt> A = MatrixX<Rt>::Random(N, N);
  MatrixX<Rt> B = MatrixX<Rt>::Zero(N, N);
  //A = S * A / A.norm();
  for (auto _ : state) {
    B = eexpokit::padm(A);
    benchmark::ClobberMemory();
  }
  state.counters["type"] = type_num<Rt>::id;
  state.SetComplexityN(state.range(0));
 }
 template <class Rt> void B_ExpExpv(benchmark::State& state)
 {
  int N = state.range(0);
  //Rt S(1);
  MatrixX<Rt> A = MatrixX<Rt>::Zero(N, N);
  MatrixX<Rt> b = MatrixX<Rt>::Ones(N, 1);
  MatrixX<Rt> c = MatrixX<Rt>::Zero(N, 1);
  //A = S * A / A.norm();
  // sparse random fill
  for (int i = 0; i < 4*N; i++)
    A((int)duals::randos::random(0.,N-1.),
      (int)duals::randos::random(0.,N-1.)) = duals::randos::random2<Rt>();
  for (auto _ : state) {
    auto ret = eexpokit::expv(1,A,b);
    if (ret.err > 1) {
      std::ofstream f("fail.m");
      f << "A=" << A.format(eexpokit::OctaveFmt) << "\n";
      break;
    }
    // c = ret.w
    benchmark::ClobberMemory();
  }
  state.counters["type"] = type_num<Rt>::id;
  state.SetComplexityN(state.range(0));
 }
 template <class Rt> void B_ExpChbv(benchmark::State& state)
 {
  int N = state.range(0);
  //Rt S(1);
  MatrixX<Rt> A = MatrixX<Rt>::Random(N, N);
  MatrixX<Rt> b = MatrixX<Rt>::Zero(N, 1);
  MatrixX<Rt> c = MatrixX<Rt>::Zero(N, 1);
  //A = S * A / A.norm();
  for (auto _ : state) {
    c = eexpokit::chbv(A,b);
    benchmark::ClobberMemory();
  }
  state.counters["type"] = type_num<Rt>::id;
  state.SetComplexityN(state.range(0));
 }
 #define MAKE_BM_SIMPLE(TYPE1,TYPE2,NF)                                  \
  BENCHMARK_TEMPLATE(B_VecVecAdd, TYPE1,TYPE2) V_RANGE(4,NF);           \
  BENCHMARK_TEMPLATE(B_VecVecSub, TYPE1,TYPE2) V_RANGE(4,NF);           \
@ -359,11 +217,7 @@ template <class Rt> void B_ExpChbv(benchmark::State& state)
  BENCHMARK_TEMPLATE(B_VecVecDiv, TYPE1,TYPE2) V_RANGE(4,NF);           \
  BENCHMARK_TEMPLATE(B_MatVec, TYPE1,TYPE2) V_RANGE(4,NF);              \
  BENCHMARK_TEMPLATE(B_MatMat, TYPE1,TYPE2) V_RANGE(1,NF);              \
-  BENCHMARK_TEMPLATE(B_MatDiv, TYPE1,TYPE2) V_RANGE(1,NF);              \
+  BENCHMARK_TEMPLATE(B_MatDiv, TYPE1,TYPE2) V_RANGE(1,NF)
  BENCHMARK_TEMPLATE(B_Expm, TYPE1) V_RANGE(1,NF);                      \
  BENCHMARK_TEMPLATE(B_ExpPadm, TYPE1) V_RANGE(1,NF);                   \
  BENCHMARK_TEMPLATE(B_ExpChbv, TYPE1) V_RANGE(1,NF);                   \
  BENCHMARK_TEMPLATE(B_ExpExpv, TYPE1) V_RANGE(1,NF)
 #define MAKE_BENCHMARKS(TYPE1,TYPE2,NF)                                 \
  MAKE_BM_SIMPLE(TYPE1,TYPE2,NF)
@ -389,7 +243,7 @@ MAKE_BM_SIMPLE(cduald, cduald,4);
 int main(int argc, char** argv)
 {
 #ifndef EIGEN_VECTORIZE
-  static_assert(false, "no vectorization?");
+  //static_assert(false, "no vectorization?");
 #endif
  std::cout << "OPT_FLAGS=" << QUOTE(OPT_FLAGS) << "\n";
  std::cout << "INSTRUCTIONSET=" << Eigen::SimdInstructionSetsInUse() << "\n";
--- a/src/include/cppduals/tests/bench_gemm.cpp
+++ b/src/include/cppduals/tests/bench_gemm.cpp
@ -6,12 +6,14 @@
 //
 // (c)2019 Michael Tesch. tesch1@gmail.com
 //
 #include <duals/dual>
 #if defined(__APPLE__) && defined(__clang__)
 #include <Accelerate/Accelerate.h>
 #else
-#ifdef EIGEN_LAPACKE
+#if defined(EIGEN_LAPACKE) || defined(__APPLE__)
 #include <Eigen/src/misc/lapacke.h>
 #else
 #include <lapacke.h>
@ -24,13 +26,13 @@ extern "C" {
 }
 #endif // defined(__APPLE__) && defined(__clang__)
 #include <duals/dual_eigen>
 #include <iostream>
 #include <fstream>
 #include <complex>
 #include <memory>
 #include "type_name.hpp"
 #include <duals/dual_eigen>
 #include <Eigen/Core>
 #include <Eigen/Dense>
@ -80,6 +82,59 @@ template <class T, class U> void B_MatMat(benchmark::State& state) {
  state.SetComplexityN(state.range(0));
 }
 // Helper templates to select correct BLAS routine
 template <typename U>
 struct blas_gemm;
 #define TRANSPOSE(X) ((X) ? CblasTrans : CblasNoTrans)
 template <>
 struct blas_gemm<float> {
  static void call(bool tA, bool tB, int m, int n, int k, float alpha, float *A,
                   int lda, float *B, int ldb, float beta, float *C, int ldc)
  {
    cblas_sgemm(CblasColMajor, TRANSPOSE(tA), TRANSPOSE(tB), m, n, k, alpha, A,
                lda, B, ldb, beta, C, ldc);
  }
 };
 template <>
 struct blas_gemm<double> {
  static void call(bool tA, bool tB, int m, int n, int k, double alpha,
                   double *A, int lda, double *B, int ldb, double beta,
                   double *C, int ldc)
  {
    cblas_dgemm(CblasColMajor, TRANSPOSE(tA), TRANSPOSE(tB), m, n, k, alpha, A,
                lda, B, ldb, beta, C, ldc);
  }
 };
 template <>
 struct blas_gemm<std::complex<float>> {
  static void call(bool tA, bool tB, int m, int n, int k,
                   std::complex<float> alpha, std::complex<float> *A, int lda,
                   std::complex<float> *B, int ldb, std::complex<float> beta,
                   std::complex<float> *C, int ldc)
  {
    cblas_cgemm(CblasColMajor, TRANSPOSE(tA), TRANSPOSE(tB), m, n, k, &alpha, A,
                lda, B, ldb, &beta, C, ldc);
  }
 };
 template <>
 struct blas_gemm<std::complex<double>> {
  static void call(bool tA, bool tB, int m, int n, int k,
                   std::complex<double> alpha, std::complex<double> *A, int lda,
                   std::complex<double> *B, int ldb, std::complex<double> beta,
                   std::complex<double> *C, int ldc)
  {
    cblas_zgemm(CblasColMajor, TRANSPOSE(tA), TRANSPOSE(tB), m, n, k, &alpha, A,
                lda, B, ldb, &beta, C, ldc);
  }
 };
 #undef TRANSPOSE
 template <class T, typename std::enable_if<!duals::is_dual<T>::value>::type* = nullptr>
 void matrix_multiplcation(T *A, int Awidth, int Aheight,
                          T *B, int Bwidth, int Bheight,
@ -99,31 +154,12 @@ void matrix_multiplcation(T *A, int Awidth, int Aheight,
  assert(A_width == B_height);
  int lda = tA ? m : k;
  int ldb = tB ? k : n;
-#define TRANSPOSE(X) ((X) ? CblasTrans : CblasNoTrans)
+
-  // http://www.netlib.org/lapack/explore-html/d7/d2b/dgemm_8f.html
+  // Call the appropriate BLAS routine based on type T
-  if (!is_complex<T>::value) {
+  blas_gemm<T>::call(tA, tB, m, n, k, T(1), A, lda, B, ldb, beta, AB, n);
    if (sizeof(T) == sizeof(float))
      cblas_sgemm(CblasColMajor, TRANSPOSE(tA), TRANSPOSE(tB),
                  m, n, k, 1.0, (float *)A, lda, (float *)B, ldb,
                  std::real(beta), (float *)AB, n);
    else
      cblas_dgemm(CblasColMajor, TRANSPOSE(tA), TRANSPOSE(tB),
                  m, n, k, 1.0, (double *)A, lda, (double *)B, ldb,
                  std::real(beta), (double *)AB, n);
  }
  else {
    std::complex<float> alphaf(1,0);
    std::complex<double> alpha(1,0);
    if (Eigen::NumTraits<T>::digits10() < 10)
      cblas_cgemm(CblasColMajor, TRANSPOSE(tA), TRANSPOSE(tB),
                  m, n, k, &alphaf, A, lda, B, ldb, &beta, AB, n);
    else
      cblas_zgemm(CblasColMajor, TRANSPOSE(tA), TRANSPOSE(tB),
                  m, n, k, &alpha, A, lda, B, ldb, &beta, AB, n);
  }
 #undef TRANSPOSE
 }
 template <class T, typename std::enable_if<duals::is_dual<T>::value>::type* = nullptr>
 void matrix_multiplcation(T *A, int Awidth, int Aheight,
                          T *B, int Bwidth, int Bheight,
@ -232,10 +268,10 @@ MAKE_BM_SIMPLE(cduald, cduald,4);
 int main(int argc, char** argv)
 {
 #ifndef EIGEN_VECTORIZE
-  static_assert(false, "no vectorization?");
+  //static_assert(false, "no vectorization?");
 #endif
 #ifndef NDEBUG
-  static_assert(false, "NDEBUG to benchmark?");
+  //static_assert(false, "NDEBUG to benchmark?");
 #endif
  std::cout << "OPT_FLAGS=" << QUOTE(OPT_FLAGS) << "\n";
  std::cout << "INSTRUCTIONSET=" << Eigen::SimdInstructionSetsInUse() << "\n";
--- a/src/include/cppduals/tests/sandbox.cpp
+++ b/src/include/cppduals/tests/sandbox.cpp
@ -12,12 +12,12 @@
 * (c)2019 Michael Tesch. tesch1@gmail.com
 */
 #include <duals/dual_eigen>
 #include <math.h>
 #include <iostream>
 #include <iomanip>
 #include "type_name.hpp"
 #include <duals/dual_eigen>
 #include <Eigen/Dense>
 #include <Eigen/Sparse>
 #include <unsupported/Eigen/MatrixFunctions>
@ -60,13 +60,45 @@ template <class T>
 struct common_type<Rando<T>,T> { typedef Rando<T> type; };
 }
-#if 0
+#if EIGEN_ARCH_ARM64
 int main(int argc, char * argv[])
 {
  emtx<cduald,50> A,B,C;
  //emtx<complexd,50> A,B,C;
  C = A * B;
  //
  using namespace Eigen::internal;
  float32x4_t a, b, v3;
  float32x2_t f2;
  float32x4x2_t yy;
  float ff[] = {1., 2., 3., 4.};
  duals::dual<double> x;
  const uint64x2_t maskq = {0, (uint64_t)-1};
  const uint64x2_t imaskq = {(uint64_t)-1, 0};
  std::cout << "a:" << std::hex << maskq[0] << "\n";
  std::cout << "b:" << maskq[1] << "\n";
  std::cout << "v3=" << v3[0] << ", " << v3[1] << ", " << v3[2] << ", " << v3[3] << "\n";
  std::cout << "f2=" << f2[0] << ", " << f2[1] << "\n";
  std::cout << ff[0] << ", " << ff[1] << ", " << ff[2] << ", " << ff[3] << "\n";
 #ifdef __ARM_NEON
  std::cout << "__ARM_NEON:" << __ARM_NEON << "\n";
 #endif
 #ifdef __ARM_NEON_FP
  std::cout << "__ARM_NEON_FP:" << __ARM_NEON_FP << "\n";
 #endif
 #ifdef __ARM_ARCH
  std::cout << "__ARM_ARCH:" << __ARM_ARCH << "\n";
 #endif
 #ifdef __ARM_BIG_ENDIAN
  std::cout << "__ARM_BIG_ENDIAN:" << __ARM_BIG_ENDIAN << "\n";
 #endif
  return 0;
 }
 #elif 0
--- a/src/include/cppduals/tests/test_1.cpp
+++ b/src/include/cppduals/tests/test_1.cpp
@ -81,9 +81,7 @@ expm4(const Eigen::EigenBase<DerivedA> & A_,
  R.setIdentity();
  R += B;
  ReturnT S = B;
  int ni = 0;
  for (int ii = 2; ii < maxt; ii++) {
    ni++;
    S = Real(1.0/ii) * S * B;
    R += S;
    auto Sn = S.norm();
--- a/src/include/cppduals/tests/test_dual.cpp
+++ b/src/include/cppduals/tests/test_dual.cpp
@ -18,6 +18,7 @@
 #include <duals/dual>
 #include <complex>
 #include <iomanip>
 #include "gtest/gtest.h"
 using duals::dualf;
@ -58,7 +59,7 @@ TEST(template_, dual_traits) {
  // depth
  EXPECT_EQ(dual_traits<float>::depth, 0);
  EXPECT_EQ(dual_traits<complexf>::depth, 0);
-  EXPECT_EQ(dual_traits<cdualf>::depth, 0);
+  EXPECT_EQ(dual_traits<cdualf>::depth, 1);
  EXPECT_EQ(dual_traits<dualf>::depth, 1);
  EXPECT_EQ(dual_traits<hyperdualf>::depth, 2);
 }
@ -186,6 +187,8 @@ TEST(template_, common_type) {
  _EXPECT_TRUE(std::is_same<decltype(cd*cd), cdualf>);
  _EXPECT_TRUE(std::is_same<decltype(cd*d), cdualf>);
  _EXPECT_TRUE(std::is_same<decltype(d*cd), cdualf>);
  _EXPECT_TRUE(std::is_same<decltype(cd*a), cdualf>);
  _EXPECT_TRUE(std::is_same<decltype(a*cd), cdualf>);
  _EXPECT_TRUE(std::is_same<decltype(cd*1), cdualf>);
  _EXPECT_TRUE(std::is_same<decltype(1*cd), cdualf>);
@ -328,6 +331,7 @@ TEST(members, rpart) {
  EXPECT_EQ(z.rpart(), 4);
  EXPECT_EQ(z.dpart(), 3);
 }
 TEST(members, dpart) {
  EXPECT_EQ(dpart(3), 0);
  dualf z(2,3);
@ -555,6 +559,15 @@ TEST(comparison, ge) {
  EXPECT_FALSE(1 >= a);
 }
 #if 0
 TEST(simple_ops, add) {
  // https://gitlab.com/tesch1/cppduals/-/issues/11
  duals::dual<std::complex<double>> x;
  x = x+std::complex<double>(1., 2.);
  duals::dual<std::complex<double>> y = sqrt(x+std::complex<double>(1., 2.));
 }
 #endif
 TEST(compound_assign, same_type) {
  // OP=
  dualf x = 2 + 4_e;
@ -977,25 +990,6 @@ TEST(non_class, random2) {
  EXPECT_NE(c2.dpart(), 0);
  EXPECT_NE(c1.rpart(), c2.rpart());
  EXPECT_NE(c1.dpart(), c2.dpart());
  // cdualf
  cdualf d1 = duals::randos::random2<cdualf>();
  cdualf d2 = duals::randos::random2<cdualf>();
  EXPECT_NE(d1.real().rpart(), 0);
  EXPECT_NE(d1.real().dpart(), 0);
  EXPECT_NE(d1.imag().rpart(), 0);
  EXPECT_NE(d1.imag().dpart(), 0);
  EXPECT_NE(d2.real().rpart(), 0);
  EXPECT_NE(d2.real().dpart(), 0);
  EXPECT_NE(d2.imag().rpart(), 0);
  EXPECT_NE(d2.imag().dpart(), 0);
  EXPECT_NE(d1.real().rpart(), d2.real().rpart());
  EXPECT_NE(d1.real().dpart(), d2.real().dpart());
  EXPECT_NE(d1.imag().rpart(), d2.imag().rpart());
  EXPECT_NE(d1.imag().dpart(), d2.imag().dpart());
 }
 TEST(smoke, funcs) {
@ -1097,6 +1091,17 @@ TEST(complex, mixing) {
  // complex<dual> * real -> complex<dual>
  // complex<real> * dual -> complex<dual> ?
  // pow_dual_complex
  B = pow(1_ef, C);
  // pow_complex_dual
  B = pow(C, 1_ef);
  // pow_dual_dual
  B = pow(1_ef, 2_ef);
  // pow_dual_scalar
  B = pow(1_ef, 2);
  // pow_scalar_dual
  B = pow(2, 1_ef);
 }
 int main(int argc, char **argv)
--- a/src/include/cppduals/tests/test_eigen.cpp
+++ b/src/include/cppduals/tests/test_eigen.cpp
@ -16,11 +16,12 @@
 * (c)2019 Michael Tesch. tesch1@gmail.com
 */
 #include "type_name.hpp"
 #include <duals/dual_eigen>
 #include "type_name.hpp"
 #include <Eigen/Dense>
 #include <Eigen/Sparse>
 #include <Eigen/StdVector>
 #include <unsupported/Eigen/MatrixFunctions>
 //#include <unsupported/Eigen/AutoDiff>
 #include "eexpokit/padm.hpp"
@ -39,10 +40,13 @@ using duals::is_complex;
 using duals::dual_traits;
 using namespace duals::literals;
-typedef std::complex<double> complexd;
+typedef long double ldouble;
 typedef std::complex<float> complexf;
-typedef std::complex<duald> cduald;
+typedef std::complex<double> complexd;
 typedef std::complex<long double> complexld;
 typedef std::complex<dualf> cdualf;
 typedef std::complex<duald> cduald;
 typedef std::complex<dualld> cdualld;
 template <class eT, int N=Eigen::Dynamic, int K = N> using emtx = Eigen::Matrix<eT, N, K>;
 template <class eT> using smtx = Eigen::SparseMatrix<eT>;
@ -60,49 +64,6 @@ template <int N=2, int K = N> using ecdf = Eigen::Matrix<cdualf, N, K> ;
  EXPECT_NEAR(rpart(A), rpart(B),tol);             \
  EXPECT_NEAR(dpart(A), dpart(B),tol)
 /// Simple taylor series, truncated when |S| is "small enough"
 template <class DerivedA, typename ReturnT = typename DerivedA::PlainObject>
 ReturnT
 expm4(const Eigen::EigenBase<DerivedA> & A_,
      typename DerivedA::RealScalar mn = std::numeric_limits<typename DerivedA::RealScalar>::epsilon() * 1000)
 {
  //std::cerr << "do    PO:" << type_name<typename DerivedA::PlainObject>() << "\n";
  typedef typename DerivedA::RealScalar Real;
  using std::isfinite;
  const DerivedA & A = A_.derived();
  int maxt = std::numeric_limits<Real>::digits;
  int s = log2(rpart(A.derived().norm())) + 1;
  s = std::max(0, s);
  auto B = A * pow(Real(2), -s);
  ReturnT R(A.rows(), A.cols());
  R.setIdentity();
  R += B;
  ReturnT S = B;
  for (int ii = 2; ii < maxt; ii++) {
    S = S * B * Real(1.0/ii);
    R += S;
    auto Sn = S.norm();
    if (!isfinite(Sn)) {
      std::cout << "expm() non-finite norm:" << Sn << " at " << ii << "\n";
      std::cout << " |R| = " << R.norm() << " s=" << s << "\n"
                << " |A| = " << rpart(A.real().norm()) << "\n"
                << " |A/2^s|=" << rpart(A.real().norm()/pow(2,s)) << "\n";
      break;
    }
    // converged yet?
    if (Sn < mn)
      break;
    if (ii == maxt - 1) {
      std::cout << "expm() didn't converge in " << maxt << " |S| = " << Sn << "\n";
      throw std::invalid_argument("no converge");
    }
  }
  for (; s; s--)
    R = R * R;
  return R;
 }
 TEST(Eigen, NumTraits) {
  //::testing::StaticAssertTypeEq<Eigen::NumTraits<duals::dual<float>>::Real, float>();
@ -448,11 +409,12 @@ TEST(measure, norm) {
  b = 3;
  Rt d(1);
  MatrixD x;
-  x << d;
+  x.array() = d;
  MatrixD a = (MatrixD() << 1,2,3, 4,5+5_ef,6, 7,8,9).finished();
  //typename MatrixD::Index index;
  EXPECT_EQ(a.sum(), 45 + 5_ef);
  EXPECT_EQ(x.sum(), 9);
  EXPECT_NEAR(rpart(a.norm()), 16.8819430161341337282, 1e-5);
  EXPECT_NEAR(rpart(a.mean()), 5, 1e-5);
  EXPECT_NEAR(dpart(a.mean()), 0.555555555555555, 1e-5);
@ -510,31 +472,49 @@ TEST(dpart, matrix) {
  EXPECT_EQ((dpart(AA) - CC).norm(),0);
 }
-TEST(func, expm) {
+TEST(eigen, exp_typechecks) {
-  typedef float T;
+  typedef emtx<dualf,3> Mat;
-  typedef dual<T> dualt;
+  typedef emtx<duald,3> Matd;
-  typedef std::complex<dual<T>> cdualt;
+  EXPECT_FALSE(Eigen::internal::is_exp_known_type<typename Mat::Scalar>::value);
-#define NN 3
+  EXPECT_FALSE(Eigen::internal::is_exp_known_type<typename Matd::Scalar>::value);
 #define N2 6
  emtx<dualt, NN> a,b;
  a = emtx<dualt, NN>::Random();
  a.array() += 1.1 + 2.2_ef;
  a.setZero();
  a = eexpokit::padm(a);
  EXPECT_LT((a - emtx<dualt, NN>::Identity()).norm(), 1e-6) << "a=" << a << "\n";
  a *= 1+2_e;
  EXPECT_LT((a - emtx<dualt, NN>::Identity()).norm(), 1e-6) << "a=" << a << "\n";
-  emtx<cdualt, NN> c;
+  typedef emtx<cdualf,3> Matc;
-  //b = a + 1_e * emtx<cdualf, 3>::Random();
+  typedef emtx<cduald,3> Matcd;
-  c.setZero();
+  EXPECT_FALSE(Eigen::internal::is_exp_known_type<typename Matc::Scalar>::value);
-  c = c.exp();
+  EXPECT_FALSE(Eigen::internal::is_exp_known_type<typename Matcd::Scalar>::value);
  //c = expm4(c);
  EXPECT_LT((c - emtx<cdualf, NN>::Identity()).norm(), 1e-6) << "b=" << b << "\n";
  #undef NN
  #undef N2
 }
 const bool _exp = true;
 const bool _padm = false;
 #define TEST_EXP(SCALAR_T, SIDE, EXP_OR_PADM) \
  TEST(func, exp_##SCALAR_T##_##SIDE##EXP_OR_PADM) { \
    emtx<SCALAR_T, SIDE> a,b; \
    a = emtx<SCALAR_T, SIDE>::Random(); \
    a.setZero(); \
    if (EXP_OR_PADM == _exp) a = a.exp(); \
    else a = eexpokit::padm(a); \
    EXPECT_LT((a - emtx<SCALAR_T, SIDE>::Identity()).norm(), 1e-6) << "a=" << a << "\n"; \
 }
 #define TEST_EXP_SIZES(SCALAR_T, EXP_OR_PADM) \
  TEST_EXP(SCALAR_T, 3, EXP_OR_PADM) \
  TEST_EXP(SCALAR_T, 4, EXP_OR_PADM) \
  TEST_EXP(SCALAR_T, 17, EXP_OR_PADM)
 #define TEST_EXP_SIZES_EOP(SCALAR_T) \
  TEST_EXP_SIZES(SCALAR_T, _padm) \
  TEST_EXP_SIZES(SCALAR_T, _exp)
 // just make sure padm is working
 TEST_EXP_SIZES_EOP(float)
 TEST_EXP_SIZES_EOP(double)
 TEST_EXP_SIZES_EOP(ldouble)
 TEST_EXP_SIZES_EOP(complexf)
 TEST_EXP_SIZES_EOP(complexd)
 TEST_EXP_SIZES_EOP(complexld)
 /* testing engine */
 #define QUOTE(...) STRFY(__VA_ARGS__)
 #define STRFY(...) #__VA_ARGS__
--- a/src/include/cppduals/tests/test_expm.cpp
+++ b/src/include/cppduals/tests/test_expm.cpp
@ -16,8 +16,8 @@
 * (c)2019 Michael Tesch. tesch1@gmail.com
 */
 #include "type_name.hpp"
 #include <duals/dual_eigen>
 #include "type_name.hpp"
 #include <Eigen/Dense>
 #include <Eigen/Sparse>
 #include <Eigen/StdVector>
@ -104,7 +104,7 @@ expm4(const Eigen::EigenBase<DerivedA> & A_,
  return R;
 }
-template <class T, int NN = 30, class DT = dual<T> >
+template <class T, int NN = 30, class DT = dual<T>, int N2 = 2*NN>
 void dexpm() {
  //typedef std::complex<float> T;
  //typedef std::complex<dual<T>> dualt;
@ -167,8 +167,6 @@ void dexpm() {
                                     << "eA2=" << eA2.block(0,0,std::min(4,NN),std::min(4,NN)) << "\n";
  EXPECT_LT((dA1 - dA2).norm(), tol) << "dA1=" << dA1.block(0,0,std::min(4,NN),std::min(4,NN)) << "\n"
                                     << "dA2=" << dA2.block(0,0,std::min(4,NN),std::min(4,NN)) << "\n";
 #undef NN
 #undef N2
 }
 #if defined(PHASE_1)
--- a/src/include/cppduals/tests/test_fmt.cpp
+++ b/src/include/cppduals/tests/test_fmt.cpp
@ -14,6 +14,7 @@
 * (c)2019 Michael Tesch. tesch1@gmail.com
 */
 #include <fmt/format.h>
 //#include <fmt/std.h>
 #define CPPDUALS_LIBFMT
 #define CPPDUALS_LIBFMT_COMPLEX
 #include <duals/dual>
@ -29,79 +30,131 @@ typedef std::complex<dualf> cdualf;
 using namespace duals::literals;
 using namespace std::complex_literals;
 TEST(libfmt, float_) {
  std::string s = fmt::format("{:.1f}", 2.f);
  EXPECT_EQ(s, "2.0");
 }
 TEST(libfmt, double_) {
  std::string s = fmt::format("{:.1f}", 2.);
  EXPECT_EQ(s, "2.0");
 }
 TEST(libfmt, complex_) {
-  std::string s = fmt::format("{}", 2.f + 3if);
+  std::string s = fmt::format("{:.1f}", 2.f + 3if);
  EXPECT_EQ(s, "(2.0+3.0if)");
  s = fmt::format("{:.2f}", 2.f + 3if);
  EXPECT_EQ(s, "(2.00+3.00if)");
 }
 TEST(libfmt, complex_el) {
  std::string s = fmt::format("{}", 2.l + 3il);
-  EXPECT_EQ(s, "(2.0+3.0il)");
+  EXPECT_EQ(s, "(2+3il)");
 }
 TEST(libfmt, complex_flags) {
  std::string s;
  s = fmt::format("{}", 2. + 3i);
-  EXPECT_EQ(s, "(2.0+3.0i)");
+  EXPECT_EQ(s, "(2+3i)");
-
+}
 TEST(libfmt, complex_flags_dollar) {
  std::string s;
  s = fmt::format("{:$}", 2. + 3i);
-  EXPECT_EQ(s, "(2.0+3.0i)");
+  EXPECT_EQ(s, "(2+3i)");
-
+}
 TEST(libfmt, complex_flags_star) {
  std::string s;
  s = fmt::format("{:*}", 2. + 3i);
  EXPECT_EQ(s, "(2+3*i)");
 }
 TEST(libfmt, complex_flags_star_f) {
  std::string s;
  s = fmt::format("{:*.1f}", 2. + 3i);
  EXPECT_EQ(s, "(2.0+3.0*i)");
 }
 TEST(libfmt, complex_flags_comma) {
  std::string s;
  s = fmt::format("{:,.1f}", 2. + 3i);
  EXPECT_EQ(s, "(2.0,3.0)");
  s = fmt::format("{:,}", 2. + 3i);
-  EXPECT_EQ(s, "(2.0,3.0)");
+  EXPECT_EQ(s, "(2,3)");
 }
 TEST(libfmt, complex_flags_a) {
  std::string s;
  // + +
  s = fmt::format("{:$+}", 2. + 3i);
  EXPECT_EQ(s, "(+2+3i)");
  s = fmt::format("{:$+.1f}", 2. + 3i);
  EXPECT_EQ(s, "(+2.0+3.0i)");
  s = fmt::format("{:*+}", 2. + 3i);
-  EXPECT_EQ(s, "(+2.0+3.0*i)");
+  EXPECT_EQ(s, "(+2+3*i)");
  s = fmt::format("{:,+}", 2. + 3i);
-  EXPECT_EQ(s, "(+2.0,+3.0)");
+  EXPECT_EQ(s, "(+2,+3)");
 }
 TEST(libfmt, complex_flags_b) {
  std::string s;
  // + -
  s = fmt::format("{:$+}", 2. - 3i);
-  EXPECT_EQ(s, "(+2.0-3.0i)");
+  EXPECT_EQ(s, "(+2-3i)");
  s = fmt::format("{:*+}", 2. - 3i);
-  EXPECT_EQ(s, "(+2.0-3.0*i)");
+  EXPECT_EQ(s, "(+2-3*i)");
  s = fmt::format("{:,+}", 2. - 3i);
  EXPECT_EQ(s, "(+2,-3)");
  s = fmt::format("{:,+.1f}", 2. - 3i);
  EXPECT_EQ(s, "(+2.0,-3.0)");
 }
 TEST(libfmt, complex_all_real)  {
  std::string s;
  s = fmt::format("{}", 2. + 0i);
  EXPECT_EQ(s, "(2)");
  s = fmt::format("{:.1f}", 2. + 0i);
  EXPECT_EQ(s, "(2.0)");
  s = fmt::format("{:*}", 2. + 0i);
-  EXPECT_EQ(s, "(2.0)");
+  EXPECT_EQ(s, "(2)");
  s = fmt::format("{:,}", 2. + 0i);
  EXPECT_EQ(s, "(2,0)");
  s = fmt::format("{:,.1f}", 2. + 0i);
  EXPECT_EQ(s, "(2.0,0.0)");
 }
 TEST(libfmt, complex_all_imag)  {
  std::string s;
  s = fmt::format("{}", 0. + 3i);
  EXPECT_EQ(s, "(3i)");
  s = fmt::format("{:.1f}", 0. + 3i);
  EXPECT_EQ(s, "(3.0i)");
  s = fmt::format("{:*}", 0. + 3i);
-  EXPECT_EQ(s, "(3.0*i)");
+  EXPECT_EQ(s, "(3*i)");
  s = fmt::format("{:,}", 0. + 3i);
-  EXPECT_EQ(s, "(0.0,3.0)");
+  EXPECT_EQ(s, "(0,3)");
 }
 TEST(libfmt, complex_plus)  {
  std::string s = fmt::format("{:+}", 1. + 3i);
  EXPECT_EQ(s, "(+1+3i)");
  s = fmt::format("{:+.1f}", 1. + 3i);
  EXPECT_EQ(s, "(+1.0+3.0i)");
  s = fmt::format("{:+}", 1. - 3i);
-  EXPECT_EQ(s, "(+1.0-3.0i)");
+  EXPECT_EQ(s, "(+1-3i)");
 }
 TEST(libfmt, complex_g) {
  std::string s = fmt::format("{:g}", 2.f + 3if);
@ -118,67 +171,102 @@ TEST(libfmt, complex_gel) {
 TEST(libfmt, dual_) {
  std::string s = fmt::format("{}", 2 + 3_ef);
  EXPECT_EQ(s, "(2+3_ef)");
  s = fmt::format("{:.1f}", 2 + 3_ef);
  EXPECT_EQ(s, "(2.0+3.0_ef)");
 }
 TEST(libfmt, dual_el) {
  std::string s = fmt::format("{}", 2 + 3_el);
  EXPECT_EQ(s, "(2+3_el)");
  s = fmt::format("{:.1f}", 2 + 3_el);
  EXPECT_EQ(s, "(2.0+3.0_el)");
 }
 TEST(libfmt, dual_flags) {
  std::string s;
  s = fmt::format("{}", 2. + 3_e);
-  EXPECT_EQ(s, "(2.0+3.0_e)");
+  EXPECT_EQ(s, "(2+3_e)");
  s = fmt::format("{:$}", 2. + 3_e);
  EXPECT_EQ(s, "(2+3_e)");
  s = fmt::format("{:$.1f}", 2. + 3_e);
  EXPECT_EQ(s, "(2.0+3.0_e)");
  s = fmt::format("{:*}", 2. + 3_e);
-  EXPECT_EQ(s, "(2.0+3.0*e)");
+  EXPECT_EQ(s, "(2+3*e)");
  s = fmt::format("{:,}", 2. + 3_e);
  EXPECT_EQ(s, "(2,3)");
  s = fmt::format("{:,.1f}", 2. + 3_e);
  EXPECT_EQ(s, "(2.0,3.0)");
  // + +
  s = fmt::format("{:$+}", 2. + 3_e);
-  EXPECT_EQ(s, "(+2.0+3.0_e)");
+  EXPECT_EQ(s, "(+2+3_e)");
  s = fmt::format("{:*+}", 2. + 3_e);
  EXPECT_EQ(s, "(+2+3*e)");
  s = fmt::format("{:*+.1f}", 2. + 3_e);
  EXPECT_EQ(s, "(+2.0+3.0*e)");
  s = fmt::format("{:,+}", 2. + 3_e);
-  EXPECT_EQ(s, "(+2.0,+3.0)");
+  EXPECT_EQ(s, "(+2,+3)");
  // + -
  s = fmt::format("{:$+}", 2. - 3_e);
-  EXPECT_EQ(s, "(+2.0-3.0_e)");
+  EXPECT_EQ(s, "(+2-3_e)");
  s = fmt::format("{:*+}", 2. - 3_e);
  EXPECT_EQ(s, "(+2-3*e)");
  s = fmt::format("{:*+.1f}", 2. - 3_e);
  EXPECT_EQ(s, "(+2.0-3.0*e)");
  s = fmt::format("{:,+}", 2. - 3_e);
-  EXPECT_EQ(s, "(+2.0,-3.0)");
+  EXPECT_EQ(s, "(+2,-3)");
  s = fmt::format("{:,+}", 2. + 3_e);
  EXPECT_EQ(s, "(+2,+3)");
 }
 TEST(libfmt, dual_all_real)  {
  std::string s;
  s = fmt::format("{}", 2 + 0_e);
-  EXPECT_EQ(s, "(2.0)");
+  EXPECT_EQ(s, "(2)");
  s = fmt::format("{:*}", 2 + 0_e);
-  EXPECT_EQ(s, "(2.0)");
+  EXPECT_EQ(s, "(2)");
  s = fmt::format("{:,}", 2 + 0_e);
-  EXPECT_EQ(s, "(2.0,0.0)");
+  EXPECT_EQ(s, "(2,0)");
 }
 TEST(libfmt, dual_all_dual)  {
  std::string s;
  s = fmt::format("a{}b", 0 + 3_e);
  EXPECT_EQ(s, "a(3_e)b");
  s = fmt::format("a{:.1f}b", 0 + 3_e);
  EXPECT_EQ(s, "a(3.0_e)b");
  s = fmt::format("a{:*}b", 0 + 3_e);
-  EXPECT_EQ(s, "a(3.0*e)b");
+  EXPECT_EQ(s, "a(3*e)b");
  s = fmt::format("a{:,}b", 0 + 3_e);
  EXPECT_EQ(s, "a(0,3)b");
  s = fmt::format("a{:,.1f}b", 0 + 3_e);
  EXPECT_EQ(s, "a(0.0,3.0)b");
 }
 TEST(libfmt, dual_plus)  {
  std::string s = fmt::format("{:+}", 1. + 3_e);
  EXPECT_EQ(s, "(+1+3_e)");
  s = fmt::format("{:+.1f}", 1. + 3_e);
  EXPECT_EQ(s, "(+1.0+3.0_e)");
  s = fmt::format("{:+}", 1. - 3_e);
  EXPECT_EQ(s, "(+1-3_e)");
 }
 TEST(libfmt, dual_g) {
  std::string s = fmt::format("{:g}", 2 + 3_ef);
  EXPECT_EQ(s, "(2+3_ef)");
@ -194,6 +282,12 @@ TEST(libfmt, dual_gel) {
  std::string s = fmt::format("{:g}", 2 + 3_el);
  EXPECT_EQ(s, "(2+3_el)");
 }
 TEST(libfmt, dual_cd) {
  std::string s = fmt::format("{}", cdualf(2 + 3_ef, 4 + 5_ef));
  EXPECT_EQ(s, "((2+3_ef)+(4+5_ef)i)");
 }
 TEST(libfmt, dual_cgt) {
  std::string s = fmt::format("{:g}", cdualf(2 + 3_ef, 4 + 5_ef));
  EXPECT_EQ(s, "((2+3_ef)+(4+5_ef)i)");
@ -208,6 +302,10 @@ TEST(libfmt, dual_cgts) {
  s = fmt::format("{:,*g}", cdualf(2 + 3_ef, 4 + 5_ef));
  EXPECT_EQ(s, "((2+3*ef),(4+5*ef))");
  // nonsense but should still work
  s = fmt::format("{:*,g}", cdualf(2 + 3_ef, 4 + 5_ef));
  EXPECT_EQ(s, "((2,3)+(4,5)*i)");
  s = fmt::format("{:,,g}", cdualf(2 + 3_ef, 4 + 5_ef));
  EXPECT_EQ(s, "((2,3),(4,5))");
 }
--- a/src/include/cppduals/tests/test_funcs.cpp
+++ b/src/include/cppduals/tests/test_funcs.cpp
@ -32,25 +32,38 @@ typedef std::complex<double> complexd;
 // not verify precision.
 #define FD_CHECK(T, F, ...)                                             \
  TEST(func##_##T, F) {                                                 \
    using std::isnan;                                                   \
    for (T x : __VA_ARGS__) {                                           \
      T prec = 100 * std::sqrt(std::numeric_limits<T>::epsilon());      \
-      T dd = dpart(F(x + dual<T>(0,1)));                                \
+      T dd = duals::dpart(F(x + dual<T>(0,1)));                         \
      /*T dx = std::numeric_limits<T>::epsilon() * (T)1000000;      */  \
      T dx = T(1)/ (1ull << (std::numeric_limits<T>::digits / 3));      \
      T fd = (F(x + dx) - F(x - dx)) / (2*dx);                          \
-      EXPECT_CNEAR(dd, fd, prec * std::abs(std::max(std::max(dd,fd),T(1)))) \
+      if (!isnan(dd) && !isnan(fd)) {                                   \
        EXPECT_CNEAR(dd, fd,                                            \
                     prec * std::abs(std::max(std::max(dd,fd),T(1))))   \
          << "dd=" << dd << " fd=" << fd << " x=" << x << " dx=" << dx; \
      }                                                                 \
    }                                                                   \
  }
 #define powL(x) pow(x,2)
 #define powR(x) pow(2,x)
 #define powLR(x) pow(x,x)
 FD_CHECK(double, exp, {-1,0,1})
-FD_CHECK(double, log, {1})
+FD_CHECK(double, log, {1, 3})
-//FD_CHECK(complexd, log, {-1,1}) TODO
+// FD_CHECK(complexd, log, {-1,1}) TODO
-FD_CHECK(double, log10, {1})
+FD_CHECK(double, log10, {1, 3})
-//FD_CHECK(complexd, log10, {-1,0,1}) TODO
+
 FD_CHECK(double, powL, {-3.,-1.,-0.4,0.,0.6,1.,2.5})
 FD_CHECK(double, powR, {-3.,-1.,-0.4,0.,0.6,1.,2.5})
 FD_CHECK(double, powLR, {-3.,-1.,-0.4,0.,0.6,1.,2.5})
 // FD_CHECK(complexd, log10, {-1,0,1}) TODO
 FD_CHECK(double, sqrt, {0.5,1.0})
 FD_CHECK(double, cbrt, {-10.,-0.01,0.01,1.0,10.})
-//FD_CHECK(complexd, sqrt, {0,1}) TODO
+// FD_CHECK(complexd, sqrt, {0,1}) TODO
 FD_CHECK(double, sin, {-1,0,1})
 FD_CHECK(double, cos, {-1,0,1})
 FD_CHECK(double, tan, {-1,0,1})
@ -59,12 +72,27 @@ FD_CHECK(double, acos, {-.9,0.,.9})
 FD_CHECK(double, atan, {-10,-1,0,1,10})
 // TODO:
-//FD_CHECK(double, sinh, {0})
+#define atan2L(x) atan2(x,2)
-//FD_CHECK(double, cosh, {0})
+#define atan2R(x) atan2(2,x)
-//FD_CHECK(double, tanh, {0})
+#define atan2LR(x) atan2(x,x)
-//FD_CHECK(double, asinh, {0})
+FD_CHECK(double, atan2L, {-10.,-1.,0.,1.,10.})
-//FD_CHECK(double, acosh, {0})
+FD_CHECK(double, atan2R, {-10.,-1.,0.01,1.,10.})
-//FD_CHECK(double, atanh, {0})
+FD_CHECK(double, atan2LR, {-10.,-1.,0.01,1.,10.})
 #define hypot2LR(x) hypot(x,x)
 FD_CHECK(double, hypot2LR, {-10.,-1.,0.01,1.,10.})
 #define scalbnL(x) scalbn(x,2)
 FD_CHECK(double, scalbnL, {-10.,-1.,0.01,1.,10.})
 FD_CHECK(double, sinh, {-0.1, 0.1})
 FD_CHECK(double, cosh, {-0.1, 0.1})
 FD_CHECK(double, tanh, {-0.1, 0.1})
 FD_CHECK(double, asinh, {-0.1, 0.1})
 FD_CHECK(double, acosh, {-1.1, 1.1})
 FD_CHECK(double, atanh, {-0.1, 0.1})
 FD_CHECK(double, log1p, {-0.1, 0.1})
 FD_CHECK(double, expm1, {-0.1, 0.1})
 FD_CHECK(double, erf, {-1,0,1})
 FD_CHECK(double, erfc, {-1,0,1})
@ -97,6 +125,8 @@ TEST(func, tgamma) {
  //EXPECT_EQ(tgamma(x).rpart(), 362880); "interestingly", compiling without optimization (-O0) causes this to fail
  EXPECT_NEAR(tgamma(x).rpart(), 362880, 362880 * 100 * std::numeric_limits<double>::epsilon());
 }
 // part selection functions
 TEST(func, rpart) {
  dualf x = 10 + 4_e;
  EXPECT_EQ(rpart(x), 10);
@ -105,16 +135,258 @@ TEST(func, dpart) {
  dualf x = 2 + 4_e;
  EXPECT_EQ(dpart(x), 4);
 }
 // non-differentiable operations on the real part.
 TEST(func, abs) {
  dualf x = -10 + 4_e;
  EXPECT_EQ(abs(x), 10);
 }
-TEST(func, arg) {
+TEST(func, fabs) {
  dualf x = -10 + 4_e;
  EXPECT_EQ(fabs(x), 10);
 }
 TEST(func, fmax) {
  dualf x = -10 + 4_e;
  dualf y = 10 + 4_e;
  EXPECT_EQ(fmax(x, y), 10);
 }
 TEST(func, fmin) {
  dualf x = -10 + 4_e;
  dualf y = 10 + 4_e;
  EXPECT_EQ(fmin(x, y), -10);
 }
 TEST(func, frexp) {
  dualf x = 6 + 4_e;
  int exp = 0;
  EXPECT_EQ(frexp(x, &exp), 0.75);
  EXPECT_EQ(exp, 3);
 }
 TEST(func, ldexp) {
  dualf x = 0.5 + 4_e;
  int exp = 1;
  EXPECT_EQ(ldexp(x, exp), 1);
 }
 TEST(func, trunc) {
  dualf x = 1.5 + 4_e;
  EXPECT_EQ(trunc(x), 1);
 }
 TEST(func, floor) {
  dualf x = 1.5 + 4_e;
  EXPECT_EQ(floor(x), 1);
 }
 TEST(func, ceil) {
  dualf x = 1.5 + 4_e;
  EXPECT_EQ(ceil(x), 2);
 }
 TEST(func, round) {
  dualf x = 1.5 + 4_e;
  EXPECT_EQ(round(x), 2);
 }
 // floating point functions
 TEST(func, fpclassify) {
  dualf x = 1.5 + 4_e;
  EXPECT_EQ(fpclassify(x), FP_NORMAL);
  EXPECT_EQ(fpclassify(std::numeric_limits<dualf>::infinity()), FP_INFINITE);
  EXPECT_EQ(fpclassify(std::numeric_limits<dualf>::quiet_NaN()), FP_NAN);
  if (std::numeric_limits<dualf>::has_denorm != std::denorm_absent) {
    EXPECT_EQ(fpclassify(std::numeric_limits<dualf>::denorm_min()), FP_SUBNORMAL);
  }
  EXPECT_EQ(fpclassify(x+std::numeric_limits<dualf>::min()), FP_NORMAL);
  EXPECT_EQ(fpclassify(2*std::numeric_limits<dualf>::max()), FP_INFINITE);
  EXPECT_EQ(fpclassify(x+std::numeric_limits<dualf>::epsilon()), FP_NORMAL);
 }
 TEST(func, isfinite) {
  dualf x = 1.5 + 4_e;
  EXPECT_EQ(isfinite(x), true);
 }
 TEST(func, isnormal) {
  dualf x = 1.5 + 4_e;
  EXPECT_EQ(isnormal(x), true);
 }
 TEST(func, isinf) {
  dualf x = 1.5 + 4_e;
  EXPECT_EQ(isinf(x), false);
 }
 TEST(func, isnan) {
  dualf x = 1.5 + 4_e;
  EXPECT_EQ(isnan(x), false);
 }
 TEST(func, signbit) {
  dualf x = 1.5 + 4_e;
  EXPECT_EQ(signbit(x), false);
 }
 TEST(func, copysign) {
  dualf x = 1.5 + 4_e;
  dualf y = -1.3 + 2_e;
  EXPECT_EQ(copysign(x, y), -1.5);
 }
 // Utility functions
 TEST(func, random) {
  dualf x = random(0.001 + 0.001_e, 1 + 1_e);
  EXPECT_GE(rpart(x), 0.001);
  EXPECT_LE(rpart(x), 1);
  EXPECT_GE(dpart(x), 0.001);
  EXPECT_LE(dpart(x), 1);
 }
 TEST(func, random2) {
  dualf x = duals::randos::random2(0.001 + 0.001_e, 1 + 1_e);
  EXPECT_GE(rpart(x), 0.001);
  EXPECT_LE(rpart(x), 1);
  EXPECT_GE(dpart(x), 0.001);
  EXPECT_LE(dpart(x), 1);
 }
 // more tests
 TEST(func, logb) {
  dualf x = 4 + 1_e;
  EXPECT_EQ(rpart(logb(.5 + 1_e)), -1.);
  EXPECT_EQ(rpart(logb(1 + 1_e)), 0.);
  EXPECT_EQ(rpart(logb(2 + 1_e)), 1.);
  EXPECT_EQ(rpart(logb(3 + 1_e)), 1.);
  EXPECT_EQ(rpart(logb(4 + 1_e)), 2.);
  EXPECT_EQ(rpart(logb(x * x)), 4.);
  EXPECT_EQ(dpart(logb(4 + 8_e)), std::numeric_limits<dualf>::infinity());
  EXPECT_EQ(dpart(logb(4.01 + 8_e)), 0.);
  EXPECT_EQ(dpart(logb(x * x)), std::numeric_limits<dualf>::infinity());
  EXPECT_EQ(dpart(logb(3 * x)), 0.);
  x += 0.01;
  EXPECT_EQ(dpart(logb(3 * x)), 0.);
  EXPECT_EQ(dpart(logb(x * x)), 0.);
 }
 TEST(func, pow) {
  dualf x = pow(0 + 0_e, 0.);
  EXPECT_EQ(rpart(x), 1);
  EXPECT_EQ(dpart(x), 0);
  dualf y = pow(0 + 0_e, 0. + 0_e);
  EXPECT_EQ(rpart(y), 1);
  EXPECT_EQ(dpart(y), 0);
  dualf z = pow(0, 0. + 0_e);
  EXPECT_EQ(rpart(z), 1);
  EXPECT_EQ(dpart(z), 0);
 }
 TEST(func, pow_complex) {
  std::complex<dualf> C(3+4_ef, 5+6_ef);
  std::complex<dualf> x = std::pow(C, 2+1_ef);
  std::complex<float> ref = std::complex<float>(std::pow(std::complex<float>(3,5), 2));
  EXPECT_NEAR(rpart(x.real()), ref.real(), 1e-5);
  EXPECT_NEAR(rpart(x.imag()), ref.imag(), 1e-6);
  //EXPECT_EQ(dpart(x), std::pow(std::complex<float>(3, 5), 2) * (2+1_ef));
 }
 //----------------------------------------------------------------------
 // Test for pow_dual_complex(const dual<T>& realBase, const std::complex<dual<T>>& complexExponent)
 //----------------------------------------------------------------------
 TEST(func, PowDualComplexTest)
 {
    // 1) Prepare inputs:
    //    realBase = 2.0 (with zero derivative for this example)
    duald realBase(2.0f);
    //    complexExponent = (1.5 + 0.7i),
    //    each part a dual with .rpart() = 1.5 or 0.7, derivative 0
    std::complex<duald> complexExponent(duald(1.5f), duald(0.7f));
    // 2) Call the function we want to test:
    //    x^y  =>  pow_dual_complex(realBase, complexExponent)
    std::complex<duald> result = std::pow(realBase, complexExponent);
    // 3) Reference: use standard pow in double
    double dblBase = 2.0;
    std::complex<double> dblExponent(1.5, 0.7);
    std::complex<double> reference = std::pow(dblBase, dblExponent);
    // 4) Compare the .rpart() of the dual’s real/imag with reference
    EXPECT_NEAR(result.real().rpart(), reference.real(), 1e-6);
    EXPECT_NEAR(result.imag().rpart(), reference.imag(), 1e-6);
 }
 //----------------------------------------------------------------------
 // Test for pow_complex_dual(const std::complex<dual<T>>& complexBase, const dual<T>& realExponent)
 //----------------------------------------------------------------------
 TEST(func, PowComplexDualTest)
 {
    // 1) Prepare inputs:
    //    complexBase = (3 + 5i), each part dual with zero derivative
    std::complex<duald> complexBase(duald(3.0f), duald(5.0f));
    //    realExponent = 2.0 as a dual
    duald realExponent(2.0f);
    // 2) Call the function we want to test:
    //    x^y => pow_complex_dual(complexBase, realExponent)
    std::complex<duald> result = pow(complexBase, realExponent);
    // 3) Reference: again, standard pow in double
    std::complex<double> dblBase(3.0, 5.0);
    double               dblExponent = 2.0;
    std::complex<double> reference   = std::pow(dblBase, dblExponent);
    // 4) Compare the .rpart() of the dual’s real/imag with reference
    EXPECT_NEAR(result.real().rpart(), reference.real(), 1e-6);
    EXPECT_NEAR(result.imag().rpart(), reference.imag(), 1e-6);
 }
 TEST(func, norm) {
  // TODO
 }
 TEST(func, conj) {
  // TODO
 }
 TEST(func, polar) {
  // TODO
 }
 TEST(func, atan) {
  EXPECT_EQ(rpart(atan(0 + 1_e)), atan(0));
  EXPECT_EQ(dpart(atan(1_e)), 1);
  EXPECT_EQ(dpart(atan(1 + 1_e)), 0.5);  // = 1 / (1 + x^2)
  EXPECT_EQ(dpart(atan(-2 + 1_e)), 1. / 5.);  // = 1 / (1 + x^2)
 }
 TEST(func, atan2) {
  // TODO
  //EXPECT_EQ(dpart(atan2(1_e, 1)), 1);
  //EXPECT_EQ(dpart(atan2(1 + 1_e, 1)), 0.5);
  //EXPECT_EQ(dpart(atan2(-2 + 1_e, 1)), (1. / 5.));
  duald y = 1 + 1_e;
  duald x = 1 + 0_e;
  auto z = atan2(y, x);
  z = atan2(y, x);
  EXPECT_EQ(rpart(z), atan2(rpart(y), rpart(x)));
  EXPECT_EQ(dpart(z), 0.5);
  y = -2 + 1_e;
  x = 1 + 0_e;
  z = atan2(y, x);
  EXPECT_EQ(rpart(z), atan2(rpart(y), rpart(x)));
  EXPECT_EQ(dpart(z), 1. / 5.);
  y = 1 + 0_e;
  x = -2 + 1_e;
  z = atan2(y, x);
  EXPECT_EQ(rpart(z), atan2(rpart(y), rpart(x)));
  EXPECT_EQ(dpart(z), -1. / 5.);
 }
 TEST(func, atan2a) {
  // TODO
  duald y = 1 + 1_e;
  EXPECT_EQ(rpart(atan2(y, 1)), atan2(rpart(y), 1));
  EXPECT_EQ(rpart(atan2(y, 2)), atan2(rpart(y), 2));
  EXPECT_EQ(dpart(atan2(y, 1)), 0.5);
  y = 2 + 1_e;
  EXPECT_EQ(dpart(atan2(y, -2)), -0.25);
 }
 TEST(func, atan2b) {
  // TODO
  duald x = 10 + 1_e;
  EXPECT_EQ(rpart(atan2(2, x)), atan2(2, rpart(x)));
  EXPECT_EQ(dpart(atan2(2, x)), -1./52);
 }
 struct pike_f1 {
  // function
--- a/src/include/cppduals/tests/test_packets.cpp
+++ b/src/include/cppduals/tests/test_packets.cpp
@ -16,8 +16,8 @@
 * (c)2019 Michael Tesch. tesch1@gmail.com
 */
 #include "type_name.hpp"
 #include <duals/dual_eigen>
 #include "type_name.hpp"
 #include <Eigen/Dense>
 #include <Eigen/Sparse>
 #include <Eigen/StdVector>
@ -60,6 +60,7 @@ template <int N=2, int K = N> using ecdf = Eigen::Matrix<cdualf, N, K> ;
 #if !defined(CPPDUALS_DONT_VECTORIZE) && !defined(EIGEN_DONT_VECTORIZE)
 #ifdef EIGEN_VECTORIZE_SSE
 TEST(Packet1cdf, pload_pstore) {
  using namespace Eigen::internal;
  cdualf cd1 = cdualf(1+2_ef,3+4_ef);
@ -68,6 +69,7 @@ TEST(Packet1cdf, pload_pstore) {
  pstore(&cd2, p1);
  EXPECT_DEQ(cd1, cd2);
 }
 #endif
 using duals::randos::random2;
@ -89,7 +91,7 @@ using duals::randos::random2;
    pstore(cd3.data(), p3);                                             \
    for (int i = 0; i < N; i++) {                                       \
      EXPECT_DNEAR(cd3[i], cd1[i] op cd2[i], tol)                       \
-        << cd1[i] << ',' << cd2[i] << " fail at " << i;                 \
+        << cd1[i] << #op << cd2[i] << "=" << cd3[i] << "!=" << (cd1[i] op cd2[i]) << " fail at " << i << "/" << N; \
    }                                                                   \
  }
@ -323,13 +325,11 @@ using duals::randos::random2;
  GEN_PACKET_TEST_UN(PTYPE,pconj,conj)
 // TODO:
 //pcplxflip
 //preduxp
 //pand
 //por
 //pxor
 //andnot
 //pbroadcast4
 //ploadquad (for packets w/ size==8)
 //pgather
 //pscatter
@ -345,7 +345,7 @@ GEN_CPACKET_TESTS(Packet2cf)
 GEN_CPACKET_TESTS(Packet4cf)
 #endif
-#ifdef EIGEN_VECTORIZE_SSE
+#if defined(EIGEN_VECTORIZE_SSE) || defined(EIGEN_VECTORIZE_NEON)
 GEN_PACKET_TESTS(Packet2df)
 GEN_PACKET_TESTS(Packet1dd)
 GEN_CPACKET_TESTS(Packet1cdf)
@ -368,20 +368,28 @@ TEST(compile, VECTORIZE) {
 #endif
 }
 TEST(compile, SSE) {
 #ifdef EIGEN_VECTORIZE_SSE
  EXPECT_TRUE(std::string(Eigen::SimdInstructionSetsInUse()).find("SSE") != std::string::npos)
-    << "Not using SSE instructions:" << Eigen::SimdInstructionSetsInUse();
+#else
-#ifndef EIGEN_VECTORIZE_SSE
+  EXPECT_TRUE(true)
  EXPECT_TRUE(false)
    << "Not using EIGEN_VECTORIZE_SSE:" << Eigen::SimdInstructionSetsInUse();
 #endif
    << "Not using EIGEN_VECTORIZE_SSE:" << Eigen::SimdInstructionSetsInUse();
 }
 TEST(compile, AVX) {
 #ifdef EIGEN_VECTORIZE_AVX
  EXPECT_TRUE(std::string(Eigen::SimdInstructionSetsInUse()).find("AVX") != std::string::npos)
-    << "Not using AVX instructions:" << Eigen::SimdInstructionSetsInUse();
+#else
-#ifndef EIGEN_VECTORIZE_AVX
+  EXPECT_TRUE(true)
  EXPECT_TRUE(false)
    << "Not using EIGEN_VECTORIZE_AVX:" << Eigen::SimdInstructionSetsInUse();
 #endif
    << "Not using EIGEN_VECTORIZE_AVX:" << Eigen::SimdInstructionSetsInUse();
 }
 TEST(compile, NEON) {
 #ifdef EIGEN_VECTORIZE_NEON
  EXPECT_TRUE(std::string(Eigen::SimdInstructionSetsInUse()).find("NEON") != std::string::npos)
 #else
  EXPECT_TRUE(true)
 #endif
    << "Not using EIGEN_VECTORIZE_NEON:" << Eigen::SimdInstructionSetsInUse();
 }
 #define QUOTE(...) STRFY(__VA_ARGS__)
--- a/src/include/cppduals/tests/test_solve.cpp
+++ b/src/include/cppduals/tests/test_solve.cpp
@ -16,8 +16,8 @@
 * (c)2019 Michael Tesch. tesch1@gmail.com
 */
 #include "type_name.hpp"
 #include <duals/dual_eigen>
 #include "type_name.hpp"
 #include <Eigen/Dense>
 #include <Eigen/Sparse>
 #include <Eigen/StdVector>
@ -78,7 +78,10 @@ void solveLu() {
  emtx<DT,NN> B = b + DT(0,1) * emtx<T,NN>::Random();
  emtx<DT,NN> C,D,E;
  C = A * B;
-  D = A.lu().solve(C);
+  // D = A.lu().solve(C);
  // D = A.colPivHouseholderQr().solve(C);
  // D = A.bdcSvd().solve(C);
  D = A.fullPivLu().solve(C);
  EXPECT_LT(rpart(B - D).norm(), tol);
  EXPECT_LT(dpart(B - D).norm(), tol);
 }
--- a/src/include/cppduals/tests/test_vectorize.cpp
+++ b/src/include/cppduals/tests/test_vectorize.cpp
@ -16,8 +16,8 @@
 * (c)2019 Michael Tesch. tesch1@gmail.com
 */
 #include "type_name.hpp"
 #include <duals/dual_eigen>
 #include "type_name.hpp"
 #include <Eigen/Dense>
 #include <Eigen/Sparse>
 #include <Eigen/StdVector>
@ -52,7 +52,7 @@ template <int N=2, int K = N> using ecdf = Eigen::Matrix<cdualf, N, K> ;
 #define _EXPECT_FALSE(...) {typedef __VA_ARGS__ fal; EXPECT_FALSE(fal::value); static_assert(!fal::value, "sa"); }
 #define ASSERT_DEQ(A,B) ASSERT_EQ(rpart(A), rpart(B)); ASSERT_EQ(dpart(A), dpart(B))
 #define ASSERT_DNEAR(A,B,tol)                                   \
-  ASSERT_NEAR(abs(rpart((A) - (B))),0,abs(rpart(A))*(tol));     \
+  ASSERT_NEAR(abs(rpart((A) - (B))),0,abs(rpart(A))*(tol)) << "rpart " << A << " " << B << "\n";     \
  ASSERT_NEAR(abs(dpart((A) - (B))),0,abs(dpart(A))*(tol))
 #define EXPECT_DEQ(A,B) EXPECT_EQ(rpart(A), rpart(B)); EXPECT_EQ(dpart(A), dpart(B))
 #define EXPECT_DNE(A,B) EXPECT_NE(rpart(A), rpart(B)); EXPECT_NE(dpart(A), dpart(B))
@ -117,9 +117,9 @@ void elemwise(int N) {
    ASSERT_DEQ(cca[i], Cca(i)) << "ca mismatch at " << i << "\n";
    ASSERT_DEQ(ccb[i], Ccb(i)) << "cb mismatch at " << i << "\n";
  }
-  ASSERT_DNEAR(sum, A.sum(), N*tol);
+  ASSERT_DNEAR(sum, A.sum(), N*tol) << "sum mismatch at " << sum << " " << A.sum() << "\n";
  ASSERT_DNEAR(sum, A.sum(), N*tol);
 }
 #if defined(PHASE_1)
 TEST(Vector, full_even_dualf) { elemwise<dualf>(512); }
 TEST(Vector, full_even_duald) { elemwise<duald>(512); }
@ -136,6 +136,7 @@ TEST(Vector, single_elem_cdualf) { elemwise<cdualf>(1); }
 TEST(Vector, two_elem_dualf) { elemwise<dualf>(2); }
 TEST(Vector, two_elem_duald) { elemwise<duald>(2); }
 TEST(Vector, two_elem_cdualf) { elemwise<cdualf>(2); }
 #endif
 #define DBOUT(X)
 #define MAKE_MULT_TEST(TYPE1, TYPE2, FIX, SIZE)                     \
@ -308,20 +309,30 @@ TEST(flags, VECTORIZE) {
 #endif
 }
 TEST(flags, SSE) {
 #ifdef EIGEN_VECTORIZE_SSE
  EXPECT_TRUE(std::string(Eigen::SimdInstructionSetsInUse()).find("SSE") != std::string::npos)
-    << "Not using SSE instructions:" << Eigen::SimdInstructionSetsInUse();
+#else
-#ifndef EIGEN_VECTORIZE_SSE
+  EXPECT_TRUE(true)
-  EXPECT_TRUE(false)
+#endif
    << "Not using EIGEN_VECTORIZE_SSE:" << Eigen::SimdInstructionSetsInUse();
 #endif
 }
 TEST(flags, AVX) {
 #ifdef EIGEN_VECTORIZE_AVX
  EXPECT_TRUE(std::string(Eigen::SimdInstructionSetsInUse()).find("AVX") != std::string::npos)
-    << "Not using AVX instructions:" << Eigen::SimdInstructionSetsInUse();
+#else
-#ifndef EIGEN_VECTORIZE_AVX
+  EXPECT_TRUE(true)
  EXPECT_TRUE(false)
    << "Not using EIGEN_VECTORIZE_AVX:" << Eigen::SimdInstructionSetsInUse();
 #endif
    << "Not using EIGEN_VECTORIZE_AVX:" << Eigen::SimdInstructionSetsInUse();
 }
 TEST(flags, NEON) {
 #ifdef EIGEN_VECTORIZE_NEON
  EXPECT_TRUE(std::string(Eigen::SimdInstructionSetsInUse()).find("NEON") != std::string::npos)
 #else
  EXPECT_TRUE(true)
 #endif
    << "Not using EIGEN_VECTORIZE_NEON:" << Eigen::SimdInstructionSetsInUse();
 }
 #define QUOTE(...) STRFY(__VA_ARGS__)
--- a/src/include/cppduals/thirdparty/CMakeLists.txt
+++ b/src/include/cppduals/thirdparty/CMakeLists.txt
@ -4,12 +4,19 @@
 # 3.10 adds support for "gtest_discover_tests" which enumerates the tests inside
 # of the code and adds them to ctest.
 #
-cmake_minimum_required (VERSION 3.10)
+cmake_minimum_required (VERSION 3.14)
 project (cppduals_thirdparty)
 include (ExternalProject)
 get_directory_property (hasParent PARENT_DIRECTORY)
 if (CMAKE_VERSION VERSION_GREATER_EQUAL 3.23)
  cmake_policy (SET CMP0135 NEW)
 endif ()
 if (CMAKE_XCODE_BUILD_SYSTEM VERSION_GREATER_EQUAL 12)
  cmake_policy (SET CMP0114 NEW)
 endif ()
 set (DEPS_ROOT "${CMAKE_BINARY_DIR}/root")
 if (hasParent)
  set (DEPS_ROOT "${CMAKE_BINARY_DIR}/thirdparty/root" PARENT_SCOPE)
@ -21,98 +28,109 @@ else (NOT WIN32)
  set (DOWNLOAD_DIR "C:/Downloads")
 endif (NOT WIN32)
 include(FetchContent)
 #
-# Google test (https://github.com/google/googletest/blob/master/googletest/README.md)
+# Google test
 # https://google.github.io/googletest/quickstart-cmake.html
 # 
 # Configure google-test
 FetchContent_Declare(
  googletest
  URL https://github.com/google/googletest/archive/refs/tags/v1.15.2.zip
  URL_HASH MD5=eb1c5c237d13ed12bf492d3997ca6b0d
  DOWNLOAD_NAME googletest-v1.15.2.zip
  DOWNLOAD_DIR "$ENV{HOME}/Downloads"
 )
-# Download and unpack googletest at configure time
+# For Windows: Prevent overriding the parent project's compiler/linker settings
-configure_file (CMakeLists-gt.txt.in googletest-download/CMakeLists.txt)
+set(gtest_force_shared_crt ON CACHE BOOL "" FORCE)
-execute_process (COMMAND ${CMAKE_COMMAND} -G "${CMAKE_GENERATOR}" .
+FetchContent_MakeAvailable(googletest)
  RESULT_VARIABLE result
  WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/googletest-download)
 if (result)
  message (FATAL_ERROR "CMake step for googletest failed: ${result}")
 endif ()
 execute_process (COMMAND ${CMAKE_COMMAND} --build .
  RESULT_VARIABLE result
  WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/googletest-download)
 if (result)
  message (FATAL_ERROR "Build step for googletest failed: ${result}")
 endif ()
-# Prevent overriding the parent project's compiler/linker
+# Enable testing
-# settings on Windows
+include(GoogleTest)
 set (gtest_force_shared_crt ON CACHE BOOL "" FORCE)
 # Add googletest directly to our build. This defines
 # the gtest and gtest_main targets.
 add_subdirectory (
  ${CMAKE_CURRENT_BINARY_DIR}/googletest-src
  ${CMAKE_CURRENT_BINARY_DIR}/googletest-build
  EXCLUDE_FROM_ALL)
 # The gtest/gtest_main targets carry header search path
 # dependencies automatically when using CMake 2.8.11 or
 # later. Otherwise we have to add them here ourselves.
 if (CMAKE_VERSION VERSION_LESS 2.8.11)
  include_directories ("${gtest_SOURCE_DIR}/include")
 endif ()
 # Can simply link against gtest or gtest_main as needed. Eg
 #add_executable (example example.cpp)
 #target_link_libraries (example gtest_main)
 #add_test (NAME example_test COMMAND example)
 #
 # Eigen
 #
-if (CPPDUALS_EIGEN_LATEST)
+if (TRUE)
-  set (EIGEN_URL http://bitbucket.org/eigen/eigen/get/default.tar.bz2)
+  set (EIGEN_URL https://gitlab.com/libeigen/eigen/-/archive/3.3.8/eigen-3.3.8.tar.bz2)
-  #set (EIGEN_MD5 ffc83130dcd37b694c6cf7e905099af9)
+  set (EIGEN_MD5 432ef01499d514f4606343276afa0ec3)
  set (EIGEN_MAX_CXX 17)
 else ()
-  set (EIGEN_URL http://bitbucket.org/eigen/eigen/get/3.3.7.tar.bz2)
+  set (EIGEN_URL https://gitlab.com/libeigen/eigen/-/archive/3.3.7/eigen-3.3.7.tar.bz2)
-  set (EIGEN_MD5 05b1f7511c93980c385ebe11bd3c93fa)
+  set (EIGEN_MD5 b9e98a200d2455f06db9c661c5610496)
  set (EIGEN_MAX_CXX 17)
 endif ()
-
+#
 # Eigen
 #
 ExternalProject_Add (eigenX
  PREFIX eigenX
  URL ${EIGEN_URL}
-  #URL_HASH MD5=${EIGEN_MD5}
+  URL_HASH MD5=${EIGEN_MD5}
  DOWNLOAD_DIR "$ENV{HOME}/Downloads"
  CONFIGURE_COMMAND ""
  BUILD_COMMAND ""
  INSTALL_COMMAND ""
  LOG_DOWNLOAD ON
  )
 ExternalProject_Get_Property (eigenX source_dir)
-if (hasParent AND NOT EIGEN3_INCLUDE_DIRS)
+if (true) # || hasParent
-  set (EIGEN3_INCLUDE_DIRS "${source_dir}" PARENT_SCOPE)
+  add_library (eigen INTERFACE IMPORTED GLOBAL)
  add_dependencies (eigen eigenX)
  set_property (TARGET eigen APPEND PROPERTY INTERFACE_INCLUDE_DIRECTORIES "${source_dir}")
  if (XCODE)
    set (IOFORMAT "IOFormat\(FullPrecision, DontAlignCols, \", \", \"\\$<SEMICOLON>\\\\n\", \"\", \"\", \"[\", \"]\"\)")
  else ()
    set (IOFORMAT "IOFormat\(FullPrecision, DontAlignCols, \", \", \"\\$<SEMICOLON>\\n\", \"\", \"\", \"[\", \"]\"\)")
  endif ()
  #target_compile_definitions (eigen INTERFACE EIGEN_DEFAULT_IO_FORMAT=${IOFORMAT})
  #target_compile_definitions (eigen INTERFACE EIGEN_DEFAULT_IO_FORMAT=EIGEN_IO_FORMAT_OCTAVE)
  set_property (TARGET eigen APPEND PROPERTY
    INTERFACE_COMPILE_DEFINITIONS EIGEN_DEFAULT_IO_FORMAT=${IOFORMAT})
 endif ()
 # if c++20, disable warning -Wdeprecated-enum-enum-conversion to eigen
 if (CMAKE_CXX_STANDARD GREATER_EQUAL 20)
  target_compile_options (eigen INTERFACE -Wno-deprecated-enum-enum-conversion)
 endif ()
 #
-# Eigen-Expokit
+# expokit
 #
-set (EEX_SHA 72bf6e445d5ae84218dcbd74580720491e0074db  )
+#set (EEX_SHA  ee28baa3bf29561501e17e5c68c2e54c85daae19  ) newer? used by spindropsSDL. md5=cebd15f9b5068c0e327753244ff6d394
 set (EEX_SHA c157dec0057be6e183a1ea2a5de353fac7e5e3a7 )
 set (EEX_MD5 89484e51f706398284235b96bc805515 )
 ExternalProject_Add (expokitX
  PREFIX expokitX
-  URL https://gitlab.com/api/v4/projects/tesch1%2Feigen-expokit/repository/archive.tbz2?sha=${EEX_SHA}
+  URL https://gitlab.com/tesch1/eigen-expokit/-/archive/${EEX_SHA}/eigen-expokit.tar.bz2
-  #URL_HASH MD5=96b79de1d01547f6d658865b7caa02ee
+  URL_HASH MD5=${EEX_MD5}
  DOWNLOAD_NAME eigen-expokit-${EEX_SHA}.tar.bz2
  DOWNLOAD_DIR "$ENV{HOME}/Downloads"
  UPDATE_COMMAND ""
  PATCH_COMMAND ""
  CONFIGURE_COMMAND ""
  BUILD_COMMAND ""
  INSTALL_COMMAND ""
  LOG_DOWNLOAD ON
  )
 ExternalProject_Get_Property (expokitX source_dir)
 if (hasParent)
-  set (EXPOKIT_INCLUDE_DIR "${source_dir}" PARENT_SCOPE)
+  add_library (expokit INTERFACE IMPORTED GLOBAL)
-endif()
+  add_dependencies (expokit expokitX)
  set_property (TARGET expokit APPEND PROPERTY INTERFACE_INCLUDE_DIRECTORIES "${source_dir}")
 endif (hasParent)
 #
 # fmt
 #
 ExternalProject_Add (fmtX
  PREFIX fmtX
-  URL https://github.com/fmtlib/fmt/archive/6.1.1.tar.gz
+  URL https://github.com/fmtlib/fmt/archive/11.1.4.tar.gz
-  URL_HASH MD5=acfb83d44afdca171ee26c597c931e7c
+  URL_HASH MD5=10c2ae163accd3b82e6b8b4dff877645
  DOWNLOAD_DIR ${DOWNLOAD_DIR}
  CONFIGURE_COMMAND ""
  BUILD_COMMAND ""
@ -121,7 +139,6 @@ ExternalProject_Add (fmtX
 ExternalProject_Get_Property (fmtX source_dir)
 ExternalProject_Get_Property (fmtX binary_dir)
 if (hasParent)
  message (" FMT3_INCLUDE_DIRS: ${source_dir}")
  add_subdirectory (${source_dir} ${binary_dir} EXCLUDE_FROM_ALL)
 endif ()
@ -131,19 +148,20 @@ if (CPPDUALS_BENCHMARK)
  #
  ExternalProject_Add (benchmarkX
    PREFIX benchmarkX
-    URL "http://github.com/google/benchmark/archive/v1.5.0.tar.gz"
+    URL "http://github.com/google/benchmark/archive/v1.9.1.tar.gz"
-    URL_HASH MD5=eb1466370f3ae31e74557baa29729e9e
+    URL_HASH MD5=92000ef8b4a7b1e9229972f8943070a7
    DOWNLOAD_DIR ${DOWNLOAD_DIR}
-    CMAKE_ARGS --target install -DBENCHMARK_ENABLE_GTEST_TESTS=OFF -DCMAKE_BUILD_TYPE=Release -DBENCHMARK_USE_LIBCXX=${CPPDUALS_USE_LIBCXX}
+    CONFIGURE_COMMAND ""
-    "-DCMAKE_INSTALL_PREFIX=<INSTALL_DIR>"
+    BUILD_COMMAND ""
-    INSTALL_DIR "${DEPS_ROOT}"
+    INSTALL_COMMAND ""
    )
  ExternalProject_Get_Property (benchmarkX source_dir)
-  ExternalProject_Get_Property (benchmarkX install_dir)
+  ExternalProject_Get_Property (benchmarkX binary_dir)
  if (hasParent)
-    set (BENCHMARK_SRC_DIR "${source_dir}" PARENT_SCOPE)
+    # https://github.com/google/benchmark#requirements
-    set (BENCHMARK_INC_DIR "${install_dir}/include" PARENT_SCOPE)
+    set (BENCHMARK_ENABLE_GTEST_TESTS OFF)
-    message (" BENCHMARK_SRC_DIR: ${BENCHMARK_SRC_DIR}")
+    set (BENCHMARK_USE_LIBCXX ${CPPDUALS_USE_LIBCXX})
    add_subdirectory (${source_dir} ${binary_dir} EXCLUDE_FROM_ALL)
  endif()
  if (Boost_FOUND AND NO)
@ -158,18 +176,6 @@ if (CPPDUALS_BENCHMARK)
      set (Boost_INCLUDE_DIRS ${Boost_INCLUDE_DIRS} PARENT_SCOPE)
    endif ()
    # piranha
    ExternalProject_Add (piranhaX PREFIX piranhaX
      URL https://github.com/bluescarni/piranha/archive/v0.11.tar.gz
      URL_HASH MD5=33482f719f6b8a6a5316f9bd148f5b10
      DOWNLOAD_DIR "$ENV{HOME}/Downloads"
      CONFIGURE_COMMAND "" BUILD_COMMAND "" INSTALL_COMMAND ""
      )
    ExternalProject_Get_Property (piranhaX source_dir)
    if (hasParent)
      set (PIRANHA_INCLUDE_DIR "${source_dir}/include" PARENT_SCOPE)
    endif ()
    # AuDi
    ExternalProject_Add (audiX PREFIX audiX
      URL https://github.com/darioizzo/audi/archive/v1.6.5.tar.gz