import std.math.*
import std.time.*

func eval_A(i: Int64, j: Int64): Float64 {
  return Float64(1) / Float64(((i + j) * (i + j + 1) / 2 + i + 1))
}

func eval_A_times_u(N: Int64, u: Array<Float64>, Au: Array<Float64>): Unit {
  for (i in 0..N) {
    Au[i] = Float64(0)
    for (j in 0..N) {
      Au[i] = Au[i] + eval_A(i, j) * u[j]
    }
  }
}

func eval_At_times_u(N: Int64, u: Array<Float64>, Au: Array<Float64>): Unit {
  for (i in 0..N) {
    Au[i] = Float64(0)
    for (j in 0..N) {
      Au[i] = Au[i] + eval_A(j, i) * u[j]
    }
  }
}

func eval_AtA_time_u(N: Int64, u: Array<Float64>, AtAu: Array<Float64>): Unit {
  var v = Array<Float64>(N, { i => Float64(0) })
  eval_A_times_u(N, u, v)
  eval_At_times_u(N, v, AtAu)
}
func main() {
  var start = DateTime.now()
  var N: Int64 = 5500
  var u = Array<Float64>(N, { i => Float64(1) })
  var v = Array<Float64>(N, { i => Float64(0) })
  var vBv: Float64 = 0.0 
  var vv: Float64 = 0.0
  for (i in 0..10) {
    eval_AtA_time_u(N, u, v)
    eval_AtA_time_u(N, v, u)
  }
  for (i in 0..N) {
    vBv = vBv + u[i] * v[i]
    vv = vv + v[i] * v[i]
  }
  print("${ sqrt(vBv / vv) }\n")
  var end = DateTime.now()
  print("time = ${ end - start }\n")
}