solution sketch: Put the four sequences of integers on the coordinate plane, with $A$ being the positive y axis, $B$ the negative x axis, $C$ the negative y axis, and $D$ the positive $x$ axis. Note that the condition modifies slightly when we do this; quadrants 2 and 4 are "reversed." Redefine f(x,y) to include the negative numbers, taking into account the fact that the condition changes. Now, we just need to prove that all terms in the first quadrant are integers if the terms in all other quadrants are integers. To do this, we just consider a few cases; but the basic idea is the following: given $f(x, y), f(x,y+1), f(x+1,y)$, we can prove that $f(x,y+1)f(x+1,y)\equiv \pm\frac{1}{f(x-1,y)f(x,y-1)}\equiv -1\mod f(x,y)$ where the $\pm$ is a plus or a minus depending on which quadrant we're working in, which would imply the existence of an integer $f(x+1,y+1)$ that satisfies the conditions. We have to take a little care here because $f(x,y)$ could be zero or a negative number but those cases aren't so different and lead to the same conclusion. In this fashion, we can prove that everything in the first quadrant is an integer, so $D \sim A$. Done!

Official solutions post (also at http://web.evanchen.cc/problems.html): We present two solutions. In what follows, we say $(A, B)$ form a great pair if $A \sim B$. First solution (Nikolai Beluhov) Let $k = a_0 = b_0 = c_0 = d_0$. We let $f$, $g$, $h$ be great functions for $(A,B)$, $(B,C)$, $(C,D)$ and write the following infinite array: \[ \begin{bmatrix} & \vdots & \vdots & {\color{blue} b_3} & \vdots & \vdots & \\[5pt] \cdots & g(2,2) & g(2,1) & {\color{blue}b_2} & f(1,2) & f(2,2) & \cdots \\[5pt] \cdots & g(1,2) & g(1,1) & {\color{blue}b_1} & f(1,1) & f(2,1) & \cdots \\[5pt] {\color{blue}c_3} & {\color{blue}c_2} & {\color{blue}c_1} & {\color{blue}k} & {\color{blue}a_1} & {\color{blue}a_2} & {\color{blue}a_3} \\[5pt] \cdots & h(2,1) & h(1,1) & {\color{blue}d_1} & & & \\[5pt] \cdots & h(2,2) & h(1,2) & {\color{blue}d_2} & & & \\[5pt] & \vdots & \vdots & {\color{blue}d_3} & & & \ddots \\[5pt] \end{bmatrix} \]The greatness condition is then equivalent to saying that any $2 \times 2$ sub-grid has determinant $\pm1$ (the sign is $+1$ in two quadrants and $-1$ in the other two), and we wish to fill in the lower-right quadrant. To this end, it suffices to prove the following. Lemma: Suppose we have a $3 \times 3$ sub-grid \[ \begin{bmatrix} a & b & c \\ x & y & z \\ p & q & \end{bmatrix} \]satisfying the determinant conditions. Then we can fill in the ninth entry in the lower right with an integer while retaining greatness. Proof. We consider only the case where the $3 \times 3$ is completely contained inside the bottom-right quadrant, since the other cases are exactly the same (or even by flipping the signs of the top row or left column appropriately). If $y = 0$ we have $-1 = bz = bx = xq$, hence $qz = -1$, and we can fill in the entry arbitrarily. Otherwise, we have $bx \equiv xq \equiv bz \equiv -1 \pmod{y}$. This is enough to imply $qz \equiv -1 \pmod y$, and so we can fill in the integer $\frac{qz+1}{y}$. $\blacksquare$ Remark: In this case (of all $+1$ determinants), I think it turns out the bottom entry is exactly equal to $qza - cyp - c - p$, which is obviously an integer. Second solution (Ankan Bhattacharya) We will give an explicit classification of great sequences: Lemma: The pair $(A,B)$ is great if and only if $a_0 = b_0$, $a_0 \mid a_1b_1 + 1$, and $a_n \mid a_{n-1} + a_{n+1}$ and $b_n \mid b_{n-1} + b_{n+1}$ for all $n$. Proof. [Proof of necessity] It is clear that $a_0 = b_0$. Then $a_0f(1, 1) - a_1b_1 = 1$, i.e. $a_0 \mid a_1b_1 + 1$. Now, focus on six entries $f(x, y)$ with $x \in \{n-1, n, n+1\}$ and $y \in \{0, 1\}$. Let $f(n-1, 1) = u$, $f(n, 1) = v$, and $f(n+1, 1) = w$, so \begin{align*} v a_{n-1} - u a_n & = 1,\\ w a_n - v a_{n+1} & = 1. \end{align*}Then \[u + w = \frac{v(a_{n-1} + a_{n+1})}{a_n}\]and from above $\gcd(v, a_n) = 1$, so $a_n \mid a_{n-1} + a_{n+1}$; similarly for $b_n$. (If $a_n = 0$, we have $va_{n-1} = 1$ and $va_{n+1} = -1$, so this is OK.) $\blacksquare$ Proof. [Proof of sufficiency] Now consider two sequences $a_0, a_1, \dots$ and $b_0, b_1, \dots$ satisfying our criteria. We build a great function $f$ by induction on $(x, y)$. More strongly, we will assume as part of the inductive hypothesis that any two adjacent entries of $f$ are relatively prime and that for any three consecutive entries horizontally or vertically, the middle one divides the sum of the other two. First we set $f(1, 1)$ so that $a_0 f(1, 1) = a_1 b_1 + 1$, which is possible. Consider an uninitialized $f(s, t)$; without loss of generality suppose $s \ge 2$. Then we know five values of $f$ and wish to set a sixth one $z$, as in the matrix below: \[ \begin{matrix} u & x\\ v & y\\ w & z \end{matrix} \](We imagine $a$-indices to increase southwards and $b$-indices to increase eastwards.) If $v \neq 0$, then the choice $y \cdot \frac{u + w}{v} - x$ works as $uy - vx = 1$. If $v = 0$, it easily follows that $\{u, w\} = \{1, -1\}$ and $y = w$ as $yw = 1$. Then we set the uninitialized entry to anything. Now we verify that this is compatible with the inductive hypothesis. From the determinant 1 condition, it easily follows that $\gcd(w, z) = \gcd(v, z) = 1$. The proof that $y \mid x + z$ is almost identical to a step performed in the ``necessary'' part of the lemma and we do not repeat it here. By induction, a desired great function $f$ exists. $\blacksquare$ We complete the solution. Let $A$, $B$, $C$, and $D$ be integer sequences for which $(A, B)$, $(B, C)$, and $(C, D)$ are great. Then $a_0 = b_0 = c_0 = d_0$, and each term in each sequence (after the zeroth term) divides the sum of its neighbors. Since $a_0$ divides all three of $a_1b_1 + 1$, $b_1c_1 + 1$, and $c_1d_1 + 1$, it follows $a_0$ divides $d_1a_1 + 1$, and thus $(D, A)$ is great as desired. Remark: To simplify the problem, we may restrict the codomain of great functions and elements in great pairs of sequences to $\mathbb{Z}_{> 0}$. This allows the parts of the solution dealing with zero entries to be ignored. Remark: Of course, this solution also shows that any odd path (in the graph induced by the great relation on sequences) completes to an even cycle. If we stipulate that great functions must have $f(0, 0) = \pm 1$, then even paths complete to cycles as well. Alternatively, we could change the great functional equation to \[f(x + 1, y + 1) f(x, y) - f(x + 1, y) f(x, y + 1) = -1.\] A quick counterexample to transitivity of $\sim$ as is without the condition $f(0,0)=1$, for concreteness: let $a_n = c_n = 3+n$ and $b_n = 3+2n$ for $n \ge 0$.

We write $a\mid b$ if there exists an integer $\ell$ such that $b=a\ell$. Note that $\ell$ may not be unique if $a=b=0$. Suppose $\gcd(a,b)=1$. Then, we have the following true statements (for nonzero elements they are trivial, check the $0$ cases by hand): \[[ab\mid c\iff(a\mid c\text{ and }b\mid c)]\text{ and }[a\mid bc\implies a\mid c].\]We'll be using these in our proof. Write $A\approx B$ if $a_0\mid 1+a_1b_1$ and $a_k\mid a_{k-1}+a_{k+1}$, $b_k\mid b_{k-1}+b_{k+1}$ for all $k\ge 1$. We claim that $A\sim B\iff A\approx B$ (this is the key part of the solution). First suppose $A\approx B$. Let $P(m,n)$ be the FE. By $P(0,0)$, we see that \[f(1,1)a_0=1+a_1b_1,\]so since $a_0\mid 1+a_1b_1$, there exists some choice of $f(1,1)\in\mathbb{Z}$ that satisfies $P(0,0)$ (note that it may not be unique). We claim that we can choose $f(1,1),\ldots,f(n,1)\in\mathbb{Z}$ such that $P(0,0),P(1,0),\ldots,P(n-1,0)$ are all satisfied. We proceed by induction on $n$, with the case $n=1$ being solved above. So suppose the claim is true for $n-1$. We see that \[P(n-1,0)\iff f(n,1)a_{n-1}=1+a_n f(n-1,1)\]and \[P(n-2,0)\implies f(n-1,1)a_{n-2}=1+a_{n-1}f(n-2,1).\]By the converse of Bezout's lemma, we have $\gcd(a_{n-2},a_{n-1})=1$ from $P(n-2,0)$. If $a_{n-2}=0$, then we have from $P(n-2,0)$ (which we know is satisfied by the inductive hypothesis) that $a_{n-1}=\pm 1$. Thus, we can choose $f(n,1)=\pm(1+a_n f(n-1,1))$ and satisfy $P(n-1,0)$, which proves the inductive step for $a_{n-2}=0$ case. So suppose $a_{n-2}\ne 0$. Multiplying $P(n-1,0)$ by $a_{n-2}$ is then equivalent to $P(n-1,0)$, so we have \[P(n-1,0)\iff f(n,1)a_{n-1}a_{n-2}=a_{n-2}+a_n a_{n-2}f(n-1,1),\]or \[P(n-1,0)\iff f(n,1)a_{n-1}a_{n-2}=a_{n-2}+a_n+a_na_{n-1}f(n-2,1)\]where we used the true $P(n-2,0)$. Note that $a_{n-1}$ divides the RHS since $a_{n-1}\mid a_{n-2}+a_n$, and $a_{n-2}$ divides the RHS since we got the RHS as $a_{n-2}$ times some integer. Thus, since $\gcd(a_{n-1},a_{n-2})$, we see that $a_{n-1}a_{n-2}$ divides the RHS, which directly means that we can choose inegral $f(n,1)$ to satisfy this latest modified version of $P(n-1,0)$. This completes the inductive step, so the claim as well. The claim now implies we can choose $f(1,1),f(2,1),\ldots$ and $f(1,1),f(1,2),\ldots$ such that $P(m,n)$ is satisfied for all $(m,n)$ such that $\min\{m,n\}=0$. Let $A'=(f(1,1),f(2,1),\ldots)$ and $B'=(f(1,1),f(1,2),\ldots)$. We show that $A'\approx B'$. We have \[P(1,0)\implies f(2,1)a_1=a_2f(1,1)+1\]and \[P(0,1)\implies f(1,2)b_1=b_2f(1,1)+1.\]Multiplying the two, we see that there is some $\ell_1\in\mathbb{Z}$ such that \[f(1,2)f(2,1)a_1b_1=1+f(1,1)\ell_1,\]so \[(f(1,2)f(2,1)+1)a_1b_1=(a_1b_1+1)+f(1,1)\ell_1.\]By $P(0,0)$, we have $a_1b_1+1=f(1,1)a_0$, so there is some $\ell_2$ such that \[(f(1,2)f(2,1)+1)a_1b_1=f(1,1)\ell_2.\]Since $\gcd(a_1,f(1,1))=\gcd(b_1,f(1,1))=1$ by the converse of Bezout on $P(0,0)$, we see that $f(1,1)\mid f(1,2)f(2,1)+1$. All we need to show now is that $f(k,1)\mid f(k-1,1)+f(k+1,1)$ and the one where we switch the variables, but that follows by symmetry. We see that \[P(k,0)\implies f(k+1,1)a_k=1+a_{k+1}f(k,1)\]and \[P(k-1,0)\implies f(k-1,1)a_k=-1+a_{k-1}f(k,1).\]Adding the two gives \[a_k(f(k+1,1)+f(k-1,1))=f(k,1)(a_{k+1}+a_{k-1})=\ell a_kf(k,1)\]since $a_k\mid a_{k+1}+a_{k-1}$. Thus, if $a_k\ne 0$, we're set, so suppose $a_k=0$. We see then from $P(k,0)$ that $f(k,1)=\pm 1$, so certainly $f(k,1)\mid f(k-1,1)+f(k+1,1)$. Thus, $A'\approx B'$, so by the same logic, we can choose $f(2,2),f(2,3),\ldots$ and $f(2,2),f(3,2),\ldots$ such that $P(m,n)$ where $\min\{m,n\}=1$ are all satisfied. We can continue in this manner to show that we can choose $f$ such that all $P(m,n)$ are satisfied by inducting on $\min\{m,n\}$. Thus, $A\approx B\implies A\sim B$. We now show that $A\sim B\implies A\approx B$. Fortunately, this is a lot easier. We have \[P(0,0)\implies f(1,1)a_0=1+a_1b_1\]which directly implies $a_0\mid 1+a_1b_1$. We showed before that \[P(n-1,0),P(n-2,0)\implies f(n,1)a_{n-1}a_{n-2}=a_{n-2}+a_n+a_na_{n-1}f(n-2,1).\]Since $a_{n-1}$ divides the LHS, it must divide the RHS. But the last term of the RHS is clearly divisible by $a_{n-1}$, so we have \[a_{n-1}\mid a_n+a_{n-2},\]and similarly for $b_n$. This shows that $A\approx B$, as desired. Thus, \[\boxed{A\sim B\iff A\approx B}.\]Showing that $A\approx B$, $B\approx C$, and $C\approx D$ implies $D\approx A$ is surprisingly trivial. Note that if $x=a_0=b_0=c_0=d_0$ is zero, then \[a_1b_1=b_1c_1=c_1d_1=-1.\]It's easy to see that this implies $a_1d_1=-1$. If $x=1$ then $x\mid 1+a_1d_1$, so suppose $|x|\ge 2$. Then, \[a_1b_1\equiv b_1c_1\equiv c_1d_1\equiv-1\pmod{x}.\]Thus, \[-1\equiv \frac{a_1b_1}{b_1c_1}\cdot(c_1d_1)\equiv a_1d_1\pmod{x},\]where the divisions are valid since $x\ge 2$ (this implies that if $b_1$ or $c_1$ was $0$, then $b_1c_1\not\equiv -1\pmod{x}$), so $x\mid 1+a_1d_1$. We obviously have $a_k\mid a_{k-1}+a_{k+1}$ and $d_k\mid d_{k-1}+d_{k+1}$ from $A\approx B$ and $C\approx D$, so we have $D\approx A$, as desired. $\blacksquare$

Here's my solution. Denote $P(m,n)$ as the assertion of $m$ and $n$ to the given functional equation. $\textbf{Lemma 01.}$ We have $\frac{f(m+2,k) + f(m,k)}{f(m+1,k)}$ being a constant for as $k$ vary while $m$ remains fixed. $\textit{Proof.}$ Notice that $P(m,n)$ and $P(m+1,n)$ gives us \[ f(m,n+1) f(m+1,n) = -1 + f(m,n) f(m+1,n+1) \]\[ f(m+1,n) f(m+2,n+1) = 1 + f(m+1,n+1) f(m+2,n) \]Add them, we will have $f(m+1,n+1) ( f(m,n) + f(m+2,n) ) = f(m+1,n) ( f(m,n+1) + f(m+2,n+1) )$. Therefore, we have \[ \frac{f(m,n) + f(m+2,n)}{f(m+1,n)} = \frac{f(m,n+1) + f(m+2,n+1)}{f(m+1,n+1)} \]By induction, one can prove that \[ \frac{f(m,k) + f(m+2,k)}{f(m+1,k)} = \frac{f(m,0) + f(m+2,0)}{f(m+1,0)} = \frac{a_m+a_{m+2}}{a_{m+1}} \] $\textbf{Lemma 02.}$ Each constant mentioned earlier in $\textbf{Lemma 01}$ must be an integer. $\textit{Proof.}$ Suppose for the sake of contradiction, for some $m$, the constant is not an integer. As the constant is in the form of $\frac{a_m + a_{m+2}}{a_{m+1}}$ where $a_i$ are integers, then if it is not an integer, then it must be expressed as $\frac{p}{q}$ where $q > 1$, and $p,q \in \mathbb{N}$. Notice that this results $q | f(m+1,k)$ for a fixed $m$ and $k$ vary. Now, notice that $q | f(m+1,n+1) f(m,n) - f(m+1,n) f(m,n+1) = 1$. So we must have $q | 1$, which is a contradiction. We now consider the relation $A \sim B$ in a table, where each cell $(i,j)$ denote $f(i,j)$. Now, we'll prove that the whole table is uniquely defined by knowing all of the values of $f(m,n)$ where $\min \{ m,n \} \le 2$. Using $\textbf{Lemma 01.}$, we can extend the above given condition to identify all the other $f(m,n)$ for any $m,n \in \mathbb{N}$. Notice that the condition $f(1,1)f(0,0)- f(0,1)f(1,0) = 1$ and $\frac{f(m,k) + f(m+2,k)}{f(m+1,k)} \in \mathbb{Z}$ being constant is sufficient to satisfy all of the condition, and in fact this is necessary. This is because: \begin{align*} a \ b \\ c \ d \\ e \ f \\ \end{align*}Consider the first $2 \times 3$ table, we must have $da - bc = 1$ and $cf - de = 1$. Now, since $\frac{a+e}{c} = \frac{b+f}{d} = k \in \mathbb{N}$, this means that $cf - de = c(kd - b) - d(kc - e) = ad - bc = 1$. Continuing doing this, we get that this two condition is sufficient. This two condition has been proven to be necessary. Now, this gives us that if $A \sim B$, then we must have $\frac{a_1 b_1 + 1}{a_0} \in \mathbb{N}$. But, similarly, since $B \sim C, C \sim D$, we must have \[ \frac{b_1 c_1 + 1}{b_0} , \frac{c_1 d_1 + 1}{c_0} \in \mathbb{N} \]Therefore, since $a_0 = b_0 = c_0 = d_0$, we have $\gcd(c_0,c_1) = 1$ and therefore, $c_0 | d_1 - b_1$, Thus, we must have $d_0 | (a_1 b_1 + 1) + (d_1 - b_1) a_1 = a_1 d_1 + 1 $. Therefore, $\frac{a_1d_1 + 1}{d_0} \in \mathbb{N}$, which therefore result $D \sim A$ as well.

Write up took forever, but problem wasn't that hard Let $f$ be the great function for $A\sim B$, let $g$ be the great function for $B\sim C$, and let $h$ be the great function for $C\sim D$. Also let $x=a_0=b_0=c_0=d_0$. Note that $f(n,0)=a_n, f(0,n)=b_n, g(n,0)=b_n, g(0,n)=c_n, h(n,0)=c_n, h(0,n)=d_n$. Consider the following matrix\[\begin{bmatrix} & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \ \\ \ldots & g(3,3) & g(3,2) & g(3,1) & b_3 & f(1,3) & f(2,3) & f(3,3) & \ldots \\ \ldots & g(2,3) & g(2,2) & g(2,1) & b_2 & f(1,2) & f(2,2) & f(3,2) & \ldots \\ \ldots & g(1,3) & g(1,2) & g(1,1) & b_1 & f(1,1) & f(2,1) & f(3,1) & \ldots \\ \ldots & c_3 & c_2 & c_1 & x & a_1 & a_2 & a_3 & \ldots \\ \ldots & h(3,1) & h(2,1) & h(1,1) & d_1 & \ & \ & \ & \ \\ \ldots & h(3,2) & h(2,2) & h(1,2) & d_2 & \ & \ & \ & \ \\ \ldots & h(3,3) & h(2,3) & h(1,3) & d_3 & \ & \ & \ & \ \\ \ & \vdots & \vdots & \vdots & \vdots & \ & \ & \ & \ \\ \end{bmatrix}\]It suffices to show that we can fill in the empty cells so that all $2\times 2$ sub-matrices have determinant $\pm1$ ($+1$ for the $g$ and empty quadrants and $-1$ for the $f$ and $h$ quadrants). In fact all we have to do is to prove the following lemma: Lemma: Suppose we have the matrix\[\begin{bmatrix} i & j & k \\ l & m & n \\ o & p & \ \\ \end{bmatrix},\]where all $2\times 2$ sub-matrices that don't include the empty cell have determinant $\pm1$. Then we can fill the empty cell with an integer $q$ so that $\begin{bmatrix} m & n \\ p & q \\ \end{bmatrix}$ has determinant $1$ (since all empty cells that need to be filled are in the empty quadrant). Proof: It suffices to show that $np\equiv- 1\pmod m$. Note that $lp\equiv \pm 1\pmod m$, $lj\equiv \pm 1\pmod m$, and $jn\equiv \pm 1\pmod m$. Thus, $l^2j^2np=(lj)^2np\equiv \pm 1\pmod m\implies np\equiv \pm 1\pmod m$. Note that since $(lj)^2\equiv 1\pmod m$, our lemma is true if $lp\cdot lj\cdot jn\equiv -1\pmod m$. We will do casework on the quadrant of $\begin{bmatrix} i & j \\ l & m \\ \end{bmatrix}$. Case 1: $\begin{bmatrix} i & j \\ l & m \\ \end{bmatrix}$ is in quadrant $g$ (in other words, quadrant II). Then $im-jl=1\implies jl\equiv -1\pmod m$ Note that $\begin{bmatrix} l & m \\ o & p \\ \end{bmatrix}$ is in quadrant $h$ and $\begin{bmatrix} j & k \\ m & n \\ \end{bmatrix}$ is in quadrant $f$. Thus, $lp\equiv jn\equiv -1\pmod m$. Then $lp\cdot lj\cdot jn\equiv -1\pmod m$, as desired. Case 2: $\begin{bmatrix} i & j \\ l & m \\ \end{bmatrix}$ is in quadrant $h$ (in other words quadrant III) Then $im-jl=-1\implies jl\equiv 1\pmod m$. Note that $\begin{bmatrix} l & m \\ o & p \\ \end{bmatrix}$ is in quadrant $h$ and $\begin{bmatrix} j & k \\ m & n \\ \end{bmatrix}$ is in the empty quadrant. Thus, $lp\equiv -1\pmod m$ and $jn\equiv 1\pmod m$. So $lp\cdot jl\cdot jn\equiv -1\pmod m$, as desired. Case 3: $\begin{bmatrix} i & j \\ l & m \\ \end{bmatrix}$ is in quadrant $f$ (in other words quadrant I) Then $im-jl=-1\implies jl\equiv 1\pmod m$. Note that $\begin{bmatrix} l & m \\ o & p \\ \end{bmatrix}$ is in the empty quadrant and $\begin{bmatrix} j & k \\ m & n \\ \end{bmatrix}$ is in quadrant $f$. Thus, $lp\equiv 1\pmod m$ and $jn\equiv -1\pmod m$. So $lp\cdot jl\cdot jn\equiv -1\pmod m$, as desired. Case 4: $\begin{bmatrix} i & j \\ l & m \\ \end{bmatrix}$ is in the empty quadrant (in other words quadrant IV). Then $im-jl=1\implies jl\equiv -1\pmod m$. Note that both $\begin{bmatrix} l & m \\ o & p \\ \end{bmatrix}$ and $\begin{bmatrix} j & k \\ m & n \\ \end{bmatrix}$ are also in the empty quadrant. So $lp\equiv jn\equiv 1\pmod m$. Thus, $lp\cdot jl\cdot jn\equiv -1\pmod m$, as desired. We have exhausted all cases, so we are done. $\blacksquare$

Notice that by taking a difference of the desired expression, we get that $\frac{f(m+1, n+1)}{f(m+1, n)} = \frac{f(m+2, n+1)+f(m, n+1)}{f(m+2, n)+f(m, n)}$. Furthermore, we have $gcd(f(m+1, n+1), f(m+1, n)) = 1$ so for all $m, n$ we have $f(m+1, 0) \mid f(m, 0)+f(m+2, 0)$ and $f(0, n+1) \mid f(0, n)+f(0, n+2)$. This suffices by induction if $f(1, 1)$ can be chosen. Thus we are done since $D$ satisfies the division condition and $f(1, 1)$ when constructing the function between $A, D$ is an integer by modulus.

Main Claim: If $a,b,c,d,e,f,g,h$ are integers such that $$ae-bd=1,bf-ce=1,dh-eg=1,$$then there exists an integer $k$ such that $$ek-hf=1.$$The condition that $k$ exists is just $hf\equiv -1\pmod e$. From the three equalities taken mod $e$, $$db\equiv -1\pmod e, bf\equiv 1\pmod e, dh\equiv 1\pmod e.$$Multiplying the second and third, $$bfdh\equiv 1\pmod e.$$However, since $db\equiv -1\pmod e$, $hf\equiv -1\pmod e$, as desired. Now, make an infinite grid, with $a_0=b_0=c_0=d_0$ in the center, and the x-axis has $a_i$ and $c_i$, and the y-axis has $b_i$ and $d_i$. Fill in the first three quadrants using the the fact that $A\sim B$, $B\sim C$, and $C\sim D$. Then, by our main claim, if a 3x3 grid containing at least one cell in the 4th quadrant has everything except the bottom right filled, we can fill the bottom right (If the top row crosses into the first quadrant or left column crosses into the third, flip the signs of the top row or left column respectively). Thus, we can fill up the fourth quadrant, so we are done.

Cute problem. Sleek statement and interesting ideas. Sequences $A,B$ are called great if $A\sim B$. We give a classification of all such great sequences. Key Claim: Any pair of sequences $(A, B)$ is great if and only if $a_0=b_0$, $a_0\mid a_1b_1+1$, and $a_n\mid a_{n+1}+a_{n-1}$ and $b_n\mid b_{n-1}+b_{n+1}$ for all non-negative integers $n$. Proof. We first show that the above condition is necessary. As given in the problem, $a_0=b_0$. Note, $f(1,1)=\frac{a_1b_1+1}{a_0}$, which means $a_0\mid a_1b_1+1$ must be true as well. For any positive integer $n$, we have \[f(n+1,1)=\frac{f(n,1)f(n+1,0)+1}{f(n,0)}=\frac{f(n,1)a_{n+1}+1}{a_n} \ \ \ \text{and} \ \ \ f(1,n+1)=\frac{f(1,n)f(0,n+1)+1}{f(0,n)}=\frac{f(1,n)b_{n+1}+1}{b_n}\]Hence, clearly, $\gcd(a_n,a_{n+1})=\gcd(b_n,b_{n+1})=1$. So, we may write $f(n+1,1)\equiv\frac1{a_n}\pmod{a_{n+1}}$, which means \[f(n,1)\equiv\frac1{a_{n-1}}\pmod{a_n}\Longrightarrow 0\equiv a_nf(n+1,1)\equiv f(n,1)a_{n+1}+1\equiv\frac{a_{n-1}+a_{n+1}}{a_{n-1}}\pmod{a_n}\]Since $\gcd(a_{n-1},a_n)=1$, we have $a_n\mid a_{n-1}+a_{n+1}$ and similarly we have $b_n\mid b_{n-1}+b_{n+1}$, as claimed. So, the conditions outlined in the claim are indeed necessary. $\blacksquare$ We now show that these conditions are sufficient. We first show that $f(n,1)$ and $f(1,n)$ are ensured to be integers. Note that since $a_0\mid a_1b_1+1$, we have $\gcd(a_0,a_1)=1$. Since $a_n\mid a_{n-1}+a_{n+1}$, we see that $\gcd(a_n,a_{n+1})=1$ for all $n$ holds. So, our proof above works backwards as well, and hence ensures that $f(n,1)$ and $f(1,n)$ are integers. Let $k\geq1$ be an integer. We show that $f(n,k)$ and $f(k,n)$ are integers, by inducting on $k$. The base case, i.e. $k=1$, has been established. In our induction hypothesis, we assume that $f(n,k)$ and $f(k,n)$ are integers for all $n$, and we will show that $f(n,k+1)$ and $f(k+1,n)$ are integers for all $n$. Note that to ensure that $f(n,k+1)$ and $f(k+1,n)$ are integers, the sufficient conditions are: $f(0,k)\mid f(0,k+1)f(1,k)+1$ (and $f(k,0)\mid f(k+1,0)f(k,1)+1$). $\gcd(f(n,k),f(n+1,k))=1$ for all integers $n$ (and $\gcd(f(k,n),f(k,n+1))=1$). $f(n,k)\mid f(n-1,k)+f(n+1,k)$ for all integers $n$ (and $f(k,n)\mid f(k,n-1)+f(k,n+1)$ for all integers $n$). Note that the conditions above in brackets are similar to their counterparts, so we will prove only prove them. Note that $f(1,k+1)=\frac{f(0,k+1)f(1,k)+1}{f(0,k)}$ is an integer, which means $\gcd(f(0,k+1),f(1,k+1))=1$. We will perform induction on $n$ to show the second bullet is always true. Base case, $n=0$ has been established. Note that $f(n-1,k)=\frac{f(n,k)f(n-1,k-1)-1}{f(n,k-1)}$ and $f(n+1,k)=\frac{f(n+1,k-1)f(n,k)+1}{f(n,k-1)}$. Adding them, we get \[f(n-1,k)+f(n+1,k)=\frac{f(n,k)(f(n-1,k-1)+f(n+1,k-1))}{f(n,k-1)}\]\[\Longrightarrow\frac{f(n-1,k)+f(n+1,k)}{f(n,k)}=\frac{f(n-1,k-1)+f(n+1,k-1)}{f(n,k-1)}\]Since all $f(n,k)$ are integers, it means $f(n,k-1)\mid f(n-1,k-1)+f(n+1,k-1)$. Subsequently, we see $f(n,k)\mid f(n-1,k)+f(n+1,k)$, which proves the third bullet. Moreover, by our induction hypothesis, we have $\gcd(f(n,k),f(n-1,k))=1$ and since $f(n+1,k)\equiv-f(n-1,k)\pmod{f(n,k)}$, we have $\gcd(f(n,k),f(n+1,k))=1$ as well. So, the induction for second bullet is complete as well. Hence, our "key claim" is true. $\blacksquare$ Finally, note that since $B,C$ and $C,D$ are great pairs as well, we have $a_0=b_0=c_0=d_0$, and $b_0\mid b_1c_1+1$, $b_0=c_0\mid c_1d_1+1$. Subtracting, we see $b_1\equiv d_1\pmod{b_0}$ and similarly, $a_1\equiv c_1\pmod{b_0}$. It follows that $d_1a_1+1\equiv b_1c_1+1\equiv 0\pmod{b_0}$. Finally, since $C,D$ are great pairs, we have $d_n\mid d_{n-1}+d_{n+1}$ as well. So, by our "key claim", $D,A$ are great pairs. $\blacksquare$ Remark. Any even cycle works, but no odd cycles work.

Place $A$ on the positive $y$ axis, $B$ on the positive $x$ axis, $C$ on the negative $y$ axis, and $D$ on the negative $x$ axis. This visualization means that there's a great function propagating out from each quadrant. We say a lattice point is solved if its value in the great function is an integer. Thus all of quadrants $1,3,4$ are solved and we want to solve quadrant $2$. Claim: $(-1,i)$ and $(j,1)$ for $i>0$ and $j<0$ is solved. Proof. First, we show $(-1,1)$ is solved. Note that if $x=a_0,b_0,c_0,d_0$ then the conditions give that $a_1b_1, b_1c_1, c_1d_1 \equiv -1 \pmod x$ so $a_1d_1 \equiv -1 \pmod x$ so $(-1,1)$ is solved. To show $(-1,i)$ is solved we show this inductively. If the value of $(-1,i)$ is $n$, then $(-1,i+1)$ has value equal to $n_1=\frac{1+na_{i+1}}{a_i}$. Then $n_1\equiv \frac{ 1 }{a_i} \pmod{a_{i+1}}$. Thus to solve $(-1,i+2)$ we need \[0 \equiv a_{i+2}n_1 + 1 \equiv \frac{a_{i+2}}{a_{i}} + 1 \pmod {a_{i+1}}\]which is independent of $D$ and thus true by $A \sim B$. The other claim is true by symmetry. $\blacksquare$ To finish, consider this matrix: \[\begin{bmatrix} & \frac 1y & \\ \frac 1x & n & y\\ & x & \end{bmatrix}\]which is a subset of quadrant $2$ taken modulo $n$. Note that $n=\frac{xy+1}{k}$ for some $k$, so $\gcd(n,x)=\gcd(n,y)=1$. Note the $\frac 1x$ and $\frac 1y$ in the grid which need to be there by the condition on the top right and bottom left $2\times 2$. We can solve the top left cell iff $\frac{1x}\cdot \frac 1y + 1 \equiv 0 \pmod n$ or equivalently $n \mid xy+1$, But this is trivial since $n=\frac{xy+1}{k}$. This shows that if we have a solved $3 \times 3$ square missing its top left square then this square is solved, which finsihes along with the claim.

It turns out that the given conditions are mostly useless, as we only need the edge cases. We show that we can construct $f(m, n)$ explicitly by induction on the value of $m+n$. Let $g, h, i$ be the functions associated with $A \sim B$, $B \sim C$, $C \sim D$. The philosophy is to ring around the rosie. There are a few cases we have to consider: $(m, n) = (1, 1)$. Then observe the congruences $b_1a_1 \equiv -1 \pmod {a_0}$, $b_1c_1 \equiv -1 \pmod {a_0}$, $c_1d_1 \equiv -1 \pmod {a_0}$, which show that $a_1d_1 \equiv -1 \pmod {a_0}$. Thus $f(1, 1) = \frac{a_1d_1+1}{a_0}$ is an integer. $m = 1$ and $n > 1$, or the analog. In this case, observe the congruences $d_{n-2} f(1, n-1) \equiv 1 \pmod {d_1}$, $d_n i(1, n-1) \equiv -1 \pmod {d_{n-1}}$, $d_{n-2} i(1, n-1) \equiv 1 \pmod {d_{n-1}}$, so $f(1, n-1) d_n \equiv -1 \pmod {d_{n-1}}$. Thus $f(1, n) = \frac{d_n f(1, n-1) + 1}{d_{n-1}}$ is an integer. $m > 1$ and $n > 1$. Then \begin{align*} f(m, n-1) f(m-1, n-2) &\equiv 1 \pmod {f(m-1, n-1)} \\ f(m-1, n-2) f(m-2, n-1) &\equiv -1 \pmod {f(m-1, n-1)} \\ f(m-1, n) f(m-2, n-1) &\equiv 1 \pmod {f(m-1, n-1)} \end{align*}which imply $f(m, n-1)f(m-1, n) \equiv -1 \pmod {f(m-1, n-1)}$, i.e. $f(m, n)$ is an integer. This exhausts all cases and completes the induction.