fix tangent interval example and add centered form example

examples/interval/CMakeLists.txt

@@ -1,5 +1,3 @@
-add_executable(interval interval.cu)
-
 include(FetchContent)
 FetchContent_Declare(
   cuinterval
@@ -8,8 +6,14 @@ FetchContent_Declare(
 )
 FetchContent_MakeAvailable(cuinterval)
 
-target_link_libraries(interval PUBLIC cuinterval)
+add_executable(natural natural.cu)
+target_link_libraries(natural PUBLIC cuinterval)
+target_compile_features(natural PRIVATE cxx_std_20 cuda_std_20)
+set_target_properties(natural PROPERTIES CUDA_ARCHITECTURES native)
+target_link_libraries(natural PRIVATE ${PROJECT_NAME})
 
-target_compile_features(interval PRIVATE cxx_std_20 cuda_std_20)
-set_target_properties(interval PROPERTIES CUDA_ARCHITECTURES native)
-target_link_libraries(interval PRIVATE ${PROJECT_NAME})
+add_executable(centered centered.cu)
+target_link_libraries(centered PUBLIC cuinterval)
+target_compile_features(centered PRIVATE cxx_std_20 cuda_std_20)
+set_target_properties(centered PROPERTIES CUDA_ARCHITECTURES native)
+target_link_libraries(centered PRIVATE ${PROJECT_NAME})

examples/interval/centered.cu

@@ -0,0 +1,76 @@
+#include "../common.h"
+
+#include <cuda_runtime.h>
+
+#include <cuinterval/cuinterval.h>
+#include <cuinterval/format.h>
+
+#include <cutangent/cutangent.cuh>
+#include <cutangent/format.h>
+
+#include <iostream>
+
+using cu::tangent;
+
+using I = cu::interval<double>;
+
+using T = tangent<I>;
+
+constexpr auto f(auto x)
+{
+    using std::pow;
+
+    return 3.0 * pow(x, 3) + pow(x, 2) - 5.0 * x - 1.0;
+}
+
+constexpr auto g(auto x, auto y)
+{
+    return x * y;
+}
+
+__global__ void centered_form(T *xs, I *res, int n)
+{
+    int i = blockIdx.x * blockDim.x + threadIdx.x;
+    if (i < n) {
+        auto r = f(xs[i]);
+
+        auto L = derivative(r);
+        auto X = value(xs[i]);
+        auto c = mid(value(xs[i]));
+
+        auto centered_form = f(c) + L * (X - c);
+
+        res[i] = centered_form;
+    }
+}
+
+int main()
+{
+    constexpr int n = 1;
+    T xs[n];
+    I res[n];
+
+    value(xs[0])      = { -1.0, 1.0 };
+    derivative(xs[0]) = { 1.0, 1.0 };
+
+    T *d_xs;
+    I *d_res;
+    CUDA_CHECK(cudaMalloc(&d_xs, n * sizeof(*xs)));
+    CUDA_CHECK(cudaMalloc(&d_res, n * sizeof(*res)));
+
+    CUDA_CHECK(cudaMemcpy(d_xs, xs, n * sizeof(*xs), cudaMemcpyHostToDevice));
+
+    centered_form<<<n, 1>>>(d_xs, d_res, n);
+
+    CUDA_CHECK(cudaMemcpy(res, d_res, n * sizeof(*res), cudaMemcpyDeviceToHost));
+
+    auto r = res[0];
+    std::cout << "f(x) = x^3 + x^2 - 5x - 1" << std::endl;
+    std::cout << "X = " << value(xs[0]) << std::endl;
+    std::cout << "centered form for f at X: " << r << std::endl;
+
+    CUDA_CHECK(cudaFree(d_xs));
+    CUDA_CHECK(cudaFree(d_res));
+
+    return 0;
+}

examples/interval/interval.cu → examples/interval/natural.cu

@@ -15,6 +15,8 @@ using cu::tangent;
 template<typename T>
 using I = cu::interval<T>;
 
+using T = tangent<I<double>>;
+
 constexpr auto f(auto x, auto y)
 {
     // Currently supported functions:
@@ -27,19 +29,19 @@ constexpr auto f(auto x, auto y)
     // auto a = x + y;
     // auto a = x - y;
     // auto a = x / y;
-    // auto a = x * y;
+    auto a = x * y;
     // auto a = sqr(x);
     // auto a = sqrt(x);
     // auto a = abs(x);
     // auto a = exp(x);
     // auto a = log(x);
     // auto a = recip(x);
     // auto a = cos(x);
-    // auto a = pown(x, 3);
+    // auto a = pow(x, 3);
     // auto a = pown(x, 4.0);
     // auto a = pow(x, 4.0);
     // auto a = pow(x, y);
-    auto a = atan2(y, x);
+    // auto a = atan2(y, x);
     // auto a = atan2(y, 2.0);
     // auto a = atan2(2.0, x);
     // auto a = max(x, y);
@@ -48,7 +50,7 @@ constexpr auto f(auto x, auto y)
     return a;
 }
 
-__global__ void kernel(I<tangent<double>> *xs, I<tangent<double>> *ys, I<tangent<double>> *res, int n)
+__global__ void kernel(T *xs, T *ys, T *res, int n)
 {
     int i = blockIdx.x * blockDim.x + threadIdx.x;
     if (i < n) {
@@ -59,14 +61,13 @@ __global__ void kernel(I<tangent<double>> *xs, I<tangent<double>> *ys, I<tangent
 int main()
 {
     constexpr int n = 1;
-    using T         = I<tangent<double>>;
     T xs[n], ys[n], res[n];
 
     // generate dummy data
-
-    xs[0] = { .lb = { 0.5, 1.0 }, .ub = { 3.0, 1.0 } } ;
-
-    ys[0] = { .lb = { 2.0, 0.0 }, .ub = { 5.0, 0.0 } };
+    value(xs[0])      = { 0.5, 3.0 };
+    derivative(xs[0]) = { 1.0, 1.0 };
+    value(ys[0])      = { 2.0, 5.0 };
+    derivative(ys[0]) = { 0.0, 0.0 };
 
     std::cout << xs[0] << std::endl;
     std::cout << ys[0] << std::endl;

include/cutangent/arithmetic/basic.cuh

@@ -62,9 +62,9 @@ fn tangent<T> operator-(tangent<T> a, T b)
 }
 
 template<typename T>
-fn tangent<T> operator-(tangent<T> a, std::integral auto b)
+fn tangent<T> operator-(tangent<T> a, auto b)
 {
-    return { a.v - static_cast<T>(b), a.d };
+    return { a.v - b, a.d };
 }
 
 template<typename T>
@@ -74,9 +74,9 @@ fn tangent<T> operator-(T a, tangent<T> b)
 }
 
 template<typename T>
-fn tangent<T> operator-(std::integral auto a, tangent<T> b)
+fn tangent<T> operator-(auto a, tangent<T> b)
 {
-    return { static_cast<T>(a) - b.v, -b.d };
+    return { a - b.v, -b.d };
 }
 
 template<typename T>
@@ -100,7 +100,7 @@ fn tangent<T> operator*(T a, tangent<T> b)
 template<typename T>
 fn tangent<T> operator*(tangent<T> a, std::integral auto b)
 {
-    return { a.v * static_cast<T>(b), a.d * static_cast<T>(b) };
+    return { a.v * b, a.d * b };
 }
 
 template<typename T>
@@ -109,6 +109,18 @@ fn tangent<T> operator*(std::integral auto a, tangent<T> b)
     return b * a;
 }
 
+template<typename T>
+fn tangent<T> operator*(tangent<T> a, auto b)
+{
+    return { a.v * b, a.d * b };
+}
+
+template<typename T>
+fn tangent<T> operator*(auto a, tangent<T> b)
+{
+    return b * a;
+}
+
 template<typename T>
 fn tangent<T> operator/(tangent<T> a, tangent<T> b)
 {

include/cutangent/arithmetic/intrinsic.cuh

@@ -75,10 +75,10 @@ namespace cu::intrinsic
     // template<> inline __device__ double exp10     (double x)           { return ::exp10(x); }
     // template<> inline __device__ double exp2      (double x)           { return ::exp2(x); }
     template<> inline __device__ __host__ tangent<double> nan() { return { ::nan(""), ::nan("") }; }
-    template<> inline __device__ tangent<double> neg_inf() { return { __longlong_as_double(0xfff0000000000000ull), 0.0 }; }
-    template<> inline __device__ tangent<double> pos_inf() { return { __longlong_as_double(0x7ff0000000000000ull), 0.0 }; }
-    template<> inline __device__ tangent<double> next_floating(tangent<double> x) { return { nextafter(x.v, pos_inf<tangent<double>>().v), 0.0 }; }
-    template<> inline __device__ tangent<double> prev_floating(tangent<double> x) { return { nextafter(x.v, neg_inf<tangent<double>>().v), 0.0 }; }
+    template<> inline __device__ tangent<double> neg_inf() { return { __longlong_as_double(0xfff0000000000000ull), __longlong_as_double(0xfff0000000000000ull) }; }
+    template<> inline __device__ tangent<double> pos_inf() { return { __longlong_as_double(0x7ff0000000000000ull), __longlong_as_double(0x7ff0000000000000ull) }; }
+    template<> inline __device__ tangent<double> next_floating(tangent<double> x) { return { nextafter(x.v, pos_inf<tangent<double>>().v), nextafter(x.d, pos_inf<tangent<double>>().d) }; }
+    template<> inline __device__ tangent<double> prev_floating(tangent<double> x) { return { nextafter(x.v, neg_inf<tangent<double>>().v), nextafter(x.d, neg_inf<tangent<double>>().d) }; }
 
 // clang-format on
 } // namespace cu::intrinsic

include/cutangent/tangent.h

@@ -28,6 +28,18 @@ struct tangent
     constexpr bool operator==(const tangent &other) const noexcept { return v == other.v; }
 };
 
+template<typename T>
+constexpr T &value(tangent<T> &x)
+{
+    return x.v;
+}
+
+template<typename T>
+constexpr T &derivative(tangent<T> &x)
+{
+    return x.d;
+}
+
 } // namespace cu
 
 #endif // CUTANGENT_TANGENT_H