S(1) = x , S(2) = S(3) = 1 , S(4) = y , S(n) = x*y*z*(S(n-1)*S(n-3)+S(n-2)^2)/S(n-4) if n > 4 .
The first observation is that this definition does produce polynomials with positive integer coefficients in the three indeterminates x , y , z. The proof is by induction. Another observation is that S(n) can be easily extended to n < 1 so that it still satisfies the non-linear equation. When S(1) = S(2) = S(3) = S(4) = 1, then S(n) becomes the Somos-4 sequence of integers. For more information about these sequences, see article "SomosSequence" in MathWorld. So far, I have not found a way to extend this to other Somos sequences such as the Somos 6 sequence. Finally, this polynomial sequence is universal in the sense that it can be used to express any solution of the general Somos-4 recursion with parameters p1,p2:
a(n) = (p1*a(n-1)*a(n-3) + p2*a(n-2)^2)/a(n-4) .
A brief table of the polynomials is:
S(1) = x
S(2) = 1
S(3) = 1
S(4) = y
S(5) = y^2*z + y*z
S(6) = x*y^3*z^2 + x*y^3*z + x*y^2*z^2
S(7) = x^2*y^5*z^3 + x^2*y^5*z^2 + x^2*y^4*z^3 + x*y^5*z^3 +
2*x*y^4*z^3 + x*y^3*z^3
S(8) = x^3*y^7*z^5 + x^3*y^7*z^4 + 3*x^3*y^6*z^5 + 3*x^3*y^6*z^4 +
x^3*y^6*z^3 + 3*x^3*y^5*z^5 + 2*x^3*y^5*z^4 + x^3*y^4*z^5 +
x^2*y^7*z^5 + 3*x^2*y^6*z^5 + 3*x^2*y^5*z^5 + x^2*y^4*z^5
If the maximum exponent of x in all terms of S(n) is taken, the sequence is
n 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 e(n) 1 0 0 0 0 1 2 3 5 7 9 12 15 18 22 26 30 35 40 45 51 57 63which is given by the following simple recursion:
e(n) = 1 + e(n-1) + e(n-3) - e(n-4) .
The maximum exponent sequence of y is e(n+2), and of z is e(n+1). The sum of the maximum exponents of x, y, and z is (n-2)*(n-3)/2, not coincidently a triangular number. The e(n) sequence is listed as A058937 in Sloane's Online Encyclopedia of Integer Sequences (OEIS). It has a simple rational generating funcion.
The polynomials factor into powers of x, y, z, and an irreducible factor F(n) which has a total degree f(n), the first few terms of which are
n 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 f(n) 0 0 0 0 1 2 4 6 8 11 14 18 22 26 31 36 42 48 54 61 68 76 84
A shifted version of the f(n) sequence is listed as A011858 in the OEIS and also has a simple generating function.
Based on earlier work, around December 1992 I came up with this sequence of polynomials with special properties. An explanation of the way in which this sequence was discovered is a long story, but it started with a study of elliptic theta functions. There is a lot more to the story which I will write up if people are interested in it.
Early on, I observed that with a single extra factor, the S(n) polynomials are made homogenous by adding adding a factor of t to the recursion to get S(n) = x*y*z*(S(n-1)*S(n-3) + t*S(n-2)^2)/S(n-4) . Thus, the new values are S(5) = y^2*z + y*z*t, S(6) = x*y^3*z^2 + x*y^3*z*t + x*y^2*z^2*t, etc.
On Oct 13, 2023 after reading arXiv:1501.02879 written by Guoce Xin, et. al., "Hankel Determinant Solutions to Several Discrete Integrable Systems and the Laurent Property" I was able to generalize their generalization in Theorem 4.4 of S(n) one step further. Let R(0) = W, R(1) = X, R(2) = Y, R(3) = Z with recursion R(n) = U*W*X*Y*Z*(R(n-1)*R(n-3) + t*R(n-2)^2)/R(n-4) . These R(n) polynomials are related to the S(n) as R(n) = S(n)*X^(2-n)*Y^(n-1) after performing the substituion of (x, y, z) with (W*Y/X^2, X*Z/Y^2, U*X^2*Y^2).