import { assert } from '../../../../../../common/util/util.js';
import { FP, FPInterval } from '../../../../../util/floating_point.js';
import { PRNG } from '../../../../../util/prng.js';
import { makeCaseCache } from '../../case_cache.js';

// We want each case to avoid losing accuracy. We want to prevent bits from
// falling off the bottom of the significand during intermediate
// calculations. We want this to be true for reasonable implementations of
// the determinant. Let's compute a budget B for the number of significant
// bits allowed for each element in the matrix.
//
// Example computation method 1:  Use lower-diagonal form:
//  A matrix is in lower-diagonal form if all the entries above and to
//  the right of the diagonal are zero.
//  1. Convert the matrix to lower-diagonal form using elementary row
//  operations.
//  2. Then the determinant is the product of the elements on the diagonal.
//
//  Analysis:
//   - Let M be the number of significand bits for the floating point type.
//     Remember to include the hidden/implicit 1, because we will assume
//     all numbers are normal.
//   - It takes N-1 elementary row operations to convert a general NxN
//     matrix to lower diagonal form.  Let's ignore the fact this uses
//     division. Let's allocate a budget 1 bit of precision per row operation.
//   - If the diagonal elements are integers bounded by D, then their product
//     is bounded by D**N.  So we split the remaining significand bit budget
//     bits N ways.
//
//  So we have B ~ (M - (N-1)) / N
//
//  For f32, M = 24:
//     N = 4 :  B =~ ( 24 - 3 ) / 4 = 21 / 4
//              The bound on matrix elements should be 2**(21/4) or roughly 38.
//     N = 3 :  B =~ ( 24 - 2 ) / 3 = 22 / 3
//              The bound on matrix elements should be 2**(22/3) or roughly 161.
//     N = 2 :  B =~ ( 24 - 1 ) / 2 = 23 / 2
//              The bound on matrix elements should be 2**(23/2) or roughly 2896.
//
//  For f16, M = 11:
//     N = 4 :  B =~ ( 11 - 3 ) / 4 = 8 / 4 = 2
//              The bound on matrix elements should be 2**2 = 4.
//     N = 3 :  B =~ ( 11 - 2 ) / 3 = 9 / 3 = 3
//              The bound on matrix elements should be 2**3 = 8.
//     N = 2 :  B =~ ( 11 - 1 ) / 2 = 10 / 2 = 5
//              The bound on matrix elements should be 2**5 = 32.
//
// Example computation method 2: Naive method:
//   The naive calculation of a determinant of an NxN matrix is the alternating
//   sum of N terms, where each term is an entry from the matrix multiplied
//   by the determinant of an (N-1)x(N-1) submatrix.
//   The recurrence is:   T(N) = (N-1) + N + T(N-1), with T(1) = 0.
//   The first N is for the adds and subtracts, the second N is for the multiplies,
//   and the T(N-1) is the determinant of the submatrix.
//   As N grows, this blows up quickly, even when accounting for reusing results.
//   It's unlikely any GPU would actually calculate it this way.
//
// Accuracy for determinant is only defined for e, where e is an integer and
// |e| < quadroot(2**21) [~38],
// due to computational complexity of calculating the general solution for 4x4,
// so custom matrices are used.

type Dim = 2 | 3 | 4;
type FPWidth = 16 | 32;

// Number of random matrices to test, per configuration.
const numSamples = 20;

// Returns a random element in a range suitable for a square matrix
// that can have an accurately computed determinant.
function randomMatrixEntry(p: PRNG, dim: Dim, fpwidth: FPWidth): number {
  // See above for the justification.
  const rangeTable = {
    16: {
      2: 32,
      3: 8,
      4: 4,
    },
    32: {
      2: 2896,
      3: 161,
      4: 38,
    },
  };
  const N = rangeTable[fpwidth][dim];
  // Centre the distribution roughly at zero.
  const balanced = p.uniformInt(N) - Math.floor(N / 2);
  return balanced;
}

// Returns a random square matrix that should have an exactly computed
// determinant for the given floating point width.
// At least some of the matrices returned by this function should have
// a non-zero determinant. This will be checked later by the nonTrivialMatrices
// function.
function randomSquareMatrix(p: PRNG, dim: Dim, fpwidth: FPWidth): number[][] {
  const result: number[][] = [...Array(dim)].map(_ => [...Array(dim)]);
  // Scale each element by a simple power of two. This should only affect
  // the exponent of the result.
  const multiplier = [1, 2, 0.25][p.uniformInt(3)];
  for (let c = 0; c < dim; c++) {
    for (let r = 0; r < dim; r++) {
      result[c][r] = multiplier * randomMatrixEntry(p, dim, fpwidth);
    }
  }
  return result;
}

type detFnType = (m: number[][]) => FPInterval;

// Returns true if at least one of the matrices has a non-zero determinant.
function nonTrivialMatrices(matrices: number[][][], detFn: detFnType): boolean {
  const detInterval = (m: number[][]) => detFn(m);
  const sumBegin = matrices.reduce((accum, m) => accum + Math.abs(detInterval(m).begin), 0);
  const sumEnd = matrices.reduce((accum, m) => accum + Math.abs(detInterval(m).end), 0);
  return sumBegin > 0 && sumEnd >= sumBegin;
}

// Cases: f32_matDxD_[non_]const
const f32_cases = ([2, 3, 4] as const)
  .flatMap(dim =>
    ([true, false] as const).map(nonConst => ({
      [`f32_mat${dim}x${dim}_${nonConst ? 'non_const' : 'const'}`]: () => {
        const p = new PRNG(dim + 32);
        const matrices: number[][][] = [...Array(numSamples)].map(_ =>
          randomSquareMatrix(p, dim, 32)
        );
        assert(nonTrivialMatrices(matrices, FP.f32.determinantInterval));
        return FP.f32.generateMatrixToScalarCases(
          matrices,
          nonConst ? 'unfiltered' : 'finite',
          FP.f32.determinantInterval
        );
      },
    }))
  )
  .reduce((a, b) => ({ ...a, ...b }), {});

// Cases: f16_matDxD_[non_]const
const f16_cases = ([2, 3, 4] as const)
  .flatMap(dim =>
    ([true, false] as const).map(nonConst => ({
      [`f16_mat${dim}x${dim}_${nonConst ? 'non_const' : 'const'}`]: () => {
        const p = new PRNG(dim + 16);
        const matrices: number[][][] = [...Array(numSamples)].map(_ =>
          randomSquareMatrix(p, dim, 16)
        );
        assert(nonTrivialMatrices(matrices, FP.f16.determinantInterval));
        return FP.f16.generateMatrixToScalarCases(
          matrices,
          nonConst ? 'unfiltered' : 'finite',
          FP.f16.determinantInterval
        );
      },
    }))
  )
  .reduce((a, b) => ({ ...a, ...b }), {});

// Cases: abstract_matDxD
const abstract_cases = ([2, 3, 4] as const)
  .map(dim => ({
    [`abstract_mat${dim}x${dim}`]: () => {
      const p = new PRNG(dim + 64);
      // Use f32 values range for abstract float.
      const matrices: number[][][] = [...Array(numSamples)].map(_ =>
        randomSquareMatrix(p, dim, 32)
      );
      assert(nonTrivialMatrices(matrices, FP.f32.determinantInterval));
      return FP.abstract.generateMatrixToScalarCases(
        matrices,
        'finite',
        // determinant has an inherited accuracy, so abstract is only expected to be as accurate as f32
        FP.f32.determinantInterval
      );
    },
  }))
  .reduce((a, b) => ({ ...a, ...b }), {});

export const d = makeCaseCache('determinant', {
  ...f32_cases,
  ...f16_cases,
  ...abstract_cases,
});
