Repository: darioizzo/audi (language: python)

Issue:
Herratic behavior of vectorized type
computations made by the vectorized type seem to be weirdly off in some cases:

```python
from pyaudi import gdual_double as gdual
from pyaudi import gdual_vdouble as gdualv
def f(x,y):
    return x*y/(y-x)
x = gdual(1.1, "x", 2)
y = gdual(2.1, "y", 2)
xv = gdualv([1, 1.1], "x", 2)
yv = gdualv([1.1, 2.1], "y", 2)
```
note that the vectorized type has the x=1.1, y=2.1 point and thus it should be computing the same values of the non vectorized type in its second component.

then:
```python
In [2]: f(x,y)                                                                                                                                                                   
Out[2]: -4.62*dx*dy-1.21*dy+2.31+4.41*dx+4.41*dx**2+1.21*dy**2
```
note that the mixed coefficient ``dx*dy`` is ``4.62``
but ...
```python
In [4]: f(xv, yv)                                                                                                                                                                
Out[4]: [11, 2.31]+[121, 4.41]*dx+[-100, -1.21]*dy+[1210, 4.41]*dx**2+[1000, 1.21]*dy**2+[-2200, -220]*dx*dy
```

note that the mixed coefficient ``dx*dy`` is ``[-2200, -220]``, while its second component should be`` 4.62``.
the mistake is not always there, it depends also from the value of the other coefficients!


Discussion / hints:
Note that if, instead, we write:
```python
xv = gdualv([1.1, 1.1], "x", 2) 
yv = gdualv([2.1, 2.1], "y", 2)   
```
the result stands correct:
```python
In [4]: f(xv, yv)                                                                                                                                                                
Out[4]: [-1.21, -1.21]*dy+[2.31, 2.31]+[1.21, 1.21]*dy**2+[4.41, 4.41]*dx**2+[-4.62, -4.62]*dx*dy+[4.41, 4.41]*dx
```

Fix patch:
diff --git a/.travis.yml b/.travis.yml
index 106fd26c..f49e689d 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -23,10 +23,11 @@ matrix:
     - docker
   - env: AUDI_BUILD="OSXDebug" AUDI_BUILD_TYPE="Debug"
     os: osx
-    osx_image: xcode10.3
+    osx_image: xcode9.4
   - env: AUDI_BUILD="OSXRelease" AUDI_BUILD_TYPE="Release"
     os: osx
-    osx_image: xcode10.3
+    osx_image: xcode9.4
+
 
 script:
   - if [[ "${AUDI_BUILD}" == OSX* ]]; then bash tools/travis_osx.sh; fi
diff --git a/include/audi/vectorized.hpp b/include/audi/vectorized.hpp
index ab5b2c2a..fcf44ca6 100644
--- a/include/audi/vectorized.hpp
+++ b/include/audi/vectorized.hpp
@@ -433,14 +433,12 @@ inline std::size_t byte_size(const vectorized<T> &c) noexcept
 template <typename T>
 inline void fma3(vectorized<T> &ret, const vectorized<T> &x, const vectorized<T> &y)
 {
-    const auto x_size = x.size(), y_size = y.size();
-
+    const auto x_size = x.size(), y_size = y.size(), ret_size = ret.size();
     constexpr auto use_fma = obake::is_mult_addable_v<T &, const T &, const T &>;
 
-    if (x_size == y_size) {
-        ret.resize(x_size);
+    // Most frequent case: all dimensions are equal case: n,n,n
+    if ((x_size == y_size) && (x_size == ret_size)) {
         auto ret_it = ret.begin();
-
         for (decltype(x.size()) i = 0; i < x_size; ++i, ++ret_it) {
             if constexpr (use_fma) {
                 obake::fma3(*ret_it, x[i], y[i]);
@@ -448,32 +446,83 @@ inline void fma3(vectorized<T> &ret, const vectorized<T> &x, const vectorized<T>
                 *ret_it += x[i] * y[i];
             }
         }
-    } else if (x_size == 1u) {
-        ret.resize(y_size);
-        auto ret_it = ret.begin();
-
-        const auto x0 = x[0];
-        for (decltype(x.size()) i = 0; i < y_size; ++i, ++ret_it) {
-            if constexpr (use_fma) {
-                obake::fma3(*ret_it, x0, y[i]);
-            } else {
-                *ret_it += x0 * y[i];
+    } else if (ret_size == 1u) { // We are in the case 1, .., ..
+        if (x_size == 1u) {      // We are in the case 1, 1, n with n > 1
+            ret.resize(y_size, ret[0]);
+            auto ret_it = ret.begin();
+            const auto x0 = x[0];
+            for (decltype(y.size()) i = 0; i < y_size; ++i, ++ret_it) {
+                if constexpr (use_fma) {
+                    obake::fma3(*ret_it, x0, y[i]);
+                } else {
+                    *ret_it += x0 * y[i];
+                }
+            }
+        } else if (y_size == 1u) { // We are in the case 1, n, 1 with n > 1
+            ret.resize(x_size, ret[0]);
+            auto ret_it = ret.begin();
+            const auto y0 = y[0];
+            for (decltype(x.size()) i = 0; i < x_size; ++i, ++ret_it) {
+                if constexpr (use_fma) {
+                    obake::fma3(*ret_it, x[i], y0);
+                } else {
+                    *ret_it += x[i] * y0;
+                }
+            }
+        } // We are in the case 1, n, n with n > 1
+        else if (y_size == x_size) {
+            ret.resize(x_size, ret[0]);
+            auto ret_it = ret.begin();
+            for (decltype(x.size()) i = 0; i < x_size; ++i, ++ret_it) {
+                if constexpr (use_fma) {
+                    obake::fma3(*ret_it, x[i], y[i]);
+                } else {
+                    *ret_it += x[i] * y[i];
+                }
+            }
+        } else {
+            throw std::invalid_argument("Coefficients of different sizes in fma3");
+        }
+    } else if (x_size == 1u) { // We are in the case n, 1
+        if (y_size == 1u) {    // We are in the case n, 1, 1
+            auto ret_it = ret.begin();
+            const auto x0 = x[0];
+            const auto y0 = y[0];
+            for (decltype(ret.size()) i = 0u; i < ret_size; ++i, ++ret_it) {
+                if constexpr (use_fma) {
+                    obake::fma3(*ret_it, x0, y0);
+                } else {
+                    *ret_it += x0 * y0;
+                }
+            }
+        } else if (y_size == ret_size) { // We are in the case n, 1, n
+            auto ret_it = ret.begin();
+            const auto x0 = x[0];
+            for (decltype(y.size()) i = 0; i < y_size; ++i, ++ret_it) {
+                if constexpr (use_fma) {
+                    obake::fma3(*ret_it, x0, y[i]);
+                } else {
+                    *ret_it += x0 * y[i];
+                }
             }
+        } else {
+            throw std::invalid_argument("Coefficients of different sizes in fma3");
         }
-    } else if (y_size == 1u) {
-        ret.resize(x_size);
-        auto ret_it = ret.begin();
 
-        const auto y0 = y[0];
-        for (decltype(x.size()) i = 0; i < x_size; ++i, ++ret_it) {
-            if constexpr (use_fma) {
-                obake::fma3(*ret_it, x[i], y0);
-            } else {
-                *ret_it += x[i] * y0;
+    } else if (y_size == 1u) {    // We are in the case n, n, 1
+        if (ret_size == x_size) { // We are in the case n, n, 1
+            auto ret_it = ret.begin();
+            const auto y0 = y[0];
+            for (decltype(x.size()) i = 0; i < x_size; ++i, ++ret_it) {
+                if constexpr (use_fma) {
+                    obake::fma3(*ret_it, x[i], y0);
+                } else {
+                    *ret_it += x[i] * y0;
+                }
             }
+        } else {
+            throw std::invalid_argument("Coefficients of different sizes in fma3");
         }
-    } else {
-        throw std::invalid_argument("Coefficients of different sizes in fma3");
     }
 }
 
diff --git a/pyaudi/common_utils.hpp b/pyaudi/common_utils.hpp
index 0d43c53e..f60d76f6 100644
--- a/pyaudi/common_utils.hpp
+++ b/pyaudi/common_utils.hpp
@@ -18,7 +18,7 @@ inline ::obake::symbol_map<T> py_dict_to_obake_sm(const py::dict &d)
     typename ::obake::symbol_map<T>::sequence_type seq;
     seq.reserve(::obake::safe_cast<decltype(seq.size())>(py::len_hint(d)));
 
-    for (const auto &[k, v] : d) {
+    for (const auto [k, v] : d) {
         seq.emplace_back(k.template cast<::std::string>(), v.template cast<T>());
     }
 
diff --git a/tools/travis_osx.sh b/tools/travis_osx.sh
index 2a22b1d5..e99a0d94 100644
--- a/tools/travis_osx.sh
+++ b/tools/travis_osx.sh
@@ -6,45 +6,28 @@ set -x
 # Exit on error.
 set -e
 
-wget https://repo.anaconda.com/miniconda/Miniconda3-latest-MacOSX-x86_64.sh -O miniconda.sh;
+wget https://repo.continuum.io/miniconda/Miniconda2-latest-MacOSX-x86_64.sh -O miniconda.sh;
+
 export deps_dir=$HOME/local
 export PATH="$HOME/miniconda/bin:$PATH"
+export PATH="$deps_dir/bin:$PATH"
+
 bash miniconda.sh -b -p $HOME/miniconda
-conda config --add channels conda-forge --force
+conda config --add channels conda-forge
+conda config --set channel_priority strict
 
-conda_pkgs="cmake>=3.3 clang clangdev eigen obake-devel mppp boost boost-cpp python=3.7"
+conda_pkgs="cmake eigen obake-devel mppp boost-cpp python=3.7 pybind11 numpy"
 
-conda create -q -p $deps_dir -y $conda_pkgs
+conda create -q -p $deps_dir -y
 source activate $deps_dir
+conda install $conda_pkgs -y
 
-export deps_dir=$HOME/local
-export PATH="$HOME/miniconda/bin:$PATH"
-export PATH="$deps_dir/bin:$PATH"
-
-export CXX=clang++
-export CC=clang
 
 mkdir build
 
-# Install Pybind11 (making sure its the same used in our pipeline)
-export PYAUDI_BUILD_DIR=`pwd`
-git clone https://github.com/pybind/pybind11.git
-cd pybind11
-git checkout 4f72ef846fe8453596230ac285eeaa0ce3278bb4
-mkdir build
-cd build
-pwd
-cmake \
-    -DPYBIND11_TEST=NO \
-    -DCMAKE_INSTALL_PREFIX=$PYAUDI_BUILD_DIR \
-    -DCMAKE_PREFIX_PATH=$PYAUDI_BUILD_DIR \
-    ..
-make install
-cd ../..
-
 # Install audi
 cd build
-cmake -DBoost_NO_BOOST_CMAKE=ON \
+CXX=clang++ CC=clang cmake -DBoost_NO_BOOST_CMAKE=ON \
     -DCMAKE_INSTALL_PREFIX=$deps_dir \
     -DCMAKE_PREFIX_PATH=$deps_dir \
     -DCMAKE_BUILD_TYPE=${AUDI_BUILD_TYPE} \


Test patch:
diff --git a/tests/gdual_v.cpp b/tests/gdual_v.cpp
index f4a5eba5..81e30c33 100644
--- a/tests/gdual_v.cpp
+++ b/tests/gdual_v.cpp
@@ -20,229 +20,31 @@ using namespace audi;
 
 using gdual_v = audi::gdual<audi::vectorized<double>>;
 
-BOOST_AUTO_TEST_CASE(construction)
+template <typename T>
+T f(const T &x, const T &y)
 {
-    // Constructing a vectorized constant
-    { // 1 - from initializer_list
-        BOOST_CHECK_THROW(gdual_v(std::vector<double>{}), std::invalid_argument);
-        BOOST_CHECK_EQUAL(gdual_v{}, gdual_v(0.));
-        gdual_v x{0.};
-        gdual_v y(std::vector<double>{1.2, 2.2});
-        BOOST_CHECK_EQUAL(x.get_order(), 0);
-        BOOST_CHECK_EQUAL(y.get_order(), 0);
-        auto x0 = x.constant_cf();
-        auto y0 = y.constant_cf();
-        BOOST_CHECK_EQUAL(x0[0], 0.);
-        BOOST_CHECK_EQUAL(y0[0], 1.2);
-        BOOST_CHECK_EQUAL(y0[1], 2.2);
-        BOOST_CHECK_EQUAL(x.get_symbol_set_size(), 0u);
-        BOOST_CHECK_EQUAL(y.get_symbol_set_size(), 0u);
-    }
-    { // 2 - from an std::vector
-        BOOST_CHECK_THROW(gdual_v(std::vector<double>{}), std::invalid_argument);
-        gdual_v x(std::vector<double>(1, 0.));
-        gdual_v y(std::vector<double>{1.2, 2.2});
-        BOOST_CHECK_EQUAL(x.get_order(), 0);
-        BOOST_CHECK_EQUAL(y.get_order(), 0);
-        auto x0 = x.constant_cf();
-        auto y0 = y.constant_cf();
-        BOOST_CHECK_EQUAL(x0[0], 0.);
-        BOOST_CHECK_EQUAL(y0[0], 1.2);
-        BOOST_CHECK_EQUAL(y0[1], 2.2);
-        BOOST_CHECK_EQUAL(x.get_symbol_set_size(), 0u);
-        BOOST_CHECK_EQUAL(y.get_symbol_set_size(), 0u);
-    }
-    // Constructing a "full" vectorized gdual_v
-    { // 1 - from initializer_list
-        BOOST_CHECK_THROW(gdual_v(std::vector<double>{}, "x", 3), std::invalid_argument);
-        BOOST_CHECK_THROW(gdual_v(std::vector<double>{2, 3}, "dx", 3), std::invalid_argument);
-        BOOST_CHECK_THROW(gdual_v(std::vector<double>{2, 3}, "x", std::numeric_limits<unsigned int>::max() - 8u),
-                          std::invalid_argument);
-        gdual_v x(std::vector<double>{0.}, "x", 3);
-        gdual_v y(std::vector<double>{1.2, 2.2}, "y", 3);
-        BOOST_CHECK_EQUAL(x.get_order(), 3);
-        BOOST_CHECK_EQUAL(y.get_order(), 3);
-        auto x0 = x.constant_cf();
-        auto y0 = y.constant_cf();
-        BOOST_CHECK_EQUAL(x0[0], 0.);
-        BOOST_CHECK_EQUAL(y0[0], 1.2);
-        BOOST_CHECK_EQUAL(y0[1], 2.2);
-        BOOST_CHECK_EQUAL(x.get_symbol_set_size(), 1u);
-        BOOST_CHECK_EQUAL(y.get_symbol_set_size(), 1u);
-        auto x1 = x.get_derivative({1});
-        auto y1 = y.get_derivative({1});
-        BOOST_CHECK_EQUAL(x1[0], 1.);
-        BOOST_CHECK_EQUAL(y1[0], 1.);
-        BOOST_CHECK_EQUAL(y1[1], 1.);
-    }
-    { // 2 - from an std::vector
-        BOOST_CHECK_THROW(gdual_v(std::vector<double>{}, "x", 3), std::invalid_argument);
-        BOOST_CHECK_THROW(gdual_v(std::vector<double>{2, 3}, "dx", 3), std::invalid_argument);
-        BOOST_CHECK_THROW(gdual_v(std::vector<double>{2, 3}, "x", std::numeric_limits<unsigned int>::max() - 8u),
-                          std::invalid_argument);
-        gdual_v x(std::vector<double>{0.}, "x", 3);
-        gdual_v y(std::vector<double>{1.2, 2.2}, "y", 3);
-        BOOST_CHECK_EQUAL(x.get_order(), 3);
-        BOOST_CHECK_EQUAL(y.get_order(), 3);
-        auto x0 = x.constant_cf();
-        auto y0 = y.constant_cf();
-        BOOST_CHECK_EQUAL(x0[0], 0.);
-        BOOST_CHECK_EQUAL(y0[0], 1.2);
-        BOOST_CHECK_EQUAL(y0[1], 2.2);
-        BOOST_CHECK_EQUAL(x.get_symbol_set_size(), 1u);
-        BOOST_CHECK_EQUAL(y.get_symbol_set_size(), 1u);
-        auto x1 = x.get_derivative({1});
-        auto y1 = y.get_derivative({1});
-        BOOST_CHECK_EQUAL(x1[0], 1.);
-        BOOST_CHECK_EQUAL(y1[0], 1.);
-        BOOST_CHECK_EQUAL(y1[1], 1.);
-    }
+    return (x * y) / (y - x);
 }
-BOOST_AUTO_TEST_CASE(arithmetic_plus)
-{
-    { // nominal case
-        gdual_v x(std::vector<double>{1., 1.}, "x", 3);
-        gdual_v y(std::vector<double>{-1., 1}, "y", 3);
-        gdual_v z(std::vector<double>{1.2, 1.3, 1.4}, "z", 3);
-        BOOST_CHECK_THROW(x + z, std::invalid_argument);
-        auto sum = x + y;
-        BOOST_CHECK_EQUAL(sum.get_symbol_set_size(), 2u);
-        auto sum0 = sum.constant_cf();
-        auto sumdx = sum.get_derivative({1, 0});
-        auto sumdy = sum.get_derivative({0, 1});
-        BOOST_CHECK_EQUAL(sum0[0], 0.);
-        BOOST_CHECK_EQUAL(sum0[1], 2.);
-        BOOST_CHECK_EQUAL(sumdx[0], 1.);
-        BOOST_CHECK_EQUAL(sumdx[1], 1.);
-        BOOST_CHECK_EQUAL(sumdy[0], 1.);
-        BOOST_CHECK_EQUAL(sumdy[1], 1.);
-    }
-    { // scalar case
-        auto N = 10u;
-        gdual_v x(std::vector<double>(N, 123.), "x", 3);
-        auto sum = x + 5.;
-        auto sum2 = x + gdual_v(std::vector<double>{5.});
-        auto sum3 = 5. + x;
-        BOOST_CHECK_EQUAL(sum, sum2);
-        BOOST_CHECK_EQUAL(sum2, sum3);
-        for (auto i = 0u; i < N; ++i) {
-            BOOST_CHECK_EQUAL(sum.constant_cf()[i], 128.);
-            BOOST_CHECK_EQUAL(sum.get_derivative({1})[i], 1.);
-        }
-    }
-}
-BOOST_AUTO_TEST_CASE(arithmetic_minus)
-{
-    {
-        gdual_v x(std::vector<double>{1., 1.}, "x", 3);
-        gdual_v y(std::vector<double>{-1., 1}, "y", 3);
-        gdual_v z(std::vector<double>{1.2, 1.3, 1.4}, "z", 3);
-        BOOST_CHECK_THROW(x - z, std::invalid_argument);
-        auto diff = x - y;
-        BOOST_CHECK_EQUAL(diff.get_symbol_set_size(), 2u);
-        auto diff0 = diff.constant_cf();
-        auto diffdx = diff.get_derivative({1, 0});
-        auto diffdy = diff.get_derivative({0, 1});
-        BOOST_CHECK_EQUAL(diff0[0], 2.);
-        BOOST_CHECK_EQUAL(diff0[1], 0.);
-        BOOST_CHECK_EQUAL(diffdx[0], 1.);
-        BOOST_CHECK_EQUAL(diffdx[1], 1.);
-        BOOST_CHECK_EQUAL(diffdy[0], -1.);
-        BOOST_CHECK_EQUAL(diffdy[1], -1.);
-    }
-    { // scalar case
-        auto N = 10u;
-        gdual_v x(std::vector<double>(N, 123.), "x", 3);
-        auto diff = x - 5.;
-        auto diff2 = x - gdual_v(std::vector<double>{5.});
-        auto diff3 = 5. - x;
-        BOOST_CHECK_EQUAL(diff, diff2);
-        BOOST_CHECK_EQUAL(diff2, -diff3);
-        for (auto i = 0u; i < N; ++i) {
-            BOOST_CHECK_EQUAL(diff.constant_cf()[i], 118.);
-            BOOST_CHECK_EQUAL(diff.get_derivative({1})[i], 1.);
-        }
-    }
-}
-BOOST_AUTO_TEST_CASE(arithmetic_mul)
-{
-    {
-        gdual_v x(std::vector<double>{1., 1.}, "x", 3);
-        gdual_v y(std::vector<double>{-1., 1}, "y", 3);
-        gdual_v z(std::vector<double>{1.2, 1.3, 1.4}, "z", 3);
-        BOOST_CHECK_THROW(x * z, std::invalid_argument);
-        auto mul = x * y;
-        BOOST_CHECK_EQUAL(mul.get_symbol_set_size(), 2u);
-        auto mul0 = mul.constant_cf();
-        auto muldx = mul.get_derivative({1, 0});
-        auto muldy = mul.get_derivative({0, 1});
-        BOOST_CHECK_EQUAL(mul0[0], -1.);
-        BOOST_CHECK_EQUAL(mul0[1], 1.);
-        BOOST_CHECK_EQUAL(muldx[0], -1.);
-        BOOST_CHECK_EQUAL(muldx[1], 1.);
-        BOOST_CHECK_EQUAL(muldy[0], 1.);
-        BOOST_CHECK_EQUAL(muldy[1], 1.);
-    }
-    { // scalar case
-        auto N = 10u;
-        gdual_v x(std::vector<double>(N, 123.), "x", 3);
-        auto mul = x * 5.;
-        auto mul2 = x * gdual_v(std::vector<double>{5.});
-        auto mul3 = 5. * x;
-        BOOST_CHECK_EQUAL(mul, mul2);
-        BOOST_CHECK_EQUAL(mul3, mul2);
-        for (auto i = 0u; i < N; ++i) {
-            BOOST_CHECK_EQUAL(mul.constant_cf()[i], 123. * 5.);
-            BOOST_CHECK_EQUAL(mul.get_derivative({1})[i], 5.);
-        }
-    }
-}
-BOOST_AUTO_TEST_CASE(arithmetic_div)
+
+template <typename T>
+T dfdxy(const T &x, const T &y)
 {
-    {
-        gdual_v x(std::vector<double>{1., 1.}, "x", 3);
-        gdual_v y(std::vector<double>{-1., 1}, "y", 3);
-        gdual_v z(std::vector<double>{1.2, 1.3, 1.4}, "z", 3);
-        BOOST_CHECK_THROW(x / z, std::invalid_argument);
-        auto div = x / y;
-        BOOST_CHECK_EQUAL(div.get_symbol_set_size(), 2u);
-        auto div0 = div.constant_cf();
-        auto divdx = div.get_derivative({1, 0});
-        auto divdy = div.get_derivative({0, 1});
-        BOOST_CHECK_EQUAL(div0[0], -1.);
-        BOOST_CHECK_EQUAL(div0[1], 1.);
-        BOOST_CHECK_EQUAL(divdx[0], -1.);
-        BOOST_CHECK_EQUAL(divdx[1], 1.);
-        BOOST_CHECK_EQUAL(divdy[0], -1.);
-        BOOST_CHECK_EQUAL(divdy[1], -1.);
-    }
-    { // scalar case
-        auto N = 10u;
-        gdual_v x(std::vector<double>(N, 50.), "x", 3);
-        auto div = x / 5.;
-        auto div2 = x / gdual_v(std::vector<double>{5.});
-        auto div3 = 5. / x;
-        BOOST_CHECK_EQUAL(div, div2);
-        BOOST_CHECK_EQUAL(div2, 1. / div3); // this will not work if floats are not representable
-        for (auto i = 0u; i < N; ++i) {
-            BOOST_CHECK_EQUAL(div.constant_cf()[i], 50. / 5.);
-            BOOST_CHECK_EQUAL(div.get_derivative({1})[i], 1 / 5.);
-        }
-    }
+    return -2 * x * y / (y - x) / (y - x) / (y - x);
 }
 
-BOOST_AUTO_TEST_CASE(comparisons)
+BOOST_AUTO_TEST_CASE(bug_issue_57)
 {
-    gdual_v x(std::vector<double>{1.1, 1.3}, "x", 4);
-    gdual_v y(std::vector<double>{3.1, 3.2}, "y", 4);
-    BOOST_CHECK(x < y);
-    BOOST_CHECK(x * x < y * y);
-    BOOST_CHECK(audi::sin(x) < y * y);
-    BOOST_CHECK(audi::sin(x) < exp(y));
-    BOOST_CHECK(y > x);
-    BOOST_CHECK(y * y > x * x);
-    BOOST_CHECK(y * y > audi::sin(x));
-    BOOST_CHECK(exp(y) > audi::sin(x));
-    BOOST_CHECK(!(y > y));
-    BOOST_CHECK(!(y < y));
+    double a = 0.;
+    double b = 1.;
+
+    gdual_v x(std::vector<double>{1, a}, "x", 2);
+    gdual_v y(std::vector<double>{2, b}, "y", 2);
+    // BOOST_CHECK_EQUAL(f(x, y).find_cf({1, 1})[1], dfdxy(a, b));
+    std::cout << "\nWRONG: " << std::endl;
+    auto wrong = x * y / (y - x);
+    std::cout << "Wrong: " << wrong << std::endl;
+    std::cout << "\nRIGHT: " << std::endl;
+    auto right = (1. / (y - x)) * (x * y) ;
+    std::cout << "Right: " << right << std::endl;
+
 }
diff --git a/tests/vectorized.cpp b/tests/vectorized.cpp
index 0710c7a0..0ecffde0 100644
--- a/tests/vectorized.cpp
+++ b/tests/vectorized.cpp
@@ -103,6 +103,135 @@ void test_math()
         BOOST_CHECK(x1 / x2 == (vectorized<T>{T(1.) / T(3.), T(2.) / T(3.)}));
         BOOST_CHECK(x1 * x2 == (vectorized<T>{T(3.), T(6.)}));
     }
+    {
+        // testing the behaviour of vectorized size promotion
+        vectorized<T> x1{T(10.), T(20.), T(30.)};
+        vectorized<T> x2{T(2.)};
+        vectorized<T> x3{T(4.)};
+        BOOST_CHECK(x1 + x2 == (vectorized<T>{T(12.), T(22.), T(32.)}));
+        BOOST_CHECK(x1 - x2 == (vectorized<T>{T(8.), T(18.), T(28.)}));
+        BOOST_CHECK(x1 * x2 == (vectorized<T>{T(20.), T(40.), T(60.)}));
+        BOOST_CHECK(x1 / x2 == (vectorized<T>{T(5.), T(10.), T(15.)}));
+
+        {
+            auto retval = x1;
+            fma3(retval, x2, x3);
+            std::cout << retval << std::endl;
+            BOOST_CHECK(retval == (vectorized<T>{T(18.), T(28.), T(38.)}));
+
+            retval = x1;
+            fma3(retval, x3, x2);
+            std::cout << retval << std::endl;
+            BOOST_CHECK(retval == (vectorized<T>{T(18.), T(28.), T(38.)}));
+
+            retval = x2;
+            fma3(retval, x1, x3);
+            std::cout << retval << std::endl;
+            BOOST_CHECK(retval == (vectorized<T>{T(42.), T(82.), T(122.)}));
+
+            retval = x2;
+            fma3(retval, x3, x1);
+            std::cout << retval << std::endl;
+            BOOST_CHECK(retval == (vectorized<T>{T(42.), T(82.), T(122.)}));
+
+            retval = x3;
+            fma3(retval, x1, x2);
+            std::cout << retval << std::endl;
+            BOOST_CHECK(retval == (vectorized<T>{T(24.), T(44.), T(64.)}));
+
+            retval = x3;
+            fma3(retval, x2, x1);
+            std::cout << retval << std::endl;
+            BOOST_CHECK(retval == (vectorized<T>{T(24.), T(44.), T(64.)}));
+        }
+    }
+
+    {
+        // testing the behaviour of vectorized size promotion
+        vectorized<T> x1{T(10.), T(20.), T(30.)};
+        vectorized<T> x2{T(2.), T(2.), T(2.)};
+        vectorized<T> x3{T(4.)};
+        BOOST_CHECK(x1 + x2 == (vectorized<T>{T(12.), T(22.), T(32.)}));
+        BOOST_CHECK(x1 - x2 == (vectorized<T>{T(8.), T(18.), T(28.)}));
+        BOOST_CHECK(x1 * x2 == (vectorized<T>{T(20.), T(40.), T(60.)}));
+        BOOST_CHECK(x1 / x2 == (vectorized<T>{T(5.), T(10.), T(15.)}));
+
+        {
+            auto retval = x1;
+            fma3(retval, x2, x3);
+            std::cout << retval << std::endl;
+            BOOST_CHECK(retval == (vectorized<T>{T(18.), T(28.), T(38.)}));
+
+            retval = x1;
+            fma3(retval, x3, x2);
+            std::cout << retval << std::endl;
+            BOOST_CHECK(retval == (vectorized<T>{T(18.), T(28.), T(38.)}));
+
+            retval = x2;
+            fma3(retval, x1, x3);
+            std::cout << retval << std::endl;
+            BOOST_CHECK(retval == (vectorized<T>{T(42.), T(82.), T(122.)}));
+
+            retval = x2;
+            fma3(retval, x3, x1);
+            std::cout << retval << std::endl;
+            BOOST_CHECK(retval == (vectorized<T>{T(42.), T(82.), T(122.)}));
+
+            retval = x3;
+            fma3(retval, x1, x2);
+            std::cout << retval << std::endl;
+            BOOST_CHECK(retval == (vectorized<T>{T(24.), T(44.), T(64.)}));
+
+            retval = x3;
+            fma3(retval, x2, x1);
+            std::cout << retval << std::endl;
+            BOOST_CHECK(retval == (vectorized<T>{T(24.), T(44.), T(64.)}));
+        }
+    }
+
+{
+        // testing the behaviour of vectorized size promotion
+        vectorized<T> x1{T(10.), T(20.), T(30.)};
+        vectorized<T> x2{T(2.), T(2.), T(2.)};
+        vectorized<T> x3{T(4.), T(4.), T(4.)};
+        BOOST_CHECK(x1 + x2 == (vectorized<T>{T(12.), T(22.), T(32.)}));
+        BOOST_CHECK(x1 - x2 == (vectorized<T>{T(8.), T(18.), T(28.)}));
+        BOOST_CHECK(x1 * x2 == (vectorized<T>{T(20.), T(40.), T(60.)}));
+        BOOST_CHECK(x1 / x2 == (vectorized<T>{T(5.), T(10.), T(15.)}));
+
+        {
+            auto retval = x1;
+            fma3(retval, x2, x3);
+            std::cout << retval << std::endl;
+            BOOST_CHECK(retval == (vectorized<T>{T(18.), T(28.), T(38.)}));
+
+            retval = x1;
+            fma3(retval, x3, x2);
+            std::cout << retval << std::endl;
+            BOOST_CHECK(retval == (vectorized<T>{T(18.), T(28.), T(38.)}));
+
+            retval = x2;
+            fma3(retval, x1, x3);
+            std::cout << retval << std::endl;
+            BOOST_CHECK(retval == (vectorized<T>{T(42.), T(82.), T(122.)}));
+
+            retval = x2;
+            fma3(retval, x3, x1);
+            std::cout << retval << std::endl;
+            BOOST_CHECK(retval == (vectorized<T>{T(42.), T(82.), T(122.)}));
+
+            retval = x3;
+            fma3(retval, x1, x2);
+            std::cout << retval << std::endl;
+            BOOST_CHECK(retval == (vectorized<T>{T(24.), T(44.), T(64.)}));
+
+            retval = x3;
+            fma3(retval, x2, x1);
+            std::cout << retval << std::endl;
+            BOOST_CHECK(retval == (vectorized<T>{T(24.), T(44.), T(64.)}));
+        }
+    }
+
 }
 
 BOOST_AUTO_TEST_CASE(construction)