#include <range/v3/view.hpp>
#include <range/v3/numeric/inner_product.hpp>
#include <iostream>

using namespace ranges::views;

int main(int argc, char * argv[]) {
    size_t N = 300000000;
    if(argc > 1) N = std::atol(argv[1]);

    auto for_ = [=](auto && f) {
        return iota(0ul, N) | transform([=](auto i) {
            return iota(0ul, N) | transform([=](auto j) {
                return f(i, j);
            });
        });
    };

    auto matrix = [=](size_t n) {
        return for_([tmp = 1. / n / n](ssize_t i, ssize_t j) {
            return tmp * (i - j) * (i + j);
        });
    };

    auto a = matrix(N), b = matrix(N);

    auto transpose = [=](const auto & m) {
        return for_([=](size_t i, size_t j) { return m[j][i]; });
    };

    auto mul = [=](const auto & a, const auto & b) {
        return for_([=](size_t i, size_t j) {
           return ranges::inner_product(a[i], b[j], 0.);            
        });
    };


    auto c = mul(a, transpose(b));

    fprintf(stderr, "%f\n", c[N/2][N/2]);
}

#[derive(Clone)]
struct Transform<I, F> {
    iter: I,
    f: F,
}

impl<B, I: Iterator, F> Iterator for Transform<I, F>
where
    F: FnMut(I::Item) -> B,
{
    type Item = B;

    fn next(&mut self) -> Option<Self::Item> {
        self.iter.next().map(&mut self.f)
    }

    fn nth(&mut self, i: usize) -> Option<Self::Item> {
        self.iter.nth(i).map(&mut self.f)
    }

    fn size_hint(&self) -> (usize, Option<usize>) {
        self.iter.size_hint()
    }
}

impl<B, I: ExactSizeIterator, F: FnMut(I::Item) -> B> ExactSizeIterator for Transform<I, F> {}

fn transform<I, O>(
    iter: impl ExactSizeIterator<Item = I> + Clone,
    f: impl Fn(I) -> O + Clone,
) -> impl ExactSizeIterator<Item = O> + Clone {
    Transform { iter, f }
}

fn for_<F, T>(
    n: usize,
    f: F,
) -> impl ExactSizeIterator<Item = impl ExactSizeIterator<Item = T>> + Clone
where
    F: Clone + Fn(usize, usize) -> T,
{
    transform(0..n, move |i| {
        let f = f.clone();
        transform(0..n, move |j| f(i, j))
    })
}

fn matrix(n: usize) -> impl ExactSizeIterator<Item = impl ExactSizeIterator<Item = f64>> + Clone {
    let tmp = 1.0 / n as f64 / n as f64;
    for_(n, move |i, j| {
        tmp * (i as f64 - j as f64) * (i as f64 + j as f64)
    })
}

fn transpose(
    m: impl ExactSizeIterator<Item = impl ExactSizeIterator<Item = f64>> + Clone,
) -> impl ExactSizeIterator<Item = impl ExactSizeIterator<Item = f64>> + Clone {
    for_(m.len(), move |i, j| {
        let mut m = m.clone();
        m.nth(j).unwrap().nth(i).unwrap()
    })
}

fn mul(
    a: impl ExactSizeIterator<Item = impl ExactSizeIterator<Item = f64>> + Clone,
    b: impl ExactSizeIterator<Item = impl ExactSizeIterator<Item = f64>> + Clone,
) -> impl ExactSizeIterator<Item = impl ExactSizeIterator<Item = f64>> + Clone {
    for_(a.len(), move |i, j| {
        let mut a = a.clone();
        let mut b = b.clone();
        a.nth(i)
            .unwrap()
            .zip(b.nth(j).unwrap())
            .map(|(l, r)| l * r)
            .sum::<f64>()
    })
}

fn main() {
    let n: usize = match std::env::args().nth(1) {
        Some(n) => n.parse().unwrap(),
        None => 3000,
    };

    let a = matrix(n);
    let b = matrix(n);
    let mut c = mul(a, transpose(b));
    println!("{}", c.nth(n / 2).unwrap().nth(n / 2).unwrap());
}