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”).

A214551
Reed Kelly's sequence: a(n) = (a(n-1) + a(n-3))/gcd(a(n-1), a(n-3)) with a(0) = a(1) = a(2) = 1.
32
1, 1, 1, 2, 3, 4, 3, 2, 3, 2, 2, 5, 7, 9, 14, 3, 4, 9, 4, 2, 11, 15, 17, 28, 43, 60, 22, 65, 25, 47, 112, 137, 184, 37, 174, 179, 216, 65, 244, 115, 36, 70, 37, 73, 143, 180, 253, 36, 6, 259, 295, 301, 80, 75, 376, 57, 44, 105, 54, 49, 22, 38, 87, 109, 147
OFFSET
0,4
COMMENTS
Like Narayana's Cows sequence A000930, except that the sums are divided by the greatest common divisor (gcd) of the prior terms.
It is a strong conjecture that 8 and 10 are missing from this sequence, but it would be nice to have a proof! See A214321 for the conjectured values. [I have often referred to this as "Reed Kelly's sequence" in talks.] - N. J. A. Sloane, Feb 18 2017
LINKS
T. D. Noe and N. J. A. Sloane, Table of n, a(n) for n = 0..10000
N. J. A. Sloane, Exciting Number Sequences (video of talk), Mar 05 2021.
FORMULA
It appears that, very roughly, a(n) ~ constant*exp(0.123...*n). - N. J. A. Sloane, Sep 07 2012. See next comment for more precise estimate.
If a(n)^(1/n) converges the limit should be near 1.126 (see link). - Benoit Cloitre, Nov 08 2015
Robert G. Wilson v reports that at around 10^7 terms a(n)^(1/n) is about exp(1/8.4). - N. J. A. Sloane, May 05 2021
EXAMPLE
a(14)=9, a(16)=3, therefore a(17)=(9+3)/gcd(9,3) = 12/3 = 4.
a(24)=28, a(26)=60, therefore a(27)=(28+60)/gcd(28,60) = 88/4 = 22.
MAPLE
a:= proc(n) a(n):= `if`(n<3, 1, (a(n-1)+a(n-3))/igcd(a(n-1), a(n-3))) end:
seq(a(n), n=0..100); # Alois P. Heinz, Oct 18 2012
MATHEMATICA
t = {1, 1, 1}; Do[AppendTo[t, (t[[-1]] + t[[-3]])/GCD[t[[-1]], t[[-3]]]], {100}]
f[l_List] := Append[l, (l[[-1]] + l[[-3]])/GCD[l[[-1]], l[[-3]]]]; Nest[f, {1, 1, 1}, 62] (* Robert G. Wilson v, Jul 23 2012 *)
RecurrenceTable[{a[0]==a[1]==a[2]==1, a[n]==(a[n-1]+a[n-3])/GCD[ a[n-1], a[n-3]]}, a, {n, 70}] (* Harvey P. Dale, May 06 2014 *)
PROG
(Perl)
use bignum;
my @seq = (1, 1, 1);
print "1 1\n2 1\n3 1\n";
for ( my $i = 3; $i < 400; $i++ )
{
my $next = ( $seq[$i-1] + $seq[$i-3] ) /
gcd( $seq[$i-1], $seq[$i-3] );
my $ind = $i+1;
print "$ind $next\n";
push( @seq, $next );
}
sub gcd {
my ($x, $y) = @_;
($x, $y) = ($y, $x % $y) while $y;
return $x;
}
(Haskell)
a214551 n = a214551_list !! n
a214551_list = 1 : 1 : 1 : zipWith f a214551_list (drop 2 a214551_list)
where f u v = (u + v) `div` gcd u v
-- Reinhard Zumkeller, Jul 23 2012
(Sage)
def A214551Rec():
x, y, z = 1, 1, 1
yield x
while True:
x, y, z = y, z, (z + x)//gcd(z, x)
yield x
A214551 = A214551Rec();
print([next(A214551) for _ in range(65)]) # Peter Luschny, Oct 18 2012
(PARI) first(n)=my(v=vector(n+1)); for(i=1, min(n, 3), v[i]=1); for(i=4, #v, v[i]=(v[i-1]+v[i-3])/gcd(v[n-1], v[i-3])); v \\ Charles R Greathouse IV, Jun 21 2017
(Python)
from math import gcd
def aupton(nn):
alst = [1, 1, 1]
for n in range(3, nn+1):
alst.append((alst[n-1] + alst[n-3])//gcd(alst[n-1], alst[n-3]))
return alst
print(aupton(64)) # Michael S. Branicky, Mar 28 2022
CROSSREFS
Similar to A000930. Cf. A341312, A341313, which are also similar.
Starting with a(2) = 3 gives A214626. - Reinhard Zumkeller, Jul 23 2012
Sequence in context: A122453 A017849 A134536 * A371769 A343435 A211010
KEYWORD
nonn,nice
AUTHOR
Reed Kelly, Jul 20 2012
STATUS
approved