using System;

namespace GlinerPoc.Neural;

/// <summary>
/// Phase 4 FP32 numerical kernels (narrow, correctness-first, scalar).
/// All kernels are OUT-OF-PLACE: inputs are never mutated (Phase 4.28) so
/// Phase 5 stage-by-stage comparisons can trust intermediate buffers.
/// Tolerance policy: abs(err) ≤ 1e-5 + 1e-4·abs(ref) unless a stage
/// documents a justified bound (P01 proposal).
/// </summary>
public static class GlinerMath
{
	public const float Atol = 1e-5f;
	public const float Rtol = 1e-4f;

	public sealed record Metrics(int Count, double MaxAbs, double MeanAbs, double MaxRel, int WorstIndex, int Violations)
	{
		/// <summary>Elementwise policy: abs(err) ≤ Atol + Rtol·abs(ref) for every element.</summary>
		public bool Pass => Count > 0 && Violations == 0;
	}

	// ---- byte → FP32 decoding (Phase 4.6) ------------------------------------

	/// <summary>Little-endian FP32 read from a byte payload (BitConverter — whitelist-safe).</summary>
	public static float ReadF32( byte[] data, int byteOffset ) =>
		BitConverter.ToSingle( data, byteOffset );

	/// <summary>Decode a full little-endian FP32 byte range to float[].</summary>
	public static float[] DecodeF32( byte[] data, int byteOffset, int floatCount )
	{
		var result = new float[floatCount];
		for ( int i = 0; i < floatCount; i++ )
		{
			result[i] = BitConverter.ToSingle( data, byteOffset + i * 4 );
		}
		return result;
	}

	// ---- Linear: Y = XWᵀ (+ b), X [S,in], W [out,in], b [out] -----------------

	public static float[] Linear( float[] x, int rows, float[] weight, float[] bias, int inDim, int outDim )
	{
		if ( x.Length != rows * inDim )
		{
			throw new ArgumentException( $"[GLI:ERROR] Linear input has {x.Length} values, expected {rows}x{inDim}." );
		}
		if ( weight.Length != inDim * outDim )
		{
			throw new ArgumentException( $"[GLI:ERROR] Linear weight has {weight.Length} values, expected {outDim}x{inDim}." );
		}
		if ( bias is not null && bias.Length != outDim )
		{
			throw new ArgumentException( $"[GLI:ERROR] Linear bias has {bias.Length} values, expected {outDim}." );
		}

		var y = new float[rows * outDim];
		for ( int r = 0; r < rows; r++ )
		{
			int xBase = r * inDim;
			int yBase = r * outDim;
			for ( int o = 0; o < outDim; o++ )
			{
				int wBase = o * inDim;
				float sum = bias is not null ? bias[o] : 0f;
				for ( int i = 0; i < inDim; i++ )
				{
					sum += x[xBase + i] * weight[wBase + i];
				}
				y[yBase + o] = sum;
			}
		}
		return y;
	}

	// ---- Residual: out = a + b (out-of-place) ---------------------------------

	public static float[] Add( float[] a, float[] b )
	{
		if ( a.Length != b.Length )
		{
			throw new ArgumentException( $"[GLI:ERROR] Residual shape mismatch {a.Length} vs {b.Length}." );
		}
		var y = new float[a.Length];
		for ( int i = 0; i < a.Length; i++ )
		{
			y[i] = a[i] + b[i];
		}
		return y;
	}

	// ---- LayerNorm (last dimension, biased variance, eps inside sqrt) ----------

	public static float[] LayerNorm( float[] x, int rows, int dim, float[] weight, float[] bias, float eps )
	{
		if ( x.Length != rows * dim )
		{
			throw new ArgumentException( $"[GLI:ERROR] LayerNorm input has {x.Length} values, expected {rows}x{dim}." );
		}
		var y = new float[x.Length];
		for ( int r = 0; r < rows; r++ )
		{
			int b = r * dim;
			// PyTorch CPU layer_norm accumulates in float32; double-precision
			// stats round differently and the 1/√(var+eps) amplification makes
			// that visible on near-constant rows. Match the oracle: float32
			// statistics and float32 normalization arithmetic.
			float mean = 0f;
			for ( int i = 0; i < dim; i++ )
			{
				mean += x[b + i];
			}
			mean /= dim;
			float variance = 0f;
			for ( int i = 0; i < dim; i++ )
			{
				float d = x[b + i] - mean;
				variance += d * d;
			}
			variance /= dim;
			float invStd = 1f / MathF.Sqrt( variance + eps );
			for ( int i = 0; i < dim; i++ )
			{
				y[b + i] = (x[b + i] - mean) * invStd * weight[i] + bias[i];
			}
		}
		return y;
	}

	// ---- GELU: exact erf form, 0.5·x·(1 + erf(x/√2)) ---------------------------

	public static float[] GeluErf( float[] x )
	{
		var y = new float[x.Length];
		for ( int i = 0; i < x.Length; i++ )
		{
			double v = x[i];
			y[i] = (float)(0.5 * v * (1.0 + Erf( v / Math.Sqrt( 2.0 ) )));
		}
		return y;
	}

	/// <summary>
	/// erf via Abramowitz &amp; Stegun 7.1.26 (max abs error 1.5e-7), double
	/// precision. The runtime BCL has no Math.Erf; this error is far inside the
	/// 1e-5 kernel tolerance for GELU (scaled by 0.5|x| it stays ≤ ~7.6e-7 at
	/// |x|=10).
	/// </summary>
	public static double Erf( double x )
	{
		const double p = 0.3275911;
		const double a1 = 0.254829592;
		const double a2 = -0.284496736;
		const double a3 = 1.421413741;
		const double a4 = -1.453152027;
		const double a5 = 1.061405429;
		double sign = x < 0 ? -1.0 : 1.0;
		double ax = Math.Abs( x );
		double t = 1.0 / (1.0 + p * ax);
		double y = 1.0 - ((((a5 * t + a4) * t + a3) * t + a2) * t + a1) * t * Math.Exp( -ax * ax );
		return sign * y;
	}

	// ---- ReLU (classifier activation, Phase 6.15) -------------------------------

	public static float[] Relu( float[] x )
	{
		var y = new float[x.Length];
		for ( int i = 0; i < x.Length; i++ )
		{
			y[i] = x[i] > 0f ? x[i] : 0f;
		}
		return y;
	}

	// ---- Stable row softmax ([rows, dim], dim = last) ---------------------------

	public static float[] SoftmaxRows( float[] x, int rows, int dim )
	{
		if ( x.Length != rows * dim )
		{
			throw new ArgumentException( $"[GLI:ERROR] Softmax input has {x.Length} values, expected {rows}x{dim}." );
		}
		var y = new float[x.Length];
		for ( int r = 0; r < rows; r++ )
		{
			int b = r * dim;
			float max = x[b];
			for ( int i = 1; i < dim; i++ )
			{
				if ( x[b + i] > max )
				{
					max = x[b + i];
				}
			}
			float sum = 0f;
			for ( int i = 0; i < dim; i++ )
			{
				float e = MathF.Exp( x[b + i] - max );
				y[b + i] = e;
				sum += e;
			}
			float inv = 1f / sum;
			for ( int i = 0; i < dim; i++ )
			{
				y[b + i] *= inv;
			}
		}
		return y;
	}

	// ---- comparison --------------------------------------------------------------

	public static Metrics Compare( float[] actual, float[] reference )
	{
		if ( actual.Length != reference.Length )
		{
			throw new ArgumentException( $"[GLI:ERROR] Compare length mismatch {actual.Length} vs {reference.Length}." );
		}
		double maxAbs = 0.0;
		double sumAbs = 0.0;
		double maxRel = 0.0;
		int worst = 0;
		int violations = 0;
		for ( int i = 0; i < reference.Length; i++ )
		{
			double d = Math.Abs( actual[i] - reference[i] );
			sumAbs += d;
			if ( d > maxAbs )
			{
				maxAbs = d;
				worst = i;
			}
			double r = Math.Abs( reference[i] );
			if ( d > Atol + Rtol * r )
			{
				violations++;
			}
			if ( r > 1e-12 )
			{
				double rel = d / r;
				if ( rel > maxRel )
				{
					maxRel = rel;
				}
			}
		}
		return new Metrics( actual.Length, maxAbs, sumAbs / actual.Length, maxRel, worst, violations );
	}
}