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.

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

Author

Reed Kelly, Jul 20 2012

Keywords

• nonn
• nice

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

Examples

			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.
		

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.

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

Programs

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

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 *)

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

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;
}

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

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

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