Any solution

A sketch: We have $(b_n - 1)(b_n + 1) = \prod_{i=1}^{n} a_i$ for each $n$. Use pigeonhole to get a long subsequence of $b_i$'s, say $b_{k_i}$ such that the set of prime divisors of $b_{k_i} - 1$ is a subset for $b_{k_j} - 1$ and similarly for $b_{k_j} + 1, b_{k_i} + 1$ for all $i < j$. This will in turn imply that $b_{k_i} - 1 | 2(b_{k_j} - 1)$ and $b_{k_i} + 1 | 2(b_{k_j} + 1)$.Then you will get the relation $y-1|2(x-1), y+1|2(x+1)$ for some $x,y$. If $x = y$, then $a_i = 1$ so let $2x = k(y-1) + 1$ for some $k > 2$. Then $y+1| k(y-1)+4$ which gives us $y+1 | 2k-4$, implying that $k \geq \frac{y}{2}$ and hence $b_{k_{i+1}} > \frac{(b_{k_i} - 1)^2 }{4} > b_{k_i}^{\frac{3}{2}}$ for $b_{k_i} > 20$. However, we can get a pretty long sequence since for each prime $\leq 2016$, it is either not in both $b_i + 1, b_i-1$, or in exactly one of them except where they can share a possible factor of $2$. This gives us $2*3^{\pi(2016)}$ possibilities of the set of prime factors of $b_i-1$ and $b_i + 1$ since if both do not share a possible factor of $2$, we have $3^{\pi(2016)}$ possibilities and if they do, it is still upper bounded by $3^{\pi(2016)}$ where $\pi(n)$ is the prime number counting function. Since $\pi(2016) = 305$, by pigeonhole this gives us a suitable subsequence of length $\frac{2^{2015}}{3^{305}}$. Now the final term of this sequence is at least $20^{\frac{3}{2}^{\frac{2^{2015}}{3^{305}} - 10}}$ since the tenth term of the sequence is at least $2^{10} > 20^2$. A simple calculation will show that this is larger than $2016^{2^{2015}}$ which is an upper bound for $b_n^2$ for all $n$ and hence we must have $a_i = 1$ somewhere.

Official solution

First we derive the following important results (well known in literature):

Now, finally:

Here is my solution: https://calimath.org/pdf/SerbiaMO2016-6.pdf And I uploaded the solution with motivation to: https://youtu.be/FpprIyA0VW0

Assume the contrary and let $a_{i}>1$ for all $i$. The key observation is that if $a_{k}a_{k+1}\dots a_{n}$ is a perfect square, then we can bound some $a_{i}$. Lemma. If $a,b>1$ are positive integers such that $a+1$ and $ab^2+1$ are both perfect squares, then $b^2>4a$. Proof: Let $a+1=x^2$, $ab^2+1=y^2$. Then we have that \[y^2 = ab^2+1 = (x^2-1)b^2+1 < b^2x^2 \Longrightarrow y\leq bx-1\]\[\Longrightarrow ab^2+1\leq (bx-1)^2 = b^2x^2-2bx+1 = (a+1)b^2 - 2b\sqrt{a+1}+1\Longrightarrow b\geq 2\sqrt{a+1}>2\sqrt{a}\] We utilize the Lemma in the following way. Construct the tuplets $b_{k} = (c_{k,p_{1}}, c_{k,p_{2}}, \dots ,c_{k,p_{\pi(2016)}})$ for every $1\leq k\leq 2^{2016}$ where \[c_{k, p_{j}} = \begin{cases} 0, \text{if } \nu_{p_{j}}\left(\prod\limits_{i=1}^{k} a_{i}\right)\equiv 0\pmod 2\\ 1, \text{if } \nu_{p_{j}}\left(\prod\limits_{i=1}^{k} a_{i}\right)\equiv 1\pmod 2\\ \end{cases}\]accounts for the parity of the power of each prime in the factorization of $\prod\limits_{i=1}^{k} a_{i}$. Denote $p=\pi(2016)$ going forwards. Clearly $\pi(2016)<2016\cdot \left(1-\frac{1}{2}\right)\left(1-\frac{1}{3}\right)\left(1-\frac{1}{5}\right)\left(1-\frac{1}{7}\right)<461$ for a rough estimate. Therefore, by the PHP, we have that at least $\frac{2^{2016}}{2^{p}}=2^{2016-p}$ of the tuplets are the same. That means that if $b_{i} = b_{j}$, then $c_{i, p_{k}} = c_{j, p_{k}}$ for all $1\leq k\leq p$ and hence $\prod\limits_{s = i+1}^{j} a_{s}$ is a perfect square, so we can assume that there exist indices $i_{1}<i_{2}<\dots <i_{2^{2016-p}}$ such that $\prod\limits_{s = i_{a}+1}^{i_{b}} a_{s}$ is a perfect square for all $1\leq a<b\leq 2^{2016-p}$. This, however, by the Lemma implies that: \[\frac{\prod\limits_{s=1}^{i_{k+1}} a_{s}}{\prod\limits_{s=1}^{i_{k}} a_{s}}\geq 4\prod\limits_{s=1}^{i_{k}}a_{s}\]If we denote $d_{k} = \prod\limits_{s=1}^{i_{k}}a_{s}$, the above inequality implies that $d_{k+1}\geq 4d_{k}^3>d_{k}^3$, which by induction gives the bound \[d_{2^{2016-p}}>d_{1}^{3^{2^{2016-p}-1}}\geq a_{1}^{3^{2^{2016-p}-1}}\geq 3^{3^{2^{2016-p}-1}}\]\[\Longrightarrow \prod\limits_{i=1}^{2^{2016}} a_{i}\geq d_{2^{2016-p}}\geq 3^{3^{2^{2016-p}-1}}> 3^{3^{2^{2016-460}-1}} > 3^{3^{2^{1555}}}>2016^{2^{2016}}\geq \prod\limits_{i=1}^{2^{2016}} a_{i}\]which is a contradiction and hence our assumption is wrong and some $a_{i}$ must be equal to $1$. The only inequality we haven't proved is that $3^{3^{2^{1555}}}>2016^{2^{2016}}$ which we'll do in a fairly straightforward way: \[3^{3^{2^{1555}}}>2016^{2^{2016}}\iff 3^{2^{1555}}\ln(3)>2^{2016}\ln(2016)\iff 2^{1555}\ln(3)+\ln\ln(3)>2016\ln(2)+\ln\ln(2016)\]But $2^{1555}\ln(3)+\ln\ln(3)-2016\ln(2)-\ln\ln(2016)>(2^{1555}-2016)\ln 3-\ln \ln 2016>0$, so we're done.