// SPDX-License-Identifier: BSD-3-Clause
/*
 * victory_math.cu - mixed int/float/double math over 64K elements.
 *
 * Each thread computes a float lane (sqrtf + fmaf), a double lane (sqrt +
 * multiply), and an integer lane (wrapping multiply/add, xor, shift, popc).
 * The host recomputes all three lanes and verifies the float/double lanes
 * within a relative tolerance and the integer lane exactly.
 *
 * Output lines are stable and greppable:
 *
 *   DEVICE <name> cc=<major>.<minor> sm=<count> mem=<bytes>
 *   MATH n=<count> badf=<float_mismatches> badd=<double_mismatches> badi=<int_mismatches> csumf=<sum> csumd=<sum>
 *   STATUS OK|FAIL
 *
 * Build with a static CUDA runtime (--cudart=static) so the guest only needs
 * the driver's libcuda.so, not a CUDA toolkit installation:
 *
 *   nvcc -O2 --cudart=static -o victory_math victory_math.cu
 *
 * Then, from the guest:  /tmp/victory_math
 */

#include <cstdio>
#include <cstdlib>
#include <cstdint>
#include <cmath>
#include <cuda_runtime.h>

#define CHECK(x)                                                         \
	do {                                                             \
		cudaError_t e_ = (x);                                    \
		if (e_ != cudaSuccess) {                                 \
			printf("CUDA_ERROR %s at line %d\n",             \
			       cudaGetErrorString(e_), __LINE__);        \
			printf("STATUS FAIL\n");                         \
			return 1;                                        \
		}                                                        \
	} while (0)

/* Portable popcount: device intrinsic under nvcc, host builtin elsewhere. */
__host__ __device__ static inline unsigned int popc(unsigned int v)
{
#if defined(__CUDA_ARCH__)
	return (unsigned int)__popc(v);
#else
	return (unsigned int)__builtin_popcount(v);
#endif
}

/* Integer lane, replicated bit-for-bit on the host. */
__host__ __device__ static inline int int_lane(unsigned int u)
{
	unsigned int w = u * 31u + 7u;

	return (int)((w ^ (u >> 3)) + popc(u & 0xFFu));
}

__global__ void mathmix(const unsigned int *in, float *of, double *od,
			int *oi, int n)
{
	int i = blockIdx.x * blockDim.x + threadIdx.x;

	if (i < n) {
		unsigned int u = in[i];

		of[i] = sqrtf((float)u) +
			fmaf((float)(u % 7u), 1.5f, -2.0f);
		od[i] = sqrt((double)u) * 1.0000001;
		oi[i] = int_lane(u);
	}
}

int main(void)
{
	cudaDeviceProp prop;
	const int n = 1 << 16; /* 64K elements */
	size_t ib = (size_t)n * sizeof(unsigned int);
	size_t fb = (size_t)n * sizeof(float);
	size_t db = (size_t)n * sizeof(double);
	unsigned int *h_in;
	float *h_of, *h_rf;
	double *h_od, *h_rd;
	int *h_oi, *h_ri;
	unsigned int *d_in;
	float *d_of;
	double *d_od;
	int *d_oi;
	double csumf = 0.0, csumd = 0.0;
	int badf = 0, badd = 0, badi = 0, i;

	setvbuf(stdout, NULL, _IONBF, 0);

	CHECK(cudaGetDeviceProperties(&prop, 0));
	printf("DEVICE %s cc=%d.%d sm=%d mem=%llu\n", prop.name,
	       prop.major, prop.minor, prop.multiProcessorCount,
	       (unsigned long long)prop.totalGlobalMem);

	h_in = (unsigned int *)malloc(ib);
	h_of = (float *)malloc(fb);
	h_rf = (float *)malloc(fb);
	h_od = (double *)malloc(db);
	h_rd = (double *)malloc(db);
	h_oi = (int *)malloc(ib);
	h_ri = (int *)malloc(ib);
	if (!h_in || !h_of || !h_rf || !h_od || !h_rd || !h_oi || !h_ri)
		return 1;

	for (i = 0; i < n; i++)
		h_in[i] = (unsigned int)(i % 100000) + 1u;

	CHECK(cudaMalloc(&d_in, ib));
	CHECK(cudaMalloc(&d_of, fb));
	CHECK(cudaMalloc(&d_od, db));
	CHECK(cudaMalloc(&d_oi, ib));
	CHECK(cudaMemcpy(d_in, h_in, ib, cudaMemcpyHostToDevice));

	mathmix<<<(n + 255) / 256, 256>>>(d_in, d_of, d_od, d_oi, n);
	CHECK(cudaGetLastError());
	CHECK(cudaDeviceSynchronize());
	CHECK(cudaMemcpy(h_of, d_of, fb, cudaMemcpyDeviceToHost));
	CHECK(cudaMemcpy(h_od, d_od, db, cudaMemcpyDeviceToHost));
	CHECK(cudaMemcpy(h_oi, d_oi, ib, cudaMemcpyDeviceToHost));

	for (i = 0; i < n; i++) {
		unsigned int u = h_in[i];

		h_rf[i] = sqrtf((float)u) + fmaf((float)(u % 7u), 1.5f, -2.0f);
		h_rd[i] = sqrt((double)u) * 1.0000001;
		h_ri[i] = int_lane(u);

		if (fabsf(h_of[i] - h_rf[i]) > 1e-4f * (1.0f + fabsf(h_rf[i])))
			badf++;
		if (fabs(h_od[i] - h_rd[i]) > 1e-12 * (1.0 + fabs(h_rd[i])))
			badd++;
		if (h_oi[i] != h_ri[i])
			badi++;
		csumf += (double)h_of[i];
		csumd += h_od[i];
	}

	printf("MATH n=%d badf=%d badd=%d badi=%d csumf=%.6f csumd=%.6f\n", n,
	       badf, badd, badi, csumf, csumd);

	if (badf != 0 || badd != 0 || badi != 0) {
		printf("STATUS FAIL\n");
		return 1;
	}

	printf("STATUS OK\n");
	CHECK(cudaFree(d_in));
	CHECK(cudaFree(d_of));
	CHECK(cudaFree(d_od));
	CHECK(cudaFree(d_oi));
	free(h_in);
	free(h_of);
	free(h_rf);
	free(h_od);
	free(h_rd);
	free(h_oi);
	free(h_ri);
	return 0;
}