First salvo of eve::complex math functions

Implements the same subset of complex function than provided by libstd

examples/tutorial/CMakeLists.txt

@@ -22,3 +22,4 @@ make_unit( "examples.tutorial" load_between.cpp     )
 make_unit( "examples.tutorial" load_ignore.cpp      )
 make_unit( "examples.tutorial" load_ignore_both.cpp )
 make_unit( "examples.tutorial" square_or_diff.cpp   )
+make_unit( "examples.tutorial" mandelbrot.cpp   )

examples/tutorial/mandelbrot.cpp

@@ -0,0 +1,137 @@
+//! [mandelbrot-all]
+#include <chrono>
+#include <iostream>
+#include <numeric>
+#include <vector>
+#include <eve/module/core.hpp>
+#include <eve/module/complex.hpp>
+
+typedef float T;
+using wide_t = eve::wide<T>;
+using wide_i = eve::wide<int>;
+using wide_l = eve::as_logical_t<wide_t>;
+using cwide_t= eve::wide<eve::complex<T>>;
+
+struct mandelbrot
+{
+  static constexpr auto algt = eve::alignment_v<T>;
+  int size, max_iter;
+  T x_min, x_max, y_min, y_max, x_range, y_range;
+  std::vector<int> iterations;
+
+  mandelbrot(int size_, int max_iter_)
+    : size(size_)
+    , max_iter(max_iter_)
+  {
+    x_min   = -2.5;
+    x_max   = 1;
+    y_min   = -1;
+    y_max   = 1;
+    x_range = x_max - x_min;
+    y_range = y_max - y_min;
+    iterations.resize(size * size);
+  }
+
+  //! [mandelbrot-complex-std]
+  void evaluate_complex_std()
+  {
+    using c_t =  std::complex<T>;
+    for (int i = 0; i < size; ++i) {
+      T  x0(T(i) / T(size) * x_range + x_min);
+      for (int j = 0; j < size; ++j) {
+        int iteration = 0;
+        c_t z0(T(j) / T(size) * y_range + y_min, x0);
+        c_t z(0);
+        T n2 =  std::norm(z);
+        while (n2 < 4 && iteration < max_iter) {
+          n2 = std::norm(z);
+          z =  z*z + z0;
+          ++iteration;
+        }
+        iterations[j + i * size] = iteration;
+      }
+    }
+  }
+  //! [mandelbrot-complex-std]
+
+  //! [mandelbrot-complex-scalar]
+  void evaluate_complex_scalar()
+  {
+    using c_t =  eve::complex<T>;
+    for (int i = 0; i < size; ++i) {
+      c_t z0(T(i) / T(size) * x_range + x_min);
+      for (int j = 0; j < size; ++j) {
+        int iteration = 0;
+        eve::imag(z0) = T(j) / T(size) * y_range + y_min;
+        c_t z(0);
+        T n2 =  eve::sqr_abs(z);
+        while (n2 < 4 && iteration < max_iter) {
+          n2 = eve::sqr_abs(z);
+          z =  eve::sqr(z) + z0;
+          ++iteration;
+        }
+        iterations[j + i * size] = iteration;
+      }
+    }
+  }
+  //! [mandelbrot-complex-scalar]
+
+
+  //! [mandelbrot-complex-simd]
+  void evaluate_complex_simd()
+  {
+    wide_t step([](auto n, auto){return n; }); // produce a vector containing {0, 1, ..., wide_t::static_size-1}
+    for (int i = 0; i < size; ++i) {
+      cwide_t z0{T(i) / T(size) * x_range + x_min};
+      wide_t fac{y_range / T(size)};
+      wide_t y_min_t{y_min};
+      for (int j = 0; j < size; j += eve::cardinal_v<wide_t>) {
+        int iteration = 0;
+        eve::imag(z0) = eve::fma(step + j, fac, y_min_t);
+
+        cwide_t z{0};
+        wide_i iter{0};
+        wide_l mask;
+        do {
+          wide_t n2 = eve::sqr_abs(z);
+          z = eve::sqr(z) + z0;
+          mask = n2 < 4;
+          ++iteration;
+          iter = eve::inc[mask](iter);
+        } while (eve::any(mask) && iteration < max_iter);
+        eve::store(iter, &iterations[j + i * size]);
+      }
+    }
+  }
+  //! [mandelbrot-complex-simd]
+};
+
+int main()
+{
+  namespace chr = std::chrono;
+  using hrc     = chr::high_resolution_clock;
+
+  int size          = 256;
+  int max_iteration = 100;
+  mandelbrot image(size, max_iteration);
+
+  auto t0 = hrc::now();
+  image.evaluate_complex_std();
+  auto t1 = hrc::now();
+  std::cout << " complex std " << chr::duration_cast<chr::milliseconds>(t1 - t0).count() << std::endl;
+
+  t0 = hrc::now();
+  image.evaluate_complex_scalar();
+  t1 = hrc::now();
+  std::cout << " complex scalar " << chr::duration_cast<chr::milliseconds>(t1 - t0).count() << std::endl;
+
+  t0 = hrc::now();
+  image.evaluate_complex_simd();
+  t1 = hrc::now();
+  std::cout << " complex simd " << chr::duration_cast<chr::milliseconds>(t1 - t0).count() << std::endl;
+
+}
+//! [mandelbrot-all]
+
+// This code can be compiled using
+// g++ mandelbrot.cpp -mavx2 -mfma -std=c++20 -O3 -DNDEBUG -o mandelbrot -I/path_to_eve/include

include/eve/module/complex/regular/complex.hpp

@@ -11,3 +11,6 @@
 #include <eve/module/complex/regular/real.hpp>
 #include <eve/module/complex/regular/i.hpp>
 #include <eve/module/complex/regular/imag.hpp>
+#include <eve/module/complex/regular/exp_i.hpp>
+#include <eve/module/complex/regular/exp_ipi.hpp>
+#include <eve/module/complex/regular/polar.hpp>

include/eve/module/complex/regular/detail/acos.hpp

@@ -0,0 +1,163 @@
+//==================================================================================================
+/*
+  EVE - Expressive Vector Engine
+  Copyright : EVE Contributors & Maintainers
+  SPDX-License-Identifier: MIT
+*/
+//==================================================================================================
+#pragma once
+
+#include <eve/detail/overload.hpp>
+#include <eve/module/core.hpp>
+#include <eve/module/math.hpp>
+#include <eve/module/complex.hpp>
+#include <eve/module/complex/regular/traits.hpp>
+
+namespace eve
+{
+
+  namespace detail
+  {
+    template<typename Z>
+    EVE_FORCEINLINE auto complex_unary_dispatch( eve::tag::acos_, Z const& a0) noexcept
+    {
+      // This implementation is a simd transcription and adaptation of the boost_math code
+      // which itself is a transcription of the pseudo-code in:
+      //
+      // "Implementing the complex Arcsine and Arccosine Functions using Exception Handling."
+      // T E Hull, Thomas F Fairgrieve and Ping Tak Peter Tang.
+      // ACM Transactions on Mathematical Software, Vol 23, No 3, Sept 1997.
+      //
+      auto [a0r, a0i] = a0;
+      using rtype = decltype(a0r);
+      const rtype a_crossover(1.5);
+      const rtype b_crossover(0.6417);
+      auto  ltzra0 = is_ltz(a0r);
+      auto  gtzia0 = is_gtz(a0i);
+      //
+      // Begin by insuring a0r >= 0 and imag(a0) >= 0 :
+      //
+      rtype x = eve::abs(a0r);
+      rtype y = eve::abs(a0i);
+      rtype proper_real = eve::acos(x);
+      auto lexone = (x <= one(as(x)));
+      auto is_proper_real = logical_and(is_real(a0), lexone);
+
+      auto s_min = eve::sqrtsmallestposval(as(x))*4;
+      auto s_max = eve::sqrtvalmax(as(x))/8;
+      rtype xp1 = eve::inc(x);
+      rtype xm1 = eve::dec(x);
+      auto not_in_safe_zone = (((x > s_max) || (x < s_min)) || ((y > s_max) || (y < s_min)));
+      //compute for safe zone
+      rtype r, i;
+      rtype yy = eve::sqr(y);
+      rtype tr = eve::sqrt(sqr(xp1) + yy); //hypot for pedantic ?
+      rtype ts = eve::sqrt(sqr(xm1) + yy); //hypot for pedantic ?
+      rtype a = eve::average(tr, ts);
+      rtype b = x/a;
+      //compute r for b >  b_crossover
+      rtype apx = a + x;
+      r = if_else(lexone,
+                  eve::atan(eve::sqrt(half(as(x))*apx*(yy/(tr+xp1)+(ts-xm1)))/x),
+                  eve::atan((y*eve::sqrt(half(as(x))*(apx/(tr+xp1)+apx/(ts+xm1))))/x)
+        );
+      // r is computed
+      r = if_else((b <= b_crossover), eve::acos(b), r);
+      //compute am1 temporary for i for a <= a_crossover
+      rtype tmp = yy/(tr+xp1);
+      rtype am1 = if_else(lexone,
+                          eve::average(tmp, yy/(ts-xm1)),
+                          eve::average(tmp, (ts+xm1)));
+      i =  if_else((a <= a_crossover),
+                   eve::log1p(am1 + eve::sqrt(am1*(eve::inc(a)))),
+                   eve::log(a + eve::sqrt(eve::dec(eve::sqr(a))))
+        );
+      // i is computed
+      //compute for exception zone
+      if (eve::any(not_in_safe_zone))
+      {
+        auto zone1 =  (y <= eve::eps(as(x))*eve::abs(xm1));
+        if (eve::any(logical_and(zone1, not_in_safe_zone)))
+        {
+          rtype rr = if_else(lexone, proper_real, zero);
+          rtype ii = if_else(lexone, y/eve::sqrt(-xp1*xm1),
+                             if_else((valmax(as(x))/xp1 > xm1),
+                                     eve::log1p(xm1 + eve::sqrt(xp1*xm1)),
+                                     log_2(as(x)) + eve::log(x)
+                               )
+            );
+          r = if_else(zone1, rr, r);
+          i = if_else(zone1, ii, i);
+        }
+        auto zone2 = (y <= s_min);
+        auto not_treated = logical_notand(zone1, not_in_safe_zone);
+        if (eve::any(logical_and(zone2, not_treated)))
+        {
+          rtype sqrty =  eve::sqrt(y);
+          r = if_else(zone2, sqrty, r);
+          i = if_else(zone2, sqrty, i);
+        }
+        auto zone3 = (dec(eps(as(x))*y) >= x);
+        not_treated = logical_notand(zone2, not_treated);
+        if (eve::any(logical_and(zone3, not_treated)))
+        {
+          r = if_else(zone3,  pio_2(as(x)), r);
+          i = if_else(zone3,  log_2(as(x)) + eve::log(y), i);
+        }
+        auto zone4 = (x > one(as(x)));
+        not_treated = logical_notand(zone3, not_treated);
+        if (eve::any(logical_and(zone4, not_treated)))
+        {
+          r = if_else(zone4, eve::atan(y/x), r);
+          i = if_else(zone4, log_2(as(x)) + eve::log(y) + half(as(x))*eve::log1p(sqr(x/y)), i);
+        }
+        not_treated = logical_notand(zone4, not_treated);
+        if (eve::any(not_treated))
+        {
+          rtype aa = eve::sqrt(eve::inc(sqr(y)));
+          r = if_else(not_treated, pio_2(as(x)), r);
+          i = if_else(not_treated, half(as(x))*eve::log1p(2*y*(y+aa)), i);
+        }
+      }
+      if (eve::any(eve::is_not_finite(a0)))
+      {
+        auto nanx = eve::is_nan(x);
+        auto nany = eve::is_nan(y);
+        auto infx = (x == inf(as(x))) ;
+        auto infy = (y == inf(as(x))) ;
+        if (eve::any(infx))
+        {
+          r = eve::if_else(infx, zero, r);
+          i = eve::if_else(infx, inf(as(x)), i);
+          r = eve::if_else(logical_and(infx, infy), pio_4(as(x)), r);
+          i = eve::if_else(logical_and(infx, infy), inf(as(x)), i);
+
+          r =  if_else(logical_and(infx, nany), y, r);
+          i =  if_else(logical_and(infx, nany), minf(as(x)), i);
+        }
+        if (eve::any(nanx))
+        {
+          r =  if_else(nanx, x, r);
+          i =  if_else(nanx, x, i);
+          i =  if_else(logical_and(nanx, infy), y, i);
+        }
+        auto test = logical_notand(logical_or(infx, nanx), infy);
+        if (eve::any(test))
+        {
+          r = if_else(logical_and(infy, test), pio_2(as(x)), r);
+          i = if_else(logical_and(infy, test), y, i);
+        }
+        test = logical_notand(logical_or(infx, nanx), nany);
+        r = if_else(test,if_else(eve::is_imag(a0), pio_2(as(x)), y), r);
+        i = if_else(test,y,i);
+      }
+      // use proper real results
+      r = if_else(is_proper_real, proper_real,   r);
+      i = if_else(is_proper_real, zero, i);
+      // restore signs
+      r = if_else(ltzra0, pi(as(x))-r, r);
+      i = if_else(gtzia0, -i, i);
+      return Z{r, i};
+    }
+  }
+}

include/eve/module/complex/regular/detail/arithmetic.hpp

@@ -15,6 +15,15 @@
 
 namespace eve::detail
 {
+
+  //==============================================================================================
+  //  Unary functions
+  //==============================================================================================
+  EVE_FORCEINLINE auto complex_unary_dispatch(eve::tag::abs_, auto const& z) noexcept
+  {
+    return eve::hypot(real(z), imag(z));
+  }
+
   EVE_FORCEINLINE auto complex_unary_dispatch(eve::tag::arg_, auto const& z) noexcept
   {
     return eve::atan2(imag(z), real(z) );
@@ -33,21 +42,29 @@ namespace eve::detail
     return if_else(is_infinite(z), Z(inf(real_t{}), copysign(zero(real_t{}), imag(z))), z);
   }
 
-  EVE_FORCEINLINE auto complex_unary_dispatch(eve::tag::abs_, auto const& z) noexcept
-  {
-    return eve::hypot(real(z), imag(z));
-  }
-
   template<typename Z>
   EVE_FORCEINLINE auto complex_unary_dispatch(eve::tag::sqr_, Z const& z) noexcept
   {
     auto [zr, zi] = z;
     return Z{diff_of_prod(zr, zr, zi, zi), 2*zr*zi};
   }
 
-  template<typename Z>
-  EVE_FORCEINLINE auto complex_binary_dispatch(eve::tag::dist_, Z const& x, Z const& y) noexcept
+  EVE_FORCEINLINE auto complex_unary_dispatch(eve::tag::sqr_abs_, auto const& z) noexcept
+  {
+    auto [zr, zi] = z;
+    return sum_of_prod(zr, zr, zi, zi);
+  }
+
+
+  //==============================================================================================
+  //  Binary functions
+  //==============================================================================================
+
+  EVE_FORCEINLINE auto complex_binary_dispatch( eve::tag::average_
+                                              , auto const& z1, auto const& z2
+                                              ) noexcept
   {
-    return eve::abs(x-y);
+    using z_t = decltype(z1+z2);
+    return z_t{eve::average(real(z1), real(z2)), eve::average(imag(z1), imag(z2))};
   }
 }