A245770 - OEIS

%I #36 Dec 07 2019 12:18:27

%S 1,9,6,5,93,84626,3,2,502884,69399,5105820,4944,2,816406286,986,

%T 4825342,70,9,480865,2,66,93,55058,25,4081284,174,270,85,555964462,

%U 4895,9644288,7566,344,28475648,8,65,201,45648,23,4543,393,2602,412,724,660,55,74881,20962,25,1715364,892,600

%N Prime sieve of Pi.

%H Robert G. Wilson v, <a href="/A245770/b245770.txt">Table of n, a(n) for n = 1..1029</a> (first 422 terms from Manfred Scheucher)

%H Manfred Scheucher, <a href="/A245770/a245770_1.sage.txt">Sage Script</a> (Note: should also run in pure python with a few modifications)

%e Find the first occurrence of prime 2 in the digits of Pi (only 25 digits in this illustration):

%e 3141592653589793238462643..., and replace it with a space:

%e 314159 653589793238462643...  Repeat the process with 3:

%e 14159 653589793238462643...  Then 5:

%e 141 9 653589793238462643...  Then 7:

%e 141 9 653589 93238462643...  Then 11, 13, 17, etc., until the first occurrence of every prime is eliminated from the digits of Pi.

%e 1   9 6  5   93  84626  ...  Then consolidate gaps between the remaining digits into a single comma:

%e 1,9,6,5,93,84626,...         to produce the first terms in the prime sieve of Pi.

%o (Python)

%o def arccot(x, unity):

%o ....sum = xpower = unity // x

%o ....n = 3

%o ....sign = -1

%o ....while 1:

%o ........xpower = xpower // (x*x)

%o ........term = xpower // n

%o ........if not term:

%o ............break

%o ........sum += sign * term

%o ........sign = -sign

%o ........n += 2

%o ....return sum

%o def pi(digits):

%o ....unity = 10**(digits + 10)

%o ....pi = 4 * (4*arccot(5, unity) - arccot(239, unity))

%o ....return pi // 10**10

%o def primes(n):

%o ....""" Returns  a list of primes < n """

%o ....sieve = [True] * n

%o ....for i in range(3,int(n**0.5)+1,2):

%o ........if sieve[i]:

%o ............sieve[i*i::2*i]=[False]*((n-i*i-1)/(2*i)+1)

%o ....return [2] + [i for i in range(3,n,2) if sieve[i]]

%o a = pi(370)

%o b = primes(100000000)

%o y = str(a)

%o for x in b:

%o ....if str(x) in y:

%o ........y = y.replace(str(x)," ",1)#replace first occurrence only

%o ........

%o while "  " in y:

%o ....y = y.replace("  "," ")#replace long chains of spaces with a single space

%o z = y.split(" ")#split terms into a list

%o z = filter(None, z)#remove null terms

%o f = map(int,z)#convert to integers

%o print(f[0:-1])

%o # _David Consiglio, Jr._, Jan 03 2015

%Y Cf. A000796, A246022, A248804 (Prime sieve of e), A248831 (Prime sieve of sqrt(2)).

%K nonn,base

%O 1,2

%A _Gil Broussard_, Aug 01 2014

%E a(14) corrected by _Manfred Scheucher_, May 25 2015