using System.Buffers; using System.Collections; using System.Collections.Concurrent; using System.Numerics; using System.Text; using System.Text.Json; namespace HyperMath.Core { public partial class Tensor<T> : IDisposable, ICloneable, IEquatable<Tensor<T>>, IEnumerable<T> where T : unmanaged, INumber<T>, IExponentialFunctions<T>, ITrigonometricFunctions<T> { #region Fields and Properties private readonly T[] data; private readonly int[] shape; private readonly int[] strides; private readonly int totalSize; private bool disposed = false; public int[] Shape => (int[])shape.Clone(); public int TotalSize => totalSize; public static int ParallelizationThreshold { get; set; } = 1000; public int Rank => shape.Length; #endregion #region Constructors public Tensor(params int[] dimensions) { if (dimensions == null || dimensions.Length == 0) throw new ArgumentException("Tensor must have at least one dimension.", nameof(dimensions)); foreach (var dim in dimensions) { if (dim <= 0) throw new ArgumentException("Each dimension size must be positive.", nameof(dimensions)); } shape = dimensions.ToArray(); strides = CalculateStrides(shape); totalSize = strides[0] * shape[0]; data = ArrayPool<T>.Shared.Rent(totalSize); Array.Clear(data, 0, totalSize); } private Tensor(T[] data, int[] shape, int[] strides) { this.data = data; this.shape = shape; this.strides = strides; this.totalSize = strides[0] * shape[0]; } #endregion #region Indexer and Element Access public T this[params int[] indices] { get { EnsureNotDisposed(); int flatIndex = CalculateFlatIndex(indices); return data[flatIndex]; } set { EnsureNotDisposed(); int flatIndex = CalculateFlatIndex(indices); data[flatIndex] = value; } } private int CalculateFlatIndex(int[] indices) { if (indices == null) throw new ArgumentNullException(nameof(indices)); if (indices.Length != shape.Length) throw new ArgumentException("Number of indices must match number of tensor dimensions.", nameof(indices)); int flatIndex = 0; for (int i = 0; i < indices.Length; i++) { int index = indices[i]; int dim = shape[i]; if (index < 0 || index >= dim) throw new ArgumentOutOfRangeException(nameof(indices), $"Index {index} out of bounds for dimension {i} of size {dim}."); flatIndex += indices[i] * strides[i]; } return flatIndex; } #endregion #region Cloning and Equality public object Clone() { EnsureNotDisposed(); T[] newData = ArrayPool<T>.Shared.Rent(totalSize); Array.Copy(data, 0, newData, 0, totalSize); int[] newShape = (int[])shape.Clone(); int[] newStrides = (int[])strides.Clone(); return new Tensor<T>(newData, newShape, newStrides); } public bool Equals(Tensor<T>? other) { EnsureNotDisposed(); if (other == null) return false; if (this == other) return true; if (!shape.SequenceEqual(other.shape)) return false; for (int i = 0; i < totalSize; i++) { if (!data[i].Equals(other.data[i])) return false; } return true; } #endregion #region IEnumerable Implementation public IEnumerator<T> GetEnumerator() { EnsureNotDisposed(); for (int i = 0; i < totalSize; i++) { yield return data[i]; } } IEnumerator IEnumerable.GetEnumerator() { return GetEnumerator(); } #endregion #region Arithmetic Operations public static Tensor<T> operator +(Tensor<T> a, Tensor<T> b) { return a.Add(b); } public static Tensor<T> operator -(Tensor<T> a, Tensor<T> b) { return a.Subtract(b); } public static Tensor<T> operator *(Tensor<T> a, Tensor<T> b) { return a.Multiply(b); } public static Tensor<T> operator /(Tensor<T> a, Tensor<T> b) { return a.Divide(b); } public static Tensor<T> operator +(Tensor<T> a, T scalar) { return a.AddScalar(scalar); } public static Tensor<T> operator +(T scalar, Tensor<T> a) { return a.AddScalar(scalar); } public static Tensor<T> operator -(Tensor<T> a, T scalar) { return a.SubtractScalar(scalar); } public static Tensor<T> operator -(T scalar, Tensor<T> a) { return a.Negate().AddScalar(scalar); } public static Tensor<T> operator *(Tensor<T> a, T scalar) { return a.MultiplyScalar(scalar); } public static Tensor<T> operator *(T scalar, Tensor<T> a) { return a.MultiplyScalar(scalar); } public static Tensor<T> operator /(Tensor<T> a, T scalar) { return a.DivideScalar(scalar); } public static Tensor<T> operator /(T scalar, Tensor<T> a) { return a.Reciprocal().MultiplyScalar(scalar); } public Tensor<T> Add(Tensor<T> other) { return ElementWiseOperation(other, (x, y) => x + y); } public Tensor<T> Subtract(Tensor<T> other) { return ElementWiseOperation(other, (x, y) => x - y); } public Tensor<T> Multiply(Tensor<T> other) { return ElementWiseOperation(other, (x, y) => x * y); } public Tensor<T> Divide(Tensor<T> other) { return ElementWiseOperation(other, (x, y) => x / y); } public Tensor<T> AddScalar(T scalar) { return ScalarOperation((x) => x + scalar); } public Tensor<T> SubtractScalar(T scalar) { return ScalarOperation((x) => x - scalar); } public Tensor<T> MultiplyScalar(T scalar) { return ScalarOperation((x) => x * scalar); } public Tensor<T> DivideScalar(T scalar) { return ScalarOperation((x) => x / scalar); } public Tensor<T> Negate() { return ScalarOperation((x) => -x); } public Tensor<T> Reciprocal() { return ScalarOperation((x) => T.One / x); } public async Task<Tensor<T>> AddAsync(Tensor<T> other) { return await Task.Run(() => Add(other)); } public async Task<Tensor<T>> SubtractAsync(Tensor<T> other) { return await Task.Run(() => Subtract(other)); } public async Task<Tensor<T>> MultiplyAsync(Tensor<T> other) { return await Task.Run(() => Multiply(other)); } public async Task<Tensor<T>> DivideAsync(Tensor<T> other) { return await Task.Run(() => Divide(other)); } #endregion #region Reduction Operations public T Sum() { EnsureNotDisposed(); if (totalSize >= ParallelizationThreshold) { T sum = T.Zero; object lockObj = new object(); Parallel.For(0, totalSize, () => T.Zero, (i, state, localSum) => { localSum += data[i]; return localSum; }, (localSum) => { lock (lockObj) { sum += localSum; } }); return sum; } else { T sum = T.Zero; for (int i = 0; i < totalSize; i++) { sum += data[i]; } return sum; } } public T Mean() { EnsureNotDisposed(); T sum = Sum(); return sum / T.CreateChecked(totalSize); } public T Max() { EnsureNotDisposed(); if (totalSize == 0) throw new InvalidOperationException("Cannot compute Max of an empty tensor."); if (totalSize >= ParallelizationThreshold) { T max = data[0]; object lockObj = new object(); Parallel.For(0, totalSize, () => data[0], (i, state, localMax) => { if (data[i] > localMax) localMax = data[i]; return localMax; }, (localMax) => { lock (lockObj) { if (localMax > max) max = localMax; } }); return max; } else { T max = data[0]; for (int i = 1; i < totalSize; i++) { if (data[i] > max) max = data[i]; } return max; } } public T Min() { EnsureNotDisposed(); if (totalSize == 0) throw new InvalidOperationException("Cannot compute Min of an empty tensor."); if (totalSize >= ParallelizationThreshold) { T min = data[0]; object lockObj = new object(); Parallel.For(0, totalSize, () => data[0], (i, state, localMin) => { if (data[i] < localMin) localMin = data[i]; return localMin; }, (localMin) => { lock (lockObj) { if (localMin < min) min = localMin; } }); return min; } else { T min = data[0]; for (int i = 1; i < totalSize; i++) { if (data[i] < min) min = data[i]; } return min; } } public Tensor<T> Sum(int axis) { return Reduce(axis, (a, b) => a + b, T.Zero); } public Tensor<T> Mean(int axis) { var sumTensor = Sum(axis); int reduceDim = shape[axis]; T divisor = T.CreateChecked(reduceDim); return sumTensor.DivideScalar(divisor); } public Tensor<T> Max(int axis) { return Reduce(axis, (a, b) => a > b ? a : b); } public Tensor<T> Min(int axis) { return Reduce(axis, (a, b) => a < b ? a : b); } #endregion #region Advanced Operations public Tensor<T> MatMul(Tensor<T> other) { EnsureNotDisposed(); other.EnsureNotDisposed(); if (this.Rank != 2 || other.Rank != 2) throw new ArgumentException("MatMul is only supported for 2D tensors."); if (this.shape[1] != other.shape[0]) throw new ArgumentException($"Inner tensor dimensions must match for MatMul. Found {this.shape[1]} and {other.shape[0]}."); int m = this.shape[0]; int n = this.shape[1]; int p = other.shape[1]; Tensor<T> result = new Tensor<T>(m, p); if (totalSize >= ParallelizationThreshold) { Parallel.For(0, m, i => { for (int j = 0; j < p; j++) { T sum = T.Zero; for (int k = 0; k < n; k++) { sum += this[i, k] * other[k, j]; } result[i, j] = sum; } }); } else { for (int i = 0; i < m; i++) { for (int j = 0; j < p; j++) { T sum = T.Zero; for (int k = 0; k < n; k++) { sum += this[i, k] * other[k, j]; } result[i, j] = sum; } } } return result; } public Tensor<T> Transpose(int axis1, int axis2) { EnsureNotDisposed(); if (axis1 < 0 || axis1 >= shape.Length) throw new ArgumentOutOfRangeException(nameof(axis1), "Axis1 is out of bounds."); if (axis2 < 0 || axis2 >= shape.Length) throw new ArgumentOutOfRangeException(nameof(axis2), "Axis2 is out of bounds."); int[] newShape = (int[])shape.Clone(); int[] newStrides = (int[])strides.Clone(); (newShape[axis1], newShape[axis2]) = (newShape[axis2], newShape[axis1]); (newStrides[axis1], newStrides[axis2]) = (newStrides[axis2], newStrides[axis1]); T[] newData = ArrayPool<T>.Shared.Rent(totalSize); int[] indices = new int[shape.Length]; for (int i = 0; i < totalSize; i++) { int remaining = i; for (int j = 0; j < shape.Length; j++) { indices[j] = remaining / strides[j]; remaining %= strides[j]; } (indices[axis1], indices[axis2]) = (indices[axis2], indices[axis1]); int newFlatIndex = 0; for (int j = 0; j < newShape.Length; j++) { newFlatIndex += indices[j] * newStrides[j]; } newData[newFlatIndex] = data[i]; } return new Tensor<T>(newData, newShape, newStrides); } #endregion #region Caching Broadcasting private class BroadcastCacheEntry { public int[] ResultShape { get; } public int[] BroadcastStridesA { get; } public int[] BroadcastStridesB { get; } public BroadcastCacheEntry(int[] resultShape, int[] broadcastStridesA, int[] broadcastStridesB) { ResultShape = resultShape; BroadcastStridesA = broadcastStridesA; BroadcastStridesB = broadcastStridesB; } } private static readonly ConcurrentDictionary<(string, string), BroadcastCacheEntry> BroadcastCache = new ConcurrentDictionary<(string, string), BroadcastCacheEntry>(); private BroadcastCacheEntry GetOrAddBroadcastCache(int[] shapeA, int[] shapeB, int[] otherStrides) { string keyA = string.Join(",", shapeA); string keyB = string.Join(",", shapeB); var cacheKey = (keyA, keyB); return BroadcastCache.GetOrAdd(cacheKey, _ => { int[] resultShape = BroadcastShapes(shapeA, shapeB); int[] broadcastStridesA = GetBroadcastStrides(shapeA, resultShape, strides); int[] broadcastStridesB = GetBroadcastStrides(shapeB, resultShape, otherStrides); return new BroadcastCacheEntry(resultShape, broadcastStridesA, broadcastStridesB); }); } #endregion #region Helper Methods private Tensor<T> ElementWiseOperation(Tensor<T> other, Func<T, T, T> operation) { EnsureNotDisposed(); other.EnsureNotDisposed(); var broadcastInfo = GetOrAddBroadcastCache(this.shape, other.shape, other.strides); int[] resultShape = broadcastInfo.ResultShape; int resultSize = resultShape.Aggregate(1, (a, b) => a * b); Tensor<T> result = new Tensor<T>(resultShape); if (resultSize >= ParallelizationThreshold) { Parallel.For(0, resultSize, i => { int[] indices = UnflattenIndex(i, resultShape); int indexA = GetFlatIndex(indices, broadcastInfo.BroadcastStridesA); int indexB = GetFlatIndex(indices, broadcastInfo.BroadcastStridesB); result.data[i] = operation(this.data[indexA], other.data[indexB]); }); } else { for (int i = 0; i < resultSize; i++) { int[] indices = UnflattenIndex(i, resultShape); int indexA = GetFlatIndex(indices, broadcastInfo.BroadcastStridesA); int indexB = GetFlatIndex(indices, broadcastInfo.BroadcastStridesB); result.data[i] = operation(this.data[indexA], other.data[indexB]); } } return result; } private Tensor<T> ScalarOperation(Func<T, T> operation) { EnsureNotDisposed(); T[] resultData = ArrayPool<T>.Shared.Rent(totalSize); if (totalSize >= ParallelizationThreshold) { Parallel.For(0, totalSize, i => { resultData[i] = operation(data[i]); }); } else { for (int i = 0; i < totalSize; i++) { resultData[i] = operation(data[i]); } } return new Tensor<T>(resultData, (int[])shape.Clone(), (int[])strides.Clone()); } private Tensor<T> Reduce(int axis, Func<T, T, T> reduction, T initial) { EnsureNotDisposed(); if (axis < 0 || axis >= shape.Length) throw new ArgumentOutOfRangeException(nameof(axis)); int[] newShape = shape.Where((s, i) => i != axis).ToArray(); int[] newStrides = CalculateStrides(newShape); int reduceDim = shape[axis]; int newTotalSize = newShape.Aggregate(1, (a, b) => a * b); T[] resultData = ArrayPool<T>.Shared.Rent(newTotalSize); for (int i = 0; i < newTotalSize; i++) { resultData[i] = initial; } int[][] indicesMap = new int[totalSize][]; if (totalSize >= ParallelizationThreshold) { Parallel.For(0, totalSize, i => { indicesMap[i] = UnflattenIndex(i, shape); }); } else { for (int i = 0; i < totalSize; i++) { indicesMap[i] = UnflattenIndex(i, shape); } } if (totalSize >= ParallelizationThreshold) { object lockObj = new object(); Parallel.For(0, totalSize, i => { var indices = indicesMap[i]; int resultIndex = 0; int stride = 1; for (int j = newShape.Length - 1, k = shape.Length - 1; j >= 0; j--, k--) { if (k == axis) k--; resultIndex += indices[k] * stride; stride *= newShape[j]; } T value = data[i]; lock (lockObj) { resultData[resultIndex] = reduction(resultData[resultIndex], value); } }); } else { for (int i = 0; i < totalSize; i++) { var indices = indicesMap[i]; int resultIndex = 0; int stride = 1; for (int j = newShape.Length - 1, k = shape.Length - 1; j >= 0; j--, k--) { if (k == axis) k--; resultIndex += indices[k] * stride; stride *= newShape[j]; } T value = data[i]; resultData[resultIndex] = reduction(resultData[resultIndex], value); } } return new Tensor<T>(resultData, newShape, newStrides); } private Tensor<T> Reduce(int axis, Func<T, T, T> reduction) { EnsureNotDisposed(); if (axis < 0 || axis >= shape.Length) throw new ArgumentOutOfRangeException(nameof(axis)); int[] newShape = shape.Where((s, i) => i != axis).ToArray(); int[] newStrides = CalculateStrides(newShape); int newTotalSize = newShape.Aggregate(1, (a, b) => a * b); T[] resultData = ArrayPool<T>.Shared.Rent(newTotalSize); bool[] initialized = new bool[newTotalSize]; int[][] indicesMap = new int[totalSize][]; for (int i = 0; i < totalSize; i++) { indicesMap[i] = UnflattenIndex(i, shape); } for (int i = 0; i < totalSize; i++) { var indices = indicesMap[i]; int resultIndex = 0; int stride = 1; for (int j = newShape.Length - 1, k = shape.Length - 1; j >= 0; j--, k--) { if (k == axis) k--; resultIndex += indices[k] * stride; stride *= newShape[j]; } T value = data[i]; if (!initialized[resultIndex]) { resultData[resultIndex] = value; initialized[resultIndex] = true; } else { resultData[resultIndex] = reduction(resultData[resultIndex], value); } } return new Tensor<T>(resultData, newShape, newStrides); } private static int[] CalculateStrides(int[] shape) { int[] strides = new int[shape.Length]; int stride = 1; for (int i = shape.Length - 1; i >= 0; i--) { strides[i] = stride; stride *= shape[i]; } return strides; } private static int[] BroadcastShapes(int[] shapeA, int[] shapeB) { int maxRank = Math.Max(shapeA.Length, shapeB.Length); int[] resultShape = new int[maxRank]; for (int i = 0; i < maxRank; i++) { int dimA = i < (maxRank - shapeA.Length) ? 1 : shapeA[i - (maxRank - shapeA.Length)]; int dimB = i < (maxRank - shapeB.Length) ? 1 : shapeB[i - (maxRank - shapeB.Length)]; if (dimA != dimB && dimA != 1 && dimB != 1) throw new InvalidOperationException("Shapes are not broadcastable."); resultShape[i] = Math.Max(dimA, dimB); } return resultShape; } private static int[] GetBroadcastStrides(int[] originalShape, int[] resultShape, int[] originalStrides) { int originalRank = originalShape.Length; int resultRank = resultShape.Length; int[] broadcastStrides = new int[resultRank]; int offset = resultRank - originalRank; for (int i = 0; i < resultRank; i++) { if (i < offset) { broadcastStrides[i] = 0; } else { if (originalShape[i - offset] == resultShape[i]) broadcastStrides[i] = originalStrides[i - offset]; else if (originalShape[i - offset] == 1) broadcastStrides[i] = 0; else throw new InvalidOperationException("Shapes are not broadcastable."); } } return broadcastStrides; } private int[] UnflattenIndex(int flatIndex, int[] shape) { int[] indices = new int[shape.Length]; for (int i = 0; i < shape.Length; i++) { indices[i] = flatIndex / strides[i]; flatIndex %= strides[i]; } return indices; } private static int GetFlatIndex(int[] indices, int[] strides) { int flatIndex = 0; for (int i = 0; i < indices.Length; i++) { flatIndex += indices[i] * strides[i]; } return flatIndex; } private void EnsureNotDisposed() { if (disposed) throw new ObjectDisposedException(nameof(Tensor<T>)); } #endregion #region Serialization public string Display() { EnsureNotDisposed(); StringBuilder display = new StringBuilder(); if (shape.Length == 1) { display.Append("["); for (int i = 0; i < shape[0]; i++) { display.Append(this[i]); if (i < shape[0] - 1) display.Append(", "); } display.AppendLine("]"); } else if (shape.Length == 2) { for (int i = 0; i < shape[0]; i++) { display.Append("["); for (int j = 0; j < shape[1]; j++) { display.Append(this[i, j]); if (j < shape[1] - 1) display.Append(", "); } display.AppendLine("]"); } } else { display.AppendLine("Tensor data display is only implemented for 1D and 2D tensors."); } return display.ToString(); } public void ToJson(Stream stream) { EnsureNotDisposed(); using var writer = new Utf8JsonWriter(stream, new JsonWriterOptions { Indented = false }); writer.WriteStartObject(); writer.WritePropertyName("shape"); writer.WriteStartArray(); foreach (var dim in shape) { writer.WriteNumberValue(dim); } writer.WriteEndArray(); writer.WritePropertyName("data"); writer.WriteStartArray(); for (int i = 0; i < totalSize; i++) { if (typeof(T) == typeof(float)) writer.WriteNumberValue(Convert.ToSingle(data[i])); else if (typeof(T) == typeof(double)) writer.WriteNumberValue(Convert.ToDouble(data[i])); else if (typeof(T) == typeof(int)) writer.WriteNumberValue(Convert.ToInt32(data[i])); else if (typeof(T) == typeof(long)) writer.WriteNumberValue(Convert.ToInt64(data[i])); else throw new InvalidOperationException($"Unsupported type {typeof(T)}"); } writer.WriteEndArray(); writer.WriteEndObject(); } public static Tensor<T> FromJson(Stream stream) { using var jsonDoc = JsonDocument.Parse(stream); var root = jsonDoc.RootElement; var shapeElement = root.GetProperty("shape"); int[] shape = shapeElement.EnumerateArray().Select(x => x.GetInt32()).ToArray(); var tensor = new Tensor<T>(shape); var dataElement = root.GetProperty("data"); int index = 0; var data = tensor.data; foreach (var item in dataElement.EnumerateArray()) { if (typeof(T) == typeof(float)) data[index++] = (T)(object)item.GetSingle(); else if (typeof(T) == typeof(double)) data[index++] = (T)(object)item.GetDouble(); else if (typeof(T) == typeof(int)) data[index++] = (T)(object)item.GetInt32(); else if (typeof(T) == typeof(long)) data[index++] = (T)(object)item.GetInt64(); else throw new InvalidOperationException($"Unsupported type {typeof(T)}"); } return tensor; } public void ToBinary(Stream stream) { EnsureNotDisposed(); using var writer = new BinaryWriter(stream, System.Text.Encoding.UTF8, leaveOpen: true); writer.Write(shape.Length); foreach (var dim in shape) { writer.Write(dim); } for (int i = 0; i < totalSize; i++) { if (typeof(T) == typeof(float)) writer.Write(Convert.ToSingle(data[i])); else if (typeof(T) == typeof(double)) writer.Write(Convert.ToDouble(data[i])); else if (typeof(T) == typeof(int)) writer.Write(Convert.ToInt32(data[i])); else if (typeof(T) == typeof(long)) writer.Write(Convert.ToInt64(data[i])); else throw new InvalidOperationException($"Unsupported type {typeof(T)}"); } } public static Tensor<T> FromBinary(Stream stream) { using var reader = new BinaryReader(stream, System.Text.Encoding.UTF8, leaveOpen: true); int dimensions = reader.ReadInt32(); int[] shape = new int[dimensions]; for (int i = 0; i < dimensions; i++) { shape[i] = reader.ReadInt32(); } var tensor = new Tensor<T>(shape); var data = tensor.data; for (int i = 0; i < tensor.totalSize; i++) { if (typeof(T) == typeof(float)) data[i] = (T)(object)reader.ReadSingle(); else if (typeof(T) == typeof(double)) data[i] = (T)(object)reader.ReadDouble(); else if (typeof(T) == typeof(int)) data[i] = (T)(object)reader.ReadInt32(); else if (typeof(T) == typeof(long)) data[i] = (T)(object)reader.ReadInt64(); else throw new InvalidOperationException($"Unsupported type {typeof(T)}"); } return tensor; } #endregion #region Dispose Pattern public void Dispose() { if (!disposed) { ArrayPool<T>.Shared.Return(data, clearArray: true); disposed = true; } } #endregion } } 

based on best practices and academic knowledge in software engineering and numerical computing implement the following improvements

Consider caching the broadcasting strides and shape calculations if these operations are repeatedly applied to the same tensor shapes.