#include <iostream>
#include <climits>
using namespace std;

// One kernel handles all 4 operations
__global__ void reduce(int* arr, int* out_min, int* out_max, int* out_sum, int n) {
    // Shared memory: each block gets its own fast scratch space
    __shared__ int s_min[256], s_max[256], s_sum[256];

    int tid = threadIdx.x;
    int i   = blockIdx.x * blockDim.x + threadIdx.x;

    // Load from global memory into shared memory (or identity values if out of bounds)
    s_min[tid] = (i < n) ? arr[i] : INT_MAX;
    s_max[tid] = (i < n) ? arr[i] : INT_MIN;
    s_sum[tid] = (i < n) ? arr[i] : 0;
    __syncthreads(); // wait for all threads to finish loading

    // Reduction loop: stride starts at half block, halves each round
    for (int stride = blockDim.x / 2; stride > 0; stride >>= 1) {
        if (tid < stride) {
            s_min[tid] = min(s_min[tid], s_min[tid + stride]);
            s_max[tid] = max(s_max[tid], s_max[tid + stride]);
            s_sum[tid] = s_sum[tid] + s_sum[tid + stride];
        }
        __syncthreads(); // sync after every stride level
    }

    // Thread 0 of each block writes the block's result atomically
    if (tid == 0) {
        atomicMin(out_min, s_min[0]);
        atomicMax(out_max, s_max[0]);
        atomicAdd(out_sum, s_sum[0]);
    }
}

int main() {
    const int N = 1024;
    int h_arr[N];
    for (int i = 0; i < N; i++) h_arr[i] = rand() % 1000;

    // Allocate GPU memory
    int *d_arr, *d_min, *d_max, *d_sum;
    cudaMalloc(&d_arr, N * sizeof(int));
    cudaMalloc(&d_min, sizeof(int));
    cudaMalloc(&d_max, sizeof(int));
    cudaMalloc(&d_sum, sizeof(int));

    // Initialize output values to identity elements
    int init_min = INT_MAX, init_max = INT_MIN, init_sum = 0;
    cudaMemcpy(d_arr, h_arr, N * sizeof(int), cudaMemcpyHostToDevice);
    cudaMemcpy(d_min, &init_min, sizeof(int), cudaMemcpyHostToDevice);
    cudaMemcpy(d_max, &init_max, sizeof(int), cudaMemcpyHostToDevice);
    cudaMemcpy(d_sum, &init_sum, sizeof(int), cudaMemcpyHostToDevice);

    // Launch: 256 threads per block, enough blocks to cover N elements
    reduce<<<(N + 255) / 256, 256>>>(d_arr, d_min, d_max, d_sum, N);

    // Copy results back
    int h_min, h_max, h_sum;
    cudaMemcpy(&h_min, d_min, sizeof(int), cudaMemcpyDeviceToHost);
    cudaMemcpy(&h_max, d_max, sizeof(int), cudaMemcpyDeviceToHost);
    cudaMemcpy(&h_sum, d_sum, sizeof(int), cudaMemcpyDeviceToHost);

    cout << "Min: "  << h_min << endl;
    cout << "Max: "  << h_max << endl;
    cout << "Sum: "  << h_sum << endl;
    cout << "Avg: "  << (float)h_sum / N << endl;

    cudaFree(d_arr); cudaFree(d_min); cudaFree(d_max); cudaFree(d_sum);
    return 0;
}
