Metadata-Version: 2.4
Name: numerx
Version: 0.1.2
Summary: Biblioteca de arrays simple para C y Python, inspirada en Numeric/Numarray/NumPy
License: MIT
Classifier: License :: OSI Approved :: MIT License
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: C
Requires-Python: >=3.8
Description-Content-Type: text/markdown
License-File: LICENSE
Dynamic: classifier
Dynamic: description
Dynamic: description-content-type
Dynamic: license
Dynamic: license-file
Dynamic: requires-python
Dynamic: summary

# numerx

`numerx` es una biblioteca de arrays para matemáticas básicas, usable
tanto desde C puro como desde Python. Está inspirada en **Numeric**,
**Numarray** y **NumPy** — toma la simplicidad de los dos primeros, la
ergonomía del tercero, y le añade su propio estilo.

**No busca competir con NumPy.** Es un proyecto hecho por gusto por la
programación, con la meta de tener una implementación pequeña y
entendible de punta a punta (dos archivos: `numerx.c` y `numerx.h`).

## Estado: 0.1.2

Cubre lo básico y algunas operaciones más para cálculos un poco más
elevados:

- **Indexado con tupla `arr[i, j]` (0.1.0)**: lectura y escritura
  directa de un elemento de un array 2D, con la misma sintaxis con
  coma que NumPy — `arr[i, j]` y `arr[i, j] = valor`. Soportado solo
  para 2D en esta versión (`arr[i, j, k]` para 3D queda para una
  versión futura); ambos índices deben ser `int` y admiten valores
  negativos, igual que el indexado 1D ya existente. `arr[i] = valor` y
  `arr[slice] = valor` también son nuevos en 0.1.0 — antes de esta
  versión no existía asignación por índice en absoluto, solo lectura.
  Con un slice, el escalar se repite en cada posición del rango
  (`arr[1:4] = 9.0` deja `9.0` en las tres posiciones) — es
  broadcasting de un único valor, no una asignación posición a
  posición. Asignar una secuencia distinta por posición
  (`arr[0:2] = [1.0, 2.0]`) no está soportado todavía (sería una forma
  de broadcasting más general) y lanza `TypeError` en vez de aceptar
  en silencio solo el primer valor.
- **`.reshape(*dims)` (0.1.0)**: reorganiza los mismos elementos en una
  forma nueva, en el mismo orden row-major — mismo concepto que
  `numpy.reshape`. Acepta tanto argumentos sueltos (`arr.reshape(2,
  3)`) como una tupla explícita (`arr.reshape((2, 3))`), igual que
  NumPy soporta ambas formas. Es **copia, no vista**: numerx todavía no
  tiene el concepto de buffer compartido entre dos `NxArray`, así que
  mutar el resultado de `reshape()` nunca afecta al array original.
  Lanza `ValueError` si el producto de las dimensiones nuevas no
  coincide con la cantidad de elementos del array original.
- **`.T` y `.size` (0.1.0)**: `.T` es la transpuesta de un array 2D
  como propiedad (`arr.T`, sin paréntesis) — usa la misma función que
  `numerx.transpose(arr)` por debajo, expuesta también como atajo de
  atributo; lanza `ValueError` sobre un array 1D, que no tiene
  transpuesta definida. `.size` es el número total de elementos
  (producto de todas las dimensiones del `shape`), disponible para
  cualquier `ndim`.
- **Comparaciones también como método de instancia (0.1.0)**:
  `arr.lt(otro)`, `arr.gt(otro)`, `arr.le(otro)`, `arr.ge(otro)` — el
  mismo resultado que `numerx.lt(arr, otro)` y compañía, pero con
  sintaxis de método en vez de función libre. Ambas formas siguen
  disponibles; ninguna reemplaza a la otra.
- **`eye(n)` (0.1.0)**: matriz identidad de `n`×`n` — `1.0` en la
  diagonal principal, `0.0` en el resto. `eye(0)` da un array 2D de
  shape `(0, 0)`, un resultado válido y no un error (mismo criterio
  que el resto de la biblioteca desde 0.0.8 para el caso vacío). Solo
  la variante cuadrada en esta versión.
- **`arange(...)` (0.1.0)**: rango `[start, stop)` en pasos de `step`
  — mismo contrato que `range()` de Python y `numpy.arange`, exclusivo
  en `stop`. Acepta 1, 2 o 3 argumentos posicionales: `arange(stop)`
  (con `start=0.0`, `step=1.0` implícitos), `arange(start, stop)` (con
  `step=1.0`), o `arange(start, stop, step)` completo. Si la dirección
  del rango pedido no coincide con el signo de `step` (por ejemplo
  `arange(5.0, 2.0)` con el `step=1.0` por defecto), el resultado es un
  array vacío, no un error — igual que `range(5, 2)` en Python. El
  conteo interno de elementos aplica un pequeño margen para absorber
  errores de redondeo de punto flotante en la división (por ejemplo,
  `arange(0.0, 0.3, 0.1)` da exactamente 3 elementos, no 2, aunque
  `0.3/0.1` no sea un entero exacto en IEEE 754).
- **`stack(arrays)` / `concatenate(arrays)` (0.1.0)**: combinan varios
  `numerx.array` 1D ya existentes en uno nuevo — ambas reciben una
  lista o tupla de arrays, no argumentos sueltos (`numerx.stack([a, b,
  c])`, no `numerx.stack(a, b, c)`), mismo criterio que
  `numpy.stack`/`numpy.concatenate`. `concatenate` une las longitudes
  en un array 1D más largo y admite arrays de longitudes distintas;
  `stack` agrega una dimensión nueva, apilando cada array como una
  fila de una matriz 2D, y por eso **sí** exige que todos midan igual.
  Ambas solo soportan arrays 1D en esta versión (2D requeriría elegir
  un eje, una decisión de diseño que 0.1.0 deja para más adelante).
- **`exp(a)` (0.1.0)**: `e^x` elemento a elemento. A diferencia de
  `sqrt`/`abs`/`neg` (vectorizadas con SIMD desde 0.0.8) y de `pow`
  con exponente entero pequeño (SIMD nuevo en esta misma versión, ver
  bullet siguiente), `exp` usa `exp()` de libm sin ningún camino
  SIMD — SSE2/AVX2 no tienen una instrucción de exponencial directa
  (a diferencia de `sqrt`, que sí la tiene), y aproximarla en vectores
  requeriría una implementación propia con su propio análisis de
  error, que queda fuera del alcance de esta versión (ver sección de
  benchmarks para el impacto medido de esto). Nunca da `NaN` para una
  entrada finita: valores grandes desbordan a `+inf`, valores muy
  negativos decaen a `0.0`, siguiendo IEEE 754. *(0.1.2 sí agrega
  SIMD a `exp` — ver el bullet de esa versión más abajo; el resto de
  este bullet, incluyendo el comportamiento en overflow/underflow,
  sigue vigente sin cambios.)*
- **`exp(a)` con SIMD (0.1.2)**: mismo algoritmo que usan NumPy y
  Eigen para esto (`Eigen::internal::pexp`), a su vez descendiente de
  Cephes (Stephen Moshier) — reducción de rango `exp(x) = 2^n *
  exp(r)` con `r` chico, un interpolante racional de grado bajo para
  `exp(r)`, y reconstrucción de `2^n` armando el patrón de bits IEEE
  754 del exponente a mano (no hay instrucción SIMD para eso
  tampoco). Antes de escribir el código con intrínsecos, el algoritmo
  se verificó en Python puro (con la misma manipulación de bits) en
  más de 6.5 millones de puntos contra `exp()` de libm, cubriendo
  desde subnormales hasta el borde de overflow: máximo 2 ULP de
  diferencia, la gran mayoría en 0-1 ULP — mismo orden de precisión
  que NumPy/Eigen documentan para este mismo algoritmo. El `floor()`
  que necesita la reducción de rango se hace con el truco estándar de
  SSE2 puro (truncar a entero + corregir el caso negativo) en vez de
  `_mm_floor_pd`, que es SSE4.1 y este proyecto solo garantiza SSE2
  como mínimo en x86 (ver sección de SIMD más abajo). El camino
  SSE2 y el AVX2 dan resultados bit-idénticos entre sí en la
  verificación hecha (a diferencia de reducciones como `sum`, acá
  cada elemento se calcula de forma independiente, así que el ancho
  del vector no cambia el resultado de ningún elemento individual).
  `exp_inplace` (C) usa el mismo camino.
- **`log(a)` (0.1.1)**: logaritmo natural elemento a elemento, vía
  `log()` de libm. Mismo criterio que `exp` *antes* de 0.1.2: sin
  camino SIMD, misma razón de fondo (ni SSE2 ni AVX2 traen una
  instrucción de logaritmo directa, y vectorizarlo de verdad
  requeriría una aproximación propia con su propio análisis de
  error) — a diferencia de `exp`, `log` **no** recibió el mismo
  tratamiento en 0.1.2 y sigue siendo puramente escalar en esta
  versión. Dominio restringido a `x > 0`, distinto de `exp` (definida
  para todo `x` finito): `log(0.0)` da `-inf`, `log(x)` con `x < 0` da
  `NaN` — mismo criterio de "no abortar, seguir IEEE 754" que ya sigue
  `sqrt` con negativos. `log_inplace` (C) usa el mismo camino.
- **Reducción por eje: `sum`/`max`/`min`/`mean` con `axis` (0.1.1)**:
  hasta esta versión, `sum`/`max`/`min`/`mean` solo reducían un array
  completo a un único `double` (o `float` de Python). Ahora aceptan un
  `axis` opcional que colapsa solo esa dimensión, devolviendo un array
  con una dimensión menos en vez de un escalar — mismo criterio que
  `numpy.sum(a, axis=k)` y compañía. En C son cuatro funciones nuevas
  (`nx_sum_axis`, `nx_max_axis`, `nx_min_axis`, `nx_mean_axis`, ver
  tabla de reducciones más abajo); en Python es el mismo
  `numerx.sum`/`max`/`min`/`mean` de siempre, con `axis=None` como
  valor por defecto que preserva el comportamiento previo a 0.1.1 sin
  cambios. Sin `axis` (o `axis=None`), sigue devolviendo un escalar
  igual que antes; con `axis=k`, devuelve un `numerx.array`. Devuelve
  `NULL` en C (`ValueError` en Python) si `axis` está fuera de rango,
  o si el array tiene alguna dimensión de tamaño 0 — caso distinto al
  de las reducciones totales de arriba, que sí definen `sum` de un
  array vacío como `0.0`; acá `NULL`/`ValueError` es más honesto que
  inventar una forma con ceros. Reducir un array 1D con `axis=0` da un
  array 1D de shape `(1,)`, **no** un escalar de `ndim==0` como sí
  hace NumPy en este caso — divergencia documentada a propósito (ver
  `numerx.h`), no un descuido: `nx_new` rechaza `ndim<=0` desde 0.0.5
  como protección contra underflow en el malloc interno, y no valía
  la pena debilitar esa validación en todo el archivo por este único
  caso borde, cuando `shape=(1,)` representa el mismo valor con la
  misma infraestructura ya existente. No hay `keepdims` en esta
  versión: la dimensión reducida se quita del shape de salida, no se
  deja en 1.
- **`pow` con exponente entero: camino rápido con SIMD (0.1.0)**:
  `pow(a, 2.0)`/`a ** 2` y similares ya no pasan por la `pow()`
  genérica de libm cuando el exponente, convertido a `double`, resulta
  ser un entero exacto dentro de `[-32, 32]` — en ese caso se calcula
  por exponentiación binaria (square-and-multiply, `O(log n)`
  multiplicaciones) usando solo `+`/`*`, que sí vectorizan directo en
  SSE2/AVX2/NEON, con el mismo resultado matemático exacto que
  `pow()`. Esto era, antes de este cambio, la operación menos
  competitiva de toda la biblioteca (~23-27x más lento que NumPy en
  arrays grandes, por hacer dos funciones trascendentes —
  `exp(exponente * log(base))` — por elemento); ahora queda en el
  mismo rango de competitividad que `add`/`mul`/`sqrt` (ver sección de
  benchmarks). Fuera de ese rango de exponente, o con exponente no
  entero, se sigue usando `pow()` de libm sin ninguna pérdida de
  corrección — el límite de `±32` es una decisión de costo/beneficio,
  no de qué exponentes son representables. `pow(x, 0)` siempre da
  `1.0` (incluida la convención `pow(0,0)=1.0` que la propia libm ya
  sigue).
- **Anidación hasta 3D (0.0.8)**: el constructor de Python
  (`numerx.array(...)`) ahora acepta listas anidadas de hasta 3
  niveles (`[[[1,2],[3,4]], [[5,6],[7,8]]]`), no solo 1D/2D. El
  parser se reescribió para ser genuinamente recursivo (antes: una
  bandera booleana con dos ramas de código casi duplicadas, sin
  camino para un tercer nivel) — el mismo código valida y llena 1D,
  2D, y 3D, en vez de un bloque nuevo por cada dimensión. `.tolist()`
  y `repr()` ya eran recursivos sobre `ndim` desde 0.0.7, así que
  funcionan con 3D sin cambios. En C, `NX_ARRAY3D(d0, d1, d2, ...)` +
  `nx_from_array3d` construyen un array 3D en una sola expresión,
  igual que `NX_ARRAY2D` para 2D. 4D o más sigue sin soportarse desde
  listas de Python (error explícito, no un intento silencioso) — el
  struct `NxArray` en sí no tiene ese límite (`nx_new` con C directo
  acepta cualquier `ndim`), es una limitación del parser de listas.
  *(0.1.2 sube este límite del parser de listas de 3 a 4 — ver
  siguiente bullet.)*
- **Anidación hasta 4D (0.1.2)**: mismo parser recursivo de 0.0.8
  (`nx_infer_shape`/`nx_validate_shape_recursive`/`nx_fill_from_list`),
  sin cambios propios — el único tope real era la constante
  `NX_MAX_NDIM` en `numerx.c`, subida de 3 a 4. `.reshape()` usa la
  misma constante para su rango válido de dimensiones, así que también
  acepta hasta 4D desde este cambio. `reduce_axis_generic` (usada por
  `sum`/`max`/`min`/`mean` por eje) y `.tolist()` ya eran genéricas
  sobre `ndim` sin límite propio desde antes — funcionan con 4D sin
  cambios, igual que ya pasó con 3D en 0.0.8. 5D o más sigue sin
  soportarse desde listas de Python por el mismo motivo que 4D no se
  soportaba antes de esta versión (límite de `NX_MAX_NDIM`, no del
  struct `NxArray` en sí, que sigue sin límite de dimensiones
  construyendo con `nx_new` desde C directo).
- **Más álgebra/matemática básica (0.0.8)**: `norm(a)` — norma
  euclidiana (L2) de un vector 1D, reutiliza `dot(a,a)` internamente;
  `clip(a, lo, hi)` — acota cada elemento al rango `[lo, hi]`, con
  variante `_inplace` en C; `linspace(start, stop, n)` — genera `n`
  valores espaciados uniformemente entre `start` y `stop`, inclusive
  en ambos extremos (igual que `numpy.linspace` por defecto), con
  cuidado explícito de que el último valor sea exacto sin arrastre de
  error de punto flotante incluso con `n` muy grande.
- **SIMD también en sqrt/abs/neg (0.0.8)**: además de `add`/`sub`/
  `mul`/`div` (0.0.6), ahora `sqrt`, `abs`, y `-a` también se
  ejecutan vectorizados. `sqrt` tiene instrucción SIMD directa;
  `abs`/`neg` se logran manipulando el bit de signo de IEEE 754 (sin
  instrucción directa para esto en SSE2/AVX2). Medido: `sqrt` pasó de
  ser escalar puro a estar en el mismo rango de competitividad que
  `add`/`mul` frente a NumPy (ver sección de benchmarks). `pow` queda
  fuera de esta ronda (tiene un parámetro extra — el exponente — que
  no encaja en el mismo patrón de las demás unarias) — sí obtiene su
  propio camino SIMD en 0.1.0, con un algoritmo distinto (ver bullet
  más arriba).
- **Organización interna del código (0.0.8)**: sin separar
  `numerx.c`/`numerx.h` en más archivos (siguen siendo los únicos
  dos), se reorganizaron con secciones etiquetadas consistentemente
  (antes: 3 separadores genéricos tipo "PARTE 1"/"PARTE 2" en más de
  1300 líneas del `.c`; ahora: 13 secciones con nombres descriptivos
  — SIMD, ciclo de vida, aritmética, reducciones, álgebra lineal,
  etc. — más una tabla de contenidos al inicio de cada archivo para
  saltar directo con búsqueda de texto). Pensado para que el código
  sea más cómodo de modificar, no como un cambio de arquitectura.
- **Facilidad de uso (0.0.7)**: en Python, `.shape` (tupla con el
  tamaño de cada dimensión) y `.ndim` (número de dimensiones) —
  igual que en NumPy. `.tolist()` convierte a una lista de Python
  respetando la forma del array (1D → lista plana, 2D/3D → listas
  anidadas, recursivo sobre `ndim`) — a diferencia de `list(arr)`, que
  siempre aplana sin importar `ndim` (`list(arr)` se mantiene igual,
  por compatibilidad; `.tolist()` es lo nuevo para cuando se necesita
  preservar la estructura, por ejemplo en un roundtrip
  Python→numerx→Python). En C, `nx_ndim(arr)`/`nx_shape(arr, i)` son
  getters descriptivos (antes había que leer `arr->ndim`/`arr->shape[i]`
  directamente), y la macro `NX_ARRAY2D(rows, cols, ...)` construye un
  array 2D en una sola expresión, análoga a `NX_ARRAY` para 1D.
- **SIMD (0.0.6)**: `add`, `sub`, `mul`, `div` (y desde 0.0.8, también
  `sqrt`/`abs`/`neg`) ahora se ejecutan vectorizados en x86/x86_64
  (SSE2 siempre disponible como mínimo, AVX2 si la CPU lo soporta) y
  en ARM64 (NEON, incluye Android arm64-v8a). La elección de qué
  camino correr se hace en runtime (`nx_simd_level()` /
  `numerx.simd_level()` desde Python) — el mismo binario funciona
  correctamente sin importar qué CPU real lo ejecute, cayendo siempre
  a un camino escalar puro si no hay SIMD disponible o soportado.
  **Estado de verificación honesto** (no todo se verificó igual):
  x86/x86_64 se compiló, ejecutó y midió en el entorno de desarrollo
  de este proyecto; ARM64/ARMv7 (NEON) se escribió siguiendo la
  documentación oficial de ARM con mucho cuidado, pero **nunca se
  compiló ni ejecutó** en ese entorno — no hay compilador cruzado ARM
  ni NDK de Android disponibles ahí. Ver sección "SIMD y Android" más
  abajo para detalles y cómo compilar/verificar esa parte tú mismo.
- **Thread-safety (0.0.6)**: la detección de SIMD se protegió con
  `pthread_once` + atómicos de C11 tras confirmar con ThreadSanitizer
  una condición de carrera real en la primera versión (múltiples hilos
  leyendo/escribiendo la variable de nivel SIMD sin sincronización).
  Relevante para cualquier uso en producción con múltiples hilos
  (apps Android, servicios web, workers paralelos).
- Arrays de `float64`, con `shape`/`ndim`. Desde Python, el
  constructor acepta listas planas (1D), listas anidadas de 2 niveles
  (2D), o de 3 niveles (3D, desde 0.0.8) — todas las listas del mismo
  nivel deben tener el mismo largo, o se lanza un error explícito.
- Operaciones elemento a elemento: `+`, `-`, `*`, `/` (sin broadcasting
  todavía — ambos operandos deben tener la misma shape). La división
  por cero nunca lanza error: sigue IEEE 754 (`x/0.0` → `inf`,
  `0.0/0.0` → `nan`), incluso en los caminos vectorizados.
- Operaciones sobre un solo array: `sqrt`, `abs`, `pow` (exponente
  escalar), y `-a` (negación). En C, cada una tiene una variante
  que reserva un array nuevo y otra `_inplace` que modifica el array
  existente sin `malloc`. En Python se exponen como función libre y
  como operador nativo (`a ** 2`, `abs(a)`, `-a`).
- Reducciones a un solo número: `sum(a)`, `max(a)`, `min(a)`, `mean(a)`,
  `norm(a)` (0.0.8, solo para 1D). `sum` de un array vacío da `0.0`;
  `max`/`min`/`mean`/`norm` de un array vacío o con la dimensión
  incorrecta lanzan `ValueError` en vez de inventar un resultado.
- Comparaciones elemento a elemento: `lt(a,b)`, `gt(a,b)`, `le(a,b)`,
  `ge(a,b)` — devuelven un array de `0.0`/`1.0` (no hay tipo booleano
  en esta versión). Solo como funciones libres, no como operadores
  nativos (`<`, `>`, etc.) — ver nota en `numerx.h` sobre por qué.
- Álgebra: `matmul(a, b)` / `a @ b` — multiplicación de matrices entre
  dos arrays 2D (requiere columnas de `a` == filas de `b`); `dot(a, b)`
  — producto punto entre dos vectores 1D de igual longitud, devuelve
  un `float`; `transpose(a)` — invierte filas y columnas de un array
  2D; `norm(a)` (0.0.8) — norma euclidiana de un vector 1D. Todas son
  operaciones distintas de `*`, que sigue siendo elemento a elemento.
- Matemática básica: `clip(a, lo, hi)` y `linspace(start, stop, n)`
  (0.0.8) — ver arriba.
- `copy(a)` — clon independiente con memoria propia (modificar la
  copia no afecta al original).
- Igualdad (`==`, `!=`) elemento a elemento.
- Indexado con enteros (incluyendo negativos) y slices; sobre 2D
  también con tupla `arr[i, j]`, y con asignación en los tres casos
  (0.1.0 — ver arriba).
- Iteración (`for x in arr`, `list(arr)`).
- Constructores de conveniencia: `zeros(n)`, `ones(n)`.
- **Encadenamiento de operaciones**: como los operadores devuelven un
  `numerx.array` nuevo cada vez, expresiones largas tipo NumPy
  funcionan directo, sin variables intermedias — ver ejemplo abajo.
- **Robustez (0.0.5)**: `nx_new`/`nx_from_data` ahora validan sus
  argumentos (`ndim <= 0`, `shape == NULL`, `data == NULL`) y
  devuelven `NULL` en vez de comportamiento indefinido. Antes de este
  fix, un `ndim` negativo causaba un crash real (`malloc` de tamaño
  gigante por underflow al convertir el entero negativo a `size_t`) —
  encontrado y corregido durante esta versión. También se agregaron
  `nx_get_safe`/`nx_set_safe`, variantes con bounds-check explícito de
  `nx_get`/`nx_set` (que se mantienen sin chequeo, para uso en loops
  calientes donde el índice ya se sabe válido).
- **`NX_ARRAY(...)` (0.0.5)**: macro variádica para construir un array
  1D en una sola expresión en C, sin declarar un `double[]` a mano
  primero — ver sección de uso desde C.
- **`matmul` más rápido (0.0.5)**: además de la optimización de 0.0.4
  (reordenar loops), esta versión compila con `-O3` en vez del `-O2`
  por defecto, lo que habilita auto-vectorización más agresiva del
  compilador. Medido: la brecha contra NumPy en 300×300 bajó de
  ~14-27x a ~8-11x más lento (ver sección de benchmarks).

Lo que **no** tiene esta versión: broadcasting, más de un dtype
(todo es `float64`), arrays de 4 o más dimensiones desde listas de
Python (aunque el struct en sí no tiene ese límite si se construye
desde C), y ningún tipo de bloqueo de cache (tiling) o BLAS en
`matmul` — sigue siendo más lento que NumPy para matrices grandes,
aunque bastante menos que en versiones anteriores (ver sección de
benchmarks). *(0.1.1 sí agrega tiling a `matmul` — este bullet
describe el estado en 0.0.5, ver sección de benchmarks para el
estado actual. 0.1.2 también sube el límite de 3D a 4D desde listas
de Python, ver la entrada de 0.0.8 más abajo para el detalle.)* El in-place de C (`nx_sqrt_inplace`, `nx_exp_inplace`,
etc.) no se expone a Python — ahí `a ** 2` siempre crea un array
nuevo, para no romper la semántica de que las operaciones no mutan
sus operandos. `len(arr)` sobre un array 2D o 3D da el número total
de elementos, no el número de filas ni de "bloques" de la primera
dimensión. `exp` es la única operación de un solo array sin camino
SIMD (ver bullet de 0.1.0 arriba); `stack`/`concatenate`/`arange`
solo cubren 1D; `arr[i, j]` con tupla solo cubre 2D.

## Instalación

```bash
pip install numerx
```

## Uso desde Python

```python
import numerx

a = numerx.array([1.0, 2.0, 3.0])
b = numerx.array([10.0, 20.0, 30.0])

print(a + b)        # numerx.array([11, 22, 33])
print(b - a)         # numerx.array([9, 18, 27])
print(a * b)         # numerx.array([10, 40, 90])
print(b / a)         # numerx.array([10, 10, 10])
print(-a)             # numerx.array([-1, -2, -3])

print(a == a)        # True
print(len(a))        # 3
print(a[0], a[-1])   # 1.0 3.0
print(a[1:3])        # numerx.array([2, 3])

z = numerx.zeros(4)  # numerx.array([0, 0, 0, 0])
o = numerx.ones(3)   # numerx.array([1, 1, 1])

# sqrt, abs, pow: como función libre o como operador nativo
s = numerx.array([4.0, 9.0, 16.0])
print(numerx.sqrt(s))   # numerx.array([2, 3, 4])
print(s ** 0.5)          # numerx.array([2, 3, 4])  (equivalente)

neg = numerx.array([-1.0, 2.0, -3.5])
print(numerx.abs(neg))  # numerx.array([1, 2, 3.5])
print(abs(neg))          # numerx.array([1, 2, 3.5])  (equivalente)

p = numerx.array([2.0, 3.0, 4.0])
print(numerx.pow(p, 2.0))  # numerx.array([4, 9, 16])
print(p ** 2)                # numerx.array([4, 9, 16])  (equivalente)

# sum, max, min: reducen a un solo float, no a un numerx.array
print(numerx.sum(a))   # 6.0
print(numerx.max(a))   # 3.0
print(numerx.min(a))   # 1.0

# sum con axis (0.1.1): reduce solo esa dimensión, da numerx.array
m = numerx.array([[1.0, 2.0], [3.0, 4.0]])
print(numerx.sum(m, axis=0))  # numerx.array([4, 6])  -- suma por columna
print(numerx.sum(m, axis=1))  # numerx.array([3, 7])  -- suma por fila

# arrays 2D + multiplicación de matrices
A = numerx.array([[1.0, 2.0], [3.0, 4.0]])
B = numerx.array([[5.0, 6.0], [7.0, 8.0]])
print(A @ B)              # numerx.array([[19, 22], [43, 50]])
print(numerx.matmul(A, B)) # igual que arriba (equivalente)

# mean: reduce a un solo float, igual que sum/max/min
print(numerx.mean(a))  # 2.0

# dot: producto punto entre dos vectores 1D -> float
v1 = numerx.array([1.0, 2.0, 3.0])
v2 = numerx.array([4.0, 5.0, 6.0])
print(numerx.dot(v1, v2))  # 32.0  (1*4 + 2*5 + 3*6)

# transpose: invierte filas y columnas de un array 2D
M = numerx.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]])
print(numerx.transpose(M))  # numerx.array([[1, 4], [2, 5], [3, 6]])

# copy: clon independiente, modificarlo no afecta al original
original = numerx.array([1.0, 2.0, 3.0])
clon = numerx.copy(original)

# comparaciones elemento a elemento: devuelven 0.0/1.0, no un booleano
x = numerx.array([1.0, 5.0, 3.0])
y = numerx.array([2.0, 2.0, 3.0])
print(numerx.lt(x, y))  # numerx.array([1, 0, 0])
print(numerx.ge(x, y))  # numerx.array([0, 1, 1])

# shape, ndim, tolist: información descriptiva y conversión que
# respeta la estructura del array
grid = numerx.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]])
print(grid.ndim)      # 2
print(grid.shape)     # (2, 3)
print(grid.tolist())  # [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]] -- preserva las filas
print(list(grid))     # [1.0, 2.0, 3.0, 4.0, 5.0, 6.0] -- list() sigue aplanando, sin cambios

# roundtrip real: una lista anidada de Python sobrevive el viaje
# ida y vuelta a través de numerx sin perder su forma
original = [[1.5, 2.5], [3.5, 4.5]]
assert numerx.array(original).tolist() == original

# construcción 3D (0.0.8): listas anidadas de hasta 3 niveles
cube = numerx.array([
    [[1.0, 2.0], [3.0, 4.0]],
    [[5.0, 6.0], [7.0, 8.0]],
])
print(cube.shape)     # (2, 2, 2)
print(cube.tolist())  # preserva los 3 niveles de anidación

# norm, clip, linspace (0.0.8)
v = numerx.array([3.0, 4.0])
print(numerx.norm(v))  # 5.0  (norma euclidiana: sqrt(3^2 + 4^2))

vals = numerx.array([-5.0, 0.0, 3.0, 10.0, 15.0])
print(numerx.clip(vals, 0.0, 10.0))  # numerx.array([0, 0, 3, 10, 10])

print(numerx.linspace(0.0, 1.0, 5))  # numerx.array([0, 0.25, 0.5, 0.75, 1])

# indexado con tupla (0.1.0): lectura y escritura de arr[i, j] en 2D
m = numerx.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]])
print(m[0, 2])       # 3.0
m[0, 0] = 99.0
print(m.tolist())    # [[99.0, 2.0, 3.0], [4.0, 5.0, 6.0]] -- la escritura mutó "m"

# asignación por slice (0.1.0): el escalar se repite en todo el rango
a2 = numerx.array([1.0, 2.0, 3.0, 4.0, 5.0])
a2[1:4] = 9.0
print(a2)  # numerx.array([1, 9, 9, 9, 5])

# reshape (0.1.0): copia, no vista -- el original nunca se muta
v = numerx.array([1.0, 2.0, 3.0, 4.0, 5.0, 6.0])
print(v.reshape(2, 3))     # numerx.array([[1, 2, 3], [4, 5, 6]])
print(v.reshape((2, 3)))   # igual que arriba: args sueltos o tupla explícita, da lo mismo
print(v)                    # numerx.array([1, 2, 3, 4, 5, 6]) -- "v" sigue intacto

# .T y .size (0.1.0)
M = numerx.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]])
print(M.T)     # numerx.array([[1, 4], [2, 5], [3, 6]])
print(M.size)  # 6

# comparaciones también como método de instancia (0.1.0) -- mismo
# resultado que numerx.lt(a, b), sintaxis distinta
a = numerx.array([1.0, 5.0, 3.0])
b = numerx.array([2.0, 2.0, 3.0])
print(a.lt(b))  # numerx.array([1, 0, 0])

# eye, arange (0.1.0)
print(numerx.eye(3))                     # numerx.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
print(numerx.arange(0.0, 1.0, 0.25))     # numerx.array([0, 0.25, 0.5, 0.75])  -- exclusivo en 1.0

# stack, concatenate (0.1.0): combinan arrays 1D ya existentes
s1 = numerx.array([1.0, 2.0, 3.0])
s2 = numerx.array([4.0, 5.0, 6.0])
print(numerx.stack([s1, s2]))         # numerx.array([[1, 2, 3], [4, 5, 6]]) -- agrega una dimensión
print(numerx.concatenate([s1, s2]))   # numerx.array([1, 2, 3, 4, 5, 6])      -- extiende la longitud

# exp (0.1.0): e^x elemento a elemento, sin operador nativo asociado
print(numerx.exp(numerx.array([0.0, 1.0])))  # numerx.array([1, 2.71828])

# log (0.1.1): logaritmo natural elemento a elemento, dominio x > 0
print(numerx.log(numerx.array([1.0, 2.71828])))  # numerx.array([0, 0.999999])

# encadenamiento: funciona directo, sin variables intermedias
signal = numerx.array([1.2, 0.8, 2.5, 1.7, 0.9])
bias   = numerx.array([0.2, 0.1, 0.5, 0.3, 0.2])
gain   = numerx.array([2.0, 1.5, 3.0, 2.5, 1.8])
weight = numerx.array([0.5, 1.0, 0.2, 0.8, 0.6])

score = (((((signal - bias) * gain) ** 2 + weight) * gain - bias) ** 2)
print(score)
```

## Uso desde C

El mismo `numerx.h`/`numerx.c` se puede usar sin Python — basta con
compilar `numerx.c` sin definir `NUMERX_PYTHON`. Requiere enlazar la
biblioteca matemática estándar (`-lm`) desde 0.0.2, por `sqrt`/`pow`.

```c
#include "numerx.h"

int main(void) {
    NxArray *a = nx_ones1d(3);
    NxArray *b = nx_zeros1d(3);
    nx_set(b, 0, 5.0);
    nx_set(b, 1, 6.0);
    nx_set(b, 2, 7.0);

    NxArray *sum = nx_add(a, b);
    nx_print(sum);  // [6, 7, 8]

    /* sqrt/abs/pow: variante nueva (no muta) o _inplace (sin malloc) */
    NxArray *raiz = nx_sqrt(sum);       // sum queda intacto
    nx_print(raiz);                     // [2.44949, 2.64575, 2.82843]
    nx_sqrt_inplace(sum);               // ahora sum SÍ cambia
    nx_print(sum);                      // igual que raiz

    nx_free(a);
    nx_free(b);
    nx_free(sum);
    nx_free(raiz);
    return 0;
}
```

```bash
gcc -o mi_programa numerx.c mi_programa.c -lm -lpthread
```

Desde 0.0.5, construir un array con valores concretos ya no requiere
declarar un `double[]` a mano — `NX_ARRAY` hace ambos pasos en una
expresión:

```c
#include "numerx.h"

int main(void) {
    /* Antes (sigue funcionando, dos pasos): */
    double datos[3] = {1.0, 2.0, 3.0};
    NxArray *a = nx_from_array1d(datos, 3);

    /* Ahora, en una sola expresión: */
    NxArray *b = NX_ARRAY(10.0, 20.0, 30.0);

    NxArray *suma = nx_add(a, b);
    nx_print(suma);  // [11, 22, 33]

    nx_free(a);
    nx_free(b);
    nx_free(suma);
    return 0;
}
```

**Importante**: cada `NX_ARRAY(...)` reserva su propio `NxArray`, así
que sigue haciendo falta `nx_free()` — igual que si se hubiera creado
a mano. Si se usa como argumento directo de otra función, hay que
guardar el puntero en una variable para poder liberarlo después.

Desde 0.0.7, `NX_ARRAY2D` hace lo mismo pero para matrices 2D — se
indican filas y columnas primero, luego los valores en orden
row-major (fila por fila). `nx_ndim`/`nx_shape` son getters
descriptivos para consultar la forma sin acceder al struct
directamente:

```c
#include "numerx.h"
#include <stdio.h>

int main(void) {
    NxArray *m = NX_ARRAY2D(2, 3,  1.0, 2.0, 3.0,
                                   4.0, 5.0, 6.0);
    /* matriz 2x3: [[1,2,3],[4,5,6]] */

    printf("ndim=%d filas=%zu columnas=%zu\n",
           nx_ndim(m), nx_shape(m, 0), nx_shape(m, 1));

    nx_free(m);
    return 0;
}
```

Si el conteo de valores no coincide con `rows*cols`, `NX_ARRAY2D`
devuelve `NULL` en vez de leer memoria fuera de rango — error típico
al contar mal los valores a mano.

Desde 0.0.8, `NX_ARRAY3D` extiende el mismo patrón a 3 dimensiones —
se indican las tres dimensiones primero, luego los valores en orden
row-major (el bloque completo de la primera dimensión, luego el
siguiente, etc.):

```c
#include "numerx.h"

int main(void) {
    NxArray *t = NX_ARRAY3D(2, 2, 2,
        1.0, 2.0, 3.0, 4.0,
        5.0, 6.0, 7.0, 8.0);
    /* array 3D con shape (2,2,2) */

    nx_free(t);
    return 0;
}
```

`norm`, `clip`, `linspace` (0.0.8) siguen el mismo espíritu de las
demás funciones matemáticas:

```c
#include "numerx.h"
#include <stdio.h>

int main(void) {
    NxArray *v = NX_ARRAY(3.0, 4.0);
    int ok;
    printf("norm: %g\n", nx_norm(v, &ok));  // 5.0

    NxArray *clipped = nx_clip(v, 0.0, 3.5);
    nx_print(clipped);  // [3, 3.5]

    NxArray *ls = nx_linspace(0.0, 1.0, 5);
    nx_print(ls);  // [0, 0.25, 0.5, 0.75, 1]

    nx_free(v);
    nx_free(clipped);
    nx_free(ls);
    return 0;
}
```

Desde 0.1.0, `nx_reshape` da una copia independiente con una forma
nueva (mismo contenido, mismo orden row-major — no es una vista);
`nx_eye`/`nx_arange` generan arrays desde cero, igual que
`nx_zeros`/`nx_linspace`; `nx_stack`/`nx_concatenate` combinan arrays
1D ya existentes en uno nuevo; y `nx_exp`/`nx_exp_inplace` siguen el
mismo patrón de variante-que-reserva + variante-en-sitio que
`nx_sqrt`/`nx_abs`/`nx_neg`:

```c
#include "numerx.h"
#include <stdio.h>

int main(void) {
    NxArray *ident = nx_eye(3);
    nx_print(ident);  // [1, 0, 0, 0, 1, 0, 0, 0, 1]  (fila por fila)

    NxArray *r = nx_arange(0.0, 1.0, 0.25);
    nx_print(r);  // [0, 0.25, 0.5, 0.75]  -- exclusivo en 1.0

    NxArray *s1 = NX_ARRAY(1.0, 2.0, 3.0);
    NxArray *s2 = NX_ARRAY(4.0, 5.0, 6.0);
    NxArray *combinados[2] = { s1, s2 };

    NxArray *apilados = nx_stack(combinados, 2);
    nx_print(apilados);  // [1, 2, 3, 4, 5, 6]  -- 2x3, agrega una dimensión

    NxArray *unidos = nx_concatenate(combinados, 2);
    nx_print(unidos);  // [1, 2, 3, 4, 5, 6]  -- 1D de longitud 6

    NxArray *e = nx_exp(s1);       // e^x, s1 queda intacto
    nx_print(e);                    // [2.71828, 7.38906, 20.0855]
    nx_exp_inplace(s1);              // ahora s1 SÍ cambia
    nx_print(s1);                    // igual que e

    size_t nueva_forma[2] = { 3, 1 };
    NxArray *columna = nx_reshape(s2, nueva_forma, 2);
    nx_print(columna);  // [4, 5, 6]  -- mismos datos, shape (3,1)

    nx_free(ident); nx_free(r);
    nx_free(s1); nx_free(s2);
    nx_free(apilados); nx_free(unidos);
    nx_free(e); nx_free(columna);
    return 0;
}
```

### API de C

Funciones de SIMD (0.0.6):

| Función | Qué hace |
|---|---|
| `nx_simd_level()` | detecta (una sola vez, cacheado, thread-safe) qué SIMD soporta la CPU actual: `NX_SIMD_NONE`, `NX_SIMD_SSE2`, `NX_SIMD_AVX2`, o `NX_SIMD_NEON` |
| `nx_simd_force(level)` | fuerza un nivel específico, ignorando la detección — para pruebas/benchmarks; forzar un nivel que la CPU no soporta de verdad puede crashear |

`nx_add`/`nx_sub`/`nx_mul`/`nx_div` usan automáticamente el nivel de
SIMD detectado — no hace falta llamar a `nx_simd_level()` para
beneficiarse, ya se consulta internamente.

Funciones para el caso 1D (el más común, evitan declarar un
`shape[1]` a mano):

| Función | Qué hace |
|---|---|
| `nx_new1d(n)` | reserva un array de `n` elementos, sin inicializar |
| `nx_zeros1d(n)` | array de `n` ceros |
| `nx_ones1d(n)` | array de `n` unos |
| `nx_from_array1d(data, n)` | copia un `double[n]` existente |
| `NX_ARRAY(v1, v2, ...)` | macro (0.0.5): construye el array en una expresión, sin declarar `double[]` a mano |
| `nx_get(arr, i)` / `nx_set(arr, i, v)` | leer/escribir un elemento, sin bounds-check (uso en loops calientes) |
| `nx_get_safe(arr, i, &ok)` / `nx_set_safe(arr, i, v)` | variantes con bounds-check (0.0.5): `ok=0` / devuelve `0` si `i` está fuera de rango, en vez de comportamiento indefinido |
| `nx_ndim(arr)` / `nx_shape(arr, dim)` | getters descriptivos (0.0.7): número de dimensiones / tamaño de la dimensión `dim` (`0` si `dim` está fuera de rango) |
| `nx_print(arr)` | imprime `[a, b, c]` a stdout, sin salto de línea |
| `nx_free(arr)` | libera el array |

Funciones n-dimensionales (para cuando se necesita más de una
dimensión, reciben `shape` como arreglo y `ndim` explícito):

`nx_new`, `nx_zeros`, `nx_ones`, `nx_from_data`, `nx_from_array2d`,
`nx_same_shape`. Desde 0.0.5, `nx_new` valida `ndim <= 0` y
`shape == NULL` (devuelve `NULL` en vez de comportamiento indefinido),
y `nx_from_data` valida `data == NULL` de la misma forma — antes de
este fix, un `ndim` negativo causaba un crash real por underflow de
`size_t` en el cálculo del tamaño a reservar. Desde 0.0.7,
`NX_ARRAY2D(rows, cols, v1, v2, ...)` construye un array 2D en una
expresión (análoga a `NX_ARRAY` para 1D) — internamente valida que
`rows*cols` coincida con la cantidad de valores dados, devolviendo
`NULL` si no calzan. Desde 0.1.0, `nx_reshape(a, new_shape, new_ndim)`
da una **copia** independiente de `a` reorganizada en `new_shape` —
mismos datos, mismo orden row-major, nunca una vista sobre `a` (numerx
no tiene todavía el concepto de buffer compartido entre dos
`NxArray`); devuelve `NULL` si el producto de `new_shape` no coincide
con `a->size`.

Operaciones binarias (misma shape, sin broadcasting): `nx_add`, `nx_sub`,
`nx_mul`, `nx_div` — devuelven un `NxArray*` nuevo, o `NULL` si las
shapes no coinciden.

Comparaciones elemento a elemento (mismo patrón que las binarias,
devuelven un array de `0.0`/`1.0`): `nx_lt`, `nx_gt`, `nx_le`, `nx_ge`.

`nx_copy(a)` — clon independiente, memoria propia (declarada desde
0.0.1 en el header pero implementada recién en 0.0.4).

Operaciones sobre un solo array, cada una con variante que reserva
(`nx_X`) y variante en sitio (`nx_X_inplace`):

| Reserva nuevo | En sitio | Qué hace |
|---|---|---|
| `nx_sqrt(a)` | `nx_sqrt_inplace(a)` | raíz cuadrada (negativos → NaN) |
| `nx_abs(a)` | `nx_abs_inplace(a)` | valor absoluto |
| `nx_pow(a, e)` | `nx_pow_inplace(a, e)` | potencia, `e` es un `double` escalar — SIMD (0.1.0) si `e` es entero exacto en `[-32, 32]` (ver nota abajo), libm genérica fuera de ese rango |
| `nx_neg(a)` | `nx_neg_inplace(a)` | negación (`-a`) |
| `nx_exp(a)` | `nx_exp_inplace(a)` | `e^x` (0.1.0) — vía libm, con SIMD desde 0.1.2 (ver nota abajo) |
| `nx_log(a)` | `nx_log_inplace(a)` | logaritmo natural, dominio `x > 0` (0.1.1) — vía libm, **sin** SIMD (ver nota abajo) |

Desde 0.0.8, `sqrt`/`abs`/`neg` también usan SIMD internamente (mismo
dispatch en runtime que `add`/`sub`/`mul`/`div` desde 0.0.6) — sin
cambio en la firma ni el comportamiento, solo más rápido. Desde 0.1.0,
`pow` también vectoriza, pero solo para el caso común: si el
exponente, convertido a `double`, resulta ser un entero exacto dentro
de `[-32, 32]`, se calcula por exponentiación binaria
(square-and-multiply) usando solo `+`/`*` — que sí vectorizan directo
— en vez de la `pow()` genérica de libm (que internamente hace
`exp(exponente * log(base))`, dos funciones trascendentes sin
instrucción SIMD directa). Mismo resultado matemático exacto que
`pow()` en ese rango, medido ~23-27x más lento que NumPy antes del
camino rápido, ahora en el mismo orden que `add`/`mul`/`sqrt` (ver
sección de benchmarks). Fuera de `[-32, 32]`, o con exponente no
entero, sigue siendo la `pow()` genérica de libm sin ningún camino
SIMD, exactamente como `exp` (0.1.0): SSE2/AVX2 no tienen una
instrucción de exponencial/logaritmo directa (a diferencia de `sqrt`,
que sí la tiene), y aproximarlas en vectores requeriría una
implementación propia con su propio análisis de error, que queda
fuera del alcance de esta versión. Ver sección de benchmarks para el
impacto medido de ambos casos frente a NumPy.

Reducciones (de un array completo a un solo `double`):

| Función | Qué hace |
|---|---|
| `nx_sum(a)` | suma de todos los elementos (array vacío → `0.0`) |
| `nx_max(a, &ok)` | máximo elemento; `ok=0` si `a` está vacío |
| `nx_min(a, &ok)` | mínimo elemento; `ok=0` si `a` está vacío |
| `nx_mean(a, &ok)` | promedio; `ok=0` si `a` está vacío (evita dividir por 0) |
| `nx_norm(a, &ok)` | norma euclidiana de un vector 1D; `ok=0` si `a` no es 1D (0.0.8) |

Reducción por eje (0.1.1) — a diferencia de la tabla de arriba, estas
colapsan solo una dimensión y devuelven un `NxArray *`, no un
`double`:

| Función | Qué hace |
|---|---|
| `nx_sum_axis(a, axis)` | suma a lo largo de `axis`; shape de salida = shape de `a` sin esa dimensión |
| `nx_max_axis(a, axis)` | máximo a lo largo de `axis`, mismo criterio de shape |
| `nx_min_axis(a, axis)` | mínimo a lo largo de `axis`, mismo criterio de shape |
| `nx_mean_axis(a, axis)` | promedio a lo largo de `axis`, mismo criterio de shape |

`axis` se cuenta desde 0. Las cuatro devuelven `NULL` si `a->ndim < 1`,
si `axis` está fuera de rango (`axis < 0` o `axis >= a->ndim`), o si
`a` tiene alguna dimensión de tamaño 0 (un array así no tiene ningún
elemento que reducir en ese eje — caso distinto al de las reducciones
totales de arriba, que sí definen `sum` de vacío como `0.0`; acá
`NULL` es más honesto que inventar una forma con ceros). Reducir un
array 1D con `axis=0` da un `NxArray` 1D de shape `(1,)`, **no** un
array de `ndim==0` como sí da NumPy en este mismo caso — la única
divergencia de `axis` respecto a NumPy, documentada a propósito en
`numerx.h`: `nx_new` rechaza `ndim<=0` desde 0.0.5 como protección
contra underflow de `size_t` en el malloc interno, y no vale la pena
debilitar esa validación en todo el archivo por este único caso
borde, cuando `shape=(1,)` representa el mismo valor con la misma
infraestructura ya existente. Si lo que se quiere es el `double`
suelto, seguir usando `nx_sum`/`max`/`min`/`mean` sin `axis`. Del lado
de Python, no hay funciones separadas: `axis` es un kwarg opcional de
`numerx.sum`/`max`/`min`/`mean` (`axis=None` por defecto, mismo
comportamiento que antes de 0.1.1); con `axis=k` devuelven
`numerx.array` en vez de `float`, y `ValueError` en vez de `NULL`
cuando `axis` es inválido.

Álgebra:

- `nx_matmul(a, b)` — multiplicación de matrices entre dos arrays 2D
  (`a->shape[1]` debe ser igual a `b->shape[0]`). Devuelve `NULL` si
  las dimensiones no calzan o si `a`/`b` no son 2D. Distinta de
  `nx_mul`, que sigue siendo elemento a elemento. **Optimizada en
  0.0.4**: el orden de los loops internos es `i-k-j` en vez del
  `i-j-k` de libro de texto, lo que hace que el acceso a memoria de
  `b` sea contiguo fila por fila en vez de saltar de fila en fila —
  mejora medida de ~28x a ~15.5x más lento que NumPy en 300×300 (ver
  sección de benchmarks). El código no salta términos donde el
  multiplicador es `0.0`: aunque sumar `0.0*x` normalmente no cambia
  nada, si `x` es `NaN` o `Inf` sí importaría (`0.0*NaN = NaN` según
  IEEE 754) — se prefirió mantener la correctitud exacta en ese caso
  raro sobre esa optimización adicional.
- `nx_dot(a, b, &ok)` — producto punto entre dos vectores 1D de igual
  longitud, devuelve un `double`. `ok=0` si `a`/`b` no son 1D o
  difieren en longitud.
- `nx_transpose(a)` — invierte filas y columnas de un array 2D.
  Devuelve `NULL` si `a` no es 2D.
- `nx_norm(a, &ok)` (0.0.8) — norma euclidiana (L2) de un vector 1D;
  reutiliza `nx_dot(a, a, &ok)` internamente (`sqrt(sum(x_i^2))`).

Matemática básica:

- `nx_clip(a, lo, hi)` / `nx_clip_inplace(a, lo, hi)` — acota cada
  elemento al rango `[lo, hi]`. Mismo patrón que las unarias de
  arriba (variante que reserva + variante en sitio).
- `nx_linspace(start, stop, n)` — genera un array 1D de `n` valores
  espaciados uniformemente entre `start` y `stop`, inclusive en
  ambos extremos. Devuelve `NULL` si `n == 0`. Con `n == 1`, devuelve
  `[start]` (`stop` se ignora, igual que `numpy.linspace`). El
  último valor se fuerza explícitamente a `stop` en vez de dejarlo a
  la acumulación de sumas de punto flotante, así que
  `nx_linspace(0, 1, N)` siempre termina exacto en `1.0` sin importar
  qué tan grande sea `N`.
- `nx_arange(start, stop, step)` (0.1.0) — rango `[start, stop)`
  exclusivo en `stop`, mismo contrato que `range()` de Python y
  `numpy.arange` — a propósito distinto de `nx_linspace`, que sí
  incluye el extremo final y recibe "cuántos puntos quiero" en vez de
  "de cuánto es cada paso". Devuelve `NULL` si `step == 0.0`; si la
  dirección del rango no coincide con el signo de `step`, el resultado
  es un array de tamaño `0` (válido, no un error). El conteo interno
  de elementos es `ceil(raw_count - epsilon)`, con `raw_count =
  (stop-start)/step` y un `epsilon` relativo pequeño restado *antes*
  del `ceil` — absorbe el ruido de redondeo típico de IEEE 754 en la
  división sin llegar a esconder un `ceil()` genuinamente distinto
  (ver el comentario junto a la implementación).
- `nx_eye(n)` (0.1.0) — matriz identidad `n`×`n`. `n == 0` da un array
  2D de shape `(0, 0)` válido, no `NULL`. Solo la variante cuadrada en
  esta versión.

Construcción y combinación (0.1.0) — a diferencia de las funciones de
arriba (que generan un array desde cero), estas combinan `NxArray`
**ya existentes**:

- `nx_concatenate(arrays, count)` — une `count` arrays 1D en uno solo,
  en el orden dado; resultado 1D de longitud igual a la suma de todas
  las longitudes. Admite arrays de longitudes distintas. Devuelve
  `NULL` si `count == 0`, si `arrays == NULL`, si algún elemento es
  `NULL`, o si algún array no es 1D — 2D requeriría elegir un eje
  (`axis=0` vs `axis=1`), una decisión de diseño que 0.1.0 deja para
  más adelante en vez de resolver a medias.
- `nx_stack(arrays, count)` — apila `count` arrays 1D de igual
  longitud `len` como filas de una matriz 2D nueva de shape `(count,
  len)` — a diferencia de `nx_concatenate` (que extiende la longitud
  sin agregar dimensión), esto agrega una dimensión nueva. Mismo
  espíritu que `numpy.stack(axis=0)` sobre arrays 1D. Devuelve `NULL`
  con las mismas condiciones que `nx_concatenate`, más si no todos los
  arrays comparten la misma longitud — no hay una forma sensata de
  rellenar el faltante, así que se rechaza en vez de adivinar.

## Benchmarks

Números reales, no estimaciones. Corridos en el entorno de
desarrollo — tus resultados variarán según CPU y build de NumPy.

### Desde Python (`bench_python.py`)

Suma y multiplicación elemento a elemento, arrays de `float64`. Desde
0.0.6, ambos caminos usan SIMD (ver sección dedicada más arriba):

```
n=      10  add: numerx=     0.10us  numpy=     0.39us  (numerx es    0.3x numpy)   mul: numerx=     0.10us  numpy=     0.39us  (numerx es    0.3x numpy)
n=     100  add: numerx=     0.11us  numpy=     0.50us  (numerx es    0.2x numpy)   mul: numerx=     0.12us  numpy=     0.40us  (numerx es    0.3x numpy)
n=    1000  add: numerx=     0.22us  numpy=     0.55us  (numerx es    0.4x numpy)   mul: numerx=     0.34us  numpy=     0.55us  (numerx es    0.6x numpy)
n=   10000  add: numerx=     3.24us  numpy=     5.83us  (numerx es    0.6x numpy)   mul: numerx=     3.16us  numpy=     6.58us  (numerx es    0.5x numpy)
n=  100000  add: numerx=    61.91us  numpy=    66.90us  (numerx es    0.9x numpy)   mul: numerx=    59.25us  numpy=    67.50us  (numerx es    0.9x numpy)
n= 1000000  add: numerx=   910.38us  numpy=   915.78us  (numerx es    1.0x numpy)   mul: numerx=   901.97us  numpy=   906.05us  (numerx es    1.0x numpy)
```

Lectura honesta: antes de 0.0.6 (sin SIMD), en 1M elementos numerx
estaba ~1.1-1.3x más lento que NumPy en `add`. Con SIMD (AVX2 en esta
máquina), ahora está prácticamente empatado en todos los tamaños, y es
más rápido que NumPy en el rango de 10 a 100k elementos — probablemente
porque numerx no paga el overhead de despacho de dtype que NumPy sí
paga por soportar muchos tipos, y ahora ambos usan vectorización real.

`sqrt`, vectorizada desde 0.0.8 (antes era escalar puro):

```
n=      10  sqrt: numerx=     0.31us  numpy=     0.43us  (numerx es   0.7x numpy)
n=     100  sqrt: numerx=     0.40us  numpy=     0.50us  (numerx es   0.8x numpy)
n=    1000  sqrt: numerx=     1.33us  numpy=     1.46us  (numerx es   0.9x numpy)
n=   10000  sqrt: numerx=    12.70us  numpy=     9.86us  (numerx es   1.3x numpy)
n=  100000  sqrt: numerx=   143.36us  numpy=    99.97us  (numerx es   1.4x numpy)
n= 1000000  sqrt: numerx=  1696.98us  numpy=  1005.45us  (numerx es   1.7x numpy)
```

En el mismo rango de competitividad que `add`/`mul` — antes de
vectorizarla, `sqrt` habría sido dramáticamente más lenta que el
`sqrt` ya vectorizado de NumPy en arrays grandes, por ser un loop
escalar puro contra código SIMD.

`exp` (0.1.0) — a propósito **sin** vectorizar (ver la nota junto a
`nx_exp` en la tabla de API de C):

```
n=      10  exp: numerx=     0.46us  numpy=     0.43us  (numerx es    1.1x numpy)
n=     100  exp: numerx=     0.65us  numpy=     0.47us  (numerx es    1.4x numpy)
n=    1000  exp: numerx=     5.72us  numpy=     1.79us  (numerx es    3.2x numpy)
n=   10000  exp: numerx=    47.47us  numpy=    10.55us  (numerx es    4.5x numpy)
n=  100000  exp: numerx=   487.32us  numpy=   111.00us  (numerx es    4.4x numpy)
n= 1000000  exp: numerx=  5180.36us  numpy=  1080.03us  (numerx es    4.8x numpy)
```

Lectura honesta: acá sí se nota la diferencia de no vectorizar. En
arrays chicos (`n<=100`) `exp` queda cerca de NumPy, dominado por el
mismo overhead fijo por llamada que se ve en el resto de las
funciones. Pero a partir de `n=1000` la brecha crece y se estabiliza
alrededor de ~4-5x más lento — bien por encima del ~1.3-1.7x de
`sqrt`, que sí es SIMD. Es exactamente la brecha que se espera de
libm escalar puro contra un loop vectorizado: no es un bug ni una
implementación descuidada, es el costo de no haber escrito una
aproximación de exponencial vectorizada propia en esta versión (ver
la nota de diseño junto a `nx_exp` en `numerx.c`).

`pow(a, 2)` (0.1.0) — el contraste directo con `exp`: mismo tipo de
operación unaria, pero esta **sí** tiene camino SIMD porque el
exponente es un entero pequeño (ver la nota junto a `nx_pow` en la
tabla de API de C):

```
n=      10  pow(a,2): numerx=     0.20us  numpy=     0.48us  (numerx es    0.4x numpy)
n=     100  pow(a,2): numerx=     0.36us  numpy=     0.50us  (numerx es    0.7x numpy)
n=    1000  pow(a,2): numerx=     0.58us  numpy=     0.95us  (numerx es    0.6x numpy)
n=   10000  pow(a,2): numerx=     6.10us  numpy=     2.35us  (numerx es    2.6x numpy)
n=  100000  pow(a,2): numerx=    71.37us  numpy=    43.46us  (numerx es    1.6x numpy)
n= 1000000  pow(a,2): numerx=   961.05us  numpy=   564.95us  (numerx es    1.7x numpy)
```

Lectura honesta: la mejora frente al estado anterior a 0.1.0 es
dramática — `pow(a, 2.0)` con la `pow()` genérica de libm medía
~23-27x más lento que NumPy en arrays grandes (la operación menos
competitiva de toda la biblioteca en ese momento); con el camino
rápido de exponentiación binaria, ahora queda en el mismo rango que
`add`/`mul`/`sqrt`, prácticamente empatado con NumPy en arrays chicos
y medianos, ~1.6-1.7x en el millón de elementos. La comparación con
`exp` de arriba deja claro de dónde sale la diferencia: no es que
`pow` sea una operación más simple que `exp` en sí misma, es que el
caso común de `pow` (exponente entero pequeño) admite un algoritmo
distinto (`+`/`*` puros) que sí vectoriza, mientras que `exp` no tiene
ese atajo — está genuinamente atado a una función trascendente por
elemento sin instrucción SIMD directa.

### matmul: numerx vs NumPy (matrices cuadradas n×n)

La brecha se redujo considerablemente respecto a versiones previas,
por dos cambios acumulados: el reordenamiento de loops de 0.0.4
(`i-k-j` en vez de `i-j-k`, acceso a memoria contiguo), y desde 0.0.5,
compilar con `-O3` en vez del `-O2` por defecto — lo que permite al
compilador auto-vectorizar el loop interno de forma más agresiva.
NumPy sigue delegando esta operación a BLAS, una biblioteca de
décadas de madurez específicamente optimizada para multiplicación de
matrices (bloqueo de cache más agresivo, y según el build, SIMD y/o
múltiples núcleos) — por eso sigue ganando, pero por bastante menos
que antes.

```
n=    2x2      matmul: numerx=        0.20us  numpy=      1.59us  (numerx es      0.1x numpy)
n=   10x10     matmul: numerx=        0.92us  numpy=      2.87us  (numerx es      0.3x numpy)
n=   50x50     matmul: numerx=       50.07us  numpy=      7.01us  (numerx es      7.1x numpy)
n=  100x100    matmul: numerx=      377.08us  numpy=     51.63us  (numerx es      7.3x numpy)
n=  300x300    matmul: numerx=     9942.17us  numpy=   1239.10us  (numerx es      8.0x numpy)
```

Lectura honesta: en matrices muy chicas numerx sigue ganando o queda
prácticamente empatado, igual que antes. La mejora real está en
matrices medianas y grandes: **0.0.3 (sin optimizar) daba ~28x más
lento en 300×300; 0.0.4 (reordenamiento de loops) bajó eso a ~14.6x;
0.0.5 (agregando `-O3`) lo baja más, a ~8x** — sin cambiar el
algoritmo ni el resultado matemático en ningún punto de esta cadena.
Esta máquina muestra variación entre corridas (~7x-11x según el
momento, por ser un entorno compartido), pero la tendencia de mejora
acumulada es consistente. Sigue siendo más lento que BLAS.

**0.1.1 agrega bloqueo de cache (tiling) a `nx_matmul`** (ver la nota
grande junto a `nx_matmul` en `numerx.c` para el detalle completo),
pero esta tabla en particular NO lo muestra: el tiling solo aporta
una ganancia clara en matrices bastante más grandes que 300×300 (a
partir de ~512×512 en las mediciones hechas al implementarlo) — en
el rango 2×2 a 300×300 que cubre este benchmark, la diferencia contra
la versión sin tiling es marginal y se pierde dentro de la propia
variación de esta máquina compartida, así que los números de arriba
no cambiaron de forma visible tras agregarlo. Cerrar la brecha
restante contra BLAS de verdad (en cualquier rango de tamaño)
requeriría microkernels y empaquetado de datos (packing) más
elaborados que un tiling simple, que quedan fuera del alcance de
esta versión.

### Desde C puro (`bench_c.c`)

Esto mide el overhead propio de numerx (el `malloc` del struct
`NxArray` en cada operación) contra un loop crudo equivalente sin
ningún wrapping:

```
n              numerx(us)        raw(us)   overhead
10                  0.067          0.007      9.93x
100                 0.212          0.070      3.03x
1000                1.655          0.612      2.70x
10000              16.145          7.715      2.09x
100000            164.839         94.764      1.74x
1000000          2074.036       1363.173      1.52x
```

El overhead relativo baja según crece el array, porque el costo fijo
de alocar el struct se diluye entre más trabajo real.

Para reproducir:

```bash
python3 bench_python.py
gcc -O2 -o bench_c numerx.c bench_c.c -lm -lpthread && ./bench_c
```

## SIMD: arquitecturas y estado de verificación

numerx detecta en runtime qué instrucciones SIMD soporta la CPU actual
y elige el camino más rápido disponible, cayendo siempre a un camino
escalar puro (el original, sin SIMD) si nada más aplica. Esto es
consulta consciente, no un compromiso a ciegas: el mismo `.so`
compilado funciona correctamente en cualquier CPU de la misma familia
de arquitectura, aunque tenga menos extensiones de las que la máquina
donde se compiló tenía disponibles.

```python
import numerx
print(numerx.simd_level())  # 'none', 'sse2', 'avx2', o 'neon'
```

**Siendo honestos sobre qué se verificó y qué no**, porque no todo se
comprobó con el mismo rigor:

| Arquitectura | SIMD | ¿Compilado y ejecutado en este proyecto? |
|---|---|---|
| x86_64 | SSE2, AVX2 | **Sí** — compilado, ejecutado, medido con benchmarks reales, y verificado con AddressSanitizer y ThreadSanitizer |
| x86 (32-bit) | SSE2 | **Sí** — compilado y ejecutado |
| ARM64 (incluye Android arm64-v8a) | NEON | **No** — escrito siguiendo la documentación oficial de intrínsecos ARM (`<arm_neon.h>`), pero nunca compilado ni ejecutado: no hay compilador cruzado ARM64 ni NDK de Android disponibles en el entorno donde se desarrolló esta versión |
| ARMv7 (incluye Android armeabi-v7a) | — | **No implementado**: por seguridad, esta versión no asume soporte de NEON en ARMv7 (ahí SÍ es opcional, a diferencia de ARM64) y siempre usa el camino escalar. La detección real (leer `HWCAP_NEON` vía `getauxval`) queda pendiente para quien tenga acceso a hardware/toolchain ARMv7 real |

Si vas a usar numerx en ARM64/ARMv7 en producción (incluyendo
Android), **por favor compílalo y corre la suite de tests (`test_c.c`)
en tu propio entorno antes de confiar en el camino NEON** — no se dio
por sentado que "se ve bien en el código" equivale a "funciona
correctamente en hardware real".

### Compilar para Android (NDK)

No verificado en este proyecto (sin NDK disponible en el entorno de
desarrollo), pero estos son los comandos estándar documentados por
Google para compilar una biblioteca C nativa con el NDK, adaptados a
`numerx.c`:

```bash
# arm64-v8a (dispositivos ARM64 modernos)
$NDK/toolchains/llvm/prebuilt/<host-tag>/bin/aarch64-linux-android21-clang \
    -O3 -shared -fPIC \
    -o libnumerx.so numerx.c

# armeabi-v7a (dispositivos ARM de 32 bits, más antiguos)
$NDK/toolchains/llvm/prebuilt/<host-tag>/bin/armv7a-linux-androideabi21-clang \
    -O3 -shared -fPIC \
    -o libnumerx.so numerx.c
```

Donde `$NDK` es la ruta a tu instalación del Android NDK, `<host-tag>`
es la plataforma de tu máquina de desarrollo (`linux-x86_64`,
`darwin-x86_64`, etc.), y `21` es el nivel mínimo de API de Android
objetivo (ajustar según tu `minSdkVersion`). No hace falta enlazar
`-lpthread` explícitamente en Android: pthread es parte de bionic (la
libc de Android) y ya está disponible sin flags adicionales.

## Desarrollo / compilar desde fuente

```bash
# extensión de Python
python3 setup.py build_ext --inplace

# tests de Python (cubre 0.0.1 a 0.1.0, ver test_python.py)
python3 test_python.py

# tests de C puro
gcc -Wall -Wextra -o test_c numerx.c test_c.c -lm -lpthread && ./test_c
```

## Licencia

MIT — ver el archivo `LICENSE`.
