// implement 45min
#include <stdio.h>
#include <string.h>
#include <algorithm>
#include <iostream>
#include <math.h>
#include <assert.h>
#include <vector>
#include <queue>
#include <string>
#include <map>
#include <set>

using namespace std;
typedef long long ll;
typedef unsigned int uint;
typedef unsigned long long ull;
static const double EPS = 1e-9;
static const double PI = acos(-1.0);

#define REP(i, n) for (int i = 0; i < (int)(n); i++)
#define FOR(i, s, n) for (int i = (s); i < (int)(n); i++)
#define FOREQ(i, s, n) for (int i = (s); i <= (int)(n); i++)
#define FORIT(it, c) for (__typeof((c).begin())it = (c).begin(); it != (c).end(); it++)
#define MEMSET(v, h) memset((v), h, sizeof(v))

typedef vector<ll> Array;
typedef vector<Array> Matrix;
#include <complex>
typedef complex<double> Complex;

void PrintMatrix(const Matrix &mat) {
  REP(y, mat.size()) {
    REP(x, mat[y].size()) {
      cout << mat[y][x] << " ";
    }
    cout << endl;
  }
}

int psize = -1;
double ptheta = -1e+100;
vector<Complex> wtable;
// a.size must 2^x
vector<Complex> FFT(double theta, const vector<Complex> &a) {
  const int n = a.size();
  if (psize != n || ptheta != theta) {
    wtable.clear();
    REP(i, a.size() >> 1) {
      wtable.push_back(exp(i * theta * Complex(0, 1)));
    }
    psize = n;
    ptheta = theta;
  }
  int step = 1;
  vector<Complex> ret = a;
  for (int m = n; m >= 2; m >>= 1) {
    int mh = m >> 1;
    for (int i = 0; i < mh; i++) {
      Complex w = wtable[i * step];
      //Complex w = exp(i * theta * Complex(0, 1));
      for (int j = i; j < n; j += m) {
        int k = j + mh;
        Complex x = ret[j] - ret[k];
        ret[j] += ret[k];
        ret[k] = w * x;
      }
    }
    //theta *= 2;
    step *= 2;
  }
  int i = 0;
  for (int j = 1; j < n - 1; j++) {
    for (int k = n >> 1; k > (i ^= k); k >>= 1) {;}
    if (j < i) { swap(ret[i], ret[j]); }
  }
  return ret;
}

void Transpose(vector<vector<Complex> > &mat) {
  int n = mat.size();
  for (int y = 0; y < n; y++) {
    for (int x = y; x < n; x++) {
      swap(mat[y][x], mat[x][y]);
    }
  }
}
Matrix Convolution2D(const Matrix &lhs, const Matrix &rhs) {
  int n = 1;
  while (n < (int)max(lhs.size(), rhs.size()) * 2) { n <<= 1; }
  vector<vector<Complex> > temp1(n, vector<Complex>(n));
  vector<vector<Complex> > temp2(n, vector<Complex>(n));
  for (int y = 0; y < n / 2; y++) {
    for (int x = 0; x < n / 2; x++) {
      if (y < (int)lhs.size() &&
          x < (int)lhs[y].size()) {
        temp1[y][x] = Complex(lhs[y][x], 0);
      }
      if (y < (int)rhs.size() &&
          x < (int)rhs[y].size()) {
        temp2[y][x] = Complex(rhs[y][x], 0);
      }
    }
  }
  REP(y, n) {
    temp1[y] = FFT(2.0 * PI / n, temp1[y]);
    temp2[y] = FFT(2.0 * PI / n, temp2[y]);
  }
  Transpose(temp1);
  Transpose(temp2);
  REP(y, n) {
    temp1[y] = FFT(2.0 * PI / n, temp1[y]);
    temp2[y] = FFT(2.0 * PI / n, temp2[y]);
  }
  REP(y, n) {
    REP(x, n) {
      temp1[y][x] *= temp2[y][x];
    }
  }
  REP(y, n) {
    temp1[y] = FFT(-2.0 * PI / n, temp1[y]);
    REP(x, n) {
      temp1[y][x] /= n;
    }
  }
  Transpose(temp1);
  REP(y, n) {
    temp1[y] = FFT(-2.0 * PI / n, temp1[y]);
    REP(x, n) {
      temp1[y][x] /= n;
    }
  }
  Matrix ret(n, Array(n, 0));
  REP(y, n) {
    REP(x, n) {
      assert(fabs(temp1[y][x].real() - round(temp1[y][x].real())) < 1e-5);
      ret[y][x] = temp1[y][x].real() + 0.5;
      //ret[y][x] = temp1[y][x].real() / n / n + 0.5;
    }
  }
  return ret;
}


int n;
Matrix CreateCounterpartMatrix(const Matrix &lhs) {
  Matrix ret(n, Array(n, 0));
  for (int y = 0; y < n; y++) {
    for (int x = 0; x < n; x++) {
      ret[n - y - 1][n - x - 1] = lhs[y][x];
    }
  }
  return ret;
}

ll square(ll x) { return x * x; }
Matrix field;
map<ll, ll> hist;

int main() {
  while (scanf("%d", &n) > 0) {
    Matrix field = Matrix(n, Array(n, 0));
    hist.clear();
    REP(y, n) {
      REP(x, n) {
        ll c;
        int v = scanf("%lld", &c);
        assert(v == 1);
        field[y][x] = c;
        hist[0] -= (c * c + c) / 2;
      }
    }
    Matrix rhs = CreateCounterpartMatrix(field);
    Matrix ans = Convolution2D(field, rhs);
    REP(y, ans.size()) {
      REP(x, ans.size()) {
        if (ans[y][x] == 0) { continue; }
        ll d = square(n - y - 1) + square(n - x - 1);
        hist[d] += ans[y][x];
      }
    }
    long double ave = 0;
    ll cnt = 0;
    FORIT(it, hist) {
      if (it->second == 0) { continue; }
      if (it->first != 0) { it->second /= 2; }
      assert(it->second >= 0);
      ave += (long double)sqrt(it->first) * it->second;
      cnt += it->second;
    }
    ave /= cnt;
    printf("%.10Lf\n", ave);
    int index = 0;
    FORIT(it, hist) {
      if (it->second == 0) { continue; }
      printf("%lld %lld\n", it->first, it->second);
      index++;
      if (index == 10000) { break; }
    }
  }
}