|
%I #28 Oct 10 2019 04:07:53
%S 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,2,3,3,3,3,3,3,3,3,3,3,4,4,6,6,6,6,
%T 6,6,8,8,8,8,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,
%U 10,16,16,16,16,16,16,19,19,19,19,20,20,29,29,29,29,29,29,31,31,31,31,31,31,31,31
%N Number of n-smooth Carmichael numbers.
%H Wikipedia, <a href="https://en.wikipedia.org/wiki/Carmichael_number">Carmichael number</a>
%H Wikipedia, <a href="https://en.wikipedia.org/wiki/Carmichael_number#Korselt.27s_criterion">Korselt's criterion</a>
%H Wikipedia, <a href="https://en.wikipedia.org/wiki/Smooth_number">Smooth number</a>
%e A number is b-smooth whenever its prime factors are all not greater than b.
%e By Korselt's criterion, every Carmichael number is squarefree, therefore n-smooth numbers are products of subsets of primes below n.
%e For n = 19, the n-smooth Carmichael numbers are
%e 571 = 3*11*17
%e 1105 = 5*13*17
%e 1729 = 7*13*19
%e There are no other Carmichael numbers that are 19-smooth, therefore a(19)=3.
%p extend:= proc(c,p)
%p if andmap(q -> (p-1 mod q <> 0), c) then (c, c union {p}) else c fi
%p end proc:
%p carm:= proc(c) local n;
%p n:= convert(c,`*`);
%p map(t -> n-1 mod (t-1), c) = {0}
%p end proc:
%p A[1]:= 0: A[2]:= 0: Primes:= []:
%p for k from 3 to 100 do
%p A[k]:= A[k-1];
%p if isprime(k) then
%p Cands:= {{}};
%p for i from 1 to nops(Primes) do
%p if (k - 1) mod Primes[i] <> 0 then
%p Cands:= map(extend,Cands,Primes[i])
%p fi;
%p od;
%p Cands:= map(`union`,select(c -> nops(c) > 1, Cands),{k});
%p A[k]:= A[k] + nops(select(carm,Cands));
%p Primes:= [op(Primes),k];
%p fi
%p od:
%p seq(A[i],i=1..100); # _Robert Israel_, Mar 14 2017
%t a[n_] := a[n] = If[n<17, 0, a[n-1] + If[! PrimeQ[n], 0, Block[{t, k = PrimePi@n, p}, p = Prime@k; Length@ Select[ Subsets[ Prime@ Range[2, k-1], {2, k-2}], (t = Times @@ #; Mod[t-1, p-1] == 0 && And @@ IntegerQ /@ ((p t - 1)/ (#-1))) &]]]]; Array[a, 80] (* _Giovanni Resta_, Mar 14 2017 *)
%o # (Python 3)
%o from collections import Counter
%o from functools import reduce
%o from itertools import combinations, chain
%o from operator import mul
%o # http://www.sympy.org
%o import sympy as sp
%o limit = int(input())
%o # credit: http://stackoverflow.com/a/16915734
%o # as proved, there are no Carmichael numbers with
%o # less than 3 prime factors, and thus modification
%o def powerset(iterable):
%o xs = list(iterable)
%o return chain.from_iterable( combinations(xs,n) for n in range(3, len(xs)+1))
%o # all computed numbers will be limit-smooth
%o def carmichael(limit):
%o for d in powerset(sp.primerange(3, limit)):
%o n = reduce(mul, d, 1)
%o broke = False
%o for p in d:
%o if (n-1)%(p-1) != 0:
%o broke = True
%o break
%o if not broke:
%o yield(d[-1])
%o # from list of pairs of (n, number_of_integers_with_n_as_greatest_prime_factor)
%o # creates the sequence
%o def prefix(lst):
%o r = []
%o s = 0
%o rpointer = 0
%o lstpointer = 0
%o while lstpointer < len(lst):
%o while lst[lstpointer][0] > rpointer:
%o r.append(s)
%o rpointer += 1
%o s += lst[lstpointer][1]
%o lstpointer += 1
%o r.append(s+lst[-1][1])
%o return r
%o c = Counter(carmichael(limit))
%o for i, e in enumerate(prefix(sorted(c.items()))):
%o print(i, e)
%Y Cf. A002997, A081702.
%K nonn
%O 1,17
%A _Michal Radwanski_, Jan 06 2017
|
