// spectral_norm.cm — Spectral norm kernel
// Pure numeric, dense array access, mathematical computation.

var math = use('math/radians')

function eval_a(i, j) {
  return 1.0 / ((i + j) * (i + j + 1) / 2 + i + 1)
}

function eval_a_times_u(n, u, au) {
  var i = 0
  var j = 0
  var sum = 0
  for (i = 0; i < n; i++) {
    sum = 0
    for (j = 0; j < n; j++) {
      sum += eval_a(i, j) * u[j]
    }
    au[i] = sum
  }
}

function eval_at_times_u(n, u, atu) {
  var i = 0
  var j = 0
  var sum = 0
  for (i = 0; i < n; i++) {
    sum = 0
    for (j = 0; j < n; j++) {
      sum += eval_a(j, i) * u[j]
    }
    atu[i] = sum
  }
}

function eval_ata_times_u(n, u, atau) {
  var v = array(n, 0)
  eval_a_times_u(n, u, v)
  eval_at_times_u(n, v, atau)
}

function spectral_norm(n) {
  var u = array(n, 1)
  var v = array(n, 0)
  var i = 0
  var vbv = 0
  var vv = 0

  for (i = 0; i < 10; i++) {
    eval_ata_times_u(n, u, v)
    eval_ata_times_u(n, v, u)
  }

  vbv = 0
  vv = 0
  for (i = 0; i < n; i++) {
    vbv += u[i] * v[i]
    vv += v[i] * v[i]
  }

  return math.sqrt(vbv / vv)
}

return {
  spectral_100: function(n) {
    var i = 0
    var result = 0
    for (i = 0; i < n; i++) {
      result = spectral_norm(100)
    }
    return result
  },

  spectral_200: function(n) {
    var i = 0
    var result = 0
    for (i = 0; i < n; i++) {
      result = spectral_norm(200)
    }
    return result
  }
}