diff --git a/.travis.yml b/.travis.yml
index f53ea49074..5c3f6d59d9 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -124,7 +124,8 @@ matrix:
             - ubuntu-toolchain-r-test
 
     - os: linux
-      env: TOOLSET=gcc COMPILER=g++-6 CXXSTD=c++03 TEST_SUITE=quadrature
+      dist: bionic
+      env: TOOLSET=gcc COMPILER=g++-6 CXXSTD=c++03 TEST_SUITE="quadrature include=/usr/include/eigen3"
       addons:
         apt:
           packages:
@@ -132,11 +133,13 @@ matrix:
             - libgmp-dev
             - libmpfr-dev
             - libfftw3-dev
+            - libeigen3-dev
           sources:
             - ubuntu-toolchain-r-test
 
     - os: linux
-      env: TOOLSET=gcc COMPILER=g++-6 CXXSTD=c++11 TEST_SUITE=quadrature
+      dist: bionic
+      env: TOOLSET=gcc COMPILER=g++-6 CXXSTD=c++11 TEST_SUITE="quadrature include=/usr/include/eigen3"
       addons:
         apt:
           packages:
@@ -144,11 +147,13 @@ matrix:
             - libgmp-dev
             - libmpfr-dev
             - libfftw3-dev
+            - libeigen3-dev
           sources:
             - ubuntu-toolchain-r-test
 
     - os: linux
-      env: TOOLSET=gcc COMPILER=g++-6 CXXSTD=c++14 TEST_SUITE=quadrature
+      dist: bionic
+      env: TOOLSET=gcc COMPILER=g++-6 CXXSTD=c++14 TEST_SUITE="quadrature include=/usr/include/eigen3"
       addons:
         apt:
           packages:
@@ -156,11 +161,13 @@ matrix:
             - libgmp-dev
             - libmpfr-dev
             - libfftw3-dev
+            - libeigen3-dev
           sources:
             - ubuntu-toolchain-r-test
 
     - os: linux
-      env: TOOLSET=gcc COMPILER=g++-6 CXXSTD=gnu++14 TEST_SUITE=quadrature
+      dist: bionic
+      env: TOOLSET=gcc COMPILER=g++-6 CXXSTD=gnu++14 TEST_SUITE="quadrature include=/usr/include/eigen3"
       addons:
         apt:
           packages:
@@ -168,11 +175,13 @@ matrix:
             - libgmp-dev
             - libmpfr-dev
             - libfftw3-dev
+            - libeigen3-dev
           sources:
             - ubuntu-toolchain-r-test
 
     - os: linux
-      env: TOOLSET=gcc COMPILER=g++-6 CXXSTD=c++1z TEST_SUITE=quadrature
+      dist: bionic
+      env: TOOLSET=gcc COMPILER=g++-6 CXXSTD=c++1z TEST_SUITE="quadrature include=/usr/include/eigen3"
       addons:
         apt:
           packages:
@@ -180,6 +189,7 @@ matrix:
             - libgmp-dev
             - libmpfr-dev
             - libfftw3-dev
+            - libeigen3-dev
           sources:
             - ubuntu-toolchain-r-test
 
@@ -400,7 +410,8 @@ matrix:
             - ubuntu-toolchain-r-test
 
     - os: linux
-      env: TOOLSET=gcc COMPILER=g++-5 CXXSTD=c++14 TEST_SUITE=quadrature
+      dist: bionic
+      env: TOOLSET=gcc COMPILER=g++-5 CXXSTD=c++14 TEST_SUITE="quadrature include=/usr/include/eigen3"
       addons:
         apt:
           packages:
@@ -408,6 +419,7 @@ matrix:
             - libgmp-dev
             - libmpfr-dev
             - libfftw3-dev
+            - libeigen3-dev
           sources:
             - ubuntu-toolchain-r-test
 
@@ -470,7 +482,7 @@ matrix:
             - ubuntu-toolchain-r-test
 
     - os: linux
-      dist: trusty
+      dist: bionic
       compiler: g++-8
       env: TOOLSET=gcc COMPILER=g++-8 CXXSTD=c++14 TEST_SUITE=misc
       addons:
@@ -481,13 +493,14 @@ matrix:
             - ubuntu-toolchain-r-test
 
     - os: linux
-      dist: trusty
+      dist: bionic
       compiler: g++-8
-      env: TOOLSET=gcc COMPILER=g++-8 CXXSTD=c++14 TEST_SUITE=quadrature
+      env: TOOLSET=gcc COMPILER=g++-8 CXXSTD=c++14 TEST_SUITE="quadrature include=/usr/include/eigen3"
       addons:
         apt:
           packages:
             - g++-8
+            - libeigen3-dev
           sources:
             - ubuntu-toolchain-r-test
 
@@ -549,7 +562,7 @@ matrix:
             - ubuntu-toolchain-r-test
 
     - os: linux
-      dist: trusty
+      dist: bionic
       compiler: g++-8
       env: TOOLSET=gcc COMPILER=g++-8 CXXSTD=gnu++03 TEST_SUITE=misc
       addons:
@@ -560,13 +573,14 @@ matrix:
             - ubuntu-toolchain-r-test
 
     - os: linux
-      dist: trusty
+      dist: bionic
       compiler: g++-8
-      env: TOOLSET=gcc COMPILER=g++-8 CXXSTD=gnu++03 TEST_SUITE=quadrature
+      env: TOOLSET=gcc COMPILER=g++-8 CXXSTD=gnu++03 TEST_SUITE="quadrature include=/usr/include/eigen3"
       addons:
         apt:
           packages:
             - g++-8
+            - libeigen3-dev
           sources:
             - ubuntu-toolchain-r-test
 
@@ -617,7 +631,7 @@ matrix:
             - ubuntu-toolchain-r-test
 
     - os: linux
-      dist: trusty
+      dist: bionic
       compiler: g++-8
       env: TOOLSET=gcc COMPILER=g++-8 CXXSTD=gnu++17 TEST_SUITE=distribution_tests
       addons:
@@ -628,13 +642,14 @@ matrix:
             - ubuntu-toolchain-r-test
 
     - os: linux
-      dist: trusty
+      dist: bionic
       compiler: g++-8
-      env: TOOLSET=gcc COMPILER=g++-8 CXXSTD=gnu++17 TEST_SUITE=quadrature
+      env: TOOLSET=gcc COMPILER=g++-8 CXXSTD=gnu++17 TEST_SUITE="quadrature include=/usr/include/eigen3"
       addons:
         apt:
           packages:
             - g++-8
+            - libeigen3-dev
           sources:
             - ubuntu-toolchain-r-test
 
@@ -707,12 +722,14 @@ matrix:
             - llvm-toolchain-xenial-6.0
 
     - os: linux
+      dist: bionic
       compiler: clang++-6.0
-      env: TOOLSET=clang COMPILER=clang++-6.0 CXXSTD=c++11 TEST_SUITE=quadrature
+      env: TOOLSET=clang COMPILER=clang++-6.0 CXXSTD=c++11 TEST_SUITE="quadrature include=/usr/include/eigen3"
       addons:
         apt:
           packages:
             - clang-6.0
+            - libeigen3-dev
           sources:
             - ubuntu-toolchain-r-test
             - llvm-toolchain-xenial-6.0
@@ -754,8 +771,12 @@ matrix:
       osx_image: xcode11
 
     - os: osx
-      env: TOOLSET=clang COMPILER=clang++ CXXSTD=c++14 TEST_SUITE=quadrature
+      env: TOOLSET=clang COMPILER=clang++ CXXSTD=c++14 TEST_SUITE="quadrature include=/usr/local/include/eigen3"
       osx_image: xcode11
+      addons:
+        homebrew:
+          packages:
+            eigen
 
     - os: osx
       env: TOOLSET=clang COMPILER=clang++ CXXSTD=c++14 TEST_SUITE=float128_tests
diff --git a/appveyor.yml b/appveyor.yml
index bbee7ac0aa..01848bbcd2 100644
--- a/appveyor.yml
+++ b/appveyor.yml
@@ -16,97 +16,95 @@ platform:
 
 environment:
   matrix:
-    - ARGS: --toolset=msvc-9.0  address-model=32
-      TEST_SUITE: special_fun distribution_tests
-    - ARGS: --toolset=msvc-10.0 address-model=32
-      TEST_SUITE: special_fun distribution_tests
-    - ARGS: --toolset=msvc-11.0 address-model=32
-      TEST_SUITE: special_fun distribution_tests
-    - ARGS: --toolset=msvc-12.0 address-model=32
-      TEST_SUITE: special_fun distribution_tests
-    - ARGS: --toolset=msvc-14.0 address-model=32
-      TEST_SUITE: special_fun distribution_tests
-    - ARGS: --toolset=msvc-9.0  address-model=32
+    # - ARGS: --toolset=msvc-9.0  address-model=32
+    #   TEST_SUITE: special_fun distribution_tests
+    # - ARGS: --toolset=msvc-10.0 address-model=32
+    #   TEST_SUITE: special_fun distribution_tests
+    # - ARGS: --toolset=msvc-11.0 address-model=32
+    #   TEST_SUITE: special_fun distribution_tests
+    # - ARGS: --toolset=msvc-12.0 address-model=32
+    #   TEST_SUITE: special_fun distribution_tests
+    # - ARGS: --toolset=msvc-14.0 address-model=32
+    #   TEST_SUITE: special_fun distribution_tests
+    - ARGS: --include=/c/tools/vcpkg/installed/eigen3_x86-windows --toolset=msvc-9.0 address-model=32
       TEST_SUITE: misc quadrature ../example//examples
     - ARGS: --toolset=msvc-10.0 address-model=32
       TEST_SUITE: misc quadrature ../example//examples
-    - ARGS: --toolset=msvc-11.0 address-model=32
+    - ARGS: --toolset=msvc-11.0 address-model=32 include=/c/tools/vcpkg/installed/eigen3_x86-windows
       TEST_SUITE: misc quadrature ../example//examples
-    - ARGS: --toolset=msvc-12.0 address-model=32
+    - ARGS: --toolset=msvc-12.0 address-model=32 --include=/c/tools/vcpkg/installed/eigen3_x86-windows
       TEST_SUITE: misc quadrature ../example//examples
-    - ARGS: --toolset=msvc-14.0 address-model=32
+    - ARGS: --toolset=msvc-14.0 address-model=32 --include=/c/tools/vcpkg/installed/eigen3_x86-windows
       TEST_SUITE: misc quadrature ../example//examples
-    - ARGS: --toolset=msvc-12.0 address-model=64
-      TEST_SUITE: special_fun distribution_tests
-    - ARGS: --toolset=msvc-12.0 address-model=64
+    # - ARGS: --toolset=msvc-12.0 address-model=64
+    #   TEST_SUITE: special_fun distribution_tests
+    - ARGS: --toolset=msvc-12.0 address-model=64 --include=/c/tools/vcpkg/installed/eigen3_x86-windows
       TEST_SUITE: misc quadrature ../example//examples
-    - ARGS: --toolset=msvc-14.0 address-model=64
-      TEST_SUITE: special_fun
-    - ARGS: --toolset=msvc-14.0 address-model=64
-      TEST_SUITE: distribution_tests
-    - ARGS: --toolset=msvc-14.0 address-model=64
+    # - ARGS: --toolset=msvc-14.0 address-model=64
+    #   TEST_SUITE: special_fun
+    # - ARGS: --toolset=msvc-14.0 address-model=64
+    #   TEST_SUITE: distribution_tests
+    - ARGS: --toolset=msvc-14.0 address-model=64 --include=/c/tools/vcpkg/installed/eigen3_x86-windows
       TEST_SUITE: misc quadrature ../example//examples
 
-    - APPVEYOR_BUILD_WORKER_IMAGE: Visual Studio 2017
-      ARGS: --toolset=msvc-14.1 address-model=64 cxxstd=17
-      TEST_SUITE: special_fun distribution_tests
+    # - APPVEYOR_BUILD_WORKER_IMAGE: Visual Studio 2017
+    #   ARGS: --toolset=msvc-14.1 address-model=64 cxxstd=17
+    #   TEST_SUITE: special_fun distribution_tests
 
     - APPVEYOR_BUILD_WORKER_IMAGE: Visual Studio 2017
-      ARGS: --toolset=msvc-14.1 address-model=64 cxxstd=17
+      ARGS: --toolset=msvc-14.1 address-model=64 cxxstd=17 --include=/c/tools/vcpkg/installed/eigen3_x86-windows
       TEST_SUITE: misc quadrature ../example//examples
 
-    - APPVEYOR_BUILD_WORKER_IMAGE: Visual Studio 2017
-      ARGS: --toolset=msvc-14.1 address-model=32
-      TEST_SUITE: special_fun
+    # - APPVEYOR_BUILD_WORKER_IMAGE: Visual Studio 2017
+    #   ARGS: --toolset=msvc-14.1 address-model=32
+    #   TEST_SUITE: special_fun
 
-    - APPVEYOR_BUILD_WORKER_IMAGE: Visual Studio 2017
-      ARGS: --toolset=msvc-14.1 address-model=32
-      TEST_SUITE: distribution_tests
+    # - APPVEYOR_BUILD_WORKER_IMAGE: Visual Studio 2017
+    #   ARGS: --toolset=msvc-14.1 address-model=32
+    #   TEST_SUITE: distribution_tests
 
     - APPVEYOR_BUILD_WORKER_IMAGE: Visual Studio 2017
-      ARGS: --toolset=msvc-14.1 address-model=32
+      ARGS: --toolset=msvc-14.1 address-model=32 --include=/c/tools/vcpkg/installed/eigen3_x86-windows
       TEST_SUITE: misc quadrature ../example//examples
 
-    - ARGS: --toolset=gcc address-model=64
-      TEST_SUITE: float128_tests
-      PATH: C:\mingw-w64\x86_64-6.3.0-posix-seh-rt_v5-rev1\mingw64\bin;%PATH%
+    # - ARGS: --toolset=gcc address-model=64
+    #   TEST_SUITE: float128_tests
+    #   PATH: C:\mingw-w64\x86_64-6.3.0-posix-seh-rt_v5-rev1\mingw64\bin;%PATH%
 
-    - ARGS: --toolset=gcc address-model=64 cxxflags=-std=gnu++1z
-      TEST_SUITE: float128_tests
-      PATH: C:\mingw-w64\x86_64-6.3.0-posix-seh-rt_v5-rev1\mingw64\bin;%PATH%
+    # - ARGS: --toolset=gcc address-model=64 cxxflags=-std=gnu++1z
+    #   TEST_SUITE: float128_tests
+    #   PATH: C:\mingw-w64\x86_64-6.3.0-posix-seh-rt_v5-rev1\mingw64\bin;%PATH%
 
-    - ARGS: --toolset=gcc address-model=64
-      TEST_SUITE: special_fun
-      PATH: C:\mingw-w64\x86_64-6.3.0-posix-seh-rt_v5-rev1\mingw64\bin;%PATH%
+    # - ARGS: --toolset=gcc address-model=64
+    #   TEST_SUITE: special_fun
+    #   PATH: C:\mingw-w64\x86_64-6.3.0-posix-seh-rt_v5-rev1\mingw64\bin;%PATH%
 
-    - ARGS: --toolset=gcc address-model=64 cxxflags=-std=gnu++1z
-      TEST_SUITE: special_fun
-      PATH: C:\mingw-w64\x86_64-6.3.0-posix-seh-rt_v5-rev1\mingw64\bin;%PATH%
+    # - ARGS: --toolset=gcc address-model=64 cxxflags=-std=gnu++1z
+    #   TEST_SUITE: special_fun
+    #   PATH: C:\mingw-w64\x86_64-6.3.0-posix-seh-rt_v5-rev1\mingw64\bin;%PATH%
 
-    - ARGS: --toolset=gcc address-model=64
-      TEST_SUITE: distribution_tests
-      PATH: C:\mingw-w64\x86_64-6.3.0-posix-seh-rt_v5-rev1\mingw64\bin;%PATH%
+    # - ARGS: --toolset=gcc address-model=64
+    #   TEST_SUITE: distribution_tests
+    #   PATH: C:\mingw-w64\x86_64-6.3.0-posix-seh-rt_v5-rev1\mingw64\bin;%PATH%
 
-    - ARGS: --toolset=gcc address-model=64 cxxflags=-std=gnu++1z
-      TEST_SUITE: distribution_tests
-      PATH: C:\mingw-w64\x86_64-6.3.0-posix-seh-rt_v5-rev1\mingw64\bin;%PATH%
+    # - ARGS: --toolset=gcc address-model=64 cxxflags=-std=gnu++1z
+    #   TEST_SUITE: distribution_tests
+    #   PATH: C:\mingw-w64\x86_64-6.3.0-posix-seh-rt_v5-rev1\mingw64\bin;%PATH%
 
-    - ARGS: --toolset=gcc address-model=64
-      TEST_SUITE: misc
-      PATH: C:\mingw-w64\x86_64-6.3.0-posix-seh-rt_v5-rev1\mingw64\bin;%PATH%
+    # - ARGS: --toolset=gcc address-model=64
+    #   TEST_SUITE: misc
+    #   PATH: C:\mingw-w64\x86_64-6.3.0-posix-seh-rt_v5-rev1\mingw64\bin;%PATH%
 
-    - ARGS: --toolset=gcc address-model=64 cxxflags=-std=gnu++1z
+    - ARGS: --toolset=gcc address-model=64 cxxflags=-std=gnu++1z --include=/c/tools/vcpkg/installed/eigen3_x86-windows
       TEST_SUITE: quadrature
       PATH: C:\mingw-w64\x86_64-6.3.0-posix-seh-rt_v5-rev1\mingw64\bin;%PATH%
 
-    - ARGS: --toolset=gcc address-model=64
-      TEST_SUITE: ../example//examples
-      PATH: C:\mingw-w64\x86_64-6.3.0-posix-seh-rt_v5-rev1\mingw64\bin;%PATH%
-
-    - ARGS: --toolset=gcc address-model=64 cxxflags=-std=gnu++1z
-      TEST_SUITE: ../example//examples
-      PATH: C:\mingw-w64\x86_64-6.3.0-posix-seh-rt_v5-rev1\mingw64\bin;%PATH%
-
+    # - ARGS: --toolset=gcc address-model=64
+    #   TEST_SUITE: ../example//examples
+    #   PATH: C:\mingw-w64\x86_64-6.3.0-posix-seh-rt_v5-rev1\mingw64\bin;%PATH%
+    # - ARGS: --toolset=gcc address-model=64 cxxflags=-std=gnu++1z
+    #   TEST_SUITE: ../example//examples
+    #   PATH: C:\mingw-w64\x86_64-6.3.0-posix-seh-rt_v5-rev1\mingw64\bin;%PATH%
 
 install:
   - cd ..
@@ -122,6 +120,8 @@ install:
   - xcopy /s /e /q %APPVEYOR_BUILD_FOLDER% libs\math
   - git submodule update --init tools/boostdep
   - python tools/boostdep/depinst/depinst.py math
+  - vcpkg install eigen3
+  - vcpkg integrate install
   - bootstrap
   - b2 headers
 
diff --git a/include/boost/math/differentiation/autodiff.hpp b/include/boost/math/differentiation/autodiff.hpp
index a369593099..e4986321c5 100644
--- a/include/boost/math/differentiation/autodiff.hpp
+++ b/include/boost/math/differentiation/autodiff.hpp
@@ -146,15 +146,15 @@ class fvar {
   // RealType(ca) | RealType | RealType is copy constructible from the arithmetic types.
   explicit fvar(root_type const&);  // Initialize a constant. (No epsilon terms.)
 
-  template <typename RealType2>
+  template <typename RealType2, typename std::enable_if<std::is_convertible<RealType2, RealType>::value>::type* = nullptr>
   fvar(RealType2 const& ca);  // Supports any RealType2 for which static_cast<root_type>(ca) compiles.
 
   // r = cr | RealType& | Assignment operator.
   fvar& operator=(fvar const&) = default;
 
   // r = ca | RealType& | Assignment operator from the arithmetic types.
-  // Handled by constructor that takes a single parameter of generic type.
-  // fvar& operator=(root_type const&); // Set a constant.
+  template<typename RealType2>
+  typename std::enable_if<std::is_convertible<RealType2, RealType>::value, fvar&>::type& operator=(RealType2 const&);
 
   // r += cr | RealType& | Adds cr to r.
   template <typename RealType2, size_t Order2>
@@ -394,6 +394,12 @@ class fvar {
   template <typename RealType2, size_t Order2>
   friend std::ostream& operator<<(std::ostream&, fvar<RealType2, Order2> const&);
 
+  template <typename RealType2, size_t Order2>
+  friend std::istream& operator>>(std::istream&, fvar<RealType2, Order2> &);
+
+  template <typename RealType2, size_t Order2>
+  friend std::wistream& operator>>(std::wistream&, fvar<RealType2, Order2> &);
+
   // C++11 Compatibility
 #ifdef BOOST_NO_CXX17_IF_CONSTEXPR
   template <typename RootType>
@@ -679,19 +685,18 @@ fvar<RealType, Order>::fvar(root_type const& ca) : v{{static_cast<RealType>(ca)}
 
 // Can cause compiler error if RealType2 cannot be cast to root_type.
 template <typename RealType, size_t Order>
-template <typename RealType2>
+template <typename RealType2, typename std::enable_if<std::is_convertible<RealType2, RealType>::value>::type*>
 fvar<RealType, Order>::fvar(RealType2 const& ca) : v{{static_cast<RealType>(ca)}} {}
 
-/*
 template<typename RealType, size_t Order>
-fvar<RealType,Order>& fvar<RealType,Order>::operator=(root_type const& ca)
+template<typename RealType2>
+typename std::enable_if<std::is_convertible<RealType2, RealType>::value, fvar<RealType, Order>&>::type&  fvar<RealType,Order>::operator=(RealType2 const& ca)
 {
     v.front() = static_cast<RealType>(ca);
-    if constexpr (0 < Order)
+    BOOST_IF_CONSTEXPR (0 < Order)
         std::fill(v.begin()+1, v.end(), static_cast<RealType>(0));
     return *this;
 }
-*/
 
 template <typename RealType, size_t Order>
 template <typename RealType2, size_t Order2>
@@ -1741,6 +1746,22 @@ std::ostream& operator<<(std::ostream& out, fvar<RealType, Order> const& cr) {
   return out << ')';
 }
 
+template <typename RealType, size_t Order>
+std::istream& operator>>(std::istream& in, fvar<RealType, Order> & cr) {
+  in >> cr.v.front();
+  BOOST_IF_CONSTEXPR (0 < Order)
+      std::fill(cr.v.begin()+1, cr.v.end(), static_cast<RealType>(0));
+  return in;
+}
+
+template <typename RealType, size_t Order>
+std::wistream& operator>>(std::wistream& in, fvar<RealType, Order> & cr) {
+  in >> cr.v.front();
+  BOOST_IF_CONSTEXPR (0 < Order)
+      std::fill(cr.v.begin()+1, cr.v.end(), static_cast<RealType>(0));
+  return in;
+}
+
 // Additional functions
 
 template <typename RealType, size_t Order>
diff --git a/test/Jamfile.v2 b/test/Jamfile.v2
index f83b60b78c..f95bbfbb35 100644
--- a/test/Jamfile.v2
+++ b/test/Jamfile.v2
@@ -1336,6 +1336,7 @@ test-suite quadrature :
    [ run test_autodiff_6.cpp ../../test/build//boost_unit_test_framework : : : <toolset>gcc-mingw:<cxxflags>-Wa,-mbig-obj <debug-symbols>off <toolset>msvc:<cxxflags>/bigobj release [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : <linkflags>-lquadmath ] [ requires cxx11_inline_namespaces ] ]
    [ run test_autodiff_7.cpp ../../test/build//boost_unit_test_framework : : : <toolset>gcc-mingw:<cxxflags>-Wa,-mbig-obj <debug-symbols>off <toolset>msvc:<cxxflags>/bigobj release [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : <linkflags>-lquadmath ] [ requires cxx11_inline_namespaces ] ]
    [ run test_autodiff_8.cpp ../../test/build//boost_unit_test_framework : : : <toolset>gcc-mingw:<cxxflags>-Wa,-mbig-obj <debug-symbols>off <toolset>msvc:<cxxflags>/bigobj release [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : <linkflags>-lquadmath ] [ requires cxx11_inline_namespaces ] ]
+   [ run test_autodiff_9.cpp ../../test/build//boost_unit_test_framework : : : <toolset>gcc-mingw:<cxxflags>-Wa,-mbig-obj <debug-symbols>off <toolset>msvc:<cxxflags>/bigobj release [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : <linkflags>-lquadmath ] [ requires cxx11_inline_namespaces ] ]
    [ compile compile_test/autodiff_incl_test.cpp : <toolset>gcc-mingw:<cxxflags>-Wa,-mbig-obj <debug-symbols>off <toolset>msvc:<cxxflags>/bigobj release [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : <linkflags>-lquadmath ] [ requires cxx11_inline_namespaces ] ]
 ;
 
diff --git a/test/test_autodiff_9.cpp b/test/test_autodiff_9.cpp
new file mode 100644
index 0000000000..8b373c39fd
--- /dev/null
+++ b/test/test_autodiff_9.cpp
@@ -0,0 +1,312 @@
+#include <Eigen/Core>
+#include <Eigen/Dense>
+
+#include "test_autodiff.hpp"
+
+namespace Eigen {
+template <typename RealType, size_t Order>
+struct NumTraits<boost::math::differentiation::autodiff_v1::detail::
+                     template fvar<RealType, Order>> : NumTraits<RealType> {
+  using fvar =
+      boost::math::differentiation::autodiff_v1::detail::template fvar<RealType,
+                                                                       Order>;
+
+  enum {
+    RequireInitialization = 1,
+    ReadCost = 1,
+    AddCost = 16,
+    MulCost = 16,
+  };
+};
+
+#define BOOST_AUTODIFF_EIGEN_SCALAR_TRAITS(A)                                 \
+  template <class RealType, size_t Order, typename BinaryOp>                  \
+  struct ScalarBinaryOpTraits<boost::math::differentiation::autodiff_v1::     \
+                                  detail::template fvar<RealType, Order>,     \
+                              A, BinaryOp> {                                  \
+    typedef boost::math::differentiation::autodiff_v1::detail::template fvar< \
+        RealType, Order>                                                      \
+        ReturnType;                                                           \
+  };                                                                          \
+  template <class RealType, size_t Order, typename BinaryOp>                  \
+  struct ScalarBinaryOpTraits<A,                                              \
+                              boost::math::differentiation::autodiff_v1::     \
+                                  detail::template fvar<RealType, Order>,     \
+                              BinaryOp> {                                     \
+    typedef boost::math::differentiation::autodiff_v1::detail::template fvar< \
+        RealType, Order>                                                      \
+        ReturnType;                                                           \
+  };
+
+BOOST_AUTODIFF_EIGEN_SCALAR_TRAITS(float);
+BOOST_AUTODIFF_EIGEN_SCALAR_TRAITS(double);
+BOOST_AUTODIFF_EIGEN_SCALAR_TRAITS(long double);
+BOOST_AUTODIFF_EIGEN_SCALAR_TRAITS(short);
+BOOST_AUTODIFF_EIGEN_SCALAR_TRAITS(unsigned short);
+BOOST_AUTODIFF_EIGEN_SCALAR_TRAITS(int);
+BOOST_AUTODIFF_EIGEN_SCALAR_TRAITS(unsigned int);
+BOOST_AUTODIFF_EIGEN_SCALAR_TRAITS(long);
+BOOST_AUTODIFF_EIGEN_SCALAR_TRAITS(unsigned long);
+
+}  // namespace Eigen
+
+BOOST_AUTO_TEST_SUITE(test_autodiff_9)
+
+using boost::math::differentiation::autodiff_v1::detail::fvar;
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(eigen_init, T, all_float_types) {
+  constexpr int size = 5;
+
+  constexpr std::size_t n = 5;
+  typedef fvar<T, n> fTn;
+  Eigen::Matrix<fTn, size, 1> x;
+  x[0] = fTn(1.5);
+  x[1] = make_fvar<T, n - 1>(2.5);
+  x[2] = fvar<T, n - 2>(3.5);
+  x[3] = 4;
+  x[4] = 5.5;
+
+  constexpr std::size_t m = 2;
+  typedef fvar<T, m> fTm;
+  Eigen::Matrix<fTm, size, 1> y;
+  y = x.template cast<fTm>();
+  BOOST_CHECK_EQUAL(x[0].derivative(0), y[0].derivative(0));
+  BOOST_CHECK_EQUAL(x[1].derivative(0), y[1].derivative(0));
+  BOOST_CHECK_EQUAL(x[2].derivative(0), y[2].derivative(0));
+  BOOST_CHECK_EQUAL(x[3].derivative(0), y[3].derivative(0));
+  BOOST_CHECK_EQUAL(x[4].derivative(0), y[4].derivative(0));
+
+  // constexpr std::size_t p = 3;
+  // typedef fvar<T, p> fTp;
+  // Eigen::Matrix<fTp, 1, size> z =
+  //     Eigen::Matrix<fTn, size, 1>::Random().transpose().template cast<fTp>();
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(eigen_general, T, all_float_types) {
+  using std::cos;
+  using std::exp;
+  using std::log;
+  using std::pow;
+  using std::sin;
+
+  constexpr int dim = 4;
+  constexpr std::size_t n = 2;
+  constexpr double p = 3.456;
+
+  typedef fvar<T, n> fTn;
+  Eigen::Matrix<fTn, dim, 1> x;
+  x[0] = make_fvar<T, n>(-1);
+  x[1] = make_fvar<T, n>(0);
+  x[2] = make_fvar<T, n>(1);
+  x[3] = make_fvar<T, n>(5);
+
+  Eigen::Matrix<fTn, dim, 1> y1, y2, y3, y4, y5;
+  y1 = x.array().sin();
+  y2 = x.array().cos();
+  y3 = x.array().exp();
+  y4 = x.array().pow(p);
+  y5 = x.array().log();
+
+  // Check sin
+  BOOST_CHECK_EQUAL(y1[0].derivative(0), sin(x[0].derivative(0)));
+  BOOST_CHECK_EQUAL(y1[1].derivative(0), sin(x[1].derivative(0)));
+  BOOST_CHECK_EQUAL(y1[2].derivative(0), sin(x[2].derivative(0)));
+  BOOST_CHECK_EQUAL(y1[3].derivative(0), sin(x[3].derivative(0)));
+  BOOST_CHECK_EQUAL(y1[0].derivative(1), cos(x[0].derivative(0)));
+  BOOST_CHECK_EQUAL(y1[1].derivative(1), cos(x[1].derivative(0)));
+  BOOST_CHECK_EQUAL(y1[2].derivative(1), cos(x[2].derivative(0)));
+  BOOST_CHECK_EQUAL(y1[3].derivative(1), cos(x[3].derivative(0)));
+  BOOST_CHECK_EQUAL(y1[0].derivative(2), -sin(x[0].derivative(0)));
+  BOOST_CHECK_EQUAL(y1[1].derivative(2), -sin(x[1].derivative(0)));
+  BOOST_CHECK_EQUAL(y1[2].derivative(2), -sin(x[2].derivative(0)));
+  BOOST_CHECK_EQUAL(y1[3].derivative(2), -sin(x[3].derivative(0)));
+
+  // Check cos
+  BOOST_CHECK_EQUAL(y2[0].derivative(0), cos(x[0].derivative(0)));
+  BOOST_CHECK_EQUAL(y2[1].derivative(0), cos(x[1].derivative(0)));
+  BOOST_CHECK_EQUAL(y2[2].derivative(0), cos(x[2].derivative(0)));
+  BOOST_CHECK_EQUAL(y2[3].derivative(0), cos(x[3].derivative(0)));
+  BOOST_CHECK_EQUAL(y2[0].derivative(1), -sin(x[0].derivative(0)));
+  BOOST_CHECK_EQUAL(y2[1].derivative(1), -sin(x[1].derivative(0)));
+  BOOST_CHECK_EQUAL(y2[2].derivative(1), -sin(x[2].derivative(0)));
+  BOOST_CHECK_EQUAL(y2[3].derivative(1), -sin(x[3].derivative(0)));
+  BOOST_CHECK_EQUAL(y2[0].derivative(2), -cos(x[0].derivative(0)));
+  BOOST_CHECK_EQUAL(y2[1].derivative(2), -cos(x[1].derivative(0)));
+  BOOST_CHECK_EQUAL(y2[2].derivative(2), -cos(x[2].derivative(0)));
+  BOOST_CHECK_EQUAL(y2[3].derivative(2), -cos(x[3].derivative(0)));
+
+  // Check exp
+  BOOST_CHECK_EQUAL(y3[0].derivative(0), exp(x[0].derivative(0)));
+  BOOST_CHECK_EQUAL(y3[1].derivative(0), exp(x[1].derivative(0)));
+  BOOST_CHECK_EQUAL(y3[2].derivative(0), exp(x[2].derivative(0)));
+  BOOST_CHECK_EQUAL(y3[3].derivative(0), exp(x[3].derivative(0)));
+  BOOST_CHECK_EQUAL(y3[0].derivative(1), exp(x[0].derivative(0)));
+  BOOST_CHECK_EQUAL(y3[1].derivative(1), exp(x[1].derivative(0)));
+  BOOST_CHECK_EQUAL(y3[2].derivative(1), exp(x[2].derivative(0)));
+  BOOST_CHECK_EQUAL(y3[3].derivative(1), exp(x[3].derivative(0)));
+  BOOST_CHECK_EQUAL(y3[0].derivative(2), exp(x[0].derivative(0)));
+  BOOST_CHECK_EQUAL(y3[1].derivative(2), exp(x[1].derivative(0)));
+  BOOST_CHECK_EQUAL(y3[2].derivative(2), exp(x[2].derivative(0)));
+  BOOST_CHECK_EQUAL(y3[3].derivative(2), exp(x[3].derivative(0)));
+
+  // Check pow (without negative or zero)
+  T powTol = 1e-4;
+  BOOST_CHECK_CLOSE(y4[2].derivative(0), pow(x[2].derivative(0), p), powTol);
+  BOOST_CHECK_CLOSE(y4[3].derivative(0), pow(x[3].derivative(0), p), powTol);
+  BOOST_CHECK_CLOSE(y4[2].derivative(1), p * pow(x[2].derivative(0), p - 1),
+                    powTol);
+  BOOST_CHECK_CLOSE(y4[3].derivative(1), p * pow(x[3].derivative(0), p - 1),
+                    powTol);
+  BOOST_CHECK_CLOSE(y4[2].derivative(2),
+                    (p - 1) * p * pow(x[2].derivative(0), p - 2), powTol);
+  BOOST_CHECK_CLOSE(y4[3].derivative(2),
+                    (p - 1) * p * pow(x[3].derivative(0), p - 2), powTol);
+
+  // Check log (without negative or zero)
+  T logTol = 1e-5;
+  BOOST_CHECK_CLOSE(y5[2].derivative(0), log(x[2].derivative(0)), logTol);
+  BOOST_CHECK_CLOSE(y5[3].derivative(0), log(x[3].derivative(0)), logTol);
+  BOOST_CHECK_CLOSE(y5[2].derivative(1), 1 / x[2].derivative(0), logTol);
+  BOOST_CHECK_CLOSE(y5[3].derivative(1), 1 / x[3].derivative(0), logTol);
+  BOOST_CHECK_CLOSE(y5[2].derivative(2), -1 / pow(x[2].derivative(0), 2),
+                    logTol);
+  BOOST_CHECK_CLOSE(y5[3].derivative(2), -1 / pow(x[3].derivative(0), 2),
+                    logTol);
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(eigen_scalar, T, all_float_types) {
+  constexpr int dim = 4;
+  constexpr size_t n = 4;
+
+  typedef fvar<T, n> fTn;
+  fTn x = make_fvar<T, n>(4);
+  Eigen::Matrix<fTn, dim, 1> X;
+  Eigen::Matrix<fTn, dim, 1> Z;
+  Eigen::Matrix<fTn, dim, dim> I = Eigen::Matrix<fTn, dim, dim>::Identity();
+  X[0] = x;
+  Z[0] = x;
+  X[1] = 2;
+  Z[1] = x;
+  X[2] = x;
+  Z[2] = -3;
+  X[3] = 4;
+  Z[3] = 5;
+
+  // y = x*x + 2*x - 3*x + 20
+  fTn y1 = X.transpose() * Z;
+  fTn y2 = Z.transpose() * I * X;
+
+  BOOST_CHECK_EQUAL(y1.derivative(0), x * x + 2 * x - 3 * x + 20);
+  BOOST_CHECK_EQUAL(y1.derivative(1), 2 * x - 1);
+  BOOST_CHECK_EQUAL(y1.derivative(2), 2);
+  BOOST_CHECK_EQUAL(y1.derivative(3), 0);
+  BOOST_CHECK_EQUAL(y1.derivative(4), 0);
+  BOOST_CHECK_EQUAL(y2.derivative(0), x * x + 2 * x - 3 * x + 20);
+  BOOST_CHECK_EQUAL(y2.derivative(1), 2 * x - 1);
+  BOOST_CHECK_EQUAL(y2.derivative(2), 2);
+  BOOST_CHECK_EQUAL(y2.derivative(3), 0);
+  BOOST_CHECK_EQUAL(y2.derivative(4), 0);
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(eigen_vector, T, all_float_types) {
+  // Note: Can't handle taking gradient all at once. That would require the same
+  // approach as `make_ftuple`, which has upper limits.
+  using std::cos;
+  using std::sin;
+
+  constexpr int dim = 3;
+  constexpr size_t n = 4;
+
+  typedef fvar<T, n> fTn;
+  fTn x = make_fvar<T, n>(5);
+  T xD0 = x.derivative(0);
+  Eigen::Matrix<fTn, dim, 1> X;
+  X[0] = 1;
+  X[1] = x;
+  X[2] = x * x;
+  Eigen::Matrix<fTn, dim, dim> M;
+  M(0, 0) = 1;
+  M(0, 1) = 2;
+  M(0, 2) = x;
+  M(1, 0) = 1 / x;
+  M(1, 1) = 4;
+  M(1, 2) = 0;
+  M(2, 0) = 5;
+  M(2, 1) = sin(x);
+  M(2, 2) = cos(x * x);
+
+  Eigen::Matrix<fTn, dim, 1> Y = M * X;
+
+  T powTol = 1e-4;
+
+  // Y[0] = 1 + 2*x + x*x*x
+  BOOST_CHECK_CLOSE(Y[0].derivative(0), 1 + 2 * xD0 + pow(xD0, 3), powTol);
+  BOOST_CHECK_EQUAL(Y[0].derivative(1), 2 + 3 * xD0 * xD0);
+  BOOST_CHECK_EQUAL(Y[0].derivative(2), 6 * xD0);
+  BOOST_CHECK_EQUAL(Y[0].derivative(3), 6);
+  BOOST_CHECK_EQUAL(Y[0].derivative(4), 0);
+
+  // Y[1] = 1/x + 4*x + 0
+  BOOST_CHECK_EQUAL(Y[1].derivative(0), 1 / xD0 + 4 * xD0);
+  BOOST_CHECK_EQUAL(Y[1].derivative(1), -1 / (xD0 * xD0) + 4);
+  BOOST_CHECK_CLOSE(Y[1].derivative(2), 2 / (pow(xD0, 3)), powTol);
+  BOOST_CHECK_CLOSE(Y[1].derivative(3), -6 / (pow(xD0, 3) * xD0), powTol);
+  BOOST_CHECK_CLOSE(Y[1].derivative(4), 24 / (pow(xD0, 3) * xD0 * xD0), powTol);
+
+  // Y[2] = 5 + x*sin(x) + x*x*cos(x*x)
+  BOOST_CHECK_EQUAL(Y[2].derivative(0),
+                    5 + xD0 * sin(xD0) + xD0 * xD0 * cos(xD0 * xD0));
+  BOOST_CHECK_CLOSE(Y[2].derivative(1),
+                    2 * xD0 * cos(xD0 * xD0) -
+                        2 * pow(xD0, 3) * sin(xD0 * xD0) + sin(xD0) +
+                        xD0 * cos(xD0),
+                    powTol);
+  BOOST_CHECK_CLOSE(Y[2].derivative(2),
+                    -xD0 * (10 * xD0 * sin(xD0 * xD0) + sin(xD0)) +
+                        (2 - 4 * pow(xD0, 4)) * cos(xD0 * xD0) + 2 * cos(xD0),
+                    powTol);
+  BOOST_CHECK_CLOSE(
+      Y[2].derivative(3),
+      -24 * xD0 * sin(xD0 * xD0) + 8 * pow(xD0, 5) * sin(xD0 * xD0) -
+          36 * pow(xD0, 3) * cos(xD0 * xD0) - 3 * sin(xD0) - xD0 * cos(xD0),
+      powTol);
+  BOOST_CHECK_CLOSE(
+      Y[2].derivative(4),
+      -24 * sin(xD0 * xD0) + 112 * pow(xD0, 4) * sin(xD0 * xD0) +
+          4 * (4 * pow(xD0, 4) - 39) * xD0 * xD0 * cos(xD0 * xD0) +
+          xD0 * sin(xD0) - 4 * cos(xD0),
+      powTol);
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(eigen_determinant, T, all_float_types) {
+  using std::cos;
+  using std::sin;
+  constexpr int dim = 4;
+  constexpr size_t n = 1;
+
+  typedef fvar<T, n> fTn;
+  fTn x = make_fvar<T, n>(3);
+  T xD0 = x.derivative(0);
+  Eigen::Matrix<fTn, dim, dim> M = 10 * Eigen::Matrix<fTn, dim, dim>::Random();
+  M(0, 3) = x;
+  M(1, 0) = 3 * x;
+  M(1, 2) = 1 / x;
+  M(2, 2) = x * x;
+  M(3, 1) = sin(x);
+  fTn y = M.determinant();
+
+  Eigen::Matrix<T, dim, dim> dMdx = Eigen::Matrix<T, dim, dim>::Zero();
+  dMdx(0, 3) = 1;
+  dMdx(1, 0) = 3;
+  dMdx(1, 2) = -1 / (xD0 * xD0);
+  dMdx(2, 2) = 2 * xD0;
+  dMdx(3, 1) = cos(xD0);
+
+  T detTol = 1e-3;
+
+  T ans = y.derivative(0) *
+          (M.inverse() * dMdx.template cast<fTn>()).trace().derivative(0);
+  BOOST_CHECK_CLOSE(y.derivative(1), ans, detTol);
+}
+
+BOOST_AUTO_TEST_SUITE_END()