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 ops = pops;
  const n = M.length;
  const A = M.map((row, i) => row.concat(ops.unitVector(n, i)));
  const E = ops.rowEchelonForm(A);

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


const numberOfPAdicStepsNeeded = (A, b) => {
  const ops = fops;

  const lengths = M => ops.transposed(M).map(r => ops.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 ops.ceil(2 * Math.log(delta * golden) / Math.log(p));
};


export default function solve(A, b) {
  console.log(`solve(${A}, ${b})`);
  const C = invModP(A);
  console.log(`  C = ${C}`);
  console.log(`  A * C = ${pops.times(A, C)}`);

  const nrSteps = numberOfPAdicStepsNeeded(A, b);
  console.log(`  nrSteps = ${nrSteps}`);

  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);
  }
  const s = si;
  console.log(`  s <- ${s}`);
};


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

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

@@ -891,16 +891,338 @@
 );%0A%7D;%0A%0A%0A
+const rationalReconstruction = (s, h) =%3E %7B%0A  const limit = Math.sqrt(h);%0A  let u = %5Bh, s%5D;%0A  let v = %5B0, 1%5D;%0A  let sign = 1;%0A%0A  while (u%5B1%5D %3E= limit) %7B%0A    const q = Math.floor(u%5B0%5D / u%5B1%5D);%0A%0A    u = %5Bu%5B1%5D, u%5B0%5D - q * u%5B1%5D%5D;%0A    v = %5Bv%5B1%5D, v%5B0%5D + q * v%5B1%5D%5D;%0A    sign *= -1;%0A  %7D%0A%0A  return fops.div(sign * u%5B1%5D, v%5B1%5D);%0A%7D;%0A%0A%0A
 export d
@@ -1769,18 +1769,17 @@
   s 
-%3C-
+=
  $%7Bs%7D%60);
 %0A%7D;%0A
@@ -1774,16 +1774,130 @@
  $%7Bs%7D%60);
+%0A%0A  const r = s.map(row =%3E row.map(x =%3E rationalReconstruction(x, pi)));%0A  console.log(%60  r = $%7Br%7D%60);%0A%0A  return r;
 %0A%7D;%0A%0A%0Aif
@@ -2064,34 +2064,8 @@
 , -4
-,  1 %5D,%0A      %5B -4,  4,  0
  %5D,%0A
@@ -2075,20 +2075,16 @@
    %5B  1,
-  0,
   0 %5D %5D,
@@ -2110,30 +2110,8 @@
  %5D,%0A
-      %5B -1, -1, -1 %5D,%0A
     

implement rational reconstruction