Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-11-07 10:36:57 

0001 /***************************************************************************
0002  * Copyright (c) Johan Mabille, Sylvain Corlay, Wolf Vollprecht and         *
0003  * Martin Renou                                                             *
0004  * Copyright (c) QuantStack                                                 *
0005  * Copyright (c) Serge Guelton                                              *
0006  *                                                                          *
0007  * Distributed under the terms of the BSD 3-Clause License.                 *
0008  *                                                                          *
0009  * The full license is in the file LICENSE, distributed with this software. *
0010  ****************************************************************************/
0011 
0012 #ifndef XSIMD_COMMON_ARITHMETIC_HPP
0013 #define XSIMD_COMMON_ARITHMETIC_HPP
0014 
0015 #include "../../types/xsimd_batch_constant.hpp"
0016 #include "./xsimd_common_details.hpp"
0017 
0018 #include <complex>
0019 #include <limits>
0020 #include <type_traits>
0021 #include <utility>
0022 
0023 namespace xsimd
0024 {
0025 
0026     namespace kernel
0027     {
0028 
0029         using namespace types;
0030 
0031         // bitwise_lshift
0032         template <class A, class T, class /*=std::enable_if_t<std::is_integral<T>::value>*/>
0033         XSIMD_INLINE batch<T, A> bitwise_lshift(batch<T, A> const& self, batch<T, A> const& other, requires_arch<common>) noexcept
0034         {
0035             return detail::apply([](T x, T y) noexcept
0036                                  { return x << y; },
0037                                  self, other);
0038         }
0039         template <size_t shift, class A, class T, class /*=std::enable_if_t<std::is_integral<T>::value>*/>
0040         XSIMD_INLINE batch<T, A> bitwise_lshift(batch<T, A> const& self, requires_arch<common>) noexcept
0041         {
0042             constexpr auto bits = std::numeric_limits<T>::digits + std::numeric_limits<T>::is_signed;
0043             static_assert(shift < bits, "Shift must be less than the number of bits in T");
0044             return bitwise_lshift(self, shift, A {});
0045         }
0046 
0047         // bitwise_rshift
0048         template <class A, class T, class /*=std::enable_if_t<std::is_integral<T>::value>*/>
0049         XSIMD_INLINE batch<T, A> bitwise_rshift(batch<T, A> const& self, batch<T, A> const& other, requires_arch<common>) noexcept
0050         {
0051             return detail::apply([](T x, T y) noexcept
0052                                  { return x >> y; },
0053                                  self, other);
0054         }
0055         template <size_t shift, class A, class T, class /*=std::enable_if_t<std::is_integral<T>::value>*/>
0056         XSIMD_INLINE batch<T, A> bitwise_rshift(batch<T, A> const& self, requires_arch<common>) noexcept
0057         {
0058             constexpr auto bits = std::numeric_limits<T>::digits + std::numeric_limits<T>::is_signed;
0059             static_assert(shift < bits, "Shift must be less than the number of bits in T");
0060             return bitwise_rshift(self, shift, A {});
0061         }
0062 
0063         // decr
0064         template <class A, class T>
0065         XSIMD_INLINE batch<T, A> decr(batch<T, A> const& self, requires_arch<common>) noexcept
0066         {
0067             return self - T(1);
0068         }
0069 
0070         // decr_if
0071         template <class A, class T, class Mask>
0072         XSIMD_INLINE batch<T, A> decr_if(batch<T, A> const& self, Mask const& mask, requires_arch<common>) noexcept
0073         {
0074             return select(mask, decr(self), self);
0075         }
0076 
0077         // div
0078         template <class A, class T, class = std::enable_if_t<std::is_integral<T>::value>>
0079         XSIMD_INLINE batch<T, A> div(batch<T, A> const& self, batch<T, A> const& other, requires_arch<common>) noexcept
0080         {
0081             return detail::apply([](T x, T y) noexcept -> T
0082                                  { return x / y; },
0083                                  self, other);
0084         }
0085 
0086         // fma
0087         template <class A, class T>
0088         XSIMD_INLINE batch<T, A> fma(batch<T, A> const& x, batch<T, A> const& y, batch<T, A> const& z, requires_arch<common>) noexcept
0089         {
0090             return x * y + z;
0091         }
0092 
0093         template <class A, class T>
0094         XSIMD_INLINE batch<std::complex<T>, A> fma(batch<std::complex<T>, A> const& x, batch<std::complex<T>, A> const& y, batch<std::complex<T>, A> const& z, requires_arch<common>) noexcept
0095         {
0096             auto res_r = fms(x.real(), y.real(), fms(x.imag(), y.imag(), z.real()));
0097             auto res_i = fma(x.real(), y.imag(), fma(x.imag(), y.real(), z.imag()));
0098             return { res_r, res_i };
0099         }
0100 
0101         // fms
0102         template <class A, class T>
0103         XSIMD_INLINE batch<T, A> fms(batch<T, A> const& x, batch<T, A> const& y, batch<T, A> const& z, requires_arch<common>) noexcept
0104         {
0105             return x * y - z;
0106         }
0107 
0108         template <class A, class T>
0109         XSIMD_INLINE batch<std::complex<T>, A> fms(batch<std::complex<T>, A> const& x, batch<std::complex<T>, A> const& y, batch<std::complex<T>, A> const& z, requires_arch<common>) noexcept
0110         {
0111             auto res_r = fms(x.real(), y.real(), fma(x.imag(), y.imag(), z.real()));
0112             auto res_i = fma(x.real(), y.imag(), fms(x.imag(), y.real(), z.imag()));
0113             return { res_r, res_i };
0114         }
0115 
0116         // fnma
0117         template <class A, class T>
0118         XSIMD_INLINE batch<T, A> fnma(batch<T, A> const& x, batch<T, A> const& y, batch<T, A> const& z, requires_arch<common>) noexcept
0119         {
0120             return -x * y + z;
0121         }
0122 
0123         template <class A, class T>
0124         XSIMD_INLINE batch<std::complex<T>, A> fnma(batch<std::complex<T>, A> const& x, batch<std::complex<T>, A> const& y, batch<std::complex<T>, A> const& z, requires_arch<common>) noexcept
0125         {
0126             auto res_r = -fms(x.real(), y.real(), fma(x.imag(), y.imag(), z.real()));
0127             auto res_i = -fma(x.real(), y.imag(), fms(x.imag(), y.real(), z.imag()));
0128             return { res_r, res_i };
0129         }
0130 
0131         // fnms
0132         template <class A, class T>
0133         XSIMD_INLINE batch<T, A> fnms(batch<T, A> const& x, batch<T, A> const& y, batch<T, A> const& z, requires_arch<common>) noexcept
0134         {
0135             return -x * y - z;
0136         }
0137 
0138         template <class A, class T>
0139         XSIMD_INLINE batch<std::complex<T>, A> fnms(batch<std::complex<T>, A> const& x, batch<std::complex<T>, A> const& y, batch<std::complex<T>, A> const& z, requires_arch<common>) noexcept
0140         {
0141             auto res_r = -fms(x.real(), y.real(), fms(x.imag(), y.imag(), z.real()));
0142             auto res_i = -fma(x.real(), y.imag(), fma(x.imag(), y.real(), z.imag()));
0143             return { res_r, res_i };
0144         }
0145 
0146         // fmas
0147         template <class A, class T>
0148         XSIMD_INLINE batch<T, A> fmas(batch<T, A> const& x, batch<T, A> const& y, batch<T, A> const& z, requires_arch<common>) noexcept
0149         {
0150             struct even_lane
0151             {
0152                 static constexpr bool get(unsigned const i, unsigned) noexcept { return (i & 1u) == 0; }
0153             };
0154             const auto mask = make_batch_bool_constant<T, even_lane, A>();
0155             return fma(x, y, select(mask, neg(z), z));
0156         }
0157 
0158         // incr
0159         template <class A, class T>
0160         XSIMD_INLINE batch<T, A> incr(batch<T, A> const& self, requires_arch<common>) noexcept
0161         {
0162             return self + T(1);
0163         }
0164 
0165         // incr_if
0166         template <class A, class T, class Mask>
0167         XSIMD_INLINE batch<T, A> incr_if(batch<T, A> const& self, Mask const& mask, requires_arch<common>) noexcept
0168         {
0169             return select(mask, incr(self), self);
0170         }
0171 
0172         // mul
0173         template <class A, class T, class /*=std::enable_if_t<std::is_integral<T>::value>*/>
0174         XSIMD_INLINE batch<T, A> mul(batch<T, A> const& self, batch<T, A> const& other, requires_arch<common>) noexcept
0175         {
0176             return detail::apply([](T x, T y) noexcept -> T
0177                                  { return x * y; },
0178                                  self, other);
0179         }
0180 
0181         // mul_hi
0182         namespace detail
0183         {
0184             template <class T>
0185             struct mulhi_helper
0186             {
0187                 using wider = std::conditional_t<
0188                     std::is_signed<T>::value,
0189                     std::conditional_t<sizeof(T) == 1, int16_t,
0190                                        std::conditional_t<sizeof(T) == 2, int32_t, int64_t>>,
0191                     std::conditional_t<sizeof(T) == 1, uint16_t,
0192                                        std::conditional_t<sizeof(T) == 2, uint32_t, uint64_t>>>;
0193 
0194                 static XSIMD_INLINE T compute(T x, T y) noexcept
0195                 {
0196                     constexpr int shift = 8 * sizeof(T);
0197                     return static_cast<T>((static_cast<wider>(x) * static_cast<wider>(y)) >> shift);
0198                 }
0199             };
0200 
0201             // 64-bit unsigned software mul_hi via 32-bit splits
0202             XSIMD_INLINE uint64_t mulhi_u64(uint64_t x, uint64_t y) noexcept
0203             {
0204 #if defined(__SIZEOF_INT128__)
0205                 return static_cast<uint64_t>((static_cast<unsigned __int128>(x) * static_cast<unsigned __int128>(y)) >> 64);
0206 #else
0207                 uint64_t xl = x & 0xffffffffULL;
0208                 uint64_t xh = x >> 32;
0209                 uint64_t yl = y & 0xffffffffULL;
0210                 uint64_t yh = y >> 32;
0211                 uint64_t ll = xl * yl;
0212                 uint64_t lh = xl * yh;
0213                 uint64_t hl = xh * yl;
0214                 uint64_t hh = xh * yh;
0215                 uint64_t mid = (ll >> 32) + (lh & 0xffffffffULL) + (hl & 0xffffffffULL);
0216                 return hh + (lh >> 32) + (hl >> 32) + (mid >> 32);
0217 #endif
0218             }
0219 
0220             XSIMD_INLINE int64_t mulhi_i64(int64_t x, int64_t y) noexcept
0221             {
0222 #if defined(__SIZEOF_INT128__)
0223                 return static_cast<int64_t>((static_cast<__int128>(x) * static_cast<__int128>(y)) >> 64);
0224 #else
0225                 uint64_t uhi = mulhi_u64(static_cast<uint64_t>(x), static_cast<uint64_t>(y));
0226                 if (x < 0)
0227                     uhi -= static_cast<uint64_t>(y);
0228                 if (y < 0)
0229                     uhi -= static_cast<uint64_t>(x);
0230                 return static_cast<int64_t>(uhi);
0231 #endif
0232             }
0233 
0234             template <>
0235             struct mulhi_helper<uint64_t>
0236             {
0237                 static XSIMD_INLINE uint64_t compute(uint64_t x, uint64_t y) noexcept { return mulhi_u64(x, y); }
0238             };
0239 
0240             template <>
0241             struct mulhi_helper<int64_t>
0242             {
0243                 static XSIMD_INLINE int64_t compute(int64_t x, int64_t y) noexcept { return mulhi_i64(x, y); }
0244             };
0245 
0246             // 64x64 unsigned mul_hi via 32x32->64 widening mul (WMul wraps _mm*_mul_epu32).
0247             template <class A, class WMul>
0248             XSIMD_INLINE batch<uint64_t, A> mulhi_u64_core(batch<uint64_t, A> const& x,
0249                                                            batch<uint64_t, A> const& y,
0250                                                            WMul mul_epu32) noexcept
0251             {
0252                 using B = batch<uint64_t, A>;
0253                 const B mask(uint64_t(0xffffffffULL));
0254                 // mul_epu32 uses only the low 32 bits, so low operands need no mask.
0255                 B xh = x >> 32;
0256                 B yh = y >> 32;
0257                 B ll = mul_epu32(x, y);
0258                 B lh = mul_epu32(x, yh);
0259                 B hl = mul_epu32(xh, y);
0260                 B hh = mul_epu32(xh, yh);
0261                 B mid = (ll >> 32) + (lh & mask) + (hl & mask);
0262                 return hh + (lh >> 32) + (hl >> 32) + (mid >> 32);
0263             }
0264 
0265             // Signed variant: unsigned core + sign fixup via arithmetic shift-by-63.
0266             template <class A, class WMul>
0267             XSIMD_INLINE batch<int64_t, A> mulhi_i64_core(batch<int64_t, A> const& x,
0268                                                           batch<int64_t, A> const& y,
0269                                                           WMul mul_epu32) noexcept
0270             {
0271                 auto ux = ::xsimd::bitwise_cast<uint64_t>(x);
0272                 auto uy = ::xsimd::bitwise_cast<uint64_t>(y);
0273                 auto uhi = mulhi_u64_core<A>(ux, uy, mul_epu32);
0274                 auto sa = ::xsimd::bitwise_cast<uint64_t>(x >> 63);
0275                 auto sb = ::xsimd::bitwise_cast<uint64_t>(y >> 63);
0276                 return ::xsimd::bitwise_cast<int64_t>(uhi - (uy & sa) - (ux & sb));
0277             }
0278 
0279             // Fused mul_hilo: both halves from one set of partials. Returns { hi, lo }.
0280             template <class A, class WMul>
0281             XSIMD_INLINE std::pair<batch<uint64_t, A>, batch<uint64_t, A>>
0282             mulhilo_u64_core(batch<uint64_t, A> const& x,
0283                              batch<uint64_t, A> const& y,
0284                              WMul mul_epu32) noexcept
0285             {
0286                 using B = batch<uint64_t, A>;
0287                 const B mask(uint64_t(0xffffffffULL));
0288                 B xh = x >> 32;
0289                 B yh = y >> 32;
0290                 B ll = mul_epu32(x, y);
0291                 B lh = mul_epu32(x, yh);
0292                 B hl = mul_epu32(xh, y);
0293                 B hh = mul_epu32(xh, yh);
0294                 B mid = (ll >> 32) + (lh & mask) + (hl & mask);
0295                 B hi = hh + (lh >> 32) + (hl >> 32) + (mid >> 32);
0296                 B lo = (ll & mask) | (mid << 32);
0297                 return { hi, lo };
0298             }
0299         }
0300 
0301         template <class A, class T, class /*=std::enable_if_t<std::is_integral<T>::value>*/>
0302         XSIMD_INLINE batch<T, A> mul_hi(batch<T, A> const& self, batch<T, A> const& other, requires_arch<common>) noexcept
0303         {
0304             return detail::apply([](T x, T y) noexcept -> T
0305                                  { return detail::mulhi_helper<T>::compute(x, y); },
0306                                  self, other);
0307         }
0308 
0309         // mul_hilo
0310         template <class A, class T, class /*=std::enable_if_t<std::is_integral<T>::value>*/>
0311         XSIMD_INLINE std::pair<batch<T, A>, batch<T, A>>
0312         mul_hilo(batch<T, A> const& self, batch<T, A> const& other, requires_arch<common>) noexcept
0313         {
0314             return std::pair<batch<T, A>, batch<T, A>> { mul_hi<A>(self, other, A {}), self * other };
0315         }
0316 
0317         // Signed 64-bit mul_hilo: unsigned path + sign fixup on hi (lo is sign-invariant).
0318         template <class A>
0319         XSIMD_INLINE std::pair<batch<int64_t, A>, batch<int64_t, A>>
0320         mul_hilo(batch<int64_t, A> const& self, batch<int64_t, A> const& other, requires_arch<common>) noexcept
0321         {
0322             auto ux = ::xsimd::bitwise_cast<uint64_t>(self);
0323             auto uy = ::xsimd::bitwise_cast<uint64_t>(other);
0324             auto hilo = mul_hilo<A>(ux, uy, A {});
0325             auto sa = ::xsimd::bitwise_cast<uint64_t>(self >> 63);
0326             auto sb = ::xsimd::bitwise_cast<uint64_t>(other >> 63);
0327             auto hi = hilo.first - (uy & sa) - (ux & sb);
0328             return { ::xsimd::bitwise_cast<int64_t>(hi),
0329                      ::xsimd::bitwise_cast<int64_t>(hilo.second) };
0330         }
0331 
0332         // rotl
0333         template <class A, class T, class STy>
0334         XSIMD_INLINE batch<T, A> rotl(batch<T, A> const& self, STy other, requires_arch<common>) noexcept
0335         {
0336             constexpr auto bits = std::numeric_limits<T>::digits + std::numeric_limits<T>::is_signed;
0337             return (self << other) | (self >> (bits - other));
0338         }
0339         template <size_t count, class A, class T>
0340         XSIMD_INLINE batch<T, A> rotl(batch<T, A> const& self, requires_arch<common>) noexcept
0341         {
0342             constexpr auto bits = std::numeric_limits<T>::digits + std::numeric_limits<T>::is_signed;
0343             static_assert(count < bits, "Count amount must be less than the number of bits in T");
0344             return bitwise_lshift<count>(self) | bitwise_rshift<bits - count>(self);
0345         }
0346 
0347         // rotr
0348         template <class A, class T, class STy>
0349         XSIMD_INLINE batch<T, A> rotr(batch<T, A> const& self, STy other, requires_arch<common>) noexcept
0350         {
0351             constexpr auto bits = std::numeric_limits<T>::digits + std::numeric_limits<T>::is_signed;
0352             return (self >> other) | (self << (bits - other));
0353         }
0354         template <size_t count, class A, class T>
0355         XSIMD_INLINE batch<T, A> rotr(batch<T, A> const& self, requires_arch<common>) noexcept
0356         {
0357             constexpr auto bits = std::numeric_limits<T>::digits + std::numeric_limits<T>::is_signed;
0358             static_assert(count < bits, "Count must be less than the number of bits in T");
0359             return bitwise_rshift<count>(self) | bitwise_lshift<bits - count>(self);
0360         }
0361 
0362         // sadd
0363         template <class A>
0364         XSIMD_INLINE batch<float, A> sadd(batch<float, A> const& self, batch<float, A> const& other, requires_arch<common>) noexcept
0365         {
0366             return add(self, other); // no saturated arithmetic on floating point numbers
0367         }
0368         template <class A, class T, class /*=std::enable_if_t<std::is_integral<T>::value>*/>
0369         XSIMD_INLINE batch<T, A> sadd(batch<T, A> const& self, batch<T, A> const& other, requires_arch<common>) noexcept
0370         {
0371             if (std::is_signed<T>::value)
0372             {
0373                 auto self_pos_branch = min(std::numeric_limits<T>::max() - other, self);
0374                 auto self_neg_branch = max(std::numeric_limits<T>::min() - other, self);
0375                 return other + select(other >= 0, self_pos_branch, self_neg_branch);
0376             }
0377             else
0378             {
0379                 const auto diffmax = std::numeric_limits<T>::max() - self;
0380                 const auto mindiff = min(diffmax, other);
0381                 return self + mindiff;
0382             }
0383         }
0384         template <class A>
0385         XSIMD_INLINE batch<double, A> sadd(batch<double, A> const& self, batch<double, A> const& other, requires_arch<common>) noexcept
0386         {
0387             return add(self, other); // no saturated arithmetic on floating point numbers
0388         }
0389 
0390         // ssub
0391         template <class A>
0392         XSIMD_INLINE batch<float, A> ssub(batch<float, A> const& self, batch<float, A> const& other, requires_arch<common>) noexcept
0393         {
0394             return sub(self, other); // no saturated arithmetic on floating point numbers
0395         }
0396         template <class A, class T, class /*=std::enable_if_t<std::is_integral<T>::value>*/>
0397         XSIMD_INLINE batch<T, A> ssub(batch<T, A> const& self, batch<T, A> const& other, requires_arch<common>) noexcept
0398         {
0399             if (std::is_signed<T>::value)
0400             {
0401                 return sadd(self, -other);
0402             }
0403             else
0404             {
0405                 const auto diff = min(self, other);
0406                 return self - diff;
0407             }
0408         }
0409         template <class A>
0410         XSIMD_INLINE batch<double, A> ssub(batch<double, A> const& self, batch<double, A> const& other, requires_arch<common>) noexcept
0411         {
0412             return sub(self, other); // no saturated arithmetic on floating point numbers
0413         }
0414 
0415     }
0416 
0417 }
0418 
0419 #endif
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!