const FLOAT64_BYTES = Float64Array.BYTES_PER_ELEMENT;

export const DEFAULT_HERMITIAN_JACOBI_MAX_SWEEPS = 32;
export const DEFAULT_MAX_HERMITIAN_EIGEN_WORK_BYTES = 16 * 1024 * 1024;
export const DEFAULT_MAX_HERMITIAN_JACOBI_PAIR_UPDATES = 1_000_000_000;

function assertPositiveInteger(value, label) {
  if (!Number.isSafeInteger(value) || value < 1) {
    throw new TypeError(`${label} must be a positive safe integer.`);
  }
}

function assertMaximumSweeps(value) {
  if (
    !Number.isSafeInteger(value) ||
    value < 0 ||
    value > DEFAULT_HERMITIAN_JACOBI_MAX_SWEEPS
  ) {
    throw new RangeError(
      `maximumSweeps must be an integer from 0 through ${DEFAULT_HERMITIAN_JACOBI_MAX_SWEEPS}.`,
    );
  }
}

function checkedProduct(factors, label) {
  let product = 1;
  for (const factor of factors) {
    if (factor !== 0 && product > Number.MAX_SAFE_INTEGER / factor) {
      throw new RangeError(`${label} exceeds the safe integer range.`);
    }
    product *= factor;
  }
  return product;
}

function roundoffGamma(operationCount) {
  const product = operationCount * (Number.EPSILON / 2);
  if (!(product < 1)) {
    throw new RangeError("roundoff bound is undefined for this workload.");
  }
  return product / (1 - product);
}

/**
 * Estimate the bounded typed-array allocation and cyclic-Jacobi work before
 * any decomposition workspace is allocated.
 */
export function estimateHermitianEigenWorkload({
  order,
  maximumSweeps = DEFAULT_HERMITIAN_JACOBI_MAX_SWEEPS,
}) {
  assertPositiveInteger(order, "order");
  assertMaximumSweeps(maximumSweeps);

  const matrixElementCount = checkedProduct([order, order], "matrix element count");
  const pairCount = checkedProduct([order, order - 1], "Jacobi pair count") / 2;
  const worstCaseRotations = checkedProduct(
    [pairCount, maximumSweeps],
    "Jacobi rotation count",
  );
  // A rotation updates n-2 matrix row pairs and n eigenvector row pairs.
  const pairUpdatesPerRotation = Math.max(0, 2 * order - 2);
  const worstCasePairUpdates = checkedProduct(
    [worstCaseRotations, pairUpdatesPerRotation],
    "Jacobi pair-update count",
  );
  // Two split-complex matrices plus normalized and physical eigenvalue arrays.
  const typedArrayElementCount = checkedProduct(
    [4 * matrixElementCount + 2 * order, FLOAT64_BYTES],
    "Hermitian eigensolver allocation",
  ) / FLOAT64_BYTES;
  const allocationBytes = checkedProduct(
    [typedArrayElementCount, FLOAT64_BYTES],
    "Hermitian eigensolver allocation",
  );

  return Object.freeze({
    order,
    maximumSweeps,
    matrixElementCount,
    pairCount,
    pairUpdatesPerRotation,
    worstCaseRotations,
    worstCasePairUpdates,
    allocationBytes,
  });
}

function inferOrder(real, imaginary, explicitOrder) {
  if (real == null || imaginary == null) {
    throw new TypeError("real and imaginary matrix components are required.");
  }
  const order = explicitOrder ?? Math.sqrt(real.length);
  assertPositiveInteger(order, "order");
  const expectedLength = checkedProduct([order, order], "matrix element count");
  if (real.length !== expectedLength || imaginary.length !== expectedLength) {
    throw new RangeError(
      `Hermitian matrix components must each contain ${expectedLength} values.`,
    );
  }
  return order;
}

function matrixFrobeniusNorm(real, imaginary, order) {
  let sum = 0;
  for (let index = 0; index < order * order; index += 1) {
    sum += real[index] * real[index] + imaginary[index] * imaginary[index];
  }
  return Math.sqrt(sum);
}

function offDiagonalFrobeniusNorm(real, imaginary, order) {
  let upperTriangleSum = 0;
  for (let row = 0; row < order; row += 1) {
    for (let column = row + 1; column < order; column += 1) {
      const index = row * order + column;
      upperTriangleSum +=
        real[index] * real[index] + imaginary[index] * imaginary[index];
    }
  }
  return Math.sqrt(2 * upperTriangleSum);
}

function rotatePair(matrixReal, matrixImaginary, vectorsReal, vectorsImaginary, order, p, q) {
  const pq = p * order + q;
  const offReal = matrixReal[pq];
  const offImaginary = matrixImaginary[pq];
  const magnitude = Math.hypot(offReal, offImaginary);
  if (magnitude === 0) return false;

  const diagonalP = matrixReal[p * order + p];
  const diagonalQ = matrixReal[q * order + q];
  const halfDifference = (diagonalQ - diagonalP) / 2;
  const tangent = halfDifference === 0
    ? 1
    : Math.sign(halfDifference) * magnitude /
      (Math.abs(halfDifference) + Math.hypot(halfDifference, magnitude));
  const cosine = 1 / Math.hypot(1, tangent);
  const sine = tangent * cosine;
  const phaseReal = offReal / magnitude;
  const phaseImaginary = offImaginary / magnitude;

  for (let row = 0; row < order; row += 1) {
    if (row === p || row === q) continue;
    const rowP = row * order + p;
    const rowQ = row * order + q;
    const pReal = matrixReal[rowP];
    const pImaginary = matrixImaginary[rowP];
    const qReal = matrixReal[rowQ];
    const qImaginary = matrixImaginary[rowQ];

    const conjugatePhaseTimesQReal =
      phaseReal * qReal + phaseImaginary * qImaginary;
    const conjugatePhaseTimesQImaginary =
      phaseReal * qImaginary - phaseImaginary * qReal;
    const phaseTimesPReal = phaseReal * pReal - phaseImaginary * pImaginary;
    const phaseTimesPImaginary = phaseReal * pImaginary + phaseImaginary * pReal;

    const nextPReal = cosine * pReal - sine * conjugatePhaseTimesQReal;
    const nextPImaginary = cosine * pImaginary - sine * conjugatePhaseTimesQImaginary;
    const nextQReal = sine * phaseTimesPReal + cosine * qReal;
    const nextQImaginary = sine * phaseTimesPImaginary + cosine * qImaginary;

    matrixReal[rowP] = nextPReal;
    matrixImaginary[rowP] = nextPImaginary;
    matrixReal[p * order + row] = nextPReal;
    matrixImaginary[p * order + row] = -nextPImaginary;
    matrixReal[rowQ] = nextQReal;
    matrixImaginary[rowQ] = nextQImaginary;
    matrixReal[q * order + row] = nextQReal;
    matrixImaginary[q * order + row] = -nextQImaginary;
  }

  matrixReal[p * order + p] = diagonalP - tangent * magnitude;
  matrixImaginary[p * order + p] = 0;
  matrixReal[q * order + q] = diagonalQ + tangent * magnitude;
  matrixImaginary[q * order + q] = 0;
  matrixReal[pq] = 0;
  matrixImaginary[pq] = 0;
  matrixReal[q * order + p] = 0;
  matrixImaginary[q * order + p] = 0;

  for (let row = 0; row < order; row += 1) {
    const rowP = row * order + p;
    const rowQ = row * order + q;
    const pReal = vectorsReal[rowP];
    const pImaginary = vectorsImaginary[rowP];
    const qReal = vectorsReal[rowQ];
    const qImaginary = vectorsImaginary[rowQ];
    const conjugatePhaseTimesQReal =
      phaseReal * qReal + phaseImaginary * qImaginary;
    const conjugatePhaseTimesQImaginary =
      phaseReal * qImaginary - phaseImaginary * qReal;
    const phaseTimesPReal = phaseReal * pReal - phaseImaginary * pImaginary;
    const phaseTimesPImaginary = phaseReal * pImaginary + phaseImaginary * pReal;
    vectorsReal[rowP] = cosine * pReal - sine * conjugatePhaseTimesQReal;
    vectorsImaginary[rowP] = cosine * pImaginary - sine * conjugatePhaseTimesQImaginary;
    vectorsReal[rowQ] = sine * phaseTimesPReal + cosine * qReal;
    vectorsImaginary[rowQ] = sine * phaseTimesPImaginary + cosine * qImaginary;
  }
  return true;
}

function canonicalizeColumnPhases(real, imaginary, order) {
  for (let column = 0; column < order; column += 1) {
    let pivotRow = 0;
    let pivotMagnitude = -1;
    for (let row = 0; row < order; row += 1) {
      const index = row * order + column;
      const magnitude = Math.hypot(real[index], imaginary[index]);
      if (magnitude > pivotMagnitude) {
        pivotMagnitude = magnitude;
        pivotRow = row;
      }
    }
    const pivotIndex = pivotRow * order + column;
    const pivotReal = real[pivotIndex];
    const pivotImaginary = imaginary[pivotIndex];
    const normalizedPivotMagnitude = Math.hypot(pivotReal, pivotImaginary);
    if (!(normalizedPivotMagnitude > 0) || !Number.isFinite(normalizedPivotMagnitude)) {
      throw new RangeError("Jacobi eigenvector phase canonicalization failed.");
    }
    const phaseReal = pivotReal / normalizedPivotMagnitude;
    const phaseImaginary = -pivotImaginary / normalizedPivotMagnitude;
    for (let row = 0; row < order; row += 1) {
      const index = row * order + column;
      const valueReal = real[index];
      const valueImaginary = imaginary[index];
      real[index] = valueReal * phaseReal - valueImaginary * phaseImaginary;
      imaginary[index] = valueReal * phaseImaginary + valueImaginary * phaseReal;
    }
    real[pivotIndex] = normalizedPivotMagnitude;
    imaginary[pivotIndex] = 0;
  }
}

function projectedInputEntry(real, imaginary, order, row, column, scale) {
  if (row === column) return [real[row * order + row] / scale, 0];
  const forward = row * order + column;
  const reverse = column * order + row;
  return [
    (real[forward] + real[reverse]) / (2 * scale),
    (imaginary[forward] - imaginary[reverse]) / (2 * scale),
  ];
}

function residualFrobeniusNorm(
  inputReal,
  inputImaginary,
  vectorsReal,
  vectorsImaginary,
  normalizedEigenvalues,
  order,
  scale,
) {
  if (scale === 0) return 0;
  let residualSquared = 0;
  for (let column = 0; column < order; column += 1) {
    const eigenvalue = normalizedEigenvalues[column];
    for (let row = 0; row < order; row += 1) {
      let productReal = 0;
      let productImaginary = 0;
      for (let inner = 0; inner < order; inner += 1) {
        const [matrixReal, matrixImaginary] = projectedInputEntry(
          inputReal,
          inputImaginary,
          order,
          row,
          inner,
          scale,
        );
        const vectorIndex = inner * order + column;
        const vectorReal = vectorsReal[vectorIndex];
        const vectorImaginary = vectorsImaginary[vectorIndex];
        productReal += matrixReal * vectorReal - matrixImaginary * vectorImaginary;
        productImaginary += matrixReal * vectorImaginary + matrixImaginary * vectorReal;
      }
      const vectorIndex = row * order + column;
      const differenceReal = productReal - eigenvalue * vectorsReal[vectorIndex];
      const differenceImaginary =
        productImaginary - eigenvalue * vectorsImaginary[vectorIndex];
      residualSquared +=
        differenceReal * differenceReal + differenceImaginary * differenceImaginary;
    }
  }
  return Math.sqrt(residualSquared);
}

function orthogonalityError(real, imaginary, order) {
  let errorSquared = 0;
  for (let left = 0; left < order; left += 1) {
    for (let right = 0; right < order; right += 1) {
      let productReal = 0;
      let productImaginary = 0;
      for (let row = 0; row < order; row += 1) {
        const leftIndex = row * order + left;
        const rightIndex = row * order + right;
        productReal +=
          real[leftIndex] * real[rightIndex] +
          imaginary[leftIndex] * imaginary[rightIndex];
        productImaginary +=
          real[leftIndex] * imaginary[rightIndex] -
          imaginary[leftIndex] * real[rightIndex];
      }
      if (left === right) productReal -= 1;
      errorSquared += productReal * productReal + productImaginary * productImaginary;
    }
  }
  return Math.sqrt(errorSquared);
}

/**
 * Deterministic full cyclic Jacobi decomposition of a complex Hermitian
 * matrix. Split-complex matrices are row-major; eigenvectors are returned as
 * row-major columns, so component `row` of eigenvector `column` is at
 * `row * order + column`.
 *
 * The input is projected to its nearest Hermitian matrix only after its
 * anti-Hermitian discrepancy has passed the reported roundoff-scale audit.
 * No diagonal loading or eigenvalue flooring is performed.
 */
export function decomposeHermitian({
  real,
  imaginary,
  order: explicitOrder,
  maximumSweeps = DEFAULT_HERMITIAN_JACOBI_MAX_SWEEPS,
  relativeTolerance,
  maximumWorkBytes = DEFAULT_MAX_HERMITIAN_EIGEN_WORK_BYTES,
  maximumPairUpdates = DEFAULT_MAX_HERMITIAN_JACOBI_PAIR_UPDATES,
}) {
  const order = inferOrder(real, imaginary, explicitOrder);
  const workload = estimateHermitianEigenWorkload({ order, maximumSweeps });
  if (!(maximumWorkBytes >= 0) || !Number.isFinite(maximumWorkBytes)) {
    throw new TypeError("maximumWorkBytes must be finite and nonnegative.");
  }
  if (!(maximumPairUpdates >= 0) || !Number.isFinite(maximumPairUpdates)) {
    throw new TypeError("maximumPairUpdates must be finite and nonnegative.");
  }
  if (workload.allocationBytes > maximumWorkBytes) {
    throw new RangeError(
      `Hermitian eigensolver requires ${workload.allocationBytes} typed-array bytes; limit is ${maximumWorkBytes}.`,
    );
  }
  if (workload.worstCasePairUpdates > maximumPairUpdates) {
    throw new RangeError(
      `Hermitian eigensolver worst-case work is ${workload.worstCasePairUpdates} pair updates; limit is ${maximumPairUpdates}.`,
    );
  }

  let scale = 0;
  for (let index = 0; index < workload.matrixElementCount; index += 1) {
    const valueReal = real[index];
    const valueImaginary = imaginary[index];
    if (!Number.isFinite(valueReal) || !Number.isFinite(valueImaginary)) {
      throw new TypeError("Hermitian matrix components must be finite.");
    }
    // Component scaling avoids overflow in hypot() for otherwise finite
    // split-complex entries near Number.MAX_VALUE.
    scale = Math.max(scale, Math.abs(valueReal), Math.abs(valueImaginary));
  }

  const inputHermitianTolerance =
    scale === 0 ? 0 : 32 * Number.EPSILON * Math.max(1, order) * scale;
  let inputHermitianError = 0;
  for (let row = 0; row < order; row += 1) {
    inputHermitianError = Math.max(
      inputHermitianError,
      Math.abs(imaginary[row * order + row]),
    );
    for (let column = row + 1; column < order; column += 1) {
      const forward = row * order + column;
      const reverse = column * order + row;
      inputHermitianError = Math.max(
        inputHermitianError,
        Math.hypot(
          real[forward] - real[reverse],
          imaginary[forward] + imaginary[reverse],
        ),
      );
    }
  }
  if (inputHermitianError > inputHermitianTolerance) {
    throw new RangeError(
      `matrix is not Hermitian within roundoff: error ${inputHermitianError}, tolerance ${inputHermitianTolerance}.`,
    );
  }

  const matrixReal = new Float64Array(workload.matrixElementCount);
  const matrixImaginary = new Float64Array(workload.matrixElementCount);
  const vectorsReal = new Float64Array(workload.matrixElementCount);
  const vectorsImaginary = new Float64Array(workload.matrixElementCount);
  const normalizedEigenvalues = new Float64Array(order);
  const eigenvalues = new Float64Array(order);

  if (scale > 0) {
    for (let row = 0; row < order; row += 1) {
      const diagonal = row * order + row;
      matrixReal[diagonal] = real[diagonal] / scale;
      for (let column = row + 1; column < order; column += 1) {
        const forward = row * order + column;
        const reverse = column * order + row;
        const projectedReal = (real[forward] + real[reverse]) / (2 * scale);
        const projectedImaginary =
          (imaginary[forward] - imaginary[reverse]) / (2 * scale);
        matrixReal[forward] = projectedReal;
        matrixImaginary[forward] = projectedImaginary;
        matrixReal[reverse] = projectedReal;
        matrixImaginary[reverse] = -projectedImaginary;
      }
    }
  }
  for (let index = 0; index < order; index += 1) {
    vectorsReal[index * order + index] = 1;
  }

  const matrixNormNormalized = matrixFrobeniusNorm(
    matrixReal,
    matrixImaginary,
    order,
  );
  const effectiveRelativeTolerance = relativeTolerance ?? roundoffGamma(32 * order);
  if (!(effectiveRelativeTolerance > 0) || !Number.isFinite(effectiveRelativeTolerance)) {
    throw new TypeError("relativeTolerance must be finite and positive.");
  }
  const absoluteTolerance = effectiveRelativeTolerance * matrixNormNormalized;
  let offDiagonalNormalized = offDiagonalFrobeniusNorm(
    matrixReal,
    matrixImaginary,
    order,
  );
  let converged = offDiagonalNormalized <= absoluteTolerance;
  let sweeps = 0;
  let rotations = 0;
  const pairTolerance = workload.pairCount === 0
    ? 0
    : absoluteTolerance / Math.sqrt(2 * workload.pairCount);

  while (!converged && sweeps < maximumSweeps) {
    let sweepRotations = 0;
    for (let p = 0; p < order; p += 1) {
      for (let q = p + 1; q < order; q += 1) {
        const pq = p * order + q;
        if (Math.hypot(matrixReal[pq], matrixImaginary[pq]) <= pairTolerance) {
          continue;
        }
        if (
          rotatePair(
            matrixReal,
            matrixImaginary,
            vectorsReal,
            vectorsImaginary,
            order,
            p,
            q,
          )
        ) {
          sweepRotations += 1;
        }
      }
    }
    rotations += sweepRotations;
    sweeps += 1;
    offDiagonalNormalized = offDiagonalFrobeniusNorm(
      matrixReal,
      matrixImaginary,
      order,
    );
    converged = offDiagonalNormalized <= absoluteTolerance;
    if (sweepRotations === 0) break;
  }

  const stableOrder = Array.from({ length: order }, (_, index) => ({
    index,
    value: matrixReal[index * order + index],
  }));
  stableOrder.sort((left, right) => {
    if (left.value > right.value) return -1;
    if (left.value < right.value) return 1;
    return left.index - right.index;
  });
  for (let column = 0; column < order; column += 1) {
    const sourceColumn = stableOrder[column].index;
    const normalizedEigenvalue = stableOrder[column].value;
    normalizedEigenvalues[column] = Object.is(normalizedEigenvalue, -0)
      ? 0
      : normalizedEigenvalue;
    const physicalEigenvalue = normalizedEigenvalues[column] * scale;
    if (!Number.isFinite(physicalEigenvalue)) {
      throw new RangeError("Hermitian eigenvalue exceeds the finite Float64 range.");
    }
    eigenvalues[column] = Object.is(physicalEigenvalue, -0) ? 0 : physicalEigenvalue;
    for (let row = 0; row < order; row += 1) {
      const destination = row * order + column;
      const source = row * order + sourceColumn;
      matrixReal[destination] = vectorsReal[source];
      matrixImaginary[destination] = vectorsImaginary[source];
    }
  }
  // A product of unitary Jacobi rotations already has unit-norm columns.
  // Canonicalize phase only; do not hide orthogonality drift by renormalizing.
  canonicalizeColumnPhases(matrixReal, matrixImaginary, order);

  const residualNormalized = residualFrobeniusNorm(
    real,
    imaginary,
    matrixReal,
    matrixImaginary,
    normalizedEigenvalues,
    order,
    scale,
  );
  const relativeResidual = matrixNormNormalized === 0
    ? 0
    : residualNormalized / matrixNormNormalized;
  const relativeOffDiagonal = matrixNormNormalized === 0
    ? 0
    : offDiagonalNormalized / matrixNormNormalized;
  const residualTolerance = roundoffGamma(128 * order);
  const orthogonalityTolerance = residualTolerance * Math.sqrt(order);
  const measuredOrthogonalityError = orthogonalityError(
    matrixReal,
    matrixImaginary,
    order,
  );
  const numericallyVerified =
    converged &&
    relativeResidual <= residualTolerance &&
    measuredOrthogonalityError <= orthogonalityTolerance;

  return {
    order,
    eigenvalues,
    eigenvectors: {
      real: matrixReal,
      imaginary: matrixImaginary,
      layout: "row-major-columns",
    },
    converged,
    sweeps,
    rotations,
    scale,
    inputHermitianError,
    inputHermitianTolerance,
    relativeTolerance: effectiveRelativeTolerance,
    offDiagonalFrobeniusNorm: offDiagonalNormalized * scale,
    relativeOffDiagonal,
    residualFrobeniusNorm: residualNormalized * scale,
    relativeResidual,
    residualTolerance,
    orthogonalityError: measuredOrthogonalityError,
    orthogonalityTolerance,
    numericallyVerified,
    workload,
  };
}