From the 5th term on, all terms have a primitive divisor; in other words, a prime divisor that divides no earlier term in the sequence. A proof appears in the Everest-McLaren-Ward paper. - Graham Everest (g.everest(AT)uea.ac.uk), Oct 26 2005
Twelve prime terms are known, occurring at indices 4, 5, 6, 7, 8, 11, 13, 16, 43, 52, 206, 647. The last two have been checked for probable primality only. The 647th term has 18498 decimal digits. Possibly these are the only prime terms in the entire sequence. - Graham Everest (g.everest(AT)uea.ac.uk), Nov 28 2006
The density of primes dividing some term in the sequence is 11/21. - Jeremy Rouse, Sep 18 2013
a(n) is a divisor of a(n+k*(2*n-3)) for all integers n and k. - Peter H van der Kamp, May 18 2015
a(n) is a divisor of A051138(k*(2*n-3)) for all integers n and k. - Helmut Ruhland, Jan 26 2024
R. Jones and J. Rouse, Galois Theory of Iterated Endomorphisms, arXiv:0706.2384 [math.NT], 2007-2009; Proceedings of the London Mathematical Society, 100, no. 3 (2010), 763-794.
a(n) = a(3-n) = (-1)^n * A006769(2*n-3) for all n in Z.
a(n+1)/a(n) seems to be asymptotic to C^n with C = 1.226.... - Benoit Cloitre, Aug 07 2002. Confirmed by Hone - see below.
The terms of the sequence have the leading order asymptotics log a(n) ~ D n^2 with D = zeta(w1)*k^2/(2*w1) - log|sigma(k)| = 0.10222281... where zeta and sigma are the Weierstrass functions with invariants g2 = 4, g3 = -1, w1 = 1.496729323 is the real half-period of the corresponding elliptic curve, k = -1.134273216 as above. This agrees with Benoit Cloitre's numerical result with C = exp(2D) = 1.2268447... - Andrew Hone, Feb 09 2005
a(n) = (a(n-1)*a(n-3) + a(n-2)^2)/a(n-4); a(0) = a(1) = a(2) = a(3) = 1; exact formula is a(n) = A*B^n*sigma (z_0+nk)/(sigma (k))^(n^2), where sigma is the Weierstrass sigma function associated to the elliptic curve y^2 = 4*x^3-4*x+1, A = 1/sigma(z_0) = 0.112724016 - 0.824911687*i, B = sigma(k)*sigma (z_0)/sigma (z_0+k) = 0.215971963 + 0.616028193*i, k = 1.859185431, z_0 = 0.204680500 + 1.225694691*i, sigma(k) = 1.555836426, all to 9 decimal places. This is a special case of a general formula for 4th-order bilinear recurrences. The Somos-4 sequence corresponds to the sequence of points (2n-3)P on the curve, where P = (0, 1). - Andrew Hone, Oct 12 2005
a(2*n) = b(-n), a(2*n+1) = b(n-1) where b(n) = A188313(n) for all n in Z. - Michael Somos, Feb 27 2022
MAPLE
Digits:=11; f(x):=4*x^3-4*x+1; sols:=evalf(solve(f(x), x)); e1:=Re(sols[1]); e3:=Re(sols[2]); w1:=evalf(Int((f(x))^(-0.5), x=e1..infinity)); w3:=I*evalf(Int((-f(x))^(-0.5), x=-infinity..e3)); k:=2*w1-evalf(Int((f(x))^(-0.5), x=1..infinity)); z0:=w3+evalf(Int((f(x))^(-0.5), x=e3..-1)); A:=1/WeierstrassSigma(z0, 4.0, -1.0); B:=WeierstrassSigma(k, 4.0, -1.0)/WeierstrassSigma(z0+k, 4.0, -1.0)/A; for n from 0 to 10 do a[n]:=A*B^n*WeierstrassSigma(z0+n*k, 4.0, -1.0)/(WeierstrassSigma(k, 4.0, -1.0))^(n^2) od; # Andrew Hone, Oct 12 2005