/*
 * Copyright (c) 2021 James Cook <falsifian@falsifian.org>
 *
 * Permission to use, copy, modify, and/or distribute this software for any
 * purpose with or without fee is hereby granted.
 *
 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY
 * SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION
 * OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN
 * CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
 */

/*
 * Catalytic algorithm for graph reachability. This version cheats by using
 * doing recursion that uses O(log n) stack frames. reachability.c doesn't
 * cheat that way, but is more convoluted as a result.
 *
 * Throughout comments in this file, N(u, v, k) is the number of paths from
 * node u to node v with length at most 2^k.
 */

#include <cstdint>
#include <cstdlib>
#include <iostream>
#include <random>
#include <vector>

namespace {

struct algorithm_parameters {
  const std::vector<std::vector<bool>> *edges;
  int64_t num_nodes;
  /* The smallest power of 2 bigger than num_nodes. */
  int64_t rotation_size;
};

void
add_reachability(algorithm_parameters *, int, int,
		 std::vector<std::vector<uint64_t>> *,
		 std::vector<std::vector<uint64_t>> *,
		 std::vector<std::vector<std::vector<uint64_t>>> *);

/* Returns the smallest power of 2 bigger than num_nodes. */
int64_t rotation_size_needed(int64_t num_nodes) {
  int64_t result = 1;
  while (result < num_nodes + 1) {
    result *= 2;
  }
  return result;
}

void require(bool condition, const std::string &message) {
  if (!condition) {
    std::cerr << "require failed: " << message << std::endl;
    std::abort();
  }
}

/*
 * Rotates the elements of array by r.
 */
void
rotate_array(std::vector<uint64_t> *array, uint64_t r) {
  r = r % array->size();
  if (r == 0)
    return;

  /*
   * Rotation by r will form one or more cycles. We rotate one cycle at a
   * time.
   */
  int64_t lowest_nonbase_index = array->size();
  for (int64_t cycle_base = 0; cycle_base < lowest_nonbase_index; ++cycle_base) {
    uint64_t first_value = (*array)[cycle_base];
    int64_t i = cycle_base;
    for (;;) {
      int64_t next_index = (i + r) % array->size();
      if (next_index == cycle_base) {
	(*array)[i] = first_value;
	break;
      }
      (*array)[i] = (*array)[next_index];
      i = next_index;
      if (i < lowest_nonbase_index)
	lowest_nonbase_index = i;
    }
  }
}

/*
 * Recursive helper. For every (u, v), rotates (*catalyst_rotation)[u][v] by
 * increment * sum_w N(u, w, log_path_length - 1) * N(w, v, log_path_length -
 * 1). Requires log_path_length > 0.
 *
 * Makes no net change to catalyst_a or catalyst_b.
 */
void
rotate_by_reachability(algorithm_parameters *params, int log_path_length,
		       int increment,
		       std::vector<std::vector<std::vector<uint64_t>>> *catalyst_rotation,
		       std::vector<std::vector<uint64_t>> *catalyst_a,
		       std::vector<std::vector<uint64_t>> *catalyst_b) {

  for (int a_set = 0; a_set < 2; ++a_set) {
    for (int b_set = 0; b_set < 2; ++b_set) {
      /*
       * Invariants:
       *
       * If a_set == 0, then *catalyst_a has its original value (when this
       * instance of rotate_by_reachability was first called).
       *
       * If a_set == 1, for every u, v, (*catalyst_a)[u][v] has its original if
       * N(u, v, log_path_length - 1) == 0; otherwise it has its original value
       * + 1.
       *
       * The relation between b_set and catalyst_b is the same as for a_set and
       * catalyst_a.
       */
      int direction = increment * (a_set * 2 - 1) * (b_set * 2 - 1);
      for (int64_t u = 0; u < params->num_nodes; ++u) {
	for (int64_t v = 0; v < params->num_nodes; ++v) {
	  for (int64_t w = 0; w < params->num_nodes; ++w) {
	    rotate_array(&(*catalyst_rotation)[u][v],
			 direction * (*catalyst_a)[u][w] * (*catalyst_b)[w][v]);
	  }
	}
      }
      add_reachability(params, log_path_length - 1, 1 - 2 * b_set,
		       catalyst_b, catalyst_a, catalyst_rotation);
    }
    add_reachability(params, log_path_length - 1, 1 - 2 * a_set,
		     catalyst_a, catalyst_b, catalyst_rotation);
  }
}

/*
 * Recursive helper. Adds increment to *catalyst_output[u][v] for every pair (u,
 * v) such that N(u, v, log_path_length) != 0.
 *
 * Makes no net change to catalyst_aux or catalyst_rotation.
 */
void
add_reachability(algorithm_parameters *params, int log_path_length,
		 int increment,
		 std::vector<std::vector<uint64_t>> *catalyst_output,
		 std::vector<std::vector<uint64_t>> *catalyst_aux,
		 std::vector<std::vector<std::vector<uint64_t>>> *catalyst_rotation) {

  if (log_path_length == 0) {
    /* Base case. */
    for (int64_t u = 0; u < params->num_nodes; ++u)
      for (int64_t v = 0; v < params->num_nodes; ++v)
	if ((*params->edges)[u][v])
	  (*catalyst_output)[u][v] += increment;
    return;
  }

  /*
   * For each (u,v), subtract
   *    (*catalyst_rotation)[u][v][N(u, v, log_path_length)]
   * from (*catalyst_output)[u][v].
   */
  rotate_by_reachability(params, log_path_length, -1, catalyst_rotation,
			 catalyst_output, catalyst_aux);
  for (int u = 0; u < params->num_nodes; ++u)
    for (int v = 0; v < params->num_nodes; ++v)
      (*catalyst_output)[u][v] -=
	(*catalyst_rotation)[u][v][0];
  rotate_by_reachability(params, log_path_length, 1, catalyst_rotation,
			 catalyst_output, catalyst_aux);

  /*
   * Subtract increment from (*catalyst_rotation)[u][v][0].
   */
  for (int u = 0; u < params->num_nodes; ++u)
    for (int v = 0; v < params->num_nodes; ++v)
      (*catalyst_rotation)[u][v][0] -= increment;

  /*
   * For each (u,v), add
   *    increment +
   *    (*catalyst_rotation)[u][v][N(u, v, log_path_length)]
   * to (*catalyst_output)[u][v].
   */
  rotate_by_reachability(params, log_path_length, -1, catalyst_rotation,
			 catalyst_output, catalyst_aux);
  for (int u = 0; u < params->num_nodes; ++u)
    for (int v = 0; v < params->num_nodes; ++v)
      (*catalyst_output)[u][v] +=
	increment + (*catalyst_rotation)[u][v][0];
  rotate_by_reachability(params, log_path_length, 1, catalyst_rotation,
			 catalyst_output, catalyst_aux);

  /*
   * Add increment to (*catalyst_rotation)[u][v][0].
   */
  for (int u = 0; u < params->num_nodes; ++u)
    for (int v = 0; v < params->num_nodes; ++v)
      (*catalyst_rotation)[u][v][0] += increment;
}


/*
 * Given a directed graph, returns true iff there is a path from node 0 to node
 * 1.
 *
 * Nodes are numbered 0 through edge.size()-1. edge[u][v] should be true iff
 * there is an edge from u to v. edge[u].size() must equal edge.size() for all
 * nodes u.
 *
 * The elements of the three *catalyst arrays will be temporarily modified, but
 * restored to their original values by the time this function returns. The
 * first two must have length at least num_nodes^2 and the third must have
 * length at least num_nodes^3.
 */
bool reachable(const std::vector<std::vector<bool>> &edges,
	       std::vector<std::vector<uint64_t>> *catalyst_a,
	       std::vector<std::vector<uint64_t>> *catalyst_b,
	       std::vector<std::vector<std::vector<uint64_t>>> *catalyst_c) {
  struct algorithm_parameters params;
  params.edges = &edges;
  params.num_nodes = edges.size();
  params.rotation_size = rotation_size_needed(params.num_nodes);
  int max_log_path_length = 0;
  while ((1 << max_log_path_length) < params.num_nodes - 1) {
    ++max_log_path_length;
  }

  /* Make sure we have enough catalytic space. */
  require(catalyst_a->size() >= params.num_nodes, "catalyst_a size");
  require(catalyst_b->size() >= params.num_nodes, "catalyst_b size");
  require(catalyst_c->size() >= params.num_nodes, "catalyst_c size");
  for (int64_t i = 0; i < params.num_nodes; ++i) {
    require((*catalyst_a)[i].size() >= params.num_nodes, "catalyst_a element size");
    require((*catalyst_b)[i].size() >= params.num_nodes, "catalyst_b element size");
    require((*catalyst_c)[i].size() >= params.num_nodes, "catalyst_c element size");
    for (int64_t j = 0; j < params.num_nodes; ++j) {
      require((*catalyst_c)[i][j].size() == params.rotation_size, "catalyst_c subelement size");
    }
  }

  /* Change (*catalyst_a)[1] iff there's a path from 0 to 1. */
  add_reachability(&params, max_log_path_length, 1, catalyst_a, catalyst_b, catalyst_c);
  uint64_t mid_computation_value = (*catalyst_a)[0][1];
  /* Restore *catalyst_a to its starting params. */
  add_reachability(&params, max_log_path_length, -1, catalyst_a, catalyst_b, catalyst_c);
  return mid_computation_value != (*catalyst_a)[0][1];
}

}  // anonymous namespace

int main(int argc, char **argv) {
  /*
   * Read the graph.
   */
  int64_t num_nodes;
  std::cin >> num_nodes;
  std::vector<std::vector<bool>> edges{(std::size_t)num_nodes};
  for (int64_t i = 0; i < num_nodes; ++i) {
    edges[i].resize(num_nodes);
  }
  for (int64_t u = 0; u < num_nodes; ++u) {
    for (int64_t v = 0; v < num_nodes; ++v) {
      bool edge;
      std::cin >> edge;
      edges[u][v] = edge || u == v;
    }
  }

  /*
   * Allocate the catalytic space and fill it randomly.
   */
  int64_t rotation_size = rotation_size_needed(num_nodes);
  std::random_device rd;
  std::default_random_engine gen(rd());
  std::uniform_int_distribution<uint64_t> distrib;
  std::vector<std::vector<uint64_t>> catalyst_a{(std::size_t)num_nodes};
  std::vector<std::vector<uint64_t>> catalyst_b{(std::size_t)num_nodes};
  std::vector<std::vector<std::vector<uint64_t>>> catalyst_c{(std::size_t)num_nodes};
  for (int64_t i = 0; i < num_nodes; ++i) {
    catalyst_a[i].resize(num_nodes);
    catalyst_b[i].resize(num_nodes);
    catalyst_c[i].resize(num_nodes);
    for (int64_t j = 0; j < num_nodes; ++j) {
      catalyst_a[i][j] = distrib(gen);
      catalyst_b[i][j] = distrib(gen);
      catalyst_c[i][j].resize(rotation_size);
      for (int64_t k = 0; k < rotation_size; ++k) {
	catalyst_c[i][j][k] = distrib(gen);
      }
    }
  }

  /*
   * Save copies of the catalytic space, so we can verify it doesn't change.
   */
  std::vector<std::vector<uint64_t>> catalyst_a_backup = catalyst_a;
  std::vector<std::vector<uint64_t>> catalyst_b_backup = catalyst_b;
  std::vector<std::vector<std::vector<uint64_t>>> catalyst_c_backup = catalyst_c;

  std::cout
    << (reachable(edges, &catalyst_a, &catalyst_b, &catalyst_c) ? "Reachable." : "Not reachable.")
    << std::endl;

  require(catalyst_a == catalyst_a_backup, "catalyst_a changed");
  require(catalyst_b == catalyst_b_backup, "catalyst_b changed");
  require(catalyst_c == catalyst_c_backup, "catalyst_c changed");

  return 0;
}