add more basic operations

examples/basic/basic.cu

@@ -1,44 +1,91 @@
-#include <cuda.h>
-#include <cuda_runtime.h>
+#include "../common.h"
 
-#include <cumccormick/cumccormick.cuh>
+#include <cuda_runtime.h>
 
 #include <cutangent/cutangent.cuh>
+#include <cutangent/format.h>
+
+#include <iostream>
+
+using cu::tangent;
+
+constexpr auto f(auto x, auto y)
+{
+    auto print = [](auto x) { printf("{%g, %g}\n", x.v, x.d); };
+
+    auto a = x + y;
+    auto b = x - y;
+    auto c = x * y;
+    auto d = x / y;
+    auto e = max(x, y);
+    auto f = min(x, y);
+    auto g = mid(x, y, y);
+    auto h = sin(x);
+    auto i = cos(x);
+    auto j = exp(x);
+    auto k = log(x);
+    auto l = pown(x, 2);
 
-template<typename T>
-using mc = cu::mccormick<T>;
+    print(a);
+    print(b);
+    print(c);
+    print(d);
+    print(e);
+    print(f);
+    print(g);
+    print(h);
+    print(i);
+    print(j);
+    print(k);
+    print(l);
+    return a;
+}
+
+__global__ void kernel(tangent<double> *xs, tangent<double> *ys,
+                       tangent<double> *res, int n)
+{
+    int i = blockIdx.x * blockDim.x + threadIdx.x;
+    if (i < n) {
+        res[i] = f(xs[i], ys[i]);
+    }
+}
 
 int main()
 {
-    // constexpr int n = 256;
-    // using T = mc<double>;
-    // T xs[n], ys[n], res[n];
-    //
-    // // generate dummy data
-    // for (int i = 0; i < n; i++) {
-    //     double v = i;
-    //     xs[i] = { .cv = -v, .cc = v, .box = { .lb = -v, .ub = v } };
-    //     ys[i] = { .cv = -v, .cc = v, .box = { .lb = -v, .ub = v } };
+    constexpr int n = 16;
+    using T         = tangent<double>;
+    T xs[n], ys[n], res[n];
+
+    // generate dummy data
+    for (int i = 0; i < n; i++) {
+        double v = i;
+        xs[i]    = { v, 1.0 };
+        ys[i]    = { v, 0.0 };
+    }
+
+    // for (auto el : xs) {
+    //     std::cout << el << std::endl;
+    // }
+
+    T *d_xs, *d_ys, *d_res;
+    CUDA_CHECK(cudaMalloc(&d_xs, n * sizeof(*xs)));
+    CUDA_CHECK(cudaMalloc(&d_ys, n * sizeof(*ys)));
+    CUDA_CHECK(cudaMalloc(&d_res, n * sizeof(*res)));
+
+    CUDA_CHECK(cudaMemcpy(d_xs, xs, n * sizeof(*xs), cudaMemcpyHostToDevice));
+    CUDA_CHECK(cudaMemcpy(d_ys, ys, n * sizeof(*ys), cudaMemcpyHostToDevice));
+
+    kernel<<<n, 1>>>(d_xs, d_ys, d_res, n);
+
+    CUDA_CHECK(cudaMemcpy(res, d_res, n * sizeof(*res), cudaMemcpyDeviceToHost));
+
+    // for (auto el : res) {
+    //     std::cout << el << std::endl;
     // }
-    //
-    // mc<double> *d_xs, *d_ys, *d_res;
-    // CUDA_CHECK(cudaMalloc(&d_xs, n * sizeof(*xs)));
-    // CUDA_CHECK(cudaMalloc(&d_ys, n * sizeof(*ys)));
-    // CUDA_CHECK(cudaMalloc(&d_res, n * sizeof(*res)));
-    //
-    // CUDA_CHECK(cudaMemcpy(d_xs, xs, n * sizeof(*xs), cudaMemcpyHostToDevice));
-    // CUDA_CHECK(cudaMemcpy(d_ys, ys, n * sizeof(*ys), cudaMemcpyHostToDevice));
-    //
-    // kernel<<<n, 1>>>(d_xs, d_ys, d_res, n);
-    //
-    // CUDA_CHECK(cudaMemcpy(res, d_res, n * sizeof(*res), cudaMemcpyDeviceToHost));
-    //
-    // auto r = res[0];
-    // printf("beale(0, 0) = " MCCORMICK_FORMAT "\n", r.box.lb, r.cv, r.cc, r.box.ub);
-    //
-    // CUDA_CHECK(cudaFree(d_xs));
-    // CUDA_CHECK(cudaFree(d_ys));
-    // CUDA_CHECK(cudaFree(d_res));
+
+    CUDA_CHECK(cudaFree(d_xs));
+    CUDA_CHECK(cudaFree(d_ys));
+    CUDA_CHECK(cudaFree(d_res));
 
     return 0;
 }

examples/common.h

@@ -0,0 +1,21 @@
+#ifndef CUTANGENT_COMMON_H
+#define CUTANGENT_COMMON_H
+
+#include <cstddef>
+#include <cstdio>
+#include <cstdlib>
+
+#include <cuda_runtime.h>
+
+#define CUDA_CHECK(x)                                                                \
+    do {                                                                             \
+        cudaError_t err = x;                                                         \
+        if (err != cudaSuccess) {                                                    \
+            fprintf(stderr, "CUDA error in %s at %s:%d: %s (%s=%d)\n", __FUNCTION__, \
+                    __FILE__, __LINE__, cudaGetErrorString(err),                     \
+                    cudaGetErrorName(err), err);                                     \
+            abort();                                                                 \
+        }                                                                            \
+    } while (0)
+
+#endif // CUTANGENT_COMMON_H

include/cutangent/arithmetic/basic.cuh

@@ -4,28 +4,45 @@
 #include <cutangent/tangent.h>
 
 #include <algorithm>
-// #include <cmath>
-// #include <numbers>
+#include <cmath>
 
 namespace cu
 {
 
-#define cuda_fn inline constexpr __device__
+#define fn inline constexpr __device__
 
 template<typename T>
-cuda_fn tangent<T> operator+(tangent<T> a, tangent<T> b)
+fn tangent<T> operator-(tangent<T> x)
+{
+    return { -x.v, -x.d };
+}
+
+template<typename T>
+fn tangent<T> operator+(tangent<T> a, tangent<T> b)
 {
     return { .v = a.v + b.v, .d = a.d + b.d };
 }
 
 template<typename T>
-cuda_fn tangent<T> operator*(tangent<T> a, tangent<T> b)
+fn tangent<T> operator-(tangent<T> a, tangent<T> b)
+{
+    return { .v = a.v - b.v, .d = a.d - b.d };
+}
+
+template<typename T>
+fn tangent<T> operator*(tangent<T> a, tangent<T> b)
 {
     return { .v = a.v * b.v, .d = a.v * b.d + a.d * b.v };
 }
 
 template<typename T>
-cuda_fn tangent<T> max(tangent<T> a, tangent<T> b)
+fn tangent<T> operator/(tangent<T> a, tangent<T> b)
+{
+    return { .v = a.v / b.v, .d = (a.d * b.v - a.v * b.d) / (b.v * b.v) };
+}
+
+template<typename T>
+fn tangent<T> max(tangent<T> a, tangent<T> b)
 {
     using std::max;
 
@@ -34,7 +51,7 @@ cuda_fn tangent<T> max(tangent<T> a, tangent<T> b)
 }
 
 template<typename T>
-cuda_fn tangent<T> min(tangent<T> a, tangent<T> b)
+fn tangent<T> min(tangent<T> a, tangent<T> b)
 {
     using std::min;
 
@@ -43,11 +60,69 @@ cuda_fn tangent<T> min(tangent<T> a, tangent<T> b)
 }
 
 template<typename T>
-cuda_fn tangent<T> mid(tangent<T> v, tangent<T> lb, tangent<T> ub)
+fn tangent<T> clamp(tangent<T> v, tangent<T> lb, tangent<T> ub)
 {
-    return { .v = mid(v.v, lb.v, ub.v), .d = lb.d * (v.v < lb.v) + ub.d * (v.v > ub.v) };
+    using std::clamp;
+
+    return { .v = clamp(v.v, lb.v, ub.v), .d = lb.d * (v.v < lb.v) + ub.d * (v.v > ub.v) };
+}
+
+template<typename T>
+fn tangent<T> mid(tangent<T> v, tangent<T> lb, tangent<T> ub)
+{
+    return clamp(v, lb, ub);
+}
+
+template<typename T>
+fn tangent<T> sin(tangent<T> x)
+{
+    using std::cos;
+    using std::sin;
+
+    return { .v = sin(x.v), .d = cos(x.v) * x.d };
 }
 
+template<typename T>
+fn tangent<T> cos(tangent<T> x)
+{
+    using std::cos;
+    using std::sin;
+
+    return { .v = cos(x.v), .d = -sin(x.v) * x.d };
+}
+
+template<typename T>
+fn tangent<T> exp(tangent<T> x)
+{
+    using std::exp;
+
+    return { .v = exp(x.v), .d = exp(x.v) * x.d };
+}
+
+template<typename T>
+fn tangent<T> log(tangent<T> x)
+{
+    using std::log;
+
+    return { .v = log(x.v), .d = x.d / x.v };
+}
+
+template<typename T>
+fn tangent<T> sqr(tangent<T> x)
+{
+    return { .v = sqr(x.v), .d = 2.0 * x.v * x.d };
+}
+
+template<typename T>
+fn tangent<T> pown(tangent<T> x, auto n)
+{
+    using std::pow;
+
+    return { .v = pow(x.v, n), .d = n * pow(x.v, n - 1) * x.d };
+}
+
+#undef fn
+
 } // namespace cu
 
 #endif // CUTANGENT_ARITHMETIC_BASIC_CUH

include/cutangent/format.h

@@ -0,0 +1,19 @@
+#ifndef CUTANGENT_FORMAT_H
+#define CUTANGENT_FORMAT_H
+
+#include <cutangent/tangent.h>
+
+#include <ostream>
+
+namespace cu
+{
+
+template<typename T>
+std::ostream &operator<<(std::ostream &os, tangent<T> x)
+{
+    return os << "{v: " << x.v << ", d: " << x.d << "}";
+}
+
+} // namespace cu
+
+#endif // CUTANGENT_FORMAT_H

include/cutangent/tangent.h

@@ -9,12 +9,13 @@ namespace cu
 template<typename T>
 struct tangent
 {
-    T v;
-    T d;
+    T v; // value
+    T d; // derivative
 
     constexpr auto operator<=>(const tangent &) const = default;
 };
 
+
 } // namespace cu
 
 #endif // CUTANGENT_TANGENT_H