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
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
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
Comments