A006509 Cald's sequence: a(n+1) = a(n) - prime(n) if that value is positive and new, otherwise a(n) + prime(n) if new, otherwise 0; start with a(1)=1.
1, 3, 6, 11, 4, 15, 2, 19, 38, 61, 32, 63, 26, 67, 24, 71, 18, 77, 16, 83, 12, 85, 164, 81, 170, 73, 174, 277, 384, 275, 162, 35, 166, 29, 168, 317, 468, 311, 148, 315, 142, 321, 140, 331, 138, 335, 136, 347, 124, 351, 122, 355, 116, 357, 106, 363, 100, 369, 98, 375, 94, 377, 84, 391, 80, 393, 76, 407, 70, 417, 68, 421, 62, 429, 56, 435, 52, 441, 44, 445, 36, 455, 34, 465, 898, 459, 902, 453, 910, 449, 912, 1379, 900, 413, 904, 405, 908, 399, 920, 397, 938, 1485, 928, 365, 934, 1505, 2082, 1495, 2088, 1489, 888, 281, 894, 1511, 892, 261, 0, 643, 1290, 637, 1296, 635, 1308, 631, 1314, 623, 1324, 615, 1334, 607, 1340
Offset: 1
References
- F. Cald, Problem 356, Franciscan order, J. Rec. Math., 7 (No. 4, 1974), 318; 10 (No. 1, 1977-78), 62-64.
- "Cald's Sequence", Popular Computing (Calabasas, CA), Vol. 4 (No. 41, Aug 1976), pp. 16-17.
- N. J. A. Sloane and Simon Plouffe, The Encyclopedia of Integer Sequences, Academic Press, 1995 (includes this sequence).
Links
- T. D. Noe, Table of n, a(n) for n = 1..10000
- Francis Cald and others, Problem 356, Franciscan order, Solution and Guesses, J. Rec. Math., 10 (No. 1, 1977-78), 62-64: Page 62, Page 63, Page 64. [Annotated copy]
- R. G. Wilson, Letter to N. J. A. Sloane with attachment, Jan 1989
Crossrefs
Programs
-
Haskell
a006509 n = a006509_list !! (n-1) a006509_list = 1 : f [1] a000040_list where f xs'@(x:_) (p:ps) | x' > 0 && x' `notElem` xs = x' : f (x':xs) ps | x'' `notElem` xs = x'' : f (x'':xs) ps | otherwise = 0 : f (0:xs) ps where x' = x - p; x'' = x + p -- Reinhard Zumkeller, Oct 17 2011 -
Maple
M1:=500000; a:=array(0..M1); have:=array(0..M1); a[0]:=1; for n from 0 to M1 do have[n]:=0; od: have[0]:=1; have[1]:=1; M2:=2000; nmax:=M2; for n from 1 to M2 do p:=ithprime(n); i:=a[n-1]-p; j:=a[n-1]+p; if i >= 1 and have[i]=0 then a[n]:=i; have[i]:=1; elif j <= M1 and have[j]=0 then a[n]:=j; have[j]:=1; elif j <= M1 then a[n]:=0; else nmax:=n-1; break; fi; od: # To get A006509: [seq(a[n],n=0..M2)]; # To get A112877 (off by 1 because of different offset in A006509): zzz:=[]; for n from 0 to nmax do if a[n]=0 then zzz:=[op(zzz),n]; fi; od: [seq(zzz[i],i=1..nops(zzz))];
-
Mathematica
lst = {1}; f := Block[{b = Last@lst, p = Prime@ Length@lst}, If[b > p && !MemberQ[lst, b - p], AppendTo[lst, b - p], If[ !MemberQ[lst, b + p], AppendTo[lst, b + p], AppendTo[lst, 0]] ]]; Do[f, {n, 60}]; lst (* Robert G. Wilson v, Apr 25 2006 *) -
PARI
A006509_upto(N, U=0)=vector(N,i, N=if(i>1, my(p=prime(i-1)); if( N>p && !bittest(U,N-p), N-p, !bittest(U, N+p), N+p), 1); N && U += 1 << N; N) \\ M. F. Hasler, Mar 06 2024
-
Python
from sympy import primerange, prime def aupton(terms): alst = [1] for n, pn in enumerate(primerange(1, prime(terms)+1), start=1): x, y = alst[-1] - pn, alst[-1] + pn if x > 0 and x not in alst: alst.append(x) elif y > 0 and y not in alst: alst.append(y) else: alst.append(0) return alst print(aupton(130)) # Michael S. Branicky, May 30 2021 -
Python
from sympy import nextprime from itertools import islice def agen(): # generator of terms pn, an, aset = 2, 1, {1} while True: yield an an = m if (m:=an-pn) > 0 and m not in aset else p if (p:=an+pn) not in aset else 0 aset.add(an) pn = nextprime(pn) print(list(islice(agen(), 131))) # Michael S. Branicky, Mar 07 2024
Extensions
More terms from Larry Reeves (larryr(AT)acm.org), Jul 20 2001
Many more terms added by N. J. A. Sloane, Apr 20 2006, to show difference from A117128.
Entry revised by N. J. A. Sloane, Mar 06 2024
Comments