cp's OEIS Frontend

This is a front-end for the Online Encyclopedia of Integer Sequences, made by Christian Perfect. The idea is to provide OEIS entries in non-ancient HTML, and then to think about how they're presented visually. The source code is on GitHub.

A006720 Somos-4 sequence: a(0)=a(1)=a(2)=a(3)=1; for n >= 4, a(n) = (a(n-1) * a(n-3) + a(n-2)^2) / a(n-4).

Original entry on oeis.org

1, 1, 1, 1, 2, 3, 7, 23, 59, 314, 1529, 8209, 83313, 620297, 7869898, 126742987, 1687054711, 47301104551, 1123424582771, 32606721084786, 1662315215971057, 61958046554226593, 4257998884448335457, 334806306946199122193, 23385756731869683322514, 3416372868727801226636179
Offset: 0

Views

Author

N. J. A. Sloane

Keywords

Comments

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

References

  • Miklos Bona, editor, Handbook of Enumerative Combinatorics, CRC Press, 2015, page 565.
  • G. Everest, A. van der Poorten, I. Shparlinski and T. Ward, Recurrence Sequences, Amer. Math. Soc., 2003; pp. 9, 179.
  • N. J. A. Sloane and Simon Plouffe, The Encyclopedia of Integer Sequences, Academic Press, 1995 (includes this sequence).

Links

Crossrefs

Cf. A006721, A006722, A006723, A006769, A048736, A028945, A028935, A151502, A165896, A188313, A051138, A006769.
For primes see A129739, A129740, A129741.
Cf. A227199 (primes dividing some term).
Cf. A178384, A247368.

Programs

  • Haskell
    a006720 n = a006720_list !! n
a006720_list = [1,1,1,1] ++
   zipWith div (foldr1 (zipWith (+)) (map b [1..2])) a006720_list
   where b i = zipWith (*) (drop i a006720_list) (drop (4-i) a006720_list)
-- Reinhard Zumkeller, Jan 22 2012
  • Magma
    I:=[1,1,1,1]; [n le 4 select I[n] else (Self(n-1)*Self(n-3)+Self(n-2)^2)/Self(n-4): n in [1..30]]; // Vincenzo Librandi, Aug 07 2017
  • 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
A006720 := proc(n)
    option remember;
    if n <= 3 then
        1;
    else
        (procname(n-1)*procname(n-3)+procname(n-2)^2)/procname(n-4) ;
    end if;
end proc: # R. J. Mathar, Jul 12 2012
  • Mathematica
    a[0] = a[1] = a[2] = a[3] = 1; a[n_] := a[n] = (a[n - 1] a[n - 3] + a[n - 2]^2)/a[n - 4]; Array[a, 23] (* Robert G. Wilson v, Jul 04 2007 *)
RecurrenceTable[{a[0]==a[1]==a[2]==a[3]==1,a[n]==(a[n-1]a[n-3]+a[n-2]^2)/ a[n-4]},a,{n,30}] (* Harvey P. Dale, Apr 07 2018 *)
b[ n_] := If[-2<=n<=2, {2, 1, 1, 3, 23}[[n+3]], 2*a[n+2]^3*a[n+3] + a[n+1]^2*(a[n+3]*a[n+4] - a[n+2]*a[n+5])]; a[ n_] := If[OddQ[n], b[(n-3)/2], b[-n/2]]; (* Michael Somos, Feb 28 2022 *)
  • PARI
    a=vector(99);a[1]=a[2]=a[3]=a[4]=1;for(n=5,#a,a[n]=(a[n-1]*a[n-3]+a[n-2]^2)/a[n-4]); a \\ Charles R Greathouse IV, Jun 16 2011
  • Python
    from gmpy2 import divexact
A006720 = [1, 1, 1, 1]
for n in range(4, 101):
    A006720.append(divexact(A006720[n-1]*A006720[n-3]+A006720[n-2]**2,A006720[n-4]))
# Chai Wah Wu, Sep 01 2014

Formula

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

Started in 1964 by Neil J. A. Sloane | Maintained by The OEIS Foundation Inc.

Content is available under The OEIS End-User License Agreement.