spectralnorm

load core::math::sqrt;

fun a(i: int, j: int) -> float {
  imu d := (i + j) * (i + j + 1) / 2 + i + 1;

  1.0 / (d as float)
}

fun mult_av(v: []float, mut out: []float, n: int) {
  for i := 0..n {
    mut sum: float = 0.0;
    for j := 0..n {
      sum += a(i, j) * v[j];
    }

    out[i] = sum;
  }
}

fun mult_atv(v: []float, mut out: []float, n: int) {
  for i := 0..n {
    mut sum: float = 0.0;
    for j := 0..n {
      sum += a(j, i) * v[j];
    }

    out[i] = sum;
  }
}

fun main() {
  imu argv: []str = args();
  imu n: int = if argv.len > 0 {
    match argv[0].parse_int() {
      Option::Some(value) => value,
      Option::None => 100,
    }
  } else {
    100
  };

  mut u: []float = [];
  mut v: []float = [];
  mut tmp: []float = [];
  for i := 0..n {
    u.push(1.0);
    v.push(0.0);
    tmp.push(0.0);
  }

  for iter := 0..10 {
    mult_av(u, tmp, n);
    mult_atv(tmp, v, n);
    mult_av(v, tmp, n);
    mult_atv(tmp, u, n);
  }

  mut v_bv: float = 0.0;
  mut vv: float = 0.0;
  for i := 0..n {
    v_bv += u[i] * v[i];
    vv += v[i] * v[i];
  }

  showln(sqrt(v_bv / vv));
}