(draft) Unary math cuda kernels

cuda/abs.cu

@@ -0,0 +1,10 @@
+// dst[i] = abs(a[i])
+
+extern "C" __global__ void
+unary_abs(float *__restrict__ dst, float *__restrict__ a, int N){
+    int i = (blockIdx.y * gridDim.x + blockIdx.x) * blockDim.x + threadIdx.x;
+
+    if (i < N){
+        dst[i] = fabsf(a[i]);
+    }
+}

cuda/acos.cu

@@ -0,0 +1,16 @@
+// dst[i] = acos(a[i]), returns 0 for a[i] outside of domain [-1, 1]
+
+extern "C" __global__ void
+unary_acos(float *__restrict__ dst, float *__restrict__ a, int N){
+    int i = (blockIdx.y * gridDim.x + blockIdx.x) * blockDim.x + threadIdx.x;
+
+    if (i < N){
+        float x = a[i];
+        if (x >= -1.0f && x <= 1.0f){
+            dst[i] = acosf(x);
+        }
+        else{
+            dst[i] = 0.0f;
+        }
+    }
+}

cuda/acosh.cu

@@ -0,0 +1,15 @@
+//dst[i] = acosh(a[i]), returns 0 for a[i] outside of domain [1, inf)
+
+extern "C" __global__ void
+unary_acosh(float* __restrict__ dst, float* __restrict__ a, int N) {
+    int i = (blockIdx.y * gridDim.x + blockIdx.x) * blockDim.x + threadIdx.x;
+
+    if (i < N) {
+        float x = a[i];
+        if (x >= 1.0f) {
+            dst[i] = acoshf(x);
+        } else {
+            dst[i] = 0.0f;
+        }
+    }
+}

cuda/asin.cu

@@ -0,0 +1,15 @@
+// dst[i] = asin(a[i]), returns 0 for a[i] outside of domain [-1, 1]
+
+extern "C" __global__ void
+unary_asin(float* __restrict__ dst, float* __restrict__ a, int N) {
+    int i = (blockIdx.y * gridDim.x + blockIdx.x) * blockDim.x + threadIdx.x;
+
+    if (i < N) {
+        float x = a[i];
+        if (x >= -1.0f && x <= 1.0f) {
+            dst[i] = asinf(x);
+        } else {
+            dst[i] = 0.0f;
+        }
+    }
+}

cuda/asinh.cu

@@ -0,0 +1,10 @@
+//dst[i] = asinh(a[i])
+
+extern "C" __global__ void
+unary_asinh(float *__restrict__ dst, float *__restrict__ a, int N){
+    int i = (blockIdx.y * gridDim.x + blockIdx.x) * blockDim.x + threadIdx.x;
+
+    if (i < N){
+        dst[i] = asinhf(a[i]);
+    }
+}

cuda/atan.cu

@@ -0,0 +1,10 @@
+//dst[i] = atan(a[i])
+
+extern "C" __global__ void
+unary_atan(float *__restrict__ dst, float *__restrict__ a, int N){
+    int i = (blockIdx.y * gridDim.x + blockIdx.x) * blockDim.x + threadIdx.x;
+
+    if (i < N){
+        dst[i] = atanf(a[i]);
+    }
+}

cuda/atan2.cu

@@ -0,0 +1,10 @@
+//dst[i] = atan2(a[i], b[i])
+
+extern "C" __global__ void
+pw_atan2(float *__restrict__ dst, float *__restrict__ a, float *__restrict__ b, int N){
+    int i = (blockIdx.y * gridDim.x + blockIdx.x) * blockDim.x + threadIdx.x;
+
+    if (i < N){
+        dst[i] = atan2f(a[i], b[i]);
+    }
+}

cuda/atanh.cu

@@ -0,0 +1,16 @@
+// dst[i] = atanh(a[i]), returns 0 for a[i] outside of domain (-1, 1)
+
+extern "C" __global__ void
+unary_atanh(float* __restrict__ dst, float* __restrict__ a, int N) {
+
+    int i = (blockIdx.y * gridDim.x + blockIdx.x) * blockDim.x + threadIdx.x;
+
+    if (i < N) {
+        float x = a[i];
+        if (x > -1.0f && x < 1.0f) {
+            dst[i] = atanhf(x);
+        } else {
+            dst[i] = 0.0f;
+        }
+    }
+}

cuda/cos.cu

@@ -0,0 +1,9 @@
+//dst[i] = cos(a[i])
+extern "C" __global__ void
+unary_cos(float *__restrict__ dst, float *__restrict__ a, int N){
+    int i = (blockIdx.y * gridDim.x + blockIdx.x) * blockDim.x + threadIdx.x;
+
+    if (i < N){
+        dst[i] = cosf(a[i]);
+    }
+}

cuda/cosh.cu

@@ -0,0 +1,9 @@
+//dst[i] = cosh(a[i])
+extern "C" __global__ void
+unary_cosh(float *__restrict__ dst, float *__restrict__ a, int N){
+    int i = (blockIdx.y * gridDim.x + blockIdx.x) * blockDim.x + threadIdx.x;
+
+    if (i < N){
+        dst[i] = coshf(a[i]);
+    }
+}

cuda/erf.cu

@@ -0,0 +1,9 @@
+//dst[i] = erf(a[i])
+extern "C" __global__ void
+unary_erf(float *__restrict__ dst, float *__restrict__ a, int N){
+    int i = (blockIdx.y * gridDim.x + blockIdx.x) * blockDim.x + threadIdx.x;
+
+    if (i < N){
+        dst[i] = erff(a[i]);
+    }
+}

cuda/erfc.cu

@@ -0,0 +1,9 @@
+//dst[i] = erfc(a[i])
+extern "C" __global__ void
+unary_erfc(float *__restrict__ dst, float *__restrict__ a, int N){
+    int i = (blockIdx.y * gridDim.x + blockIdx.x) * blockDim.x + threadIdx.x;
+
+    if (i < N){
+        dst[i] = erfcf(a[i]);
+    }
+}

cuda/exp.cu

@@ -0,0 +1,9 @@
+//dst[i] = exp(a[i])
+extern "C" __global__ void
+unary_exp(float *__restrict__ dst, float *__restrict__ a, int N){
+    int i = (blockIdx.y * gridDim.x + blockIdx.x) * blockDim.x + threadIdx.x;
+
+    if (i < N){
+        dst[i] = expf(a[i]);
+    }
+}

cuda/gamma.cu

@@ -0,0 +1,15 @@
+// dst[i] = gamma(a[i]), returns 0 for x <= 0, returns 0 for a[i] outside of domain (0, inf) with poles at non-positive integers
+extern "C" __global__ void
+unary_gamma(float* __restrict__ dst, float* __restrict__ a, int N) {
+
+    int i = (blockIdx.y * gridDim.x + blockIdx.x) * blockDim.x + threadIdx.x;
+
+    if (i < N) {
+        float x = a[i];
+        if (x > 0.0f) {
+            dst[i] = tgammaf(x);
+        } else {
+            dst[i] = 0.0f;
+        }
+    }
+}

cuda/heaviside.cu

@@ -0,0 +1,19 @@
+//dst[i] = heaviside(a[i])
+// returns 0 for x < 0, 0.5 for x == 0, 1 for x > 0
+extern "C" __global__ void
+unary_heaviside(float* __restrict__ dst, float* __restrict__ a, int N) {
+    int i = (blockIdx.y * gridDim.x + blockIdx.x) * blockDim.x + threadIdx.x;
+
+    if (i < N) {
+        float x = a[i];
+
+        if (x > 0.0f) {
+            dst[i] = 1.0f;
+        } else if (x < 0.0f) {
+            dst[i] = 0.0f;
+        } else {
+            // x == 0
+            dst[i] = 0.5f;
+        }
+    }
+}

cuda/log.cu

@@ -0,0 +1,13 @@
+// dst[i] = log(a[i]), returns 0 for non-positive input
+extern "C" __global__ void
+unary_log(float* __restrict__ dst, float* __restrict__ a, int N) {
+    int i = (blockIdx.y * gridDim.x + blockIdx.x) * blockDim.x + threadIdx.x;
+
+    if (i < N) {
+        if (a[i] > 0.0f) {
+            dst[i] = logf(a[i]);
+        } else {
+            dst[i] = 0.0f;
+        }
+    }
+}

cuda/mod.cu

@@ -0,0 +1,10 @@
+//dst[i] = mod(a[i], b[i])
+extern "C" __global__ void
+pw_mod(float *__restrict__ dst, float *__restrict__ a, float *__restrict__ b, int N)
+{
+    int i = (blockIdx.y * gridDim.x + blockIdx.x) * blockDim.x + threadIdx.x;
+    if (i < N)
+    {
+        dst[i] = fmodf(a[i], b[i]);
+    }
+}

cuda/pow.cu

@@ -0,0 +1,28 @@
+// dst[i] = pow(a[i], b[i]), pow(a,b) for negative a and b returns -pow(-a,b) for fractional exponents, and for 0^0 returns 1
+extern "C" __global__ void
+pw_pow(float* __restrict__ dst, float* __restrict__ a, float* __restrict__ b, int N) {
+int i = (blockIdx.y * gridDim.x + blockIdx.x) * blockDim.x + threadIdx.x;
+
+    if (i < N) {
+        float x = a[i];
+        float y = b[i];
+
+        // Guard: 0^negative → return 0
+        if (x == 0.0f && y < 0.0f) {
+            dst[i] = 0.0f;
+            return;
+        }
+
+        // Check if exponent is effectively integer
+        float y_floor = floorf(y);
+        bool y_is_int = fabsf(y - y_floor) < 1e-6f;
+
+        if (x < 0.0f && !y_is_int) {
+            // Fractional exponent → move minus sign outside
+            dst[i] = -1.0f * powf(fabsf(x), y);
+        } else {
+            // Integer exponent or positive base
+            dst[i] = powf(x, y);
+        }
+    }
+}

cuda/sin.cu

@@ -0,0 +1,10 @@
+//dst[i] = sin(a[i])
+extern "C" __global__ void
+unary_sin(float *__restrict__ dst, float *__restrict__ a, int N) {
+    int i = (blockIdx.y * gridDim.x + blockIdx.x) * blockDim.x + threadIdx.x;
+
+    if (i < N)
+    {
+        dst[i] = sinf(a[i]);
+    }
+}

cuda/sinc.cu

@@ -0,0 +1,11 @@
+// dst[i] = sinc(a[i])
+extern "C" __global__ void
+unary_sinc(float *__restrict__ dst, float *__restrict__ a, int N) {
+    int i = (blockIdx.y * gridDim.x + blockIdx.x) * blockDim.x + threadIdx.x;
+
+    if (i < N)
+    {
+        float x = a[i];
+        dst[i] = (x == 0.0f) ? 1.0f : sinf(x) / x;
+    }
+}

cuda/sinh.cu

@@ -0,0 +1,10 @@
+//dst[i] = sinh(a[i])
+extern "C" __global__ void
+unary_sinh(float *__restrict__ dst, float *__restrict__ a, int N){
+    int i = (blockIdx.y * gridDim.x + blockIdx.x) * blockDim.x + threadIdx.x;
+
+    if (i < N)
+    {
+        dst[i] = sinhf(a[i]);
+    }
+}

cuda/tan.cu

@@ -0,0 +1,9 @@
+//dst[i] = tan(a[i]), for poles at π/2 + nπ, returns 0
+extern "C" __global__ void
+unary_tan(float* __restrict__ dst, float* __restrict__ a, int N) {
+    int i = (blockIdx.y * gridDim.x + blockIdx.x) * blockDim.x + threadIdx.x;
+
+    if (i < N) {
+        dst[i] = tanf(a[i]);
+    }
+}

cuda/tanh.cu

@@ -0,0 +1,9 @@
+//dst[i] = tanh(a[i])
+extern "C" __global__ void
+unary_tanh(float *__restrict__ dst, float *__restrict__ a, int N){
+    int i = (blockIdx.y * gridDim.x + blockIdx.x) * blockDim.x + threadIdx.x;
+
+    if (i < N){
+        dst[i] = tanhf(a[i]);
+    }
+}