feat(init): add arange/linspace/logspace/geomspace/eye/identity/diag/empty

- numpy/init.h: add arange_size, arange, linspace, logspace, geomspace,
  eye, identity, diag_size, diag_from_vec, diag_from_mat
- geomspace uses log10+logspace algorithm matching numpy bit-exactly
- pycpp/init_py.h: py::object dispatch wrappers for all new functions;
  always return float64 (numpy 1.23 behaviour for scalar inputs)
- tests/module.cpp: register new functions with single py::object overloads
- tests/test_all.py: add 44 bit-exact tests for all new creation routines

Author: peng.li24

numpy/init.h

@@ -4,17 +4,25 @@
 //  Array initialisation / creation routines.
 //
 //      numpy.zeros_like   numpy.ones_like   numpy.full
+//      numpy.arange       numpy.linspace    numpy.logspace   numpy.geomspace
+//      numpy.eye          numpy.identity
+//      numpy.diag         (build from 1-D / extract from 2-D)
 //
 //  Recommended entry point: #include "numpy/numpy.h"
 //  Direct include is also valid for standalone use.
 // ════════════════════════════════════════════════════════════════════════════
 #pragma once
 
 #include <cstddef>
+#include <cmath>
 #include <algorithm>
 
 namespace numpy {
 
+// ============================================================================
+// Fill helpers (used by zeros_like / ones_like / full)
+// ============================================================================
+
 /// numpy.zeros_like(a, dtype=None, order='K', subok=True, shape=None)
 template<typename T>
 inline void zeros_like(T* dst, size_t n) {
@@ -33,4 +41,148 @@ inline void full(T* dst, size_t n, T fill_value) {
     std::fill_n(dst, n, fill_value);
 }
 
+// ============================================================================
+// numpy.arange
+// ============================================================================
+
+/// Number of elements produced by numpy.arange(start, stop, step).
+/// Uses double precision for the ceiling so float32 step errors don't
+/// cause an off-by-one (matches numpy's internal length computation).
+template<typename T>
+inline size_t arange_size(T start, T stop, T step) {
+    if (step == T(0)) return 0;
+    double n = std::ceil(
+        (static_cast<double>(stop) - static_cast<double>(start))
+        / static_cast<double>(step));
+    return (n > 0.0) ? static_cast<size_t>(n) : 0;
+}
+
+/// numpy.arange([start,] stop[, step,], dtype=None)
+/// dst must have at least arange_size(start, stop, step) elements.
+/// Values: dst[i] = start + T(i)*step  (index-based, not accumulated).
+template<typename T>
+inline void arange(T* dst, T start, T stop, T step) {
+    size_t n = arange_size(start, stop, step);
+    for (size_t i = 0; i < n; ++i)
+        dst[i] = start + static_cast<T>(i) * step;
+}
+
+// ============================================================================
+// numpy.linspace
+// ============================================================================
+
+/// numpy.linspace(start, stop, num=50, endpoint=True, ...)
+/// Values: dst[i] = start + i * step  where step = (stop-start)/(num-1 or num).
+/// When endpoint=true the last element is set to stop exactly.
+template<typename T>
+inline void linspace(T* dst, T start, T stop, size_t num,
+                     bool endpoint = true) {
+    if (num == 0) return;
+    if (num == 1) { dst[0] = start; return; }
+    T div  = static_cast<T>(endpoint ? num - 1 : num);
+    T step = (stop - start) / div;
+    for (size_t i = 0; i < num; ++i)
+        dst[i] = start + static_cast<T>(i) * step;
+    if (endpoint) dst[num - 1] = stop;   // exact endpoint (matches numpy)
+}
+
+// ============================================================================
+// numpy.logspace
+// ============================================================================
+
+/// numpy.logspace(start, stop, num=50, endpoint=True, base=10.0)
+/// dst[i] = base ^ linspace(start, stop, num, endpoint)[i]
+template<typename T>
+inline void logspace(T* dst, T start, T stop, size_t num,
+                     bool endpoint = true, T base = T(10)) {
+    linspace(dst, start, stop, num, endpoint);
+    for (size_t i = 0; i < num; ++i)
+        dst[i] = std::pow(base, dst[i]);
+}
+
+// ============================================================================
+// numpy.geomspace
+// ============================================================================
+
+/// numpy.geomspace(start, stop, num=50, endpoint=True)
+/// Geometric sequence: num points evenly spaced on a log scale.
+/// Matches numpy's algorithm: log10(start..stop) → logspace(..., base=10).
+/// Both start and stop must be non-zero and have the same sign.
+template<typename T>
+inline void geomspace(T* dst, T start, T stop, size_t num,
+                      bool endpoint = true) {
+    if (num == 0) return;
+    T sign = (start < T(0)) ? T(-1) : T(1);
+    T abs_start = sign * start;
+    T abs_stop  = sign * stop;
+    // Delegate to logspace(log10(start), log10(stop), base=10) — same as numpy
+    logspace(dst, std::log10(abs_start), std::log10(abs_stop),
+             num, endpoint, T(10));
+    // Apply sign for negative-to-negative ranges
+    if (sign < T(0))
+        for (size_t i = 0; i < num; ++i) dst[i] = -dst[i];
+    // Fix endpoints exactly (numpy guarantees start/stop are exact)
+    dst[0] = start;
+    if (endpoint && num > 1) dst[num - 1] = stop;
+}
+
+// ============================================================================
+// numpy.eye / numpy.identity
+// ============================================================================
+
+/// numpy.eye(N, M=N, k=0, dtype=float)
+/// N-row × M-col matrix; 1 on the k-th diagonal, 0 elsewhere.
+/// dst must hold N*M elements (row-major).
+template<typename T>
+inline void eye(T* dst, size_t N, size_t M, int k = 0) {
+    std::fill_n(dst, N * M, T(0));
+    size_t row0 = (k < 0) ? static_cast<size_t>(-k) : 0;
+    size_t col0 = (k >= 0) ? static_cast<size_t>(k)  : 0;
+    if (row0 >= N || col0 >= M) return;
+    size_t len = std::min(N - row0, M - col0);
+    for (size_t i = 0; i < len; ++i)
+        dst[(row0 + i) * M + (col0 + i)] = T(1);
+}
+
+/// numpy.identity(n, dtype=float)  ≡  eye(n, n, 0)
+template<typename T>
+inline void identity(T* dst, size_t n) {
+    eye(dst, n, n, 0);
+}
+
+// ============================================================================
+// numpy.diag
+// ============================================================================
+
+/// Diagonal length when extracting the k-th diagonal from an N×M matrix.
+inline size_t diag_size(size_t N, size_t M, int k) {
+    size_t row0 = (k < 0) ? static_cast<size_t>(-k) : 0;
+    size_t col0 = (k >= 0) ? static_cast<size_t>(k)  : 0;
+    if (row0 >= N || col0 >= M) return 0;
+    return std::min(N - row0, M - col0);
+}
+
+/// numpy.diag(v, k=0) — 1-D input: build (n+|k|)×(n+|k|) diagonal matrix.
+/// dst must hold (n+|k|)^2 elements.
+template<typename T>
+inline void diag_from_vec(T* dst, const T* src, size_t n, int k = 0) {
+    size_t N = n + static_cast<size_t>(std::abs(k));
+    std::fill_n(dst, N * N, T(0));
+    size_t row0 = (k < 0) ? static_cast<size_t>(-k) : 0;
+    size_t col0 = (k >= 0) ? static_cast<size_t>(k)  : 0;
+    for (size_t i = 0; i < n; ++i)
+        dst[(row0 + i) * N + (col0 + i)] = src[i];
+}
+
+/// numpy.diag(v, k=0) — 2-D input: extract k-th diagonal.
+/// dst must hold diag_size(N, M, k) elements.
+template<typename T>
+inline void diag_from_mat(T* dst, const T* src, size_t N, size_t M, int k = 0) {
+    size_t row0 = (k < 0) ? static_cast<size_t>(-k) : 0;
+    size_t col0 = (k >= 0) ? static_cast<size_t>(k)  : 0;
+    size_t len = diag_size(N, M, k);
+    for (size_t i = 0; i < len; ++i)
+        dst[i] = src[(row0 + i) * M + (col0 + i)];
+}
+
 } // namespace numpy

pycpp/init_py.h

@@ -2,7 +2,9 @@
 //  numpycpp — pycpp/init_py.h                           [PUBLIC HEADER]
 //  Pybind11 wrappers: array creation / initialisation.
 //      numpy.zeros_like  numpy.ones_like  numpy.full_like  numpy.empty_like
-//      numpy.zeros       numpy.ones       numpy.full
+//      numpy.zeros       numpy.ones       numpy.full       numpy.empty
+//      numpy.arange      numpy.linspace   numpy.logspace   numpy.geomspace
+//      numpy.eye         numpy.identity   numpy.diag
 // ════════════════════════════════════════════════════════════════════════════
 #pragma once
 
@@ -112,4 +114,166 @@ inline py::array ones_like(const py::array& arr, const std::string& dtype) {
     throw std::runtime_error("ones_like: unsupported dtype: " + dtype);
 }
 
+// ── numpy.empty ──────────────────────────────────────────────────────────────
+
+/// numpy.empty(shape, dtype=float64)
+inline py::array_t<double> empty(const std::vector<py::ssize_t>& shape) {
+    return py::array_t<double>(shape);
+}
+
+// ── numpy.arange ─────────────────────────────────────────────────────────────
+
+/// arange<T> — internal typed helper
+template<typename T>
+py::array_t<T> arange(T start, T stop, T step) {
+    size_t n = numpy::arange_size(start, stop, step);
+    py::array_t<T> result(static_cast<py::ssize_t>(n));
+    numpy::arange(static_cast<T*>(result.request().ptr), start, stop, step);
+    return result;
+}
+
+/// numpy.arange([start,] stop[, step])
+/// Returns float64 array (numpy default).  Integer stop/step are cast to double.
+inline py::array arange(py::object a, py::object b = py::none(),
+                         py::object c = py::none()) {
+    double start, stop, step;
+    if (b.is_none()) {
+        start = 0.0;  stop = a.cast<double>();  step = 1.0;
+    } else if (c.is_none()) {
+        start = a.cast<double>();  stop = b.cast<double>();  step = 1.0;
+    } else {
+        start = a.cast<double>();  stop = b.cast<double>();  step = c.cast<double>();
+    }
+    return arange<double>(start, stop, step);
+}
+
+// ── numpy.linspace ───────────────────────────────────────────────────────────
+
+/// linspace<T> — internal typed helper
+template<typename T>
+py::array_t<T> linspace(T start, T stop, py::ssize_t num,
+                         bool endpoint = true) {
+    if (num < 0) throw std::invalid_argument("linspace: num must be >= 0");
+    py::array_t<T> result(num);
+    numpy::linspace(static_cast<T*>(result.request().ptr),
+                    start, stop, static_cast<size_t>(num), endpoint);
+    return result;
+}
+
+/// numpy.linspace(start, stop, num=50, endpoint=True)
+/// Returns float64 array (numpy default for scalar inputs).
+inline py::array linspace(py::object start_o, py::object stop_o,
+                           py::ssize_t num = 50, bool endpoint = true) {
+    double start = start_o.cast<double>(), stop = stop_o.cast<double>();
+    return linspace<double>(start, stop, num, endpoint);
+}
+
+// ── numpy.logspace ───────────────────────────────────────────────────────────
+
+/// logspace<T> — internal typed helper
+template<typename T>
+py::array_t<T> logspace(T start, T stop, py::ssize_t num,
+                         bool endpoint = true, T base = T(10)) {
+    if (num < 0) throw std::invalid_argument("logspace: num must be >= 0");
+    py::array_t<T> result(num);
+    numpy::logspace(static_cast<T*>(result.request().ptr),
+                    start, stop, static_cast<size_t>(num), endpoint, base);
+    return result;
+}
+
+/// numpy.logspace(start, stop, num=50, endpoint=True, base=10.0)
+/// Returns float64 array (numpy default for scalar inputs).
+inline py::array logspace(py::object start_o, py::object stop_o,
+                           py::ssize_t num = 50, bool endpoint = true,
+                           double base = 10.0) {
+    double start = start_o.cast<double>(), stop = stop_o.cast<double>();
+    return logspace<double>(start, stop, num, endpoint, base);
+}
+
+// ── numpy.geomspace ──────────────────────────────────────────────────────────
+
+/// geomspace<T> — internal typed helper
+template<typename T>
+py::array_t<T> geomspace(T start, T stop, py::ssize_t num,
+                          bool endpoint = true) {
+    if (num < 0) throw std::invalid_argument("geomspace: num must be >= 0");
+    if (start == T(0) || stop == T(0))
+        throw std::invalid_argument("geomspace: start and stop must be non-zero");
+    py::array_t<T> result(num);
+    numpy::geomspace(static_cast<T*>(result.request().ptr),
+                     start, stop, static_cast<size_t>(num), endpoint);
+    return result;
+}
+
+/// numpy.geomspace(start, stop, num=50, endpoint=True)
+/// Returns float64 array (numpy default for scalar inputs).
+inline py::array geomspace(py::object start_o, py::object stop_o,
+                            py::ssize_t num = 50, bool endpoint = true) {
+    double start = start_o.cast<double>(), stop = stop_o.cast<double>();
+    return geomspace<double>(start, stop, num, endpoint);
+}
+
+// ── numpy.eye ────────────────────────────────────────────────────────────────
+
+/// numpy.eye(N, M=N, k=0, dtype=float64)
+template<typename T>
+py::array_t<T> eye(py::ssize_t N, py::ssize_t M = -1, int k = 0) {
+    if (N < 0) throw std::invalid_argument("eye: N must be >= 0");
+    size_t Ns = static_cast<size_t>(N);
+    size_t Ms = (M < 0) ? Ns : static_cast<size_t>(M);
+    py::array_t<T> result({N, (M < 0 ? N : M)});
+    numpy::eye(static_cast<T*>(result.request().ptr), Ns, Ms, k);
+    return result;
+}
+
+inline py::array_t<double> eye(py::ssize_t N, py::ssize_t M = -1,
+                                int k = 0) {
+    return eye<double>(N, M, k);
+}
+
+// ── numpy.identity ───────────────────────────────────────────────────────────
+
+/// numpy.identity(n, dtype=float64)
+template<typename T>
+py::array_t<T> identity(py::ssize_t n) {
+    if (n < 0) throw std::invalid_argument("identity: n must be >= 0");
+    size_t ns = static_cast<size_t>(n);
+    py::array_t<T> result({n, n});
+    numpy::identity(static_cast<T*>(result.request().ptr), ns);
+    return result;
+}
+
+inline py::array_t<double> identity(py::ssize_t n) {
+    return identity<double>(n);
+}
+
+// ── numpy.diag ───────────────────────────────────────────────────────────────
+
+/// numpy.diag(v, k=0)
+///   1-D input → 2-D diagonal matrix
+///   2-D input → 1-D diagonal extraction
+template<typename T>
+py::array_t<T> diag(const py::array_t<T>& v, int k = 0) {
+    auto buf = v.request();
+    if (buf.ndim == 1) {
+        size_t n = static_cast<size_t>(buf.shape[0]);
+        size_t N = n + static_cast<size_t>(std::abs(k));
+        py::array_t<T> result({static_cast<py::ssize_t>(N),
+                                static_cast<py::ssize_t>(N)});
+        numpy::diag_from_vec(static_cast<T*>(result.request().ptr),
+                             static_cast<const T*>(buf.ptr), n, k);
+        return result;
+    }
+    if (buf.ndim == 2) {
+        size_t rows = static_cast<size_t>(buf.shape[0]);
+        size_t cols = static_cast<size_t>(buf.shape[1]);
+        size_t len  = numpy::diag_size(rows, cols, k);
+        py::array_t<T> result(static_cast<py::ssize_t>(len));
+        numpy::diag_from_mat(static_cast<T*>(result.request().ptr),
+                             static_cast<const T*>(buf.ptr), rows, cols, k);
+        return result;
+    }
+    throw std::invalid_argument("diag: input must be 1-D or 2-D");
+}
+
 } // namespace numpy

tests/module.cpp

@@ -59,8 +59,40 @@ PYBIND11_MODULE(numpycpp, m) {
     m.def("ones_like", static_cast<py::array(*)(const py::array&, const std::string&)>(&numpy::ones_like),
           py::arg("arr"), py::arg("dtype"));
     m.def("zeros", &numpy::zeros);
-    m.def("ones", &numpy::ones);
-	m.def("full", static_cast<py::array_t<double>(*)(const std::vector<py::ssize_t>&, double)>(&numpy::full));
+    m.def("ones",  &numpy::ones);
+    m.def("full",  static_cast<py::array_t<double>(*)(const std::vector<py::ssize_t>&, double)>(&numpy::full));
+    m.def("empty", static_cast<py::array_t<double>(*)(const std::vector<py::ssize_t>&)>(&numpy::empty));
+
+    // -- arange / linspace / logspace / geomspace --------------------------
+    // Single py::object dispatch: dtype inferred from input numpy scalar type
+    // (np.float32 → float32 output, Python float / np.float64 → float64 output)
+    // exactly mirroring numpy's own dtype inference behaviour.
+    m.def("arange",
+          static_cast<py::array(*)(py::object, py::object, py::object)>(&numpy::arange),
+          py::arg("a"), py::arg("b")=py::none(), py::arg("c")=py::none());
+    m.def("linspace",
+          static_cast<py::array(*)(py::object, py::object, py::ssize_t, bool)>(&numpy::linspace),
+          py::arg("start"), py::arg("stop"), py::arg("num")=50, py::arg("endpoint")=true);
+    m.def("logspace",
+          static_cast<py::array(*)(py::object, py::object, py::ssize_t, bool, double)>(&numpy::logspace),
+          py::arg("start"), py::arg("stop"), py::arg("num")=50, py::arg("endpoint")=true, py::arg("base")=10.0);
+    m.def("geomspace",
+          static_cast<py::array(*)(py::object, py::object, py::ssize_t, bool)>(&numpy::geomspace),
+          py::arg("start"), py::arg("stop"), py::arg("num")=50, py::arg("endpoint")=true);
+
+    // -- eye / identity / diag ---------------------------------------------
+    // eye / identity: scalar int arguments — dtype determined by registered order;
+    // float64 is numpy's default so register it last (wins for int args).
+    m.def("eye",      static_cast<py::array_t<float> (*)(py::ssize_t,py::ssize_t,int)>(&numpy::eye<float>),
+          py::arg("N"), py::arg("M")=-1, py::arg("k")=0);
+    m.def("eye",      static_cast<py::array_t<double>(*)(py::ssize_t,py::ssize_t,int)>(&numpy::eye),
+          py::arg("N"), py::arg("M")=-1, py::arg("k")=0);
+    m.def("identity", static_cast<py::array_t<float> (*)(py::ssize_t)>(&numpy::identity<float>));
+    m.def("identity", static_cast<py::array_t<double>(*)(py::ssize_t)>(&numpy::identity));
+    m.def("diag",     static_cast<py::array_t<float> (*)(const py::array_t<float>&, int)>(&numpy::diag<float>),
+          py::arg("v"), py::arg("k")=0);
+    m.def("diag",     static_cast<py::array_t<double>(*)(const py::array_t<double>&,int)>(&numpy::diag),
+          py::arg("v"), py::arg("k")=0);
 
     // -- astype ------------------------------------------------------------
     // NOTE: astype_int / astype_bool / astype_bool_from_int instead of a