-
Notifications
You must be signed in to change notification settings - Fork 0
/
michael_omega_palindromes.py
71 lines (59 loc) · 1.88 KB
/
michael_omega_palindromes.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
#!/usr/bin/python
# Python program for OEIS A335645
# Translation of Daniel Suteu's PARI
# by Michael S. Branicky, Feb 06 2023
# Install sympy for pypy3:
# pip --python=/usr/bin/pypy3 install sympy
# Timings with pypy3:
# 10 231645546132 77 2.6759650707244873
# 11 49795711759794 93 5.6582934856414795
# 12 2415957997595142 111 11.3309805393219
# 13 495677121121776594 131 52.19296479225159
# 14 22181673755737618122 153 135.02077960968018
# A335645 Smallest palindrome with exactly n distinct prime factors.
data = [1, 2, 6, 66, 858, 6006, 222222, 20522502, 244868442, 6172882716, 231645546132, 49795711759794, 2415957997595142, 495677121121776594, 22181673755737618122, 5521159517777159511255]
from sympy import primorial, primerange, integer_nthroot
def cond(n):
s = str(n)
return s == s[::-1]
def omega_cond(A, B, n):
A = max(A, primorial(n))
def f(m, p, j):
lst = []
for q in primerange(p, integer_nthroot(B//m, j)[0]+2):
if q == 5 and m%2 == 0:
continue
v = m*q
while v <= B:
if j == 1:
if v >= A and cond(v):
print("Found upper-bound: ", v)
lst.append(v)
elif v*(q+1) <= B:
lst += f(v, q+1, j-1)
v *= q
return lst
return sorted(f(1, 2, n))
def a(n):
if n == 0:
return 1
x = primorial(n)
y = 2*x
while True:
print("Sieving range: ", [x,y]);
v = omega_cond(x, y, n)
if len(v) > 0:
return v[0]
x = y+1
y = 2*x
print(a(13))
# ~ print([a(n) for n in range(0, 12)])
# ~ from time import time
# ~ time0 = time()
# ~ alst = []
# ~ for n in range(20):
# ~ an = a(n)
# ~ alst.append(an)
# ~ print(n, an, len(str(alst))-2, time()-time0)
# ~ print(" ", alst)
# ~ print(" ", data)