from stdint import *
import stdlib
import string
import vipermath
import stdio
import t, c
import memhub


# ============================================================
# Constants
# ============================================================
MAX_NDIM: t.CDefine = 4

float64: t.CTypedef = t.CDouble
float32: t.CTypedef = t.CFloat
int64:   t.CTypedef = t.CLong
int32:   t.CTypedef = t.CInt
uint8:   t.CTypedef = t.CUnsignedChar

pi: t.CDefine = 3.14159265358979323846
e:  t.CDefine = 2.71828182845904523536

# ============================================================
# ndarray - generic n-dimensional array (PEP 695)
#
# Element type T is a numeric type (t.CDouble / t.CFloat / t.CInt / ...).
# T is erased to i8* at the IR level but each ndarray[T] instantiation
# generates specialized code (data access uses T's load/store width).
#
# Default usage: ndarray[t.CDouble] for floating point.
# ============================================================
@t.Object
class ndarray[T]:
    data: T | t.CPtr
    shape: t.CArray[t.CSizeT, MAX_NDIM]
    strides: t.CArray[t.CSizeT, MAX_NDIM]
    ndim: t.CInt
    size: t.CSizeT
    owns_data: t.CInt
    pool: memhub.MemManager | t.CPtr

    def __new__(self, pool: memhub.MemManager | t.CPtr, n: t.CSizeT):
        self.pool = pool
        self.data = pool.alloc(n * T.__sizeof__())
        if self.data:
            memset(self.data, 0, n * T.__sizeof__())
        self.ndim = 1
        self.size = n
        self.shape[0] = n
        self.strides[0] = 1
        self.owns_data = 1

    def __add__(self, other: ndarray[T] | t.CPtr) -> ndarray[T] | t.CPtr:
        result: ndarray[T] | t.CPtr = _empty_like(self)
        if not result: return None
        i: t.CSizeT = 0
        while i < self.size:
            result.data[i] = self.data[i] + other.data[i]
            i += 1
        return result

    def __sub__(self, other: ndarray[T] | t.CPtr) -> ndarray[T] | t.CPtr:
        result: ndarray[T] | t.CPtr = _empty_like(self)
        if not result: return None
        i: t.CSizeT = 0
        while i < self.size:
            result.data[i] = self.data[i] - other.data[i]
            i += 1
        return result

    def __mul__(self, other: ndarray[T] | t.CPtr) -> ndarray[T] | t.CPtr:
        result: ndarray[T] | t.CPtr = _empty_like(self)
        if not result: return None
        i: t.CSizeT = 0
        while i < self.size:
            result.data[i] = self.data[i] * other.data[i]
            i += 1
        return result

    def __truediv__(self, other: ndarray[T] | t.CPtr) -> ndarray[T] | t.CPtr:
        result: ndarray[T] | t.CPtr = _empty_like(self)
        if not result: return None
        i: t.CSizeT = 0
        while i < self.size:
            if other.data[i] != T(0):
                result.data[i] = self.data[i] / other.data[i]
            else:
                result.data[i] = T(0)
            i += 1
        return result

    def __floordiv__(self, other: ndarray[T] | t.CPtr) -> ndarray[T] | t.CPtr:
        result: ndarray[T] | t.CPtr = _empty_like(self)
        if not result: return None
        i: t.CSizeT = 0
        while i < self.size:
            if other.data[i] != T(0):
                result.data[i] = vipermath.floor(self.data[i] / other.data[i])
            else:
                result.data[i] = T(0)
            i += 1
        return result

    def __mod__(self, other: ndarray[T] | t.CPtr) -> ndarray[T] | t.CPtr:
        result: ndarray[T] | t.CPtr = _empty_like(self)
        if not result: return None
        i: t.CSizeT = 0
        while i < self.size:
            if other.data[i] != T(0):
                result.data[i] = vipermath.fmod(self.data[i], other.data[i])
            else:
                result.data[i] = T(0)
            i += 1
        return result

    def __neg__(self) -> ndarray[T] | t.CPtr:
        result: ndarray[T] | t.CPtr = _empty_like(self)
        if not result: return None
        i: t.CSizeT = 0
        while i < self.size:
            result.data[i] = -self.data[i]
            i += 1
        return result

    def __len__(self) -> t.CInt:
        if self.ndim > 0: return t.CInt(self.shape[0])
        return 0

    def at2d(self, row: t.CSizeT, col: t.CSizeT) -> T:
        return self.data[row * self.shape[1] + col]

    def set2d(self, row: t.CSizeT, col: t.CSizeT, val: T):
        self.data[row * self.shape[1] + col] = val

    def delete(self):
        if self.owns_data and self.data:
            self.pool.free(self.data)
            self.data = None
        self.pool.free(self)

    def fill(self, val: T):
        i: t.CSizeT = 0
        while i < self.size:
            self.data[i] = val
            i += 1

    def copy(self) -> ndarray[T] | t.CPtr:
        result: ndarray[T] | t.CPtr = _empty_like(self)
        if not result: return None
        memcpy(result.data, self.data, self.size * T.__sizeof__())
        return result

    def reshape(self, new_shape: INTPTR, new_ndim: t.CInt) -> ndarray[T] | t.CPtr:
        new_size: t.CSizeT = 1
        i: t.CInt = 0
        for i in range(new_ndim):
            new_size *= t.CSizeT(new_shape[i])
        if new_size != self.size: return None
        result: ndarray[T] | t.CPtr = _alloc_ndarray(self.pool)
        if not result: return None
        result.data = self.data
        result.pool = self.pool
        result.ndim = new_ndim
        result.size = self.size
        result.owns_data = 0
        for i in range(new_ndim):
            result.shape[i] = t.CSizeT(new_shape[i])
        _compute_strides(result)
        return result

    def sum(self) -> T:
        s: T = T(0)
        i: t.CSizeT = 0
        while i < self.size:
            s += self.data[i]
            i += 1
        return s

    def mean(self) -> T:
        if self.size == 0: return T(0)
        return self.sum() / T(self.size)

    def min(self) -> T:
        if self.size == 0: return T(0)
        m: T = self.data[0]
        i: t.CSizeT = 1
        while i < self.size:
            if self.data[i] < m: m = self.data[i]
            i += 1
        return m

    def max(self) -> T:
        if self.size == 0: return T(0)
        m: T = self.data[0]
        i: t.CSizeT = 1
        while i < self.size:
            if self.data[i] > m: m = self.data[i]
            i += 1
        return m

    def argmax(self) -> t.CInt:
        if self.size == 0: return -1
        m: T = self.data[0]
        idx: t.CInt = 0
        i: t.CSizeT = 1
        while i < self.size:
            if self.data[i] > m:
                m = self.data[i]
                idx = t.CInt(i)
            i += 1
        return idx

    def argmin(self) -> t.CInt:
        if self.size == 0: return -1
        m: T = self.data[0]
        idx: t.CInt = 0
        i: t.CSizeT = 1
        while i < self.size:
            if self.data[i] < m:
                m = self.data[i]
                idx = t.CInt(i)
            i += 1
        return idx

    def dot(self, other: ndarray[T] | t.CPtr) -> T:
        s: T = T(0)
        n: t.CSizeT = self.size if (self.size < other.size) else other.size
        i: t.CSizeT = 0
        while i < n:
            s += self.data[i] * other.data[i]
            i += 1
        return s

    def T(self) -> ndarray[T] | t.CPtr:
        if self.ndim != 2: return None
        rows: t.CSizeT = self.shape[0]
        cols: t.CSizeT = self.shape[1]
        result: ndarray[T] | t.CPtr = empty2d(self.pool, rows, cols)
        if not result: return None
        # Swap shape for transpose
        result.shape[0] = cols
        result.shape[1] = rows
        _compute_strides(result)
        i: t.CSizeT = 0
        while i < rows:
            j: t.CSizeT = 0
            while j < cols:
                result.data[j * rows + i] = self.data[i * cols + j]
                j += 1
            i += 1
        return result

    def print_arr(self):
        if self.ndim == 1:
            printf("[")
            i: t.CSizeT = 0
            while i < self.size:
                if i > 0: printf(", ")
                printf("%.4f", self.data[i])
                i += 1
            printf("]\n")
        elif self.ndim == 2:
            rows: t.CSizeT = self.shape[0]
            cols: t.CSizeT = self.shape[1]
            printf("[")
            i: t.CSizeT = 0
            while i < rows:
                if i > 0: printf(" ")
                printf("[")
                j: t.CSizeT = 0
                while j < cols:
                    if j > 0: printf(", ")
                    printf("%.4f", self.data[i * cols + j])
                    j += 1
                if i < rows - 1:
                    printf("]\n")
                else:
                    printf("]]\n")
                i += 1


# ============================================================
# Convenience aliases: typed ndarray instantiations
# (must be after class ndarray[T] definition to avoid forward reference)
# ============================================================
Float64Array: t.CTypedef = ndarray[t.CDouble]
Float32Array: t.CTypedef = ndarray[t.CFloat]
Int64Array:   t.CTypedef = ndarray[t.CLong]
Int32Array:   t.CTypedef = ndarray[t.CInt]
Uint8Array:   t.CTypedef = ndarray[t.CUnsignedChar]


# ============================================================
# Internal helpers (specialized for float64 by default)
# Use _alloc_ndarray_t[T] / _empty_like_t[T] for other types
# ============================================================
def _alloc_ndarray(pool: memhub.MemManager | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    a: ndarray[t.CDouble] | t.CPtr = pool.alloc(ndarray[t.CDouble].__sizeof__())
    if a:
        memset(a, 0, ndarray[t.CDouble].__sizeof__())
    return a

def _compute_strides(a: ndarray[t.CDouble] | t.CPtr):
    if a.ndim <= 0: return
    a.strides[a.ndim - 1] = 1
    i: t.CInt = a.ndim - 2
    while i >= 0:
        a.strides[i] = a.strides[i + 1] * a.shape[i + 1]
        i -= 1

def _empty_like(a: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = _alloc_ndarray(a.pool)
    if not result: return None
    result.data = a.pool.alloc(a.size * t.CDouble.__sizeof__())
    if not result.data:
        a.pool.free(result)
        return None
    result.pool = a.pool
    result.ndim = a.ndim
    result.size = a.size
    result.owns_data = 1
    i: t.CInt = 0
    for i in range(a.ndim):
        result.shape[i] = a.shape[i]
    _compute_strides(result)
    return result


# ============================================================
# Array creation functions (float64 default)
# ============================================================
def array(pool: memhub.MemManager | t.CPtr, data: t.CDouble | t.CPtr, n: t.CSizeT) -> ndarray[t.CDouble] | t.CPtr:
    a: ndarray[t.CDouble] | t.CPtr = _alloc_ndarray(pool)
    if not a: return None
    a.data = pool.alloc(n * t.CDouble.__sizeof__())
    if not a.data:
        pool.free(a)
        return None
    memcpy(a.data, data, n * t.CDouble.__sizeof__())
    a.pool = pool
    a.ndim = 1
    a.size = n
    a.shape[0] = n
    a.strides[0] = 1
    a.owns_data = 1
    return a

def zeros(pool: memhub.MemManager | t.CPtr, n: t.CSizeT) -> ndarray[t.CDouble] | t.CPtr:
    a: ndarray[t.CDouble] | t.CPtr = _alloc_ndarray(pool)
    if not a: return None
    a.data = pool.alloc(n * t.CDouble.__sizeof__())
    if not a.data:
        pool.free(a)
        return None
    memset(a.data, 0, n * t.CDouble.__sizeof__())
    a.pool = pool
    a.ndim = 1
    a.size = n
    a.shape[0] = n
    a.strides[0] = 1
    a.owns_data = 1
    return a

def ones(pool: memhub.MemManager | t.CPtr, n: t.CSizeT) -> ndarray[t.CDouble] | t.CPtr:
    a: ndarray[t.CDouble] | t.CPtr = zeros(pool, n)
    if not a: return None
    a.fill(t.CDouble(1.0))
    return a

def full(pool: memhub.MemManager | t.CPtr, n: t.CSizeT, val: t.CDouble) -> ndarray[t.CDouble] | t.CPtr:
    a: ndarray[t.CDouble] | t.CPtr = zeros(pool, n)
    if not a: return None
    a.fill(val)
    return a

def arange(pool: memhub.MemManager | t.CPtr, start: t.CDouble, stop: t.CDouble, step: t.CDouble) -> ndarray[t.CDouble] | t.CPtr:
    if step == t.CDouble(0.0): return None
    n: t.CSizeT = 0
    if step > t.CDouble(0.0):
        v: t.CDouble = start
        while v < stop:
            n += 1
            v += step
    else:
        v: t.CDouble = start
        while v > stop:
            n += 1
            v += step
    if n == 0: return None
    a: ndarray[t.CDouble] | t.CPtr = zeros(pool, n)
    if not a: return None
    i: t.CSizeT = 0
    v = start
    while i < n:
        a.data[i] = v
        v += step
        i += 1
    return a

def linspace(pool: memhub.MemManager | t.CPtr, start: t.CDouble, stop: t.CDouble, num: t.CSizeT) -> ndarray[t.CDouble] | t.CPtr:
    if num == 0: return None
    a: ndarray[t.CDouble] | t.CPtr = zeros(pool, num)
    if not a: return None
    if num == 1:
        a.data[0] = start
        return a
    step: t.CDouble = (stop - start) / t.CDouble(num - 1)
    i: t.CSizeT = 0
    while i < num:
        a.data[i] = start + t.CDouble(i) * step
        i += 1
    return a

def empty2d(pool: memhub.MemManager | t.CPtr, rows: t.CSizeT, cols: t.CSizeT) -> ndarray[t.CDouble] | t.CPtr:
    a: ndarray[t.CDouble] | t.CPtr = _alloc_ndarray(pool)
    if not a: return None
    total: t.CSizeT = rows * cols
    a.data = pool.alloc(total * t.CDouble.__sizeof__())
    if not a.data:
        pool.free(a)
        return None
    memset(a.data, 0, total * t.CDouble.__sizeof__())
    a.pool = pool
    a.ndim = 2
    a.size = total
    a.shape[0] = rows
    a.shape[1] = cols
    a.owns_data = 1
    _compute_strides(a)
    return a

def zeros2d(pool: memhub.MemManager | t.CPtr, rows: t.CSizeT, cols: t.CSizeT) -> ndarray[t.CDouble] | t.CPtr:
    return empty2d(pool, rows, cols)

def ones2d(pool: memhub.MemManager | t.CPtr, rows: t.CSizeT, cols: t.CSizeT) -> ndarray[t.CDouble] | t.CPtr:
    a: ndarray[t.CDouble] | t.CPtr = empty2d(pool, rows, cols)
    if not a: return None
    a.fill(t.CDouble(1.0))
    return a

def eye(pool: memhub.MemManager | t.CPtr, n: t.CSizeT) -> ndarray[t.CDouble] | t.CPtr:
    a: ndarray[t.CDouble] | t.CPtr = empty2d(pool, n, n)
    if not a: return None
    i: t.CSizeT = 0
    while i < n:
        a.data[i * n + i] = t.CDouble(1.0)
        i += 1
    return a

def diag(pool: memhub.MemManager | t.CPtr, vals: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    n: t.CSizeT = vals.size
    a: ndarray[t.CDouble] | t.CPtr = empty2d(pool, n, n)
    if not a: return None
    i: t.CSizeT = 0
    while i < n:
        a.data[i * n + i] = vals.data[i]
        i += 1
    return a


# ============================================================
# Mathematical functions (element-wise, float64 default)
# ============================================================
def np_abs(a: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = _empty_like(a)
    if not result: return None
    i: t.CSizeT = 0
    while i < a.size:
        v: t.CDouble = a.data[i]
        result.data[i] = v if v >= t.CDouble(0.0) else -v
        i += 1
    return result

def np_sqrt(a: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = _empty_like(a)
    if not result: return None
    i: t.CSizeT = 0
    while i < a.size:
        result.data[i] = vipermath.sqrt(a.data[i])
        i += 1
    return result

def np_exp(a: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = _empty_like(a)
    if not result: return None
    i: t.CSizeT = 0
    while i < a.size:
        result.data[i] = vipermath.exp(a.data[i])
        i += 1
    return result

def np_log(a: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = _empty_like(a)
    if not result: return None
    i: t.CSizeT = 0
    while i < a.size:
        result.data[i] = vipermath.log(a.data[i])
        i += 1
    return result

def np_sin(a: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = _empty_like(a)
    if not result: return None
    i: t.CSizeT = 0
    while i < a.size:
        result.data[i] = vipermath.sin(a.data[i])
        i += 1
    return result

def np_cos(a: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = _empty_like(a)
    if not result: return None
    i: t.CSizeT = 0
    while i < a.size:
        result.data[i] = vipermath.cos(a.data[i])
        i += 1
    return result

def np_tan(a: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = _empty_like(a)
    if not result: return None
    i: t.CSizeT = 0
    while i < a.size:
        result.data[i] = vipermath.tan(a.data[i])
        i += 1
    return result

def np_pow(a: ndarray[t.CDouble] | t.CPtr, p: t.CDouble) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = _empty_like(a)
    if not result: return None
    i: t.CSizeT = 0
    while i < a.size:
        result.data[i] = vipermath.pow(a.data[i], p)
        i += 1
    return result


# ============================================================
# Scalar operations (float64 default)
# ============================================================
def add_scalar(a: ndarray[t.CDouble] | t.CPtr, s: t.CDouble) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = _empty_like(a)
    if not result: return None
    i: t.CSizeT = 0
    while i < a.size:
        result.data[i] = a.data[i] + s
        i += 1
    return result

def mul_scalar(a: ndarray[t.CDouble] | t.CPtr, s: t.CDouble) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = _empty_like(a)
    if not result: return None
    i: t.CSizeT = 0
    while i < a.size:
        result.data[i] = a.data[i] * s
        i += 1
    return result

def sub_scalar(a: ndarray[t.CDouble] | t.CPtr, s: t.CDouble) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = _empty_like(a)
    if not result: return None
    i: t.CSizeT = 0
    while i < a.size:
        result.data[i] = a.data[i] - s
        i += 1
    return result

def div_scalar(a: ndarray[t.CDouble] | t.CPtr, s: t.CDouble) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = _empty_like(a)
    if not result: return None
    i: t.CSizeT = 0
    while i < a.size:
        if s != t.CDouble(0.0):
            result.data[i] = a.data[i] / s
        else:
            result.data[i] = t.CDouble(0.0)
        i += 1
    return result


# ============================================================
# Matrix operations (float64 default)
# ============================================================
def matmul(a: ndarray[t.CDouble] | t.CPtr, b: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    # 2D matrix multiply: a(m,k) x b(k,n) = c(m,n)
    if a.ndim != 2 or b.ndim != 2: return None
    m: t.CSizeT = a.shape[0]
    k: t.CSizeT = a.shape[1]
    k2: t.CSizeT = b.shape[0]
    n: t.CSizeT = b.shape[1]
    if k != k2: return None
    c: ndarray[t.CDouble] | t.CPtr = empty2d(a.pool, m, n)
    if not c: return None
    i: t.CSizeT = 0
    while i < m:
        j: t.CSizeT = 0
        while j < n:
            s: t.CDouble = t.CDouble(0.0)
            p: t.CSizeT = 0
            while p < k:
                s += a.data[i * k + p] * b.data[p * n + j]
                p += 1
            c.data[i * n + j] = s
            j += 1
        i += 1
    return c

def dot_product(a: ndarray[t.CDouble] | t.CPtr, b: ndarray[t.CDouble] | t.CPtr) -> t.CDouble:
    return a.dot(b)


# ============================================================
# Reduction operations (float64 default)
# ============================================================
def np_sum(a: ndarray[t.CDouble] | t.CPtr) -> t.CDouble:
    return a.sum()

def np_mean(a: ndarray[t.CDouble] | t.CPtr) -> t.CDouble:
    return a.mean()

def np_min(a: ndarray[t.CDouble] | t.CPtr) -> t.CDouble:
    return a.min()

def np_max(a: ndarray[t.CDouble] | t.CPtr) -> t.CDouble:
    return a.max()

def np_argmax(a: ndarray[t.CDouble] | t.CPtr) -> t.CInt:
    return a.argmax()

def np_argmin(a: ndarray[t.CDouble] | t.CPtr) -> t.CInt:
    return a.argmin()


# ============================================================
# Utility (float64 default)
# ============================================================
def var(a: ndarray[t.CDouble] | t.CPtr) -> t.CDouble:
    if a.size == 0: return t.CDouble(0.0)
    m: t.CDouble = a.mean()
    s: t.CDouble = t.CDouble(0.0)
    i: t.CSizeT = 0
    while i < a.size:
        d: t.CDouble = a.data[i] - m
        s += d * d
        i += 1
    return s / t.CDouble(a.size)

def std(a: ndarray[t.CDouble] | t.CPtr) -> t.CDouble:
    return vipermath.sqrt(var(a))

def norm(a: ndarray[t.CDouble] | t.CPtr) -> t.CDouble:
    return vipermath.sqrt(a.dot(a))

def clip(a: ndarray[t.CDouble] | t.CPtr, lo: t.CDouble, hi: t.CDouble) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = _empty_like(a)
    if not result: return None
    i: t.CSizeT = 0
    while i < a.size:
        v: t.CDouble = a.data[i]
        if v < lo: v = lo
        if v > hi: v = hi
        result.data[i] = v
        i += 1
    return result

def concatenate(a: ndarray[t.CDouble] | t.CPtr, b: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    if a.ndim != 1 or b.ndim != 1: return None
    total: t.CSizeT = a.size + b.size
    result: ndarray[t.CDouble] | t.CPtr = zeros(a.pool, total)
    if not result: return None
    memcpy(result.data, a.data, a.size * t.CDouble.__sizeof__())
    memcpy(result.data + a.size, b.data, b.size * t.CDouble.__sizeof__())
    return result

def sort_arr(a: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = a.copy()
    if not result: return None
    # Simple insertion sort
    i: t.CSizeT = 1
    while i < result.size:
        key: t.CDouble = result.data[i]
        j: t.CSizeT = i - 1
        while j < result.size and result.data[j] > key:
            result.data[j + 1] = result.data[j]
            if j == 0: break
            j -= 1
        if j > 0 or result.data[0] <= key:
            result.data[j + 1] = key
        else:
            result.data[1] = result.data[0]
            result.data[0] = key
        i += 1
    return result

def reverse(a: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = _empty_like(a)
    if not result: return None
    i: t.CSizeT = 0
    while i < a.size:
        result.data[i] = a.data[a.size - 1 - i]
        i += 1
    return result


# ============================================================
# Additional math functions (element-wise, float64 default)
# ============================================================
def np_log10(a: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = _empty_like(a)
    if not result: return None
    i: t.CSizeT = 0
    while i < a.size:
        result.data[i] = vipermath.log10(a.data[i])
        i += 1
    return result

def np_log2(a: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = _empty_like(a)
    if not result: return None
    i: t.CSizeT = 0
    while i < a.size:
        result.data[i] = vipermath.log2(a.data[i])
        i += 1
    return result

def np_floor(a: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = _empty_like(a)
    if not result: return None
    i: t.CSizeT = 0
    while i < a.size:
        result.data[i] = vipermath.floor(a.data[i])
        i += 1
    return result

def np_ceil(a: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = _empty_like(a)
    if not result: return None
    i: t.CSizeT = 0
    while i < a.size:
        result.data[i] = vipermath.ceil(a.data[i])
        i += 1
    return result

def np_round(a: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = _empty_like(a)
    if not result: return None
    i: t.CSizeT = 0
    while i < a.size:
        result.data[i] = vipermath.round(a.data[i])
        i += 1
    return result

def np_sign(a: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = _empty_like(a)
    if not result: return None
    i: t.CSizeT = 0
    while i < a.size:
        if a.data[i] > t.CDouble(0.0):
            result.data[i] = t.CDouble(1.0)
        elif a.data[i] < t.CDouble(0.0):
            result.data[i] = t.CDouble(-1.0)
        else:
            result.data[i] = t.CDouble(0.0)
        i += 1
    return result

def np_tanh(a: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = _empty_like(a)
    if not result: return None
    i: t.CSizeT = 0
    while i < a.size:
        result.data[i] = vipermath.tanh(a.data[i])
        i += 1
    return result

def np_sinh(a: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = _empty_like(a)
    if not result: return None
    i: t.CSizeT = 0
    while i < a.size:
        result.data[i] = vipermath.sinh(a.data[i])
        i += 1
    return result

def np_cosh(a: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = _empty_like(a)
    if not result: return None
    i: t.CSizeT = 0
    while i < a.size:
        result.data[i] = vipermath.cosh(a.data[i])
        i += 1
    return result

def np_arcsin(a: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = _empty_like(a)
    if not result: return None
    i: t.CSizeT = 0
    while i < a.size:
        result.data[i] = vipermath.asin(a.data[i])
        i += 1
    return result

def np_arccos(a: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = _empty_like(a)
    if not result: return None
    i: t.CSizeT = 0
    while i < a.size:
        result.data[i] = vipermath.acos(a.data[i])
        i += 1
    return result

def np_arctan(a: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = _empty_like(a)
    if not result: return None
    i: t.CSizeT = 0
    while i < a.size:
        result.data[i] = vipermath.atan(a.data[i])
        i += 1
    return result

def np_arctan2(y: ndarray[t.CDouble] | t.CPtr, x: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = _empty_like(y)
    if not result: return None
    i: t.CSizeT = 0
    while i < y.size:
        result.data[i] = vipermath.atan2(y.data[i], x.data[i])
        i += 1
    return result

def np_degrees(a: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = _empty_like(a)
    if not result: return None
    i: t.CSizeT = 0
    while i < a.size:
        result.data[i] = vipermath.degrees(a.data[i])
        i += 1
    return result

def np_radians(a: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = _empty_like(a)
    if not result: return None
    i: t.CSizeT = 0
    while i < a.size:
        result.data[i] = vipermath.radians(a.data[i])
        i += 1
    return result

def np_isnan(a: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = _empty_like(a)
    if not result: return None
    i: t.CSizeT = 0
    while i < a.size:
        result.data[i] = t.CDouble(vipermath.isnan(a.data[i]))
        i += 1
    return result

def np_isinf(a: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = _empty_like(a)
    if not result: return None
    i: t.CSizeT = 0
    while i < a.size:
        result.data[i] = t.CDouble(vipermath.isinf(a.data[i]))
        i += 1
    return result

def np_maximum(a: ndarray[t.CDouble] | t.CPtr, b: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = _empty_like(a)
    if not result: return None
    i: t.CSizeT = 0
    while i < a.size:
        result.data[i] = a.data[i] if a.data[i] > b.data[i] else b.data[i]
        i += 1
    return result

def np_minimum(a: ndarray[t.CDouble] | t.CPtr, b: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = _empty_like(a)
    if not result: return None
    i: t.CSizeT = 0
    while i < a.size:
        result.data[i] = a.data[i] if a.data[i] < b.data[i] else b.data[i]
        i += 1
    return result


# ============================================================
# Additional array operations (float64 default)
# ============================================================
def cumsum(a: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = _empty_like(a)
    if not result: return None
    s: t.CDouble = t.CDouble(0.0)
    i: t.CSizeT = 0
    while i < a.size:
        s += a.data[i]
        result.data[i] = s
        i += 1
    return result

def diff(a: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    if a.size < 2: return None
    result: ndarray[t.CDouble] | t.CPtr = zeros(a.pool, a.size - 1)
    if not result: return None
    i: t.CSizeT = 0
    while i < a.size - 1:
        result.data[i] = a.data[i + 1] - a.data[i]
        i += 1
    return result

def flatten(a: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = _empty_like(a)
    if not result: return None
    result.ndim = 1
    result.shape[0] = a.size
    result.strides[0] = 1
    memcpy(result.data, a.data, a.size * t.CDouble.__sizeof__())
    return result

def trace(a: ndarray[t.CDouble] | t.CPtr) -> t.CDouble:
    if a.ndim != 2: return t.CDouble(0.0)
    n: t.CSizeT = a.shape[0] if a.shape[0] < a.shape[1] else a.shape[1]
    s: t.CDouble = t.CDouble(0.0)
    i: t.CSizeT = 0
    while i < n:
        s += a.data[i * a.shape[1] + i]
        i += 1
    return s

def outer(a: ndarray[t.CDouble] | t.CPtr, b: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = empty2d(a.pool, a.size, b.size)
    if not result: return None
    i: t.CSizeT = 0
    while i < a.size:
        j: t.CSizeT = 0
        while j < b.size:
            result.data[i * b.size + j] = a.data[i] * b.data[j]
            j += 1
        i += 1
    return result

def np_where(condition: ndarray[t.CDouble] | t.CPtr, x: ndarray[t.CDouble] | t.CPtr, y: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = _empty_like(x)
    if not result: return None
    i: t.CSizeT = 0
    while i < x.size:
        if condition.data[i] != t.CDouble(0.0):
            result.data[i] = x.data[i]
        else:
            result.data[i] = y.data[i]
        i += 1
    return result

def np_count_nonzero(a: ndarray[t.CDouble] | t.CPtr) -> t.CInt:
    count: t.CInt = 0
    i: t.CSizeT = 0
    while i < a.size:
        if a.data[i] != t.CDouble(0.0):
            count += 1
        i += 1
    return count

def np_all(a: ndarray[t.CDouble] | t.CPtr) -> t.CInt:
    i: t.CSizeT = 0
    while i < a.size:
        if a.data[i] == t.CDouble(0.0):
            return 0
        i += 1
    return 1

def np_any(a: ndarray[t.CDouble] | t.CPtr) -> t.CInt:
    i: t.CSizeT = 0
    while i < a.size:
        if a.data[i] != t.CDouble(0.0):
            return 1
        i += 1
    return 0

def np_equal(a: ndarray[t.CDouble] | t.CPtr, b: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = _empty_like(a)
    if not result: return None
    i: t.CSizeT = 0
    while i < a.size:
        result.data[i] = t.CDouble(1.0) if a.data[i] == b.data[i] else t.CDouble(0.0)
        i += 1
    return result

def np_not_equal(a: ndarray[t.CDouble] | t.CPtr, b: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = _empty_like(a)
    if not result: return None
    i: t.CSizeT = 0
    while i < a.size:
        result.data[i] = t.CDouble(1.0) if a.data[i] != b.data[i] else t.CDouble(0.0)
        i += 1
    return result

def np_less(a: ndarray[t.CDouble] | t.CPtr, b: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = _empty_like(a)
    if not result: return None
    i: t.CSizeT = 0
    while i < a.size:
        result.data[i] = t.CDouble(1.0) if a.data[i] < b.data[i] else t.CDouble(0.0)
        i += 1
    return result

def np_greater(a: ndarray[t.CDouble] | t.CPtr, b: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = _empty_like(a)
    if not result: return None
    i: t.CSizeT = 0
    while i < a.size:
        result.data[i] = t.CDouble(1.0) if a.data[i] > b.data[i] else t.CDouble(0.0)
        i += 1
    return result

def np_less_equal(a: ndarray[t.CDouble] | t.CPtr, b: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = _empty_like(a)
    if not result: return None
    i: t.CSizeT = 0
    while i < a.size:
        result.data[i] = t.CDouble(1.0) if a.data[i] <= b.data[i] else t.CDouble(0.0)
        i += 1
    return result

def np_greater_equal(a: ndarray[t.CDouble] | t.CPtr, b: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = _empty_like(a)
    if not result: return None
    i: t.CSizeT = 0
    while i < a.size:
        result.data[i] = t.CDouble(1.0) if a.data[i] >= b.data[i] else t.CDouble(0.0)
        i += 1
    return result


# ============================================================
# Additional creation functions (float64 default)
# ============================================================
def empty(pool: memhub.MemManager | t.CPtr, n: t.CSizeT) -> ndarray[t.CDouble] | t.CPtr:
    return zeros(pool, n)

def full2d(pool: memhub.MemManager | t.CPtr, rows: t.CSizeT, cols: t.CSizeT, val: t.CDouble) -> ndarray[t.CDouble] | t.CPtr:
    a: ndarray[t.CDouble] | t.CPtr = empty2d(pool, rows, cols)
    if not a: return None
    a.fill(val)
    return a

def zeros_like(a: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = _alloc_ndarray(a.pool)
    if not result: return None
    result.data = a.pool.alloc(a.size * t.CDouble.__sizeof__())
    if not result.data:
        a.pool.free(result)
        return None
    memset(result.data, 0, a.size * t.CDouble.__sizeof__())
    result.pool = a.pool
    result.ndim = a.ndim
    result.size = a.size
    result.owns_data = 1
    i: t.CInt = 0
    for i in range(a.ndim):
        result.shape[i] = a.shape[i]
    _compute_strides(result)
    return result

def ones_like(a: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    result: ndarray[t.CDouble] | t.CPtr = _empty_like(a)
    if not result: return None
    result.fill(t.CDouble(1.0))
    return result

def arange1(pool: memhub.MemManager | t.CPtr, stop: t.CDouble) -> ndarray[t.CDouble] | t.CPtr:
    return arange(pool, t.CDouble(0.0), stop, t.CDouble(1.0))

def linspace2(pool: memhub.MemManager | t.CPtr, start: t.CDouble, stop: t.CDouble) -> ndarray[t.CDouble] | t.CPtr:
    return linspace(pool, start, stop, 50)

def meshgrid(x: ndarray[t.CDouble] | t.CPtr, y: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    # Returns flattened grid pairs - simplified version
    # Returns X grid as 2D array (each row is x)
    result: ndarray[t.CDouble] | t.CPtr = empty2d(x.pool, y.size, x.size)
    if not result: return None
    i: t.CSizeT = 0
    while i < y.size:
        j: t.CSizeT = 0
        while j < x.size:
            result.data[i * x.size + j] = x.data[j]
            j += 1
        i += 1
    return result


# ============================================================
# Linear algebra helpers (float64 default)
# ============================================================
def det2x2(a: ndarray[t.CDouble] | t.CPtr) -> t.CDouble:
    # Determinant of 2x2 matrix
    if a.ndim != 2: return t.CDouble(0.0)
    if a.shape[0] != 2 or a.shape[1] != 2: return t.CDouble(0.0)
    return a.data[0] * a.data[3] - a.data[1] * a.data[2]

def inv2x2(a: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    # Inverse of 2x2 matrix
    d: t.CDouble = det2x2(a)
    if d == t.CDouble(0.0): return None
    result: ndarray[t.CDouble] | t.CPtr = empty2d(a.pool, 2, 2)
    if not result: return None
    result.data[0] = a.data[3] / d
    result.data[1] = -a.data[1] / d
    result.data[2] = -a.data[2] / d
    result.data[3] = a.data[0] / d
    return result

def cross3(a: ndarray[t.CDouble] | t.CPtr, b: ndarray[t.CDouble] | t.CPtr) -> ndarray[t.CDouble] | t.CPtr:
    # Cross product of 3D vectors
    result: ndarray[t.CDouble] | t.CPtr = zeros(a.pool, 3)
    if not result: return None
    result.data[0] = a.data[1] * b.data[2] - a.data[2] * b.data[1]
    result.data[1] = a.data[2] * b.data[0] - a.data[0] * b.data[2]
    result.data[2] = a.data[0] * b.data[1] - a.data[1] * b.data[0]
    return result

def np_interp(x: t.CDouble, xp: ndarray[t.CDouble] | t.CPtr, fp: ndarray[t.CDouble] | t.CPtr) -> t.CDouble:
    # 1D linear interpolation
    if xp.size == 0: return t.CDouble(0.0)
    if x <= xp.data[0]: return fp.data[0]
    if x >= xp.data[xp.size - 1]: return fp.data[xp.size - 1]
    i: t.CSizeT = 0
    while i < xp.size - 1:
        if x >= xp.data[i] and x <= xp.data[i + 1]:
            t_val: t.CDouble = (x - xp.data[i]) / (xp.data[i + 1] - xp.data[i])
            return fp.data[i] + t_val * (fp.data[i + 1] - fp.data[i])
        i += 1
    return fp.data[xp.size - 1]

def prod(a: ndarray[t.CDouble] | t.CPtr) -> t.CDouble:
    p: t.CDouble = t.CDouble(1.0)
    i: t.CSizeT = 0
    while i < a.size:
        p *= a.data[i]
        i += 1
    return p

def median(a: ndarray[t.CDouble] | t.CPtr) -> t.CDouble:
    if a.size == 0: return t.CDouble(0.0)
    s: ndarray[t.CDouble] | t.CPtr = sort_arr(a)
    if not s: return t.CDouble(0.0)
    mid: t.CSizeT = a.size / 2
    if a.size % 2 == 1:
        result: t.CDouble = s.data[mid]
    else:
        result = (s.data[mid - 1] + s.data[mid]) / t.CDouble(2.0)
    s.delete()
    return result

def percentile(a: ndarray[t.CDouble] | t.CPtr, q: t.CDouble) -> t.CDouble:
    if a.size == 0: return t.CDouble(0.0)
    s: ndarray[t.CDouble] | t.CPtr = sort_arr(a)
    if not s: return t.CDouble(0.0)
    idx: t.CDouble = q / t.CDouble(100.0) * t.CDouble(a.size - 1)
    lo: t.CSizeT = t.CSizeT(vipermath.floor(idx))
    hi: t.CSizeT = t.CSizeT(vipermath.ceil(idx))
    frac: t.CDouble = idx - vipermath.floor(idx)
    result: t.CDouble = s.data[lo] + frac * (s.data[hi] - s.data[lo])
    s.delete()
    return result