import * as mats from './matrices';
import { intMatrices, matrices, residueClassRing } from './types';

const iops = intMatrices;
const fops = matrices;

const p = 9999991;
const pops = mats.extend(residueClassRing(p), ['Integer'], true);


const invModP = M => {
  const n = M.length;
  const A = M.map((row, i) => row.concat(pops.unitVector(n, i)));
  const E = pops.rowEchelonForm(A);

  return E.map(row => row.slice(n));
};


const numberOfPAdicStepsNeeded = (A, b) => {
  const lengths = M => fops.transposed(M).map(r => fops.norm(r));
  const max = v => v.reduce((x, y) => x > y ? x : y);
  const product = v => v.reduce((x, y) => x * y);

  const ls = lengths(A).concat(max(lengths(b)));
  const delta = product(ls.sort().slice(A[0].length));
  const golden = (1 + Math.sqrt(5)) / 2;

  return Math.ceil(2 * Math.log(delta * golden) / Math.log(p));
};


const rationalReconstruction = (s, h) => {
  let u = [h, s];
  let v = [0, 1];
  let sign = 1;

  while (iops.gt(iops.times(u[1], u[1]), h)) {
    const q = iops.idiv(u[0], u[1]);

    u = [u[1], iops.minus(u[0], iops.times(q, u[1]))];
    v = [v[1], iops.plus(v[0], iops.times(q, v[1]))];
    sign *= -1;
  }

  return fops.div(iops.times(sign, u[1]), v[1]);
};


export default function solve(A, b) {
  const C = invModP(A);
  const nrSteps = numberOfPAdicStepsNeeded(A, b);

  let bi = b;
  let pi = 1;
  let si = 0;

  for (let i = 0; i < nrSteps; ++i) {
    const xi = pops.times(C, bi);
    bi = iops.idiv(iops.minus(bi, iops.times(A, xi)), p);
    si = iops.plus(si, iops.times(pi, xi));
    pi = iops.times(pi, p);
  }

  return si.map(row => row.map(x => rationalReconstruction(x, pi)));
};


if (require.main == module) {
  Array.prototype.toString = function() {
    return '[ ' + this.map(x => x.toString()).join(', ') + ' ]';
  };

  const testSolver = (A, b) => {
    console.log();
    console.log(`solving ${A} * x = ${b}`);
    const x = solve(A, b);
    console.log(`  x = ${x}`);
    console.log(`  A * x = ${fops.times(A, x)}`);
  };

  testSolver(
    [ [  4, -4 ],
      [  1,  0 ] ],
    [ [  1,  1,  1 ],
      [  0,  0,  0 ] ]
  );
}

@@ -382,16 +382,88 @@
 rm(A);%0A%0A
+  if (fops.eq(E.map(row =%3E row.slice(0, n)), fops.identityMatrix(n)))%0A  
   return
@@ -1354,16 +1354,51 @@
 odP(A);%0A
+  if (C == null)%0A    return null;%0A%0A
   const 

do not return a result if matrix is singular