This is the special case $(p,q,r)=(1,2,3)$ of the $3$-term Gale-Robinson recurrence:
$$x_{n+p+q+r} x_n = x_{n+p} x_{n+q+r} + x_{n+q} x_{n+p+r} + x_{n+p+q} x_{n+r}$$

Fomin and Zelevinsky proved that, treating the initial values as formal variables, all the $x_n$ are Laurent polynomials with integer coefficients in $(x_0, x_1, \cdots x_{p+q+r-1})$. Therefore, any prime which occurs in the denomintator of some $x_n$ must divide one of the initial values. In your case, that shows that the only prime which can occur in the denominator is $2$.

The sequence seems to have a lot of periodicity modulo powers of $2$, but I haven't yet found a specific statement which I can inductively prove to make sure the denominator never has a $2$ in it.

Here is a proof subject to an unproved relation. According to my laptop, for $n \leq 100$, we have
$$a x_n x_{n+48} + b x_{n+6} x_{n+42} + c x_{n+12} x_{n+36} + d x_{n+18} x_{n+30} + e x_{n+24}^2 =0$$
where $(a,b,c,d,e)$ is

```
(42872600952532756413944577, -7642585197866180463969286501605177115683023, -14777777125160439954108773045163128226074672889387272080, 148964391693661992923680078954077756067110081304751719212599407, 50595833510742832041116346653564092895564724512883353187577334991)
```

The only thing you need to see here is that $(a,b,c,d,e) \equiv (1,1,0,1,1) \mod 2$.

This shows inductively, if $x_{6n}$ is odd for $0 \leq n \leq 7$ (and it is), then $x_{6n}$ is odd for all $n$. Similarly, $x_{6n+1}$, $x_{6n+2}$ and $x_{6n+3}$ are always odd.

The remaining two residue classes are harder, but I found them by doubling my spacing:
$$f x_n x_{n+96} + g x_{n+12} x_{n+84} + h x_{n+24} x_{n+72} + i x_{n+36} x_{n+60} + j x_{n+48}^2 =0$$
with $(f,g,h,i,j)$ equal to

```
(275482421676870371359371463998680435538426963076376380974139366712401528564856518381424556318432563385462892296746457695, -1204187018838917569168734117714049614241185115821559510667269181326087284012251290418681275712354771003853097159487786429076362441774104293746291761783054990301019370687256368668064571321856, 46237948454612900518472923159014977519269442452459902128288590632859589486362269029696379416299719982558060325802294250784058152579910008058460895113571344637868831455176995514194104478127008847575580360960143927556981822085371274654718852, -25258589788169445344223776170236779359054408212184860769933304831364012583819797167161677185288240579806555892962118768351982581424384932133901083260133870630044033809364594884555608989561773335221097351699363330231907242444968619795482350223356786592137691094080623104, 102927547672248207711100989742092219264928069372556751802909240396120983835548385841382536622225530269264900739652796985313302402437855010196494257440627546507583982533406361509714125293160643210579263127894875120075097659020181644866881731502319992829790755348150678034508425757)
```

(Checked for $n \leq 200$.) Again, all that matters is that $(f,g,h,i,j)$ are $(1,0,0,0,1) \bmod 2$. This lets us show inductively that, if the first $7$ values of $x_{12n+4}$ are always $2 \bmod 4$, then they all are. Similarly, $x_{12n+5}$ is always $2 \bmod 4$, $x_{12n+10}$ and $x_{12n+11}$ are always $4 \bmod 8$.

There is obviously a lot of mystery going on with this bilinear relations. I know that they have something to do with $\theta$ functions for abelian surfaces, but the details don't really seem to be recorded anywhere.

3more comments