cp's OEIS Frontend

This is a front-end for the Online Encyclopedia of Integer Sequences, made by Christian Perfect. The idea is to provide OEIS entries in non-ancient HTML, and then to think about how they're presented visually. The source code is on GitHub.

User: Michal Radwanski

Michal Radwanski's wiki page.

Michal Radwanski has authored 1 sequences.

A280617 Number of n-smooth Carmichael numbers.

Original entry on oeis.org

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, 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, 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
Offset: 1

Author

Michal Radwanski, Jan 06 2017

Keywords

Comments

A number is b-smooth whenever its prime factors are all not greater than b.
By Korselt's criterion, every Carmichael number is squarefree, therefore n-smooth numbers are products of subsets of primes below n.

Examples

			For n = 19, the n-smooth Carmichael numbers are
  571 = 3*11*17
  1105 = 5*13*17
  1729 = 7*13*19
There are no other Carmichael numbers that are 19-smooth, therefore a(19)=3.
		

Crossrefs

Programs

  • Maple
    extend:= proc(c,p)
         if andmap(q -> (p-1 mod q <> 0), c) then (c, c union {p}) else c fi
    end proc:
    carm:= proc(c) local n;
      n:= convert(c,`*`);
      map(t -> n-1 mod (t-1), c) = {0}
    end proc:
    A[1]:= 0: A[2]:= 0: Primes:= []:
    for k from 3 to 100 do
      A[k]:= A[k-1];
      if isprime(k) then
        Cands:= {{}};
        for i from 1 to nops(Primes) do
          if (k - 1) mod Primes[i] <> 0 then
            Cands:= map(extend,Cands,Primes[i])
          fi;
        od;
        Cands:= map(`union`,select(c -> nops(c) > 1, Cands),{k});
        A[k]:= A[k] + nops(select(carm,Cands));
        Primes:= [op(Primes),k];
      fi
    od:
    seq(A[i],i=1..100); # Robert Israel, Mar 14 2017
  • Mathematica
    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 *)
  • Python
    from collections import Counter
    from functools import reduce
    from itertools import combinations, chain
    from operator import mul
    # http://www.sympy.org
    import sympy as sp
    limit = 90
    # credit: http://stackoverflow.com/a/16915734
    # as proved, there are no Carmichael numbers with
    # less than 3 prime factors, and thus modification
    def powerset(iterable):
        xs = list(iterable)
        return chain.from_iterable(combinations(xs, n) for n in range(3, len(xs) + 1))
    # all computed numbers will be limit-smooth
    def carmichael(limit):
        for d in powerset(sp.primerange(3, limit)):
            n = reduce(mul, d, 1)
            broke = False
            for p in d:
                if (n - 1) % (p - 1) != 0:
                    broke = True
                    break
            if not broke:
                yield (d[-1])
    # from list of pairs of (n, number_of_integers_with_n_as_greatest_prime_factor)
    # creates the sequence
    def prefix(lst):
        r = []
        s = 0
        rpointer = 0
        lstpointer = 0
        while lstpointer < len(lst):
            while lst[lstpointer][0] > rpointer:
                r.append(s)
                rpointer += 1
            s += lst[lstpointer][1]
            lstpointer += 1
        r.append(s + lst[-1][1])
        return r
    c = Counter(carmichael(limit))
    for i, e in enumerate(prefix(sorted(c.items()))):
        if i > 0:
            print(e, end=", ")
    
  • Python
    from math import prod
    from itertools import combinations
    from sympy import primerange
    def A280617(n):
        plist, c = list(primerange(3,n+1)), 0
        for l in range(2,len(plist)+1):
            for q in combinations(plist,l):
                k = prod(q)-1
                if not any(k%(r-1) for r in q):
                    c+=1
        return c # Chai Wah Wu, Sep 25 2024

Formula

a(n) = Sum_{k=1..primepi(n)} A283715(k). - Ondrej Kutal, Sep 27 2024