login

Year-end appeal: Please make a donation to the OEIS Foundation to support ongoing development and maintenance of the OEIS. We are now in our 61st year, we have over 378,000 sequences, and we’ve reached 11,000 citations (which often say “discovered thanks to the OEIS”).

a(n) = ( a(n-1)*a(n-7) + a(n-4)^2 ) / a(n-8); a(0) = ... = a(7) = 1.
8

%I #53 Sep 11 2024 10:07:39

%S 1,1,1,1,1,1,1,1,2,3,4,5,9,18,34,93,180,348,724,3033,9666,24986,83761,

%T 261033,1023728,3923791,26128126,105734485,381740209,1895904805,

%U 14058722881,97964968321,517832518189,4364261070929,25225712161101,181840424632390

%N a(n) = ( a(n-1)*a(n-7) + a(n-4)^2 ) / a(n-8); a(0) = ... = a(7) = 1.

%C From _Vladimir Shevelev_, Apr 04 2016: (Start)

%C For k >= 0, an infinite sequence {a(k,n)} of Somos's sequences (n>=0) is:

%C a(k,0) = a(k,1)= ... = a(k,2*k+1) = 1;

%C and then for n >= 2*k+2,

%C a(k,n) = (a(k,n-1)*a(k,n-2*k-1) + a(k,n-k-1)^2)/a(k,n-2*k-2).

%C In particular, {a(0,n)}=A006125, {a(1,n)}=A006720, {a(2,n)}=A102276, {a(3,n)}=A018896.

%C One can prove that the sequence {a(k,n)} has the first 4k+2 simple differences: 2k+1 zeros, after that k+1 1's and after that k consecutive squares, beginning with 2^2.

%C Further we have nontrivial differences. The first of them for k=0,1,2,... are 6, 16, 33, 59, 96, 146, 211, 293, 394, 516, ... that is, {k^3/3 + 5*k^2/2 + 43*k/6 + 6}.

%C (End)

%H T. D. Noe, <a href="/A018896/b018896.txt">Table of n, a(n) for n = 0..100</a>

%H Mohamed Bensaid, <a href="https://arxiv.org/abs/2409.05911">Sato tau functions and construction of Somos sequence</a>, arXiv:2409.05911 [math.NT], 2024. See p. 7.

%H David Gale, <a href="http://dx.doi.org/10.1007/BF03024306">Mathematical Entertainments</a>, Mathematical Intelligencer, volume 18, number 3, Summer 1996, page 25.

%H Eric Weisstein's World of Mathematics, <a href="http://mathworld.wolfram.com/SomosSequence.html">Somos Sequence</a>

%H <a href="/index/Tu#2wis">Index entries for two-way infinite sequences</a>

%p f:= proc(n) option remember;

%p if n <= 7 then 1 else

%p (procname(n-1)*procname(n-7)+procname(n-4)^2)/procname(n-8)

%p fi

%p end proc:

%p seq(f(n),n=0..50); # _Robert Israel_, Apr 04 2016

%t RecurrenceTable[{a[1]==a[2]==a[3]==a[4]==a[5]==a[6]==a[7]==a[8]==1, a[n]==(a[n-1]a[n-7]+ a[n-4]^2)/a[n-8]},a[n],{n,50}] (* _Harvey P. Dale_, May 02 2011 *)

%t k = 3; Set[#, 1] & /@ Map[a[k, #] &, Range[0, 2 k + 1]]; a[k_, n_] /; n >= 2 k + 2 := (a[k, n - 1] a[k, n - 2 k - 1] + a[k, n - k - 1]^2)/ a[k, n - 2 k - 2]; Table[a[k, n], {n, 0, 35}] (* _Michael De Vlieger_, Apr 04 2016 *)

%o (Haskell)

%o a018896 n = a018896_list !! n

%o a018896_list = replicate 8 1 ++ f 8 where

%o f x = ((a018896 (x - 1) * a018896 (x - 7) + a018896 (x - 4) ^ 2)

%o `div` a018896 (x - 8)) : f (x + 1)

%o -- _Reinhard Zumkeller_, Oct 01 2012

%o (Magma) [n le 8 select 1 else (Self(n-1)*Self(n-7)+Self(n-4)^2 ) / Self(n-8): n in [1..40]]; // _Vincenzo Librandi_, Dec 08 2016

%Y Cf. A006125, A006720, A102276.

%K nonn,nice

%O 0,9

%A _Michael Somos_

%E More terms from _Harvey P. Dale_, May 02 2011