Source code in the jq programming language

def sum(s): reduce s as $_ (0; . + $_);

# Sum of squares
def ss(s): sum(s|.*.);

# Create an m x n matrix
def matrix(m; n; init):
  if m == 0 then []
  elif m == 1 then [range(0;n) | init]
  elif m > 0 then
    matrix(1;n;init) as $row
    | [range(0;m) | $row ]
  else error("matrix\(m);_;_) invalid")
  end;

def dot_product(a; b):
  reduce range(0;a|length) as $i (0; . + (a[$i] * b[$i]) );

# A and B should both be numeric matrices, A being m by n, and B being n by p.
def multiply($A; $B):
  ($B[0]|length) as $p
  | ($B|transpose) as $BT
  | reduce range(0; $A|length) as $i
       ([];
       reduce range(0; $p) as $j 
         (.;
          .[$i][$j] = dot_product( $A[$i]; $BT[$j] ) ));

# $ndec decimal places
def round($ndec):
  def rpad: tostring | ($ndec - length) as $l | . + ("0" * $l);
  def abs: if . < 0 then -. else . end;
  pow(10; $ndec) as $p
  | round as $round
  | if $p * ((. - $round)|abs) < 0.1
    then ($round|tostring)  + "." + ($ndec * "0")
    else  . * $p | round / $p
    | tostring
    | capture("(?<left>[^.]*)[.](?<right>.*)")
    | .left + "." + (.right|rpad)
    end;

# pretty-print a 2-d matrix
def pp($ndec; $width):
  def pad(n): tostring | (n - length) * " " + .;
  def row: map(round($ndec) | pad($width)) | join(" ");
  reduce .[] as $row (""; . + "\n\($row|row)");

def minor($x; $d):
   ($x|length) as $nr
   | ($x[0]|length) as $nc
   | reduce range(0; $d) as $i (matrix($nr;$nc;0); .[$i][$i] = 1)
   | reduce range($d; $nr) as $i (.;
       reduce range($d;$nc) as $j (.; .[$i][$j] = $x[$i][$j] ) );

def vmadd($a; $b; $s):
  reduce range (0; $a|length) as $i ([];
    .[$i] = $a[$i] + $s * $b[$i] );

def vmul($v):
  ($v|length) as $n
  | reduce range(0;$n) as $i (null;
      reduce range(0;$n) as $j (.; .[$i][$j] = -2 * $v[$i] * $v[$j] ))
  | reduce range(0;$n) as $i (.; .[$i][$i] += 1 );

def vnorm($x):
  sum($x[] | .*.) | sqrt;

def vdiv($x; $d):
  [range (0;$x|length) | $x[.] / $d];

def mcol($m; $c):
  [range (0;$m|length) | $m[.][$c]];

def householder($m):
    ($m|length) as $nr
    | ($m[0]|length) as $nc
    | { q: [], # $nr
        z: $m,
        k: 0 }
    | until( .k >= $nc or .k >= $nr-1;
        .z = minor(.z; .k)
        | .x = mcol(.z; .k)
        | .a = vnorm(.x)
        | if ($m[.k][.k] > 0) then .a = -.a else . end
        | .e = [range (0; $nr) as $i | if ($i == .k) then 1 else 0 end]
        | .e = vmadd(.x; .e; .a)
        | .e = vdiv(.e; vnorm(.e))
        | .q[.k] = vmul(.e)
        | .z = multiply(.q[.k]; .z)
        | .k += 1 )
    | .Q = .q[0]
    | .R = multiply(.q[0]; $m)
    | .i = 1
    | until (.i >= $nc or .i >= $nr-1;
        .Q = multiply(.q[.i]; .Q)
        | .i += 1 )
    | .R = multiply(.Q; $m)
    | .Q |= transpose
    | [.Q, .R] ;

def x: [
  [12, -51,   4],
  [ 6, 167, -68],
  [-4,  24, -41],
  [-1,   1,   0],
  [ 2,   0,   3]
];

def task:
  def pp: pp(3;8);

  # Assume $a and $b are conformal
  def ssd($a; $b):
    [$a[][]] as $a
    | [$b[][]] as $b
    | ss( range(0;$a|length) | $a[.] - $b[.] );

  householder(x) as [$Q, $R]
  | multiply($Q; $R) as $m
  | "Q:",      ($Q|pp),
   "\nR:",     ($R|pp),
   "\nQ * R:", ($m|pp),
   "\nSum of squared discrepancies: \(ssd(x; $m))" 
;

task