Feat/fp16/math4

include/eve/module/core/regular/impl/rec.hpp

@@ -28,7 +28,11 @@ namespace eve::detail
   constexpr T  rec_(EVE_REQUIRES(cpu_), O const& o, T const& a) noexcept
   requires(!O::contains(mod))
   {
-    if constexpr( floating_value<T> )
+    if constexpr(O::contains(widen))
+    {
+      return add[o.drop(widen)](upgrade(a));
+    }
+    else if constexpr( floating_value<T> )
     {
       if constexpr( O::contains(lower) || O::contains(upper))
       {

include/eve/module/core/regular/manhattan.hpp

@@ -107,11 +107,14 @@ namespace eve
     EVE_FORCEINLINE constexpr auto
     manhattan_(EVE_REQUIRES(cpu_), O const & o, Ts... args) noexcept
     {
+      using r_t = common_value_t<Ts...>;
+      using e_t = element_type_t<r_t>;
       if constexpr(O::contains(widen))
         return manhattan[o.drop(widen)](upgrade(args)...);
+      else if constexpr(std::same_as<e_t, eve::float16_t>)
+        return eve::detail::apply_fp16_as_fp32(eve::manhattan[o], args...);
       else
       {
-        using r_t = common_value_t<Ts...>;
         auto l_abs = [](){
           if constexpr(integral_value<r_t> && O::contains(saturated))
           return eve::abs[saturated];

include/eve/module/core/regular/rec.hpp

@@ -15,11 +15,11 @@
 namespace eve
 {
 template<typename Options>
-struct rec_t : elementwise_callable<rec_t, Options, raw_option, pedantic_option,
+struct rec_t : elementwise_callable<rec_t, Options, raw_option, pedantic_option, widen_option,
                                     lower_option, upper_option, strict_option, mod_option>
 {
   template<eve::value T>
-  constexpr EVE_FORCEINLINE T operator()(T v) const noexcept
+  constexpr EVE_FORCEINLINE upgrade_if_t<Options,T>  operator()(T v) const noexcept
   { return EVE_DISPATCH_CALL(v); }
 
   EVE_CALLABLE_OBJECT(rec_t, rec_);
@@ -55,6 +55,7 @@ struct rec_t : elementwise_callable<rec_t, Options, raw_option, pedantic_option,
 //!      constexpr auto rec[lower](floating_value auto x)          noexcept; // 5
 //!      constexpr auto rec[upper](floating_value auto x)          noexcept; // 6
 //!      constexpr auto rec[mod = p](floating_value auto x)        noexcept; // 7
+//!      constexpr auto rec[widen](floating_value auto x)          noexcept; // 8
 //!   }
 //!   @endcode
 //!
@@ -78,7 +79,8 @@ struct rec_t : elementwise_callable<rec_t, Options, raw_option, pedantic_option,
 //!      7. Computes the result in modular arithmetic. the parameters must be flint positive
 //!         and less than the modulus. The modulus itself must be less than maxflint. Note that
 //!         mul[mod = p](a, rec[mod = p](a)) is the gcd of p and a (1 iff a and p are coprime)
-//!
+//!      8. The inverse is computed in the double sized element type (if available).
+//!         This decorator has no effect on double and  64 bits integrals.//!
 //!  @note
 //!     For [integral value](@ref eve::integral_value) `rec(x)` is equivalent to:
 //!       * If x==1 or x==-1, x is returned.

include/eve/module/math/regular/agd.hpp

@@ -83,9 +83,12 @@ namespace eve
   namespace detail
   {
     template<typename T, callable_options O>
-    constexpr EVE_FORCEINLINE T agd_(EVE_REQUIRES(cpu_), O const&, T const& x)
+    constexpr EVE_FORCEINLINE T agd_(EVE_REQUIRES(cpu_), O const& o, T const& x)
     {
-      return 2*atanh(tan(x*half(as(x))));
+      if constexpr(std::same_as<eve::element_type_t<T>, eve::float16_t>)
+        return eve::detail::apply_fp16_as_fp32(eve::agd[o], x);
+      else
+        return 2*atanh(tan(x*half(as(x))));
     }
   }
 }

include/eve/module/math/regular/cbrt.hpp

@@ -74,40 +74,45 @@ namespace eve
   namespace detail
   {
     template<eve::floating_value T, callable_options O>
-    EVE_FORCEINLINE constexpr auto cbrt_(EVE_REQUIRES(cpu_), O const& , T x) noexcept
+    EVE_FORCEINLINE constexpr auto cbrt_(EVE_REQUIRES(cpu_), O const& o, T x) noexcept
     {
-      using vt_t = element_type_t<T>;
-      auto test0 = is_eqz(x) || is_not_finite(x);
-      if constexpr( eve::scalar_value<T> )
-        if( test0 ) return x;
-      constexpr vt_t factor[5] = {0.6299605249474365823836,
-                                  0.793700525984099737376,
-                                  1.0,
-                                  1.2599210498948731647672,
-                                  1.587401051968199474751};
-      auto           ax        = eve::abs(x);
-      auto           test      = is_less(eve::abs(x), T(100) * smallestposval(eve::as<T>()));
-      ax                       = ldexp[test](ax, 54);
-      /* Reduce x.  xm now is an range  [0.5, 1.0].  */
-      auto [xm, xe] = ifrexp[raw](ax);
-      T u;
-      if constexpr( std::is_same_v<vt_t, double> )
-        u = eve::reverse_horner(xm, T(0x1.6b69cba168ff2p-2), T(0x1.8218dde9028b4p+0), T(-0x1.0eb8277cd8d5dp+1)
-                                , T(0x1.39350adad51ecp+1), T(-0x1.d5ae6cfa20f0cp+0)
-                                , T(0x1.91e2a6fe7e984p-1), T(-0x1.29801e893366dp-3));
-      else if constexpr( std::is_same_v<vt_t, float> )
-        u = eve::reverse_horner(xm, T(0x1.f87bc4p-2f), T(0x1.6527f4p-1f), T(-0x1.88324ap-3f));
-      auto t2 = sqr(u) * u;
-      u *= fma(xm, T(2), t2) / fma(T(2), t2, xm);
+      if constexpr(std::same_as<eve::element_type_t<T>, eve::float16_t>)
+        return eve::detail::apply_fp16_as_fp32(eve::cbrt[o], x);
+      else
+      {
+        using vt_t = element_type_t<T>;
+        auto test0 = is_eqz(x) || is_not_finite(x);
+        if constexpr( eve::scalar_value<T> )
+          if( test0 ) return x;
+        constexpr vt_t factor[5] = {0.6299605249474365823836,
+                                    0.793700525984099737376,
+                                    1.0,
+                                    1.2599210498948731647672,
+                                    1.587401051968199474751};
+        auto           ax        = eve::abs(x);
+        auto           test      = is_less(eve::abs(x), T(100) * smallestposval(eve::as<T>()));
+        ax                       = ldexp[test](ax, 54);
+        /* Reduce x.  xm now is an range  [0.5, 1.0].  */
+        auto [xm, xe] = ifrexp[raw](ax);
+        T u;
+        if constexpr( std::is_same_v<vt_t, double> )
+          u = eve::reverse_horner(xm, T(0x1.6b69cba168ff2p-2), T(0x1.8218dde9028b4p+0), T(-0x1.0eb8277cd8d5dp+1)
+                                 , T(0x1.39350adad51ecp+1), T(-0x1.d5ae6cfa20f0cp+0)
+                                 , T(0x1.91e2a6fe7e984p-1), T(-0x1.29801e893366dp-3));
+        else if constexpr( std::is_same_v<vt_t, float> )
+          u = eve::reverse_horner(xm, T(0x1.f87bc4p-2f), T(0x1.6527f4p-1f), T(-0x1.88324ap-3f));
+        auto t2 = sqr(u) * u;
+        u *= fma(xm, T(2), t2) / fma(T(2), t2, xm);
 
-      if constexpr( eve::scalar_value<T> ) u *= factor[2 + xe % 3];
-      else                            u *= gather(&factor[0], 2 + xe - (xe / 3) * 3);
-      u = minus[is_ltz(x)](u);
-      if constexpr( eve::scalar_value<T> ) xe = add[test](int(xe) / 3, -18);
-      else                            xe = add[test](xe / 3, -18);
-      auto z = ldexp(u, xe);
-      if constexpr( eve::scalar_value<T> ) return z;
-      else                            return if_else(test0, x, z);
+        if constexpr( eve::scalar_value<T> ) u *= factor[2 + xe % 3];
+        else                            u *= gather(&factor[0], 2 + xe - (xe / 3) * 3);
+        u = minus[is_ltz(x)](u);
+        if constexpr( eve::scalar_value<T> ) xe = add[test](int(xe) / 3, -18);
+        else                            xe = add[test](xe / 3, -18);
+        auto z = ldexp(u, xe);
+        if constexpr( eve::scalar_value<T> ) return z;
+        else                            return if_else(test0, x, z);
+      }
     }
   }
 }

include/eve/module/math/regular/cos.hpp

@@ -99,7 +99,9 @@ namespace eve
     template<typename T, callable_options O>
     constexpr EVE_FORCEINLINE T cos_(EVE_REQUIRES(cpu_), O const& o , T const& a0)
     {
-      if constexpr(O::contains(quarter_circle))
+      if constexpr(std::same_as<eve::element_type_t<T>, eve::float16_t>)
+        return eve::detail::apply_fp16_as_fp32(eve::cos[o], a0);
+      else if constexpr(O::contains(quarter_circle))
       {
         auto x2          = sqr(a0);
         auto x2nlepi2_16 = is_not_less_equal(x2, pi2o_16[upper](as(a0)));

include/eve/module/math/regular/cosd.hpp

@@ -85,9 +85,11 @@ namespace eve
   namespace detail
   {
     template<typename T, callable_options O>
-    constexpr EVE_FORCEINLINE T cosd_(EVE_REQUIRES(cpu_), O const&, T const& a0)
+    constexpr EVE_FORCEINLINE T cosd_(EVE_REQUIRES(cpu_), O const& o, T const& a0)
     {
-      if constexpr(O::contains(quarter_circle))
+      if constexpr(std::same_as<eve::element_type_t<T>, eve::float16_t>)
+        return eve::detail::apply_fp16_as_fp32(eve::cosd[o], a0);
+      else if constexpr(O::contains(quarter_circle))
       {
         return eve::cospi[quarter_circle](div_180(a0));
       }

include/eve/module/math/regular/cosh.hpp

@@ -83,47 +83,51 @@ namespace eve
   namespace detail
   {
     template<typename T, callable_options O>
-    constexpr EVE_FORCEINLINE T cosh_(EVE_REQUIRES(cpu_), O const&, T const& a0)
+    constexpr EVE_FORCEINLINE T cosh_(EVE_REQUIRES(cpu_), O const& o, T const& a0)
     {
-      //////////////////////////////////////////////////////////////////////////////
+      if constexpr( scalar_value<T> )
+      {
+        if( is_eqz(a0) ) return one(eve::as(a0));
+      }      //////////////////////////////////////////////////////////////////////////////
       // if x = abs(a0) according x < Threshold e =  exp(x) or exp(x/2) is
       // respectively computed
       // *  in the first case cosh (e+rec[pedantic](e))/2
       // *  in the second     cosh is (e/2)*e (avoiding undue overflow)
       // Threshold is maxlog - Log_2
       //////////////////////////////////////////////////////////////////////////////
-      if constexpr( scalar_value<T> )
-      {
-        if( is_eqz(a0) ) return one(eve::as(a0));
-      }
-      T ovflimitmln2 = maxlog(as(a0))-log_2(as(a0));
-      auto x = eve::abs(a0);
-      if constexpr( scalar_value<T> )
+      if constexpr(std::same_as<eve::element_type_t<T>, eve::float16_t>)
+        return eve::detail::apply_fp16_as_fp32(eve::cosh[o], a0);
+      else
       {
-        if(is_not_finite(x)) return x;
-        else if( x >= ovflimitmln2 )
+        T ovflimitmln2 = maxlog(as(a0))-log_2(as(a0));
+        auto x = eve::abs(a0);
+        if constexpr( scalar_value<T> )
+        {
+          if(is_not_finite(x)) return x;
+          else if( x >= ovflimitmln2 )
+          {
+            auto w = exp(x * half(eve::as<T>()));
+            auto t = half(eve::as<T>()) * w;
+            t *= w;
+            return t;
+          }
+          auto t = exp(x);
+          return (x > 22) ? t * half(eve::as<T>()) : average(t, rec[pedantic](t));
+        }
+        else
         {
+          auto t    = exp(x);
+          auto invt = if_else(x > 22, eve::zero, rec[pedantic](t));
+          auto c    = average(t, invt);
+          auto test = x < ovflimitmln2;
+          if( eve::all(test) ) return c;
           auto w = exp(x * half(eve::as<T>()));
-          auto t = half(eve::as<T>()) * w;
+          t      = half(eve::as<T>()) * w;
           t *= w;
-          return t;
-        }
-        auto t = exp(x);
-        return (x > 22) ? t * half(eve::as<T>()) : average(t, rec[pedantic](t));
-      }
-      else
-      {
-        auto t    = exp(x);
-        auto invt = if_else(x > 22, eve::zero, rec[pedantic](t));
-        auto c    = average(t, invt);
-        auto test = x < ovflimitmln2;
-        if( eve::all(test) ) return c;
-        auto w = exp(x * half(eve::as<T>()));
-        t      = half(eve::as<T>()) * w;
-        t *= w;
 
-        c = if_else(test, c, t);
-        return c;
+          c = if_else(test, c, t);
+          return c;
+        }
       }
     }
   }

include/eve/module/math/regular/cospi.hpp

@@ -86,9 +86,11 @@ namespace eve
   namespace detail
   {
     template<typename T, callable_options O>
-    constexpr EVE_FORCEINLINE T cospi_(EVE_REQUIRES(cpu_), O const&, T const& a0)
+    constexpr EVE_FORCEINLINE T cospi_(EVE_REQUIRES(cpu_), O const& o, T const& a0)
     {
-      if constexpr(O::contains(quarter_circle))
+      if constexpr(std::same_as<eve::element_type_t<T>, eve::float16_t>)
+        return eve::detail::apply_fp16_as_fp32(eve::cospi[o], a0);
+      else if constexpr(O::contains(quarter_circle))
       {
         return eve::cos[quarter_circle](a0*pi(eve::as<T>()));
       }