

A214551


Reed Kelly's sequence: a(n) = (a(n1) + a(n3))/gcd(a(n1), a(n3)) with a(0) = a(1) = a(2) = 1.


31



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
(list;
graph;
refs;
listen;
history;
text;
internal format)



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
Benoit Cloitre, Graph of a(n)^(1/n) for n=1 up to 381817
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(n1)+a(n3))/igcd(a(n1), a(n3))) 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[n1]+a[n3])/GCD[ a[n1], a[n3]]}, 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[$i1] + $seq[$i3] ) /
gcd( $seq[$i1], $seq[$i3] );
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();
[next(A214551) for i 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[i1]+v[i3])/gcd(v[n1], v[i3])); v \\ Charles R Greathouse IV, Jun 21 2017


CROSSREFS

Similar to A000930. Cf. A341312, A341313, which are also similar.
Cf. also A214320, A214321, A214322, A214323 (gcd's), A219898 (records), A214324, A214325, A214330, A214331, A214809, A227836, A227837.
Starting with a(2) = 3 gives A214626.  Reinhard Zumkeller, Jul 23 2012
Sequence in context: A122453 A017849 A134536 * A343435 A211010 A131731
Adjacent sequences: A214548 A214549 A214550 * A214552 A214553 A214554


KEYWORD

nonn,nice


AUTHOR

Reed Kelly, Jul 20 2012


STATUS

approved



