Fix most issues

sbastrakov · sbastrakov · commit 8ea8d96011a1 · 2021-11-18T15:53:06.000Z
diff --git a/include/alpaka/math/Complex.hpp b/include/alpaka/math/Complex.hpp
@@ -256,10 +256,9 @@ namespace alpaka
     template<typename T>
     ALPAKA_FN_HOST_ACC Complex<T> operator/(Complex<T> const& lhs, Complex<T> const& rhs)
     {
-        return lhs
-            * Complex<T>{
-                rhs.real() / (rhs.real() * rhs.real() + rhs.imag() * rhs.imag()),
-                -rhs.imag() / (rhs.real() * rhs.real() + rhs.imag() * rhs.imag())};
+        return Complex<T>{
+            (lhs.real() * rhs.real() + lhs.imag() * rhs.imag()) / (rhs.real() * rhs.real() + rhs.imag() * rhs.imag()),
+            (lhs.imag() * rhs.real() - lhs.real() * rhs.imag()) / (rhs.real() * rhs.real() + rhs.imag() * rhs.imag())};
     }
 
     //! Division of complex and a real number
@@ -273,10 +272,9 @@ namespace alpaka
     template<typename T>
     ALPAKA_FN_HOST_ACC Complex<T> operator/(T const& lhs, Complex<T> const& rhs)
     {
-        return lhs
-            * Complex<T>{
-                rhs.real() / (rhs.real() * rhs.real() + rhs.imag() * rhs.imag()),
-                -rhs.imag() / (rhs.real() * rhs.real() + rhs.imag() * rhs.imag())};
+        return Complex<T>{
+            lhs * rhs.real() / (rhs.real() * rhs.real() + rhs.imag() * rhs.imag()),
+            -lhs * rhs.imag() / (rhs.real() * rhs.real() + rhs.imag() * rhs.imag())};
     }
 
     //! Equality of two complex numbers
diff --git a/include/alpaka/math/MathStdLib.hpp b/include/alpaka/math/MathStdLib.hpp
@@ -54,6 +54,7 @@ namespace alpaka
             , public Atan2StdLib
             , public CbrtStdLib
             , public CeilStdLib
+            , public ConjStdLib
             , public CosStdLib
             , public ErfStdLib
             , public ExpStdLib
diff --git a/include/alpaka/math/MathUniformCudaHipBuiltIn.hpp b/include/alpaka/math/MathUniformCudaHipBuiltIn.hpp
@@ -56,6 +56,7 @@ namespace alpaka
             , public Atan2UniformCudaHipBuiltIn
             , public CbrtUniformCudaHipBuiltIn
             , public CeilUniformCudaHipBuiltIn
+            , public ConjUniformCudaHipBuiltIn
             , public CosUniformCudaHipBuiltIn
             , public ErfUniformCudaHipBuiltIn
             , public ExpUniformCudaHipBuiltIn
diff --git a/include/alpaka/math/arg/ArgUniformCudaHipBuiltIn.hpp b/include/alpaka/math/arg/ArgUniformCudaHipBuiltIn.hpp
@@ -30,7 +30,7 @@ namespace alpaka
 
         namespace traits
         {
-            //! The CUDA/HIP built in arg trait specialization.
+            //! The CUDA/HIP built in arg trait specialization for float and double.
             template<typename TArgument>
             struct Arg<ArgUniformCudaHipBuiltIn, TArgument, std::enable_if_t<std::is_floating_point<TArgument>::value>>
             {
diff --git a/include/alpaka/math/atan/AtanUniformCudaHipBuiltIn.hpp b/include/alpaka/math/atan/AtanUniformCudaHipBuiltIn.hpp
@@ -60,7 +60,7 @@ namespace alpaka
                 {
                     // This holds everywhere, including the branch cuts: atan(z) = -i/2 * ln((i - z) / (i + z))
                     return Complex<T>{0.0, -0.5}
-                    * log(ctx, (Complex<T>{0.0, 1.0} + arg) / (Complex<T>{0.0, 1.0} + arg));
+                    * log(ctx, (Complex<T>{0.0, 1.0} - arg) / (Complex<T>{0.0, 1.0} + arg));
                 }
             };
         } // namespace traits
diff --git a/include/alpaka/math/conj/ConjUniformCudaHipBuiltIn.hpp b/include/alpaka/math/conj/ConjUniformCudaHipBuiltIn.hpp
@@ -29,7 +29,7 @@ namespace alpaka
 
         namespace traits
         {
-            //! The CUDA/HIP built in conj trait specialization.
+            //! The CUDA/HIP built in conj trait specialization for float and double.
             template<typename TArg>
             struct Conj<ConjUniformCudaHipBuiltIn, TArg, std::enable_if_t<std::is_floating_point<TArg>::value>>
             {
diff --git a/include/alpaka/math/exp/ExpUniformCudaHipBuiltIn.hpp b/include/alpaka/math/exp/ExpUniformCudaHipBuiltIn.hpp
@@ -14,9 +14,8 @@
 #    include <alpaka/core/CudaHipMath.hpp>
 #    include <alpaka/core/Unused.hpp>
 #    include <alpaka/math/Complex.hpp>
-#    include <alpaka/math/cos/Traits.hpp>
 #    include <alpaka/math/exp/Traits.hpp>
-#    include <alpaka/math/sin/Traits.hpp>
+#    include <alpaka/math/sincos/Traits.hpp>
 
 #    include <type_traits>
 
@@ -60,7 +59,9 @@ namespace alpaka
                 __device__ auto operator()(TCtx const& ctx, Complex<T> const& arg)
                 {
                     // exp(z) = exp(x + iy) = exp(x) * (cos(y) + i * sin(y))
-                    return exp(ctx, arg.real()) * Complex<T>{cos(ctx, arg.imag()), sin(ctx, arg.imag())};
+                    auto re = T{}, im = T{};
+                    sincos(ctx, arg.imag(), im, re);
+                    return exp(ctx, arg.real()) * Complex<T>{re, im};
                 }
             };
         } // namespace traits
diff --git a/include/alpaka/math/sqrt/SqrtUniformCudaHipBuiltIn.hpp b/include/alpaka/math/sqrt/SqrtUniformCudaHipBuiltIn.hpp
@@ -63,7 +63,7 @@ namespace alpaka
                     // principal value of sqrt(z) = sqrt(|z|) * e^(i * arg(z) / 2)
                     auto const halfArg = T(0.5) * arg(ctx, argument);
                     auto re = T{}, im = T{};
-                    sincos(ctx, halfArg, re, im);
+                    sincos(ctx, halfArg, im, re);
                     return sqrt(ctx, abs(ctx, argument)) * Complex<T>(re, im);
                 }
             };
diff --git a/include/alpaka/math/tan/TanUniformCudaHipBuiltIn.hpp b/include/alpaka/math/tan/TanUniformCudaHipBuiltIn.hpp
@@ -58,10 +58,10 @@ namespace alpaka
                 template<typename TCtx>
                 __device__ auto operator()(TCtx const& ctx, Complex<T> const& arg)
                 {
-                    // tan(z) = i * (e^-iz - e^iz) / (e^-iz + e^iz)
-                    auto const exp1 = exp(ctx, Complex<T>{0.0, -1.0} * arg);
-                    auto const exp2 = exp(ctx, Complex<T>{0.0, 1.0} * arg);
-                    return Complex<T>{0.0, 1.0} * (exp1 - exp2) / (exp1 + exp2);
+                    // tan(z) = i * (e^-iz - e^iz) / (e^-iz + e^iz) = i * (1 - e^2iz) / (1 + e^2iz)
+                    // Warning: this straightforward implementation can easily result in NaN as 0/0 or inf/inf.
+                    auto const expValue = exp(ctx, Complex<T>{0.0, 2.0} * arg);
+                    return Complex<T>{0.0, 1.0} * (T{1.0} - expValue) / (T{1.0} + expValue);
                 }
             };
         } // namespace traits
diff --git a/test/unit/math/src/DataGen.hpp b/test/unit/math/src/DataGen.hpp
@@ -12,6 +12,7 @@
 #include "Defines.hpp"
 
 #include <cassert>
+#include <cmath>
 #include <limits>
 #include <random>
 
@@ -26,12 +27,12 @@ namespace alpaka
                 template<typename TData>
                 struct RngHelper
                 {
-                    static constexpr auto getMax()
+                    static auto getMax()
                     {
                         return std::numeric_limits<TData>::max();
                     }
 
-                    static constexpr auto getLowest()
+                    static auto getLowest()
                     {
                         return std::numeric_limits<TData>::lowest();
                     }
@@ -48,16 +49,16 @@ namespace alpaka
                 template<typename TData>
                 struct RngHelper<Complex<TData>>
                 {
-                    static constexpr auto getMax()
+                    static auto getMax()
                     {
-                        auto const max = std::numeric_limits<TData>::max();
-                        return Complex<TData>{max, max};
+                        /// auto const max = TData{0.01} * std::sqrt(std::numeric_limits<TData>::max());
+                        /// return Complex<TData>{max, TData{-0.7} * max};
+                        return Complex<TData>{2.0, 2.0};
                     }
 
-                    static constexpr auto getLowest()
+                    static auto getLowest()
                     {
-                        auto const lowest = std::numeric_limits<TData>::lowest();
-                        return Complex<TData>{lowest, lowest};
+                        return -getMax();
                     }
 
                     using Distribution = std::uniform_real_distribution<TData>;
@@ -98,13 +99,13 @@ namespace alpaka
                     static_assert(TArgs::capacity > 6, "Set of args must provide > 6 entries.");
                     using RngHelper = RngHelper<TData>;
                     auto rngHelper = RngHelper{};
-                    constexpr auto max = rngHelper.getMax();
-                    constexpr auto low = rngHelper.getLowest();
+                    auto const max = rngHelper.getMax();
+                    auto const low = rngHelper.getLowest();
                     std::default_random_engine eng{static_cast<std::default_random_engine::result_type>(seed)};
 
                     // These pseudo-random numbers are implementation/platform specific!
                     using Distribution = typename RngHelper::Distribution;
-                    Distribution dist(0, 1000);
+                    Distribution dist(0, /*1000*/ 10.0);
                     Distribution distOne(-1, 1);
                     for(size_t k = 0; k < TFunctor::arity_nr; ++k)
                     {
diff --git a/test/unit/math/src/Defines.hpp b/test/unit/math/src/Defines.hpp
@@ -9,6 +9,8 @@
 
 #pragma once
 
+#include <alpaka/alpaka.hpp>
+
 #include <cmath>
 #include <iomanip>
 #include <iostream>
@@ -68,6 +70,35 @@ namespace alpaka
                     return static_cast<T>(1) / sqrt(arg);
                 }
 
+                //! Stub for division expressed same way as alpaka math traits
+                template<typename TAcc, typename T>
+                ALPAKA_FN_HOST_ACC auto divides(TAcc&, T const& arg1, T const& arg2)
+                {
+                    return arg1 / arg2;
+                }
+
+                //! Stub for subtraction expressed same way as alpaka math traits
+                template<typename TAcc, typename T>
+                ALPAKA_FN_HOST_ACC auto minus(TAcc&, T const& arg1, T const& arg2)
+                {
+                    return arg1 - arg2;
+                }
+
+                //! Stub for multiplication expressed same way as alpaka math traits
+                template<typename TAcc, typename T>
+                ALPAKA_FN_HOST_ACC auto multiplies(TAcc&, T const& arg1, T const& arg2)
+                {
+                    return arg1 * arg2;
+                }
+
+                //! Stub for addition expressed same way as alpaka math traits
+                template<typename TAcc, typename T>
+                ALPAKA_FN_HOST_ACC auto plus(TAcc&, T const& arg1, T const& arg2)
+                {
+                    return arg1 + arg2;
+                }
+
+
             } // namespace math
         } // namespace unit
     } // namespace test
diff --git a/test/unit/math/src/Functor.hpp b/test/unit/math/src/Functor.hpp
@@ -133,6 +133,8 @@ namespace alpaka
 
                 ALPAKA_TEST_MATH_OP_FUNCTOR(OpCeil, Arity::Unary, std::ceil, alpaka::math::ceil, Range::Unrestricted)
 
+                ALPAKA_TEST_MATH_OP_FUNCTOR(OpConj, Arity::Unary, std::conj, alpaka::math::conj, Range::Unrestricted)
+
                 ALPAKA_TEST_MATH_OP_FUNCTOR(OpCos, Arity::Unary, std::cos, alpaka::math::cos, Range::Unrestricted)
 
                 ALPAKA_TEST_MATH_OP_FUNCTOR(OpErf, Arity::Unary, std::erf, alpaka::math::erf, Range::Unrestricted)
@@ -241,29 +243,77 @@ namespace alpaka
                     Range::Unrestricted,
                     Range::NotZero)
 
-                //! @todo add conj
-
                 // Binary functors to be used only for real types
-                //! @todo OpPow should in in RealAndComplex instead, see note there
-                using BinaryFunctorsReal = std::tuple<OpAtan2, OpFmod, OpMax, OpMin, OpRemainder, OpPow>;
-
-                // Binary functors to be used for both real and complex types
-                //! @todo OpPow should be here, something is broken for it
-                //! with NaNs in both std:: and our implementation - perhaps data is wrong for it
-                using BinaryFunctorsRealAndComplex = std::tuple</*OpPow*/>;
+                using BinaryFunctorsReal = std::tuple<OpAtan2, OpFmod, OpMax, OpMin, OpPow, OpRemainder>;
 
                 // Unary functors to be used only for real types
-                using UnaryFunctorsReal
-                    = std::tuple<OpCbrt, OpCeil, OpErf, OpFloor, OpTrunc, OpIsnan, OpIsinf, OpIsfinite>;
+                using UnaryFunctorsReal = std::tuple<
+                    OpAbs,
+                    OpAcos,
+                    OpArg,
+                    OpAsin,
+                    OpAtan,
+                    OpCbrt,
+                    OpCeil,
+                    OpConj,
+                    OpCos,
+                    OpErf,
+                    OpExp,
+                    OpFloor,
+                    OpIsnan,
+                    OpIsinf,
+                    OpIsfinite,
+                    OpLog,
+                    OpRsqrt,
+                    OpSin,
+                    OpSqrt,
+                    OpTan,
+                    OpTrunc>;
+
+                // For complex numbers also test arithmetic operations
+                ALPAKA_TEST_MATH_OP_FUNCTOR(
+                    OpDivides,
+                    Arity::Binary,
+                    std::divides<>{},
+                    alpaka::test::unit::math::divides,
+                    Range::Unrestricted,
+                    Range::NotZero)
+
+                ALPAKA_TEST_MATH_OP_FUNCTOR(
+                    OpMinus,
+                    Arity::Binary,
+                    std::minus<>{},
+                    alpaka::test::unit::math::minus,
+                    Range::Unrestricted,
+                    Range::Unrestricted)
+
+                ALPAKA_TEST_MATH_OP_FUNCTOR(
+                    OpMultiplies,
+                    Arity::Binary,
+                    std::multiplies<>{},
+                    alpaka::test::unit::math::multiplies,
+                    Range::Unrestricted,
+                    Range::Unrestricted)
+
+                ALPAKA_TEST_MATH_OP_FUNCTOR(
+                    OpPlus,
+                    Arity::Binary,
+                    std::plus<>{},
+                    alpaka::test::unit::math::plus,
+                    Range::Unrestricted,
+                    Range::Unrestricted)
+
+                // Binary functors to be used for complex types
+                using BinaryFunctorsComplex = std::tuple<OpDivides, OpMinus, OpMultiplies, OpPlus, OpPow>;
 
                 // Unary functors to be used for both real and complex types
-                using UnaryFunctorsRealAndComplex = std::tuple<
+                using UnaryFunctorsComplex = std::tuple<
                     OpAbs,
                     OpAcos,
-                    //! @todo Arg causes warning for host vs. host-device function calls, not clear why
-                    /// OpArg,
+                    OpArg,
                     OpAsin,
                     OpAtan,
+                    OpConj,
                     OpCos,
                     OpExp,
                     OpLog,
diff --git a/test/unit/math/src/math.cpp b/test/unit/math/src/math.cpp

Original file line number	Diff line number	Diff line change
`@@ -30,7 +30,7 @@ namespace alpaka`
`30`	`30`
`31`	`31`	`namespace traits`
`32`	`32`	`{`
`33`		`- //! The CUDA/HIP built in arg trait specialization.`
	`33`	`+ //! The CUDA/HIP built in arg trait specialization for float and double.`
`34`	`34`	`template<typename TArgument>`
`35`	`35`	`struct Arg<ArgUniformCudaHipBuiltIn, TArgument, std::enable_if_t<std::is_floating_point<TArgument>::value>>`
`36`	`36`	`{`
Original file line number	Diff line number	Diff line change
`@@ -60,7 +60,7 @@ namespace alpaka`
`60`	`60`	`{`
`61`	`61`	`// This holds everywhere, including the branch cuts: atan(z) = -i/2 * ln((i - z) / (i + z))`
`62`	`62`	`return Complex<T>{0.0, -0.5}`
`63`		`- * log(ctx, (Complex<T>{0.0, 1.0} + arg) / (Complex<T>{0.0, 1.0} + arg));`
	`63`	`+ * log(ctx, (Complex<T>{0.0, 1.0} - arg) / (Complex<T>{0.0, 1.0} + arg));`
`64`	`64`	`}`
`65`	`65`	`};`
`66`	`66`	`} // namespace traits`
Original file line number	Diff line number	Diff line change
`@@ -29,7 +29,7 @@ namespace alpaka`
`29`	`29`
`30`	`30`	`namespace traits`
`31`	`31`	`{`
`32`		`- //! The CUDA/HIP built in conj trait specialization.`
	`32`	`+ //! The CUDA/HIP built in conj trait specialization for float and double.`
`33`	`33`	`template<typename TArg>`
`34`	`34`	`struct Conj<ConjUniformCudaHipBuiltIn, TArg, std::enable_if_t<std::is_floating_point<TArg>::value>>`
`35`	`35`	`{`
Original file line number	Diff line number	Diff line change
`@@ -63,7 +63,7 @@ namespace alpaka`
`63`	`63`	`// principal value of sqrt(z) = sqrt(\|z\|) * e^(i * arg(z) / 2)`
`64`	`64`	`auto const halfArg = T(0.5) * arg(ctx, argument);`
`65`	`65`	`auto re = T{}, im = T{};`
`66`		`- sincos(ctx, halfArg, re, im);`
	`66`	`+ sincos(ctx, halfArg, im, re);`
`67`	`67`	`return sqrt(ctx, abs(ctx, argument)) * Complex<T>(re, im);`
`68`	`68`	`}`
`69`	`69`	`};`
Original file line number	Diff line number	Diff line change
`@@ -12,6 +12,7 @@`
`12`	`12`	`#include "Defines.hpp"`
`13`	`13`
`14`	`14`	`#include <cassert>`
	`15`	`+#include <cmath>`
`15`	`16`	`#include <limits>`
`16`	`17`	`#include <random>`
`17`	`18`
`@@ -26,12 +27,12 @@ namespace alpaka`
`26`	`27`	`template<typename TData>`
`27`	`28`	`struct RngHelper`
`28`	`29`	`{`
`29`		`- static constexpr auto getMax()`
	`30`	`+ static auto getMax()`
`30`	`31`	`{`
`31`	`32`	`return std::numeric_limits<TData>::max();`
`32`	`33`	`}`
`33`	`34`
`34`		`- static constexpr auto getLowest()`
	`35`	`+ static auto getLowest()`
`35`	`36`	`{`
`36`	`37`	`return std::numeric_limits<TData>::lowest();`
`37`	`38`	`}`
`@@ -48,16 +49,16 @@ namespace alpaka`
`48`	`49`	`template<typename TData>`
`49`	`50`	`struct RngHelper<Complex<TData>>`
`50`	`51`	`{`
`51`		`- static constexpr auto getMax()`
	`52`	`+ static auto getMax()`
`52`	`53`	`{`
`53`		`- auto const max = std::numeric_limits<TData>::max();`
`54`		`- return Complex<TData>{max, max};`
	`54`	`+ /// auto const max = TData{0.01} * std::sqrt(std::numeric_limits<TData>::max());`
	`55`	`+ /// return Complex<TData>{max, TData{-0.7} * max};`
	`56`	`+ return Complex<TData>{2.0, 2.0};`
`55`	`57`	`}`
`56`	`58`
`57`		`- static constexpr auto getLowest()`
	`59`	`+ static auto getLowest()`
`58`	`60`	`{`
`59`		`- auto const lowest = std::numeric_limits<TData>::lowest();`
`60`		`- return Complex<TData>{lowest, lowest};`
	`61`	`+ return -getMax();`
`61`	`62`	`}`
`62`	`63`
`63`	`64`	`using Distribution = std::uniform_real_distribution<TData>;`
`@@ -98,13 +99,13 @@ namespace alpaka`
`98`	`99`	`static_assert(TArgs::capacity > 6, "Set of args must provide > 6 entries.");`
`99`	`100`	`using RngHelper = RngHelper<TData>;`
`100`	`101`	`auto rngHelper = RngHelper{};`
`101`		`- constexpr auto max = rngHelper.getMax();`
`102`		`- constexpr auto low = rngHelper.getLowest();`
	`102`	`+ auto const max = rngHelper.getMax();`
	`103`	`+ auto const low = rngHelper.getLowest();`
`103`	`104`	`std::default_random_engine eng{static_cast<std::default_random_engine::result_type>(seed)};`
`104`	`105`
`105`	`106`	`// These pseudo-random numbers are implementation/platform specific!`
`106`	`107`	`using Distribution = typename RngHelper::Distribution;`
`107`		`- Distribution dist(0, 1000);`
	`108`	`+ Distribution dist(0, /1000/ 10.0);`
`108`	`109`	`Distribution distOne(-1, 1);`
`109`	`110`	`for(size_t k = 0; k < TFunctor::arity_nr; ++k)`
`110`	`111`	`{`