from stdint import *
import string
import math
import stdio
import memhub


# ============================================================
# Constants
# ============================================================
define MAX_NDIM = 4

typedef float64 = double
typedef float32 = float

define pi = 3.14159265358979323846
define e  = 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.
# ============================================================
class ndarray[T]:
    data: T*
    shape: [MAX_NDIM]sizet
    strides: [MAX_NDIM]sizet
    ndim: int
    size: sizet
    owns_data: int
    pool: memhub.MemManager*

    def __new__(self, pool: memhub.MemManager*, n: sizet):
        self.pool = pool
        self.data = pool.alloc(n * T.__sizeof__())
        if self.data:
            string.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]*) -> ndarray[T]*:
        result: ndarray[T]* = _empty_like(self)
        if not result: return None
        for i in range(int(self.size)):
            result.data[i] = self.data[i] + other.data[i]
        return result

    def __sub__(self, other: ndarray[T]*) -> ndarray[T]*:
        result: ndarray[T]* = _empty_like(self)
        if not result: return None
        for i in range(int(self.size)):
            result.data[i] = self.data[i] - other.data[i]
        return result

    def __mul__(self, other: ndarray[T]*) -> ndarray[T]*:
        result: ndarray[T]* = _empty_like(self)
        if not result: return None
        for i in range(int(self.size)):
            result.data[i] = self.data[i] * other.data[i]
        return result

    def __truediv__(self, other: ndarray[T]*) -> ndarray[T]*:
        result: ndarray[T]* = _empty_like(self)
        if not result: return None
        for i in range(int(self.size)):
            if other.data[i] != T(0):
                result.data[i] = self.data[i] / other.data[i]
            else:
                result.data[i] = T(0)
        return result

    def __floordiv__(self, other: ndarray[T]*) -> ndarray[T]*:
        result: ndarray[T]* = _empty_like(self)
        if not result: return None
        for i in range(int(self.size)):
            if other.data[i] != T(0):
                result.data[i] = math.floor(self.data[i] / other.data[i])
            else:
                result.data[i] = T(0)
        return result

    def __mod__(self, other: ndarray[T]*) -> ndarray[T]*:
        result: ndarray[T]* = _empty_like(self)
        if not result: return None
        for i in range(int(self.size)):
            if other.data[i] != T(0):
                result.data[i] = math.fmod(self.data[i], other.data[i])
            else:
                result.data[i] = T(0)
        return result

    def __neg__(self) -> ndarray[T]*:
        result: ndarray[T]* = _empty_like(self)
        if not result: return None
        for i in range(int(self.size)):
            result.data[i] = -self.data[i]
        return result

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

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

    def set2d(self, row: sizet, col: sizet, 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):
        for i in range(int(self.size)):
            self.data[i] = val

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

    def reshape(self, new_shape: INTPTR, new_ndim: int) -> ndarray[T]*:
        new_size: sizet = 1
        i: int = 0
        for i in range(int(new_ndim)):
            new_size *= sizet(new_shape[i])
        if new_size != self.size: return None
        result: ndarray[T]* = _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(int(new_ndim)):
            result.shape[i] = sizet(new_shape[i])
        _compute_strides(result)
        return result

    def sum(self) -> T:
        s: T = T(0)
        for i in range(int(self.size)):
            s += self.data[i]
        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: sizet = 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: sizet = 1
        while i < self.size:
            if self.data[i] > m: m = self.data[i]
            i += 1
        return m

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

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

    def dot(self, other: ndarray[T]*) -> T:
        s: T = T(0)
        n: sizet = self.size if (self.size < other.size) else other.size
        for i in range(int(n)):
            s += self.data[i] * other.data[i]
        return s

    def T(self) -> ndarray[T]*:
        if self.ndim != 2: return None
        rows: sizet = self.shape[0]
        cols: sizet = self.shape[1]
        result: ndarray[T]* = 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)
        for i in range(int(rows)):
            for j in range(int(cols)):
                result.data[j * rows + i] = self.data[i * cols + j]
        return result

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


# ============================================================
# Convenience aliases: typed ndarray instantiations
# (must be after class ndarray[T] definition to avoid forward reference)
# ============================================================
typedef Float64Array = ndarray[double]
typedef Float32Array = ndarray[float]
typedef Int64Array   = ndarray[int32]
typedef Int32Array   = ndarray[int]
typedef Uint8Array   = ndarray[uchar]


# ============================================================
# 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*) -> ndarray[double]*:
    a: ndarray[double]* = pool.alloc(ndarray[double].__sizeof__())
    if a: string.memset(a, 0, ndarray[double].__sizeof__())
    return a

def _compute_strides(a: ndarray[double]*):
    if a.ndim <= 0: return
    a.strides[a.ndim - 1] = 1
    i: int = a.ndim - 2
    while i >= 0:
        a.strides[i] = a.strides[i + 1] * a.shape[i + 1]
        i -= 1

def _empty_like(a: ndarray[double]*) -> ndarray[double]*:
    result: ndarray[double]* = _alloc_ndarray(a.pool)
    if not result: return None
    result.data = a.pool.alloc(a.size * double.__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: int = 0
    for i in range(int(a.ndim)):
        result.shape[i] = a.shape[i]
    _compute_strides(result)
    return result


# ============================================================
# Array creation functions (float64 default)
# ============================================================
def array(pool: memhub.MemManager*, data: double*, n: sizet) -> ndarray[double]*:
    a: ndarray[double]* = _alloc_ndarray(pool)
    if not a: return None
    a.data = pool.alloc(n * double.__sizeof__())
    if not a.data:
        pool.free(a)
        return None
    string.memcpy(a.data, data, n * double.__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*, n: sizet) -> ndarray[double]*:
    a: ndarray[double]* = _alloc_ndarray(pool)
    if not a: return None
    a.data = pool.alloc(n * double.__sizeof__())
    if not a.data:
        pool.free(a)
        return None
    string.memset(a.data, 0, n * double.__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*, n: sizet) -> ndarray[double]*:
    a: ndarray[double]* = zeros(pool, n)
    if not a: return None
    a.fill(double(1.0))
    return a

def full(pool: memhub.MemManager*, n: sizet, val: double) -> ndarray[double]*:
    a: ndarray[double]* = zeros(pool, n)
    if not a: return None
    a.fill(val)
    return a

def arange(pool: memhub.MemManager*, start: double, stop: double, step: double) -> ndarray[double]*:
    if step == double(0.0): return None
    n: sizet = 0
    if step > double(0.0):
        v: double = start
        while v < stop:
            n += 1
            v += step
    else:
        v: double = start
        while v > stop:
            n += 1
            v += step
    if n == 0: return None
    a: ndarray[double]* = zeros(pool, n)
    if not a: return None
    v = start
    for i in range(int(n)):
        a.data[i] = v
        v += step
    return a

def linspace(pool: memhub.MemManager*, start: double, stop: double, num: sizet) -> ndarray[double]*:
    if num == 0: return None
    a: ndarray[double]* = zeros(pool, num)
    if not a: return None
    if num == 1:
        a.data[0] = start
        return a
    step: double = (stop - start) / double(num - 1)
    for i in range(int(num)):
        a.data[i] = start + double(i) * step
    return a

def empty2d(pool: memhub.MemManager*, rows: sizet, cols: sizet) -> ndarray[double]*:
    a: ndarray[double]* = _alloc_ndarray(pool)
    if not a: return None
    total: sizet = rows * cols
    a.data = pool.alloc(total * double.__sizeof__())
    if not a.data:
        pool.free(a)
        return None
    string.memset(a.data, 0, total * double.__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*, rows: sizet, cols: sizet) -> ndarray[double]*:
    return empty2d(pool, rows, cols)

def ones2d(pool: memhub.MemManager*, rows: sizet, cols: sizet) -> ndarray[double]*:
    a: ndarray[double]* = empty2d(pool, rows, cols)
    if not a: return None
    a.fill(double(1.0))
    return a

def eye(pool: memhub.MemManager*, n: sizet) -> ndarray[double]*:
    a: ndarray[double]* = empty2d(pool, n, n)
    if not a: return None
    for i in range(int(n)):
        a.data[i * n + i] = double(1.0)
    return a

def diag(pool: memhub.MemManager*, vals: ndarray[double]*) -> ndarray[double]*:
    n: sizet = vals.size
    a: ndarray[double]* = empty2d(pool, n, n)
    if not a: return None
    for i in range(int(n)):
        a.data[i * n + i] = vals.data[i]
    return a


# ============================================================
# Mathematical functions (element-wise, float64 default)
# ============================================================
def np_abs(a: ndarray[double]*) -> ndarray[double]*:
    result: ndarray[double]* = _empty_like(a)
    if not result: return None
    for i in range(int(a.size)):
        v: double = a.data[i]
        result.data[i] = v if v >= double(0.0) else -v
    return result

def np_sqrt(a: ndarray[double]*) -> ndarray[double]*:
    result: ndarray[double]* = _empty_like(a)
    if not result: return None
    for i in range(int(a.size)):
        result.data[i] = math.sqrt(a.data[i])
    return result

def np_exp(a: ndarray[double]*) -> ndarray[double]*:
    result: ndarray[double]* = _empty_like(a)
    if not result: return None
    for i in range(int(a.size)):
        result.data[i] = math.exp(a.data[i])
    return result

def np_log(a: ndarray[double]*) -> ndarray[double]*:
    result: ndarray[double]* = _empty_like(a)
    if not result: return None
    for i in range(int(a.size)):
        result.data[i] = math.log(a.data[i])
    return result

def np_sin(a: ndarray[double]*) -> ndarray[double]*:
    result: ndarray[double]* = _empty_like(a)
    if not result: return None
    for i in range(int(a.size)):
        result.data[i] = math.sin(a.data[i])
    return result

def np_cos(a: ndarray[double]*) -> ndarray[double]*:
    result: ndarray[double]* = _empty_like(a)
    if not result: return None
    for i in range(int(a.size)):
        result.data[i] = math.cos(a.data[i])
    return result

def np_tan(a: ndarray[double]*) -> ndarray[double]*:
    result: ndarray[double]* = _empty_like(a)
    if not result: return None
    for i in range(int(a.size)):
        result.data[i] = math.tan(a.data[i])
    return result

def np_pow(a: ndarray[double]*, p: double) -> ndarray[double]*:
    result: ndarray[double]* = _empty_like(a)
    if not result: return None
    for i in range(int(a.size)):
        result.data[i] = math.pow(a.data[i], p)
    return result


# ============================================================
# Scalar operations (float64 default)
# ============================================================
def add_scalar(a: ndarray[double]*, s: double) -> ndarray[double]*:
    result: ndarray[double]* = _empty_like(a)
    if not result: return None
    for i in range(int(a.size)):
        result.data[i] = a.data[i] + s
    return result

def mul_scalar(a: ndarray[double]*, s: double) -> ndarray[double]*:
    result: ndarray[double]* = _empty_like(a)
    if not result: return None
    for i in range(int(a.size)):
        result.data[i] = a.data[i] * s
    return result

def sub_scalar(a: ndarray[double]*, s: double) -> ndarray[double]*:
    result: ndarray[double]* = _empty_like(a)
    if not result: return None
    for i in range(int(a.size)):
        result.data[i] = a.data[i] - s
    return result

def div_scalar(a: ndarray[double]*, s: double) -> ndarray[double]*:
    result: ndarray[double]* = _empty_like(a)
    if not result: return None
    for i in range(int(a.size)):
        if s != double(0.0):
            result.data[i] = a.data[i] / s
        else:
            result.data[i] = double(0.0)
    return result


# ============================================================
# Matrix operations (float64 default)
# ============================================================
def matmul(a: ndarray[double]*, b: ndarray[double]*) -> ndarray[double]*:
    # 2D matrix multiply: a(m,k) x b(k,n) = c(m,n)
    if a.ndim != 2 or b.ndim != 2: return None
    m: sizet = a.shape[0]
    k: sizet = a.shape[1]
    k2: sizet = b.shape[0]
    n: sizet = b.shape[1]
    if k != k2: return None
    c: ndarray[double]* = empty2d(a.pool, m, n)
    if not c: return None
    for i in range(int(m)):
        for j in range(int(n)):
            s: double = double(0.0)
            for p in range(int(k)):
                s += a.data[i * k + p] * b.data[p * n + j]
            c.data[i * n + j] = s
    return c

def dot_product(a: ndarray[double]*, b: ndarray[double]*) -> double:
    return a.dot(b)


# ============================================================
# Reduction operations (float64 default)
# ============================================================
def np_sum(a: ndarray[double]*) -> double:
    return a.sum()

def np_mean(a: ndarray[double]*) -> double:
    return a.mean()

def np_min(a: ndarray[double]*) -> double:
    return a.min()

def np_max(a: ndarray[double]*) -> double:
    return a.max()

def np_argmax(a: ndarray[double]*) -> int:
    return a.argmax()

def np_argmin(a: ndarray[double]*) -> int:
    return a.argmin()


# ============================================================
# Utility (float64 default)
# ============================================================
def var(a: ndarray[double]*) -> double:
    if a.size == 0: return double(0.0)
    m: double = a.mean()
    s: double = double(0.0)
    for i in range(int(a.size)):
        d: double = a.data[i] - m
        s += d * d
    return s / double(a.size)

def std(a: ndarray[double]*) -> double:
    return math.sqrt(var(a))

def norm(a: ndarray[double]*) -> double:
    return math.sqrt(a.dot(a))

def clip(a: ndarray[double]*, lo: double, hi: double) -> ndarray[double]*:
    result: ndarray[double]* = _empty_like(a)
    if not result: return None
    for i in range(int(a.size)):
        v: double = a.data[i]
        if v < lo: v = lo
        if v > hi: v = hi
        result.data[i] = v
    return result

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

def sort_arr(a: ndarray[double]*) -> ndarray[double]*:
    result: ndarray[double]* = a.copy()
    if not result: return None
    # Simple insertion sort
    i: sizet = 1
    while i < result.size:
        key: double = result.data[i]
        j: sizet = 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[double]*) -> ndarray[double]*:
    result: ndarray[double]* = _empty_like(a)
    if not result: return None
    for i in range(int(a.size)):
        result.data[i] = a.data[a.size - 1 - i]
    return result


# ============================================================
# Additional math functions (element-wise, float64 default)
# ============================================================
def np_log10(a: ndarray[double]*) -> ndarray[double]*:
    result: ndarray[double]* = _empty_like(a)
    if not result: return None
    for i in range(int(a.size)):
        result.data[i] = math.log10(a.data[i])
    return result

def np_log2(a: ndarray[double]*) -> ndarray[double]*:
    result: ndarray[double]* = _empty_like(a)
    if not result: return None
    for i in range(int(a.size)):
        result.data[i] = math.log2(a.data[i])
    return result

def np_floor(a: ndarray[double]*) -> ndarray[double]*:
    result: ndarray[double]* = _empty_like(a)
    if not result: return None
    for i in range(int(a.size)):
        result.data[i] = math.floor(a.data[i])
    return result

def np_ceil(a: ndarray[double]*) -> ndarray[double]*:
    result: ndarray[double]* = _empty_like(a)
    if not result: return None
    for i in range(int(a.size)):
        result.data[i] = math.ceil(a.data[i])
    return result

def np_round(a: ndarray[double]*) -> ndarray[double]*:
    result: ndarray[double]* = _empty_like(a)
    if not result: return None
    for i in range(int(a.size)):
        result.data[i] = math.round(a.data[i])
    return result

def np_sign(a: ndarray[double]*) -> ndarray[double]*:
    result: ndarray[double]* = _empty_like(a)
    if not result: return None
    for i in range(int(a.size)):
        if a.data[i] > double(0.0):
            result.data[i] = double(1.0)
        elif a.data[i] < double(0.0):
            result.data[i] = double(-1.0)
        else:
            result.data[i] = double(0.0)
    return result

def np_tanh(a: ndarray[double]*) -> ndarray[double]*:
    result: ndarray[double]* = _empty_like(a)
    if not result: return None
    for i in range(int(a.size)):
        result.data[i] = math.tanh(a.data[i])
    return result

def np_sinh(a: ndarray[double]*) -> ndarray[double]*:
    result: ndarray[double]* = _empty_like(a)
    if not result: return None
    for i in range(int(a.size)):
        result.data[i] = math.sinh(a.data[i])
    return result

def np_cosh(a: ndarray[double]*) -> ndarray[double]*:
    result: ndarray[double]* = _empty_like(a)
    if not result: return None
    for i in range(int(a.size)):
        result.data[i] = math.cosh(a.data[i])
    return result

def np_arcsin(a: ndarray[double]*) -> ndarray[double]*:
    result: ndarray[double]* = _empty_like(a)
    if not result: return None
    for i in range(int(a.size)):
        result.data[i] = math.asin(a.data[i])
    return result

def np_arccos(a: ndarray[double]*) -> ndarray[double]*:
    result: ndarray[double]* = _empty_like(a)
    if not result: return None
    for i in range(int(a.size)):
        result.data[i] = math.acos(a.data[i])
    return result

def np_arctan(a: ndarray[double]*) -> ndarray[double]*:
    result: ndarray[double]* = _empty_like(a)
    if not result: return None
    for i in range(int(a.size)):
        result.data[i] = math.atan(a.data[i])
    return result

def np_arctan2(y: ndarray[double]*, x: ndarray[double]*) -> ndarray[double]*:
    result: ndarray[double]* = _empty_like(y)
    if not result: return None
    for i in range(int(y.size)):
        result.data[i] = math.atan2(y.data[i], x.data[i])
    return result

def np_degrees(a: ndarray[double]*) -> ndarray[double]*:
    result: ndarray[double]* = _empty_like(a)
    if not result: return None
    for i in range(int(a.size)):
        result.data[i] = math.degrees(a.data[i])
    return result

def np_radians(a: ndarray[double]*) -> ndarray[double]*:
    result: ndarray[double]* = _empty_like(a)
    if not result: return None
    for i in range(int(a.size)):
        result.data[i] = math.radians(a.data[i])
    return result

def np_isnan(a: ndarray[double]*) -> ndarray[double]*:
    result: ndarray[double]* = _empty_like(a)
    if not result: return None
    for i in range(int(a.size)):
        result.data[i] = double(math.isnan(a.data[i]))
    return result

def np_isinf(a: ndarray[double]*) -> ndarray[double]*:
    result: ndarray[double]* = _empty_like(a)
    if not result: return None
    for i in range(int(a.size)):
        result.data[i] = double(math.isinf(a.data[i]))
    return result

def np_maximum(a: ndarray[double]*, b: ndarray[double]*) -> ndarray[double]*:
    result: ndarray[double]* = _empty_like(a)
    if not result: return None
    for i in range(int(a.size)):
        result.data[i] = a.data[i] if a.data[i] > b.data[i] else b.data[i]
    return result

def np_minimum(a: ndarray[double]*, b: ndarray[double]*) -> ndarray[double]*:
    result: ndarray[double]* = _empty_like(a)
    if not result: return None
    for i in range(int(a.size)):
        result.data[i] = a.data[i] if a.data[i] < b.data[i] else b.data[i]
    return result


# ============================================================
# Additional array operations (float64 default)
# ============================================================
def cumsum(a: ndarray[double]*) -> ndarray[double]*:
    result: ndarray[double]* = _empty_like(a)
    if not result: return None
    s: double = double(0.0)
    for i in range(int(a.size)):
        s += a.data[i]
        result.data[i] = s
    return result

def diff(a: ndarray[double]*) -> ndarray[double]*:
    if a.size < 2: return None
    result: ndarray[double]* = zeros(a.pool, a.size - 1)
    if not result: return None
    for i in range(int(a.size - 1)):
        result.data[i] = a.data[i + 1] - a.data[i]
    return result

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

def trace(a: ndarray[double]*) -> double:
    if a.ndim != 2: return double(0.0)
    n: sizet = a.shape[0] if a.shape[0] < a.shape[1] else a.shape[1]
    s: double = double(0.0)
    for i in range(int(n)):
        s += a.data[i * a.shape[1] + i]
    return s

def outer(a: ndarray[double]*, b: ndarray[double]*) -> ndarray[double]*:
    result: ndarray[double]* = empty2d(a.pool, a.size, b.size)
    if not result: return None
    for i in range(int(a.size)):
        for j in range(int(b.size)):
            result.data[i * b.size + j] = a.data[i] * b.data[j]
    return result

def np_where(condition: ndarray[double]*, x: ndarray[double]*, y: ndarray[double]*) -> ndarray[double]*:
    result: ndarray[double]* = _empty_like(x)
    if not result: return None
    for i in range(int(x.size)):
        if condition.data[i] != double(0.0):
            result.data[i] = x.data[i]
        else:
            result.data[i] = y.data[i]
    return result

def np_count_nonzero(a: ndarray[double]*) -> int:
    count: int = 0
    for i in range(int(a.size)):
        if a.data[i] != double(0.0):
            count += 1
    return count

def np_all(a: ndarray[double]*) -> int:
    for i in range(int(a.size)):
        if a.data[i] == double(0.0):
            return 0
    return 1

def np_any(a: ndarray[double]*) -> int:
    for i in range(int(a.size)):
        if a.data[i] != double(0.0):
            return 1
    return 0

def np_equal(a: ndarray[double]*, b: ndarray[double]*) -> ndarray[double]*:
    result: ndarray[double]* = _empty_like(a)
    if not result: return None
    for i in range(int(a.size)):
        result.data[i] = double(1.0) if a.data[i] == b.data[i] else double(0.0)
    return result

def np_not_equal(a: ndarray[double]*, b: ndarray[double]*) -> ndarray[double]*:
    result: ndarray[double]* = _empty_like(a)
    if not result: return None
    for i in range(int(a.size)):
        result.data[i] = double(1.0) if a.data[i] != b.data[i] else double(0.0)
    return result

def np_less(a: ndarray[double]*, b: ndarray[double]*) -> ndarray[double]*:
    result: ndarray[double]* = _empty_like(a)
    if not result: return None
    for i in range(int(a.size)):
        result.data[i] = double(1.0) if a.data[i] < b.data[i] else double(0.0)
    return result

def np_greater(a: ndarray[double]*, b: ndarray[double]*) -> ndarray[double]*:
    result: ndarray[double]* = _empty_like(a)
    if not result: return None
    for i in range(int(a.size)):
        result.data[i] = double(1.0) if a.data[i] > b.data[i] else double(0.0)
    return result

def np_less_equal(a: ndarray[double]*, b: ndarray[double]*) -> ndarray[double]*:
    result: ndarray[double]* = _empty_like(a)
    if not result: return None
    for i in range(int(a.size)):
        result.data[i] = double(1.0) if a.data[i] <= b.data[i] else double(0.0)
    return result

def np_greater_equal(a: ndarray[double]*, b: ndarray[double]*) -> ndarray[double]*:
    result: ndarray[double]* = _empty_like(a)
    if not result: return None
    for i in range(int(a.size)):
        result.data[i] = double(1.0) if a.data[i] >= b.data[i] else double(0.0)
    return result


# ============================================================
# Additional creation functions (float64 default)
# ============================================================
def empty(pool: memhub.MemManager*, n: sizet) -> ndarray[double]*:
    return zeros(pool, n)

def full2d(pool: memhub.MemManager*, rows: sizet, cols: sizet, val: double) -> ndarray[double]*:
    a: ndarray[double]* = empty2d(pool, rows, cols)
    if not a: return None
    a.fill(val)
    return a

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

def ones_like(a: ndarray[double]*) -> ndarray[double]*:
    result: ndarray[double]* = _empty_like(a)
    if not result: return None
    result.fill(double(1.0))
    return result

def arange1(pool: memhub.MemManager*, stop: double) -> ndarray[double]*:
    return arange(pool, double(0.0), stop, double(1.0))

def linspace2(pool: memhub.MemManager*, start: double, stop: double) -> ndarray[double]*:
    return linspace(pool, start, stop, 50)

def meshgrid(x: ndarray[double]*, y: ndarray[double]*) -> ndarray[double]*:
    # Returns flattened grid pairs - simplified version
    # Returns X grid as 2D array (each row is x)
    result: ndarray[double]* = empty2d(x.pool, y.size, x.size)
    if not result: return None
    for i in range(int(y.size)):
        for j in range(int(x.size)):
            result.data[i * x.size + j] = x.data[j]
    return result


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

def inv2x2(a: ndarray[double]*) -> ndarray[double]*:
    # Inverse of 2x2 matrix
    d: double = det2x2(a)
    if d == double(0.0): return None
    result: ndarray[double]* = 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[double]*, b: ndarray[double]*) -> ndarray[double]*:
    # Cross product of 3D vectors
    result: ndarray[double]* = 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: double, xp: ndarray[double]*, fp: ndarray[double]*) -> double:
    # 1D linear interpolation
    if xp.size == 0: return double(0.0)
    if x <= xp.data[0]: return fp.data[0]
    if x >= xp.data[xp.size - 1]: return fp.data[xp.size - 1]
    for i in range(int(xp.size - 1)):
        if x >= xp.data[i] and x <= xp.data[i + 1]:
            t_val: double = (x - xp.data[i]) / (xp.data[i + 1] - xp.data[i])
            return fp.data[i] + t_val * (fp.data[i + 1] - fp.data[i])
    return fp.data[xp.size - 1]

def prod(a: ndarray[double]*) -> double:
    p: double = double(1.0)
    for i in range(int(a.size)):
        p *= a.data[i]
    return p

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

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