Displaying 1-10 of 40 results found.
Position of 2310^n among 11-smooth numbers A051038.
+20
3
1, 283, 3847, 20996, 74228, 203084, 469053, 960396, 1797086, 3135610, 5173909, 8156188, 12377846, 18190320, 26005929, 36302854, 49629820, 66611231, 87951744, 114441450, 146960432, 186483973, 234087084, 290949702, 358361266, 437725888, 530566933, 638532124, 763398291, 907076258
COMMENTS
Also position of 2310^(n+1) in A147572.
MATHEMATICA
Table[
Sum[Floor@ Log[11, 2310^n/(2^i*3^j*5^k*7^m)] + 1,
{i, 0, Log[2, 2310^n]},
{j, 0, Log[3, 2310^n/2^i]},
{k, 0, Log[5, 2310^n/(2^i*3^j)]},
{m, 0, Log[7, 2310^n/(2^i*3^j*5^k)]}],
{n, 0, 8}]
PROG
(Python) # uses imports/function in A372401(Python)
from sympy import integer_log, prevprime
def g(x, m): return sum((x//3**i).bit_length() for i in range(integer_log(x, 3)[0]+1)) if m==3 else sum(g(x//(m**i), prevprime(m))for i in range(integer_log(x, m)[0]+1))
Numbers that cannot be written as a sum of two or fewer 11-smooth numbers ( A051038).
+20
0
479, 958, 1151, 1319, 1437, 1559, 1679, 1916, 2302, 2351, 2395, 2638, 2874, 2999, 3013, 3071, 3118, 3353, 3358, 3453, 3671, 3737, 3769, 3832, 3911, 3957, 4199, 4309, 4311, 4604, 4677, 4702, 4703, 4751, 4790, 4919, 5037, 5057, 5269, 5276, 5389, 5443, 5519, 5597, 5683
COMMENTS
Notice that A045535(4) = a(1) = 479.
MATHEMATICA
f[n_] := Union@Flatten@Table[2^a*3^b*5^c*7^d, {a, 0, Log2[n]}, {b, 0, Log[3, n/2^a]}, {c, 0, Log[5, n/(2^a*3^b)]}, {d, 0, Log[7, n/(2^a*3^b*5^c)]}];
b = Block[{nn = 3000, s}, s = f[nn]; {0, 1}~Join~
Select[Union@Flatten@Outer[Plus, s, s], # <= nn &]];
Complement[Range[3000], b]
3-smooth numbers: numbers of the form 2^i*3^j with i, j >= 0.
+10
334
1, 2, 3, 4, 6, 8, 9, 12, 16, 18, 24, 27, 32, 36, 48, 54, 64, 72, 81, 96, 108, 128, 144, 162, 192, 216, 243, 256, 288, 324, 384, 432, 486, 512, 576, 648, 729, 768, 864, 972, 1024, 1152, 1296, 1458, 1536, 1728, 1944, 2048, 2187, 2304, 2592, 2916, 3072, 3456, 3888
COMMENTS
This sequence is easily confused with A033845, which gives numbers of the form 2^i*3^j with i, j >= 1. Don't simply say "numbers of the form 2^i*3^j", but specify which sequence you mean. - N. J. A. Sloane, May 26 2024 These numbers were once called "harmonic numbers", see Lenstra links. - N. J. A. Sloane, Jul 03 2015 Successive numbers k such that phi(6k) = 2k. - Artur Jasinski, Nov 05 2008 Also numbers that are divisible by neither 6k - 1 nor 6k + 1, for all k > 0. - Robert G. Wilson v, Oct 26 2010 Also numbers m such that the rooted tree with Matula-Goebel number m has m antichains. The Matula-Goebel number of a rooted tree can be defined in the following recursive manner: to the one-vertex tree there corresponds the number 1; to a tree T with root degree 1 there corresponds the t-th prime number, where t is the Matula-Goebel number of the tree obtained from T by deleting the edge emanating from the root; to a tree T with root degree m>=2 there corresponds the product of the Matula-Goebel numbers of the m branches of T. The vertices of a rooted tree can be regarded as a partially ordered set, where u<=v holds for two vertices u and v if and only if u lies on the unique path between v and the root. An antichain is a nonempty set of mutually incomparable vertices. Example: m=4 is in the sequence because the corresponding rooted tree is \/=ARB (R is the root) having 4 antichains (A, R, B, AB). - Emeric Deutsch, Jan 30 2012 The number of terms less than or equal to n is Sum_{i=0..floor(log_2(n))} floor(log_3(n/2^i) + 1), or Sum_{i=0..floor(log_3(n))} floor(log_2(n/3^i) + 1), which requires fewer terms to compute. - Robert G. Wilson v, Aug 17 2012 In the 14th century Levi Ben Gerson proved that the only pairs of terms which differ by 1 are (1,2), (2,3), (3,4), and (8,9); see A235365, A235366, A236210. - Jonathan Sondow, Jan 20 2014 The sum of the reciprocals of the 3-smooth numbers is equal to 3. Brief proof: 1 + 1/2 + 1/3 + 1/4 + 1/6 + 1/8 + 1/9 + ... = (Sum_{k>=0} 1/2^k) * (Sum_{m>=0} 1/3^m) = (1/(1-1/2)) * (1/(1-1/3)) = (2/(2-1)) * (3/(3-1)) = 3. - Bernard Schott, Feb 19 2019 Also those integers k for which, for every prime p > 3, p^(2k) - 1 == 0 (mod 24k). - Federico Provvedi, May 23 2022 For n>1, the exponents’ parity {parity(i), parity(j)} of one out of four consecutive terms is {odd, odd}. Therefore, for n>1, at least one out of every four consecutive terms is a Zumkeller number ( A083207). If for the term whose parity is {even, odd}, even also means nonzero, then this term is also a Zumkeller number (as is the case with the last of the four consecutive terms 1296, 1458, 1536, 1728). - Ivan N. Ianakiev, Jul 10 2022 Except the initial terms 2, 3, 4, 8, 9 and 16, these are numbers k such that k^6 divides 6^k. Except the initial terms 2, 3, 4, 6, 8, 9, 16, 18 and 27, these are numbers k such that k^12 divides 12^k. - Mohammed Yaseen, Jul 21 2022
REFERENCES
J.-M. De Koninck & A. Mercier, 1001 Problèmes en Théorie Classique des Nombres, Problème 654 pp. 85, 287-8, Ellipses Paris 2004.
S. Ramanujan, Collected Papers, Ed. G. H. Hardy et al., Cambridge 1927; Chelsea, NY, 1962, p. xxiv.
R. Tijdeman, Some applications of Diophantine approximation, pp. 261-284 of Surveys in Number Theory (Urbana, May 21, 2000), ed. M. A. Bennett et al., Peters, 2003.
FORMULA
An asymptotic formula for a(n) is roughly a(n) ~ 1/sqrt(6)*exp(sqrt(2*log(2)*log(3)*n)). - Benoit Cloitre, Nov 20 2001 The characteristic function of this sequence is given by Sum_{n >= 1} x^a(n) = Sum_{n >= 1} moebius(6*n)*x^n/(1 - x^n). - Paul D. Hanna, Sep 18 2011
MAPLE
A003586 := proc(n) option remember; if n = 1 then 1; else for a from procname(n-1)+1 do numtheory[factorset](a) minus {2, 3} ; if % = {} then return a; end if; end do: end if; end proc: # R. J. Mathar, Feb 28 2011
with(numtheory): for i from 1 to 23328 do if(i/phi(i)=3)then print(i/6) fi od; # Gary Detlefs, Jun 28 2011
MATHEMATICA
a[1] = 1; j = 1; k = 1; n = 100; For[k = 2, k <= n, k++, If[2*a[k - j] < 3^j, a[k] = 2*a[k - j], {a[k] = 3^j, j++}]]; Table[a[i], {i, 1, n}] (* Hai He (hai(AT)mathteach.net) and Gilbert Traub, Dec 28 2004 *)
aa = {}; Do[If[EulerPhi[6 n] == 2 n, AppendTo[aa, n]], {n, 1, 1000}]; aa (* Artur Jasinski, Nov 05 2008 *) fQ[n_] := Union[ MemberQ[{1, 5}, # ] & /@ Union@ Mod[ Rest@ Divisors@ n, 6]] == {False}; fQ[1] = True; Select[ Range@ 4000, fQ] (* Robert G. Wilson v, Oct 26 2010 *) powerOfTwo = 12; Select[Nest[Union@Join[#, 2*#, 3*#] &, {1}, powerOfTwo-1], # < 2^powerOfTwo &] (* Robert G. Wilson v and T. D. Noe, Mar 03 2011 *) fQ[n_] := n == 3 EulerPhi@ n; Select[6 Range@ 4000, fQ]/6 (* Robert G. Wilson v, Jul 08 2011 *) mx = 4000; Sort@ Flatten@ Table[2^i*3^j, {i, 0, Log[2, mx]}, {j, 0, Log[3, mx/2^i]}] (* Robert G. Wilson v, Aug 17 2012 *) f[n_] := Block[{p2, p3 = 3^Range[0, Floor@ Log[3, n] + 1]}, p2 = 2^Floor[Log[2, n/p3] + 1]; Min[ Select[ p2*p3, IntegerQ]]]; NestList[f, 1, 54] (* Robert G. Wilson v, Aug 22 2012 *) Select[Range@4000, Last@Map[First, FactorInteger@#] <= 3 &] (* Vincenzo Librandi, Aug 25 2016 *) Select[Range[4000], Max[FactorInteger[#][[All, 1]]]<4&] (* Harvey P. Dale, Jan 11 2017 *)
PROG
(PARI) test(n)=for(p=2, 3, while(n%p==0, n/=p)); n==1;
for(n=1, 4000, if(test(n), print1(n", ")))
(PARI) list(lim)=my(v=List(), N); for(n=0, log(lim\1+.5)\log(3), N=3^n; while(N<=lim, listput(v, N); N<<=1)); vecsort(Vec(v)) \\ Charles R Greathouse IV, Jun 28 2011 (PARI) list(lim)=my(v=List(), N); for(n=0, logint(lim\=1, 3), N=3^n; while(N<=lim, listput(v, N); N<<=1)); Set(v) \\ Charles R Greathouse IV, Jan 10 2018 (Haskell)
import Data.Set (Set, singleton, insert, deleteFindMin)
smooth :: Set Integer -> [Integer]
smooth s = x : smooth (insert (3*x) $ insert (2*x) s')
where (x, s') = deleteFindMin s
a003586_list = smooth (singleton 1)
a003586 n = a003586_list !! (n-1)
(Sage)
def isA003586(n) :
return not any(d != 2 and d != 3 for d in prime_divisors(n))
@CachedFunction
if n == 1 : return 1
while not isA003586(k) : k += 1
return k
(Python)
from itertools import count, takewhile
def aupto(lim):
pows2 = list(takewhile(lambda x: x<lim, (2**i for i in count(0))))
pows3 = list(takewhile(lambda x: x<lim, (3**i for i in count(0))))
return sorted(c*d for c in pows2 for d in pows3 if c*d <= lim)
(Python)
from sympy import integer_log
def bisection(f, kmin=0, kmax=1):
while f(kmax) > kmax: kmax <<= 1
while kmax-kmin > 1:
kmid = kmax+kmin>>1
if f(kmid) <= kmid:
kmax = kmid
else:
kmin = kmid
return kmax
def f(x): return n+x-sum((x//3**i).bit_length() for i in range(integer_log(x, 3)[0]+1))
(Python) # faster for initial segment of sequence
import heapq
from itertools import islice
def A003586gen(): # generator of terms
v, oldv, h, psmooth_primes, = 1, 0, [1], [2, 3]
while True:
v = heapq.heappop(h)
if v != oldv:
yield v
oldv = v
for p in psmooth_primes:
heapq.heappush(h, v*p)
(Magma) [n: n in [1..4000] | PrimeDivisors(n) subset [2, 3]]; // Bruno Berselli, Sep 24 2012
CROSSREFS
Cf. A051037, A002473, A051038, A080197, A080681, A080682, A117221, A105420, A062051, A117222, A117220, A090184, A131096, A131097, A186711, A186712, A186771, A088468, A061987, A080683 (p-smooth numbers with other values of p), A025613 (a subsequence).
EXTENSIONS
Deleted claim that this sequence is union of 2^n ( A000079) and 3^n ( A000244) sequences -- this does not include the terms which are not pure powers. - Walter Roscello (wroscello(AT)comcast.net), Nov 16 2008
7-smooth numbers: positive numbers whose prime divisors are all <= 7.
(Formerly M0477 N0177)
+10
159
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 12, 14, 15, 16, 18, 20, 21, 24, 25, 27, 28, 30, 32, 35, 36, 40, 42, 45, 48, 49, 50, 54, 56, 60, 63, 64, 70, 72, 75, 80, 81, 84, 90, 96, 98, 100, 105, 108, 112, 120, 125, 126, 128, 135, 140, 144, 147, 150, 160, 162, 168, 175, 180, 189, 192
COMMENTS
Also called humble numbers; sometimes also called highly composite numbers, but this usually refers to A002182. Successive numbers k such that phi(210k) = 48k. - Artur Jasinski, Nov 05 2008 Numbers which are products of single-digit numbers. - N. J. A. Sloane, Jul 02 2017 Phi(a(n)) is 7-smooth. In fact, the Euler Phi function applied to p-smooth numbers, for any prime p, is p-smooth. - Richard Locke Peterson, May 09 2020 Also those integers k, such that, for every prime p > 5, p^(12k) - 1 == 0 (mod 5040k). - Federico Provvedi, Jun 06 2022 The nonprimes with this property are all terms except for 2, 3, 5 and 7, i.e.: (1, 4, 6, 8, 9, 10, 12, 14, 15, 16, 18, 20, 21, 24, 25, 27, 28, 30, 32, 35, 36, 40, 42, 45, ...); the composite terms are all but the first one of this subsequence. ["Trivial" data provided mainly for search purpose.] - M. F. Hasler, Jun 06 2023
REFERENCES
B. C. Berndt, Ramanujan's Notebooks Part IV, Springer-Verlag, see p. 52.
N. J. A. Sloane, A Handbook of Integer Sequences, Academic Press, 1973 (includes this sequence).
N. J. A. Sloane and Simon Plouffe, The Encyclopedia of Integer Sequences, Academic Press, 1995 (includes this sequence).
FORMULA
Sum_{n>=1} 1/a(n) = Product_{primes p <= 7} p/(p-1) = (2*3*5*7)/(1*2*4*6) = 35/8. - Amiram Eldar, Sep 22 2020
MATHEMATICA
Select[Range[250], Max[Transpose[FactorInteger[ # ]][[1]]]<=7&]
aa = {}; Do[If[EulerPhi[210 n] == 48 n, AppendTo[aa, n]], {n, 1, 1200}]; aa (* Artur Jasinski, Nov 05 2008 *) mxExp = 8; Select[Union[Times @@@ Flatten[Table[Tuples[{2, 3, 5, 7}, n], {n, mxExp}], 1]], # <= 2^mxExp &] (* Harvey P. Dale, Aug 13 2012 *) mx = 200; Sort@ Flatten@ Table[ 2^i*3^j*5^k*7^l, {i, 0, Log[2, mx]}, {j, 0, Log[3, mx/2^i]}, {k, 0, Log[5, mx/(2^i*3^j)]}, {l, 0, Log[7, mx/(2^i*3^j*5^k)]}] (* Robert G. Wilson v, Aug 17 2012 *)
PROG
(PARI) test(n)=m=n; forprime(p=2, 7, while(m%p==0, m=m/p)); return(m==1)
for(n=1, 200, if(test(n), print1(n", ")))
(PARI) list(lim)=my(v=List(), t); for(a=0, logint(lim\1, 7), for(b=0, logint(lim\7^a, 5), for(c=0, logint(lim\7^a\5^b, 3), t=3^c*5^b*7^a; while(t<=lim, listput(v, t); t<<=1)))); Set(v) \\ Charles R Greathouse IV, Feb 22 2017 (Haskell)
import Data.Set (singleton, deleteFindMin, fromList, union)
a002473 n = a002473_list !! (n-1)
a002473_list = f $ singleton 1 where
f s = x : f (s' `union` fromList (map (* x) [2, 3, 5, 7]))
where (x, s') = deleteFindMin s
(Magma) [n: n in [1..200] | PrimeDivisors(n) subset PrimesUpTo(7)]; // Bruno Berselli, Sep 24 2012 (Python)
import heapq
from itertools import islice
from sympy import primerange
def A002473gen(p=7): # generate all p-smooth terms
v, oldv, h, psmooth_primes, = 1, 0, [1], list(primerange(1, p+1))
while True:
v = heapq.heappop(h)
if v != oldv:
yield v
oldv = v
for p in psmooth_primes:
heapq.heappush(h, v*p)
(Python)
from sympy import integer_log
def bisection(f, kmin=0, kmax=1):
while f(kmax) > kmax: kmax <<= 1
while kmax-kmin > 1:
kmid = kmax+kmin>>1
if f(kmid) <= kmid:
kmax = kmid
else:
kmin = kmid
return kmax
def f(x):
c = n+x
for i in range(integer_log(x, 7)[0]+1):
i7 = 7**i
m = x//i7
for j in range(integer_log(m, 5)[0]+1):
j5 = 5**j
r = m//j5
for k in range(integer_log(r, 3)[0]+1):
c -= (r//3**k).bit_length()
return c
EXTENSIONS
Additional comments from Michel Lecomte, Jun 09 2007
5-smooth numbers, i.e., numbers whose prime divisors are all <= 5.
+10
114
1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36, 40, 45, 48, 50, 54, 60, 64, 72, 75, 80, 81, 90, 96, 100, 108, 120, 125, 128, 135, 144, 150, 160, 162, 180, 192, 200, 216, 225, 240, 243, 250, 256, 270, 288, 300, 320, 324, 360, 375, 384, 400, 405
COMMENTS
Sometimes called the Hamming sequence, since Hamming asked for an efficient algorithm to generate the list, in ascending order, of all numbers of the form 2^i*3^j*5^k for i,j,k >= 0. The problem was popularized by Edsger Dijkstra.
Numbers k such that 8*k = EulerPhi(30*k). - Artur Jasinski, Nov 05 2008 Also called "harmonic whole numbers", see Howard and Longair, 1982, Table I, page 121. - Hugo Pfoertner, Jul 16 2020 Also called ugly numbers, although it is not clear why. - Gus Wiseman, May 21 2021 Some woody bamboo species have extraordinarily long and stable flowering intervals that belong to this sequence. The model by Veller, Nowak & Davis justifies this observation from the evolutionary point of view. - Andrey Zabolotskiy, Jun 27 2021 Also those integers k for which, for every prime p > 5, p^(4*k) - 1 == 0 (mod 240*k). - Federico Provvedi, May 23 2022 As noted in the comments to A085152, Størmer's theorem implies that the only pairs of consecutive integers that appear as consecutive terms of this sequence are (1,2), (2,3), (3,4), (4,5), (5,6), (8,9), (9,10), (15,16), (24,25), and (80,81). These all represent significant musical intervals. - Hal M. Switkay, Dec 05 2022
FORMULA
Let s(n) = Card(k | a(k)<n) and f(n) = log(n*sqrt(30))^3/(6*log(2)*log(3)*log(5)). Then s(n) = f(n) + O(log(n)). Conjecture: s(n)=f(n) + O(log log n). For example, s(10000000) = 768 is well approximated by f(10000000) = 769.3... (see graphic given as link). - Benoit Cloitre, Dec 30 2001 The characteristic function of this sequence is given by:
Sum_{n>=1} x^a(n) = Sum_{n>=1} -Möbius(30*n)*x^n/(1-x^n). - Paul D. Hanna, Sep 18 2011 Sum_{n>=1} 1/a(n) = Product_{primes p <= 5} p/(p-1) = (2*3*5)/(1*2*4) = 15/4. - Amiram Eldar, Sep 22 2020
EXAMPLE
The sequence of terms together with their prime indices begins:
1: {} 25: {3,3}
2: {1} 27: {2,2,2}
3: {2} 30: {1,2,3}
4: {1,1} 32: {1,1,1,1,1}
5: {3} 36: {1,1,2,2}
6: {1,2} 40: {1,1,1,3}
8: {1,1,1} 45: {2,2,3}
9: {2,2} 48: {1,1,1,1,2}
10: {1,3} 50: {1,3,3}
12: {1,1,2} 54: {1,2,2,2}
15: {2,3} 60: {1,1,2,3}
16: {1,1,1,1} 64: {1,1,1,1,1,1}
18: {1,2,2} 72: {1,1,1,2,2}
20: {1,1,3} 75: {2,3,3}
24: {1,1,1,2} 80: {1,1,1,1,3}
(End)
MAPLE
option remember;
local a;
if n = 1 then
1;
else
for a from procname(n-1)+1 do
numtheory[factorset](a) minus {2, 3, 5 } ;
if % = {} then
return a;
end if;
end do:
end if;
end proc:
MATHEMATICA
mx = 405; Sort@ Flatten@ Table[ 2^a*3^b*5^c, {a, 0, Log[2, mx]}, {b, 0, Log[3, mx/2^a]}, {c, 0, Log[5, mx/(2^a*3^b)]}] (* Or *)
With[{nn=10}, Select[Union[Times@@@Flatten[Table[Tuples[{2, 3, 5}, n], {n, 0, nn}], 1]], #<=2^nn&]] (* Harvey P. Dale, Feb 28 2022 *)
PROG
(PARI) test(n)= {m=n; forprime(p=2, 5, while(m%p==0, m=m/p)); return(m==1)}
for(n=1, 500, if(test(n), print1(n", ")))
(PARI) a(n)=local(m); if(n<1, 0, n=a(n-1); until(if(m=n, forprime(p=2, 5, while(m%p==0, m/=p)); m==1), n++); n)
(PARI) list(lim)=my(v=List(), s, t); for(i=0, logint(lim\=1, 5), t=5^i; for(j=0, logint(lim\t, 3), s=t*3^j; while(s<=lim, listput(v, s); s<<=1))); Set(v) \\ Charles R Greathouse IV, Sep 21 2011; updated Sep 19 2016 (PARI) smooth(P:vec, lim)={ my(v=List([1]), nxt=vector(#P, i, 1), indx, t);
while(1, t=vecmin(vector(#P, i, v[nxt[i]]*P[i]), &indx);
if(t>lim, break); if(t>v[#v], listput(v, t)); nxt[indx]++);
Vec(v)
};
(Magma) [n: n in [1..500] | PrimeDivisors(n) subset [2, 3, 5]]; // Bruno Berselli, Sep 24 2012 (Haskell)
import Data.Set (singleton, deleteFindMin, insert)
a051037 n = a051037_list !! (n-1)
a051037_list = f $ singleton 1 where
f s = y : f (insert (5 * y) $ insert (3 * y) $ insert (2 * y) s')
where (y, s') = deleteFindMin s
(Python)
def isok(n):
while n & 1 == 0: n >>= 1
while n % 3 == 0: n //= 3
while n % 5 == 0: n //= 5
(Python)
from sympy import integer_log
def bisection(f, kmin=0, kmax=1):
while f(kmax) > kmax: kmax <<= 1
while kmax-kmin > 1:
kmid = kmax+kmin>>1
if f(kmid) <= kmid:
kmax = kmid
else:
kmin = kmid
return kmax
def f(x):
c = n+x
for i in range(integer_log(x, 5)[0]+1):
for j in range(integer_log(y:=x//5**i, 3)[0]+1):
c -= (y//3**j).bit_length()
return c
(Python) # faster for initial segment of sequence
import heapq
from itertools import islice
def A051037gen(): # generator of terms
v, oldv, h, psmooth_primes, = 1, 0, [1], [2, 3, 5]
while True:
v = heapq.heappop(h)
if v != oldv:
yield v
oldv = v
for p in psmooth_primes:
heapq.heappush(h, v*p)
CROSSREFS
The partitions with these Heinz numbers are counted by A001399. Requiring the sum of prime indices to be even gives A344297.
13-smooth numbers: numbers whose prime divisors are all <= 13.
+10
25
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 18, 20, 21, 22, 24, 25, 26, 27, 28, 30, 32, 33, 35, 36, 39, 40, 42, 44, 45, 48, 49, 50, 52, 54, 55, 56, 60, 63, 64, 65, 66, 70, 72, 75, 77, 78, 80, 81, 84, 88, 90, 91, 96, 98, 99, 100, 104, 105, 108, 110, 112, 117, 120
COMMENTS
Numbers of the form 2^r*3^s*5^t*7^u*11^v*13^w with r, s, t, u, v, w >= 0.
FORMULA
Sum_{n>=1} 1/a(n) = Product_{primes p <= 13} p/(p-1) = (2*3*5*7*11*13)/(1*2*4*6*10*12) = 1001/192. - Amiram Eldar, Sep 22 2020
EXAMPLE
33 = 3*11 and 39 = 3*13 are terms but 34 = 2*17 is not.
MATHEMATICA
mx = 120; Sort@ Flatten@ Table[ 2^i*3^j*5^k*7^l*11^m*13^n, {i, 0, Log[2, mx]}, {j, 0, Log[3, mx/2^i]}, {k, 0, Log[5, mx/(2^i*3^j)]}, {l, 0, Log[7, mx/(2^i*3^j*5^k)]}, {m, 0, Log[11, mx/(2^i*3^j*5^k*7^l)]}, {n, 0, Log[13, mx/(2^i*3^j*5^k*7^l*11^m)]}] (* Robert G. Wilson v, Aug 17 2012 *)
PROG
(PARI) test(n)=m=n; forprime(p=2, 13, while(m%p==0, m=m/p)); return(m==1)
for(n=1, 200, if(test(n), print1(n", ")))
(PARI) list(lim, p=13)=if(p==2, return(powers(2, logint(lim\1, 2)))); my(v=[], q=precprime(p-1), t=1); for(e=0, logint(lim\=1, p), v=concat(v, list(lim\t, q)*t); t*=p); Set(v) \\ Charles R Greathouse IV, Apr 16 2020 (Magma) [n: n in [1..150] | PrimeDivisors(n) subset PrimesUpTo(13)]; // Bruno Berselli, Sep 24 2012 (Python)
import heapq
from itertools import islice
from sympy import primerange
def agen(p=13): # generate all p-smooth terms
v, oldv, h, psmooth_primes, = 1, 0, [1], list(primerange(1, p+1))
while True:
v = heapq.heappop(h)
if v != oldv:
yield v
oldv = v
for p in psmooth_primes:
heapq.heappush(h, v*p)
(Python)
from sympy import integer_log, prevprime
def bisection(f, kmin=0, kmax=1):
while f(kmax) > kmax: kmax <<= 1
while kmax-kmin > 1:
kmid = kmax+kmin>>1
if f(kmid) <= kmid:
kmax = kmid
else:
kmin = kmid
return kmax
def g(x, m): return sum((x//3**i).bit_length() for i in range(integer_log(x, 3)[0]+1)) if m==3 else sum(g(x//(m**i), prevprime(m))for i in range(integer_log(x, m)[0]+1))
def f(x): return n+x-g(x, 13)
Numbers all of whose prime factors are palindromes.
+10
19
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 14, 15, 16, 18, 20, 21, 22, 24, 25, 27, 28, 30, 32, 33, 35, 36, 40, 42, 44, 45, 48, 49, 50, 54, 55, 56, 60, 63, 64, 66, 70, 72, 75, 77, 80, 81, 84, 88, 90, 96, 98, 99, 100, 101, 105, 108, 110, 112, 120, 121, 125, 126, 128, 131
EXAMPLE
10 = 2 * 5 is a term since both 2 and 5 are palindromes.
110 = 2 * 5 * 11 is a term since 2, 5 and 11 are palindromes.
MAPLE
N:= 5: # to get all terms of up to N digits
digrev:= proc(t) local L; L:= convert(t, base, 10);
add(L[-i-1]*10^i, i=0..nops(L)-1);
end proc:
PPrimes:= [2, 3, 5, 7, 11]:
for d from 3 to N by 2 do
m:= (d-1)/2;
PPrimes:= PPrimes, select(isprime, [seq(seq(n*10^(m+1)+y*10^m+digrev(n), y=0..9), n=10^(m-1)..10^m-1)]);
od:
PPrimes:= map(op, [PPrimes]):
M:= 10^N:
B:= Vector(M);
B[1]:= 1:
for p in PPrimes do
for k from 1 to floor(log[p](M)) do
R:= [$1..floor(M/p^k)];
B[p^k*R] := B[p^k*R] + B[R]
od
od:
# alternative
isA033620:= proc(n)
for d in numtheory[factorset](n) do
if not isA002113(op(1, d)) then
return false;
end if;
end do;
true ;
end proc:
for n from 1 to 300 do
if isA033620(n) then
printf("%d, ", n) ;
end if;
MATHEMATICA
palQ[n_]:=Reverse[x=IntegerDigits[n]]==x; Select[Range[131], And@@palQ/@First/@FactorInteger[#]&] (* Jayanta Basu, Jun 05 2013 *)
PROG
(Haskell)
a033620 n = a033620_list !! (n-1)
a033620_list = filter chi [1..] where
chi n = a136522 spf == 1 && (n' == 1 || chi n') where
n' = n `div` spf
(PARI) ispal(n)=n=digits(n); for(i=1, #n\2, if(n[i]!=n[#n+1-i], return(0))); 1
(Python)
from sympy import isprime, primefactors
def pal(n): s = str(n); return s == s[::-1]
def ok(n): return all(pal(f) for f in primefactors(n))
23-smooth numbers: numbers whose prime divisors are all <= 23.
+10
18
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, 30, 32, 33, 34, 35, 36, 38, 39, 40, 42, 44, 45, 46, 48, 49, 50, 51, 52, 54, 55, 56, 57, 60, 63, 64, 65, 66, 68, 69, 70, 72, 75, 76, 77, 78, 80, 81, 84, 85, 88, 90, 91, 92, 95
FORMULA
Sum_{n>=1} 1/a(n) = Product_{primes p <= 23} p/(p-1) = (2*3*5*7*11*13*17*19*23)/(1*2*4*6*10*12*16*18*22) = 676039/110592. - Amiram Eldar, Sep 22 2020
MAPLE
select(t -> max(numtheory:-factorset(t)) <= 23, [$1..1000]); # Robert Israel, Jan 22 2016
MATHEMATICA
mx = 100; Sort@ Flatten@ Table[ 2^a*3^b*5^c*7^d*11^e*13^f*17^g*19^h*23^i, {a, 0, Log[2, mx]}, {b, 0, Log[3, mx/2^a]}, {c, 0, Log[5, mx/(2^a*3^b)]}, {d, 0, Log[7, mx/(2^a*3^b*5^c)]}, {e, 0, Log[11, mx/(2^a*3^b*5^c*7^d)]}, {f, 0, Log[13, mx/(2^a*3^b*5^c*7^d*11^e)]}, {g, 0, Log[17, mx/(2^a*3^b*5^c*7^d*11^e*13^f)]}, {h, 0, Log[19, mx/(2^a*3^b*5^c*7^d*11^e*13^f*17^g)]}, {i, 0, Log[23, mx/(2^a*3^b*5^c*7^d*11^e*13^f*17^g*19^h)]}] (* Robert G. Wilson v, Jan 19 2016 *)
PROG
(PARI) test(n)=m=n; forprime(p=2, 23, while(m%p==0, m=m/p)); return(m==1)
for(n=1, 100, if(test(n), print1(n", ")))
(PARI) list(lim, p=23)=if(p==2, return(powers(2, logint(lim\1, 2)))); my(v=[], q=precprime(p-1), t=1); for(e=0, logint(lim\=1, p), v=concat(v, list(lim\t, q)*t); t*=p); Set(v) \\ Charles R Greathouse IV, Apr 16 2020 (Magma) [n: n in [1..100] | PrimeDivisors(n) subset PrimesUpTo(23)]; // Bruno Berselli, Sep 24 2012 (Python)
import heapq
from itertools import islice
from sympy import primerange
def agen(p=23): # generate all p-smooth terms
v, oldv, h, psmooth_primes, = 1, 0, [1], list(primerange(1, p+1))
while True:
v = heapq.heappop(h)
if v != oldv:
yield v
oldv = v
for p in psmooth_primes:
heapq.heappush(h, v*p)
(Python)
from sympy import integer_log, prevprime
def bisection(f, kmin=0, kmax=1):
while f(kmax) > kmax: kmax <<= 1
while kmax-kmin > 1:
kmid = kmax+kmin>>1
if f(kmid) <= kmid:
kmax = kmid
else:
kmin = kmid
return kmax
def g(x, m): return sum((x//3**i).bit_length() for i in range(integer_log(x, 3)[0]+1)) if m==3 else sum(g(x//(m**i), prevprime(m))for i in range(integer_log(x, m)[0]+1))
def f(x): return n+x-g(x, 23)
19-smooth numbers: numbers whose prime divisors are all <= 19.
+10
17
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 24, 25, 26, 27, 28, 30, 32, 33, 34, 35, 36, 38, 39, 40, 42, 44, 45, 48, 49, 50, 51, 52, 54, 55, 56, 57, 60, 63, 64, 65, 66, 68, 70, 72, 75, 76, 77, 78, 80, 81, 84, 85, 88, 90, 91, 95, 96, 98, 99, 100
FORMULA
Sum_{n>=1} 1/a(n) = Product_{primes p <= 19} p/(p-1) = (2*3*5*7*11*13*17*19)/(1*2*4*6*10*12*16*18) = 323323/55296. - Amiram Eldar, Sep 22 2020
MATHEMATICA
mx = 120; Sort@ Flatten@ Table[ 2^i*3^j*5^k*7^l*11^m*13^n*17^o*19^p, {i, 0, Log[2, mx]}, {j, 0, Log[3, mx/2^i]}, {k, 0, Log[5, mx/(2^i*3^j)]}, {l, 0, Log[7, mx/(2^i*3^j*5^k)]}, {m, 0, Log[11, mx/(2^i*3^j*5^k*7^l)]}, {n, 0, Log[13, mx/(2^i*3^j*5^k*7^l*11^m)]}, {o, 0, Log[17, mx/(2^i*3^j*5^k*7^l*11^m*13^n)]}, {p, 0, Log[19, mx/(2^i*3^j*5^k*7^l*11^m*13^n*17^o)]}] (* Robert G. Wilson v, Jan 19 2016 *) Select[Range[100], Max[FactorInteger[#][[All, 1]]]<20&] (* Harvey P. Dale, Sep 20 2018 *)
PROG
(PARI) test(n)= {m=n; forprime(p=2, 19, while(m%p==0, m=m/p)); return(m==1)}
for(n=1, 200, if(test(n), print1(n", ")))
(PARI) list(lim, p=19)=if(p==2, return(powers(2, logint(lim\1, 2)))); my(v=[], q=precprime(p-1), t=1); for(e=0, logint(lim\=1, p), v=concat(v, list(lim\t, q)*t); t*=p); Set(v) \\ Charles R Greathouse IV, Apr 16 2020 (Magma) [n: n in [1..100] | PrimeDivisors(n) subset PrimesUpTo(19)]; // Bruno Berselli, Sep 24 2012 (Python)
import heapq
from itertools import islice
from sympy import primerange
def agen(p=19): # generate all p-smooth terms
v, oldv, h, psmooth_primes, = 1, 0, [1], list(primerange(1, p+1))
while True:
v = heapq.heappop(h)
if v != oldv:
yield v
oldv = v
for p in psmooth_primes:
heapq.heappush(h, v*p)
(Python)
from sympy import integer_log
def bisection(f, kmin=0, kmax=1):
while f(kmax) > kmax: kmax <<= 1
while kmax-kmin > 1:
kmid = kmax+kmin>>1
if f(kmid) <= kmid:
kmax = kmid
else:
kmin = kmid
return kmax
def g(x, m): return sum((x//3**i).bit_length() for i in range(integer_log(x, 3)[0]+1)) if m==3 else sum(g(x//(m**i), prevprime(m))for i in range(integer_log(x, m)[0]+1))
def f(x): return n+x-g(x, 19)
17-smooth numbers: numbers whose prime divisors are all <= 17.
+10
16
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 20, 21, 22, 24, 25, 26, 27, 28, 30, 32, 33, 34, 35, 36, 39, 40, 42, 44, 45, 48, 49, 50, 51, 52, 54, 55, 56, 60, 63, 64, 65, 66, 68, 70, 72, 75, 77, 78, 80, 81, 84, 85, 88, 90, 91, 96, 98, 99, 100, 102, 104, 105
FORMULA
Sum_{n>=1} 1/a(n) = Product_{primes p <= 17} p/(p-1) = (2*3*5*7*11*13*17)/(1*2*4*6*10*12*16) = 17017/3072. - Amiram Eldar, Sep 22 2020
MATHEMATICA
mx = 120; Sort@ Flatten@ Table[ 2^i*3^j*5^k*7^l*11^m*13^n*17^o, {i, 0, Log[2, mx]}, {j, 0, Log[3, mx/2^i]}, {k, 0, Log[5, mx/(2^i*3^j)]}, {l, 0, Log[7, mx/(2^i*3^j*5^k)]}, {m, 0, Log[11, mx/(2^i*3^j*5^k*7^l)]}, {n, 0, Log[13, mx/(2^i*3^j*5^k*7^l*11^m)]}, {o, 0, Log[17, mx/(2^i*3^j*5^k*7^l*11^m*13^n)]}] (* Robert G. Wilson v, Aug 17 2012 *)
PROG
(PARI) test(n)= {m=n; forprime(p=2, 17, while(m%p==0, m=m/p)); return(m==1)}
for(n=1, 200, if(test(n), print1(n", ")))
(PARI) list(lim, p=17)=if(p==2, return(powers(2, logint(lim\1, 2)))); my(v=[], q=precprime(p-1), t=1); for(e=0, logint(lim\=1, p), v=concat(v, list(lim\t, q)*t); t*=p); Set(v) \\ Charles R Greathouse IV, Apr 16 2020 (Magma) [n: n in [1..150] | PrimeDivisors(n) subset PrimesUpTo(17)]; // Bruno Berselli, Sep 24 2012 (Python)
import heapq
from itertools import islice
from sympy import primerange
def agen(p=17): # generate all p-smooth terms
v, oldv, h, psmooth_primes, = 1, 0, [1], list(primerange(1, p+1))
while True:
v = heapq.heappop(h)
if v != oldv:
yield v
oldv = v
for p in psmooth_primes:
heapq.heappush(h, v*p)
(Python)
from sympy import integer_log, prevprime
def bisection(f, kmin=0, kmax=1):
while f(kmax) > kmax: kmax <<= 1
while kmax-kmin > 1:
kmid = kmax+kmin>>1
if f(kmid) <= kmid:
kmax = kmid
else:
kmin = kmid
return kmax
def g(x, m): return sum((x//3**i).bit_length() for i in range(integer_log(x, 3)[0]+1)) if m==3 else sum(g(x//(m**i), prevprime(m))for i in range(integer_log(x, m)[0]+1))
def f(x): return n+x-g(x, 17)
Search completed in 0.026 seconds
|