Given a prime number \$p\$ output the asymptotic density of the set of positive integers which have \$p\$ as their second-smallest distinct prime factor


Input: one of the following formats:

  • A prime number \$p\$
  • A integer \$n\$ to represent the \$n\$-th (0-index or 1-index) prime number \$p\$
  • No input; return probabilities for all prime numbers

Output: The asymptotic density of the set of positive integers which have \$p\$ as their second-smallest distinct prime factor.

  • This may be a numerator, denominator pair (in either, specified order), a fraction object, a floating point number, etc.
  • Type or machine inaccuracies are acceptable as long as the method is fundamentally sound but if these do adversely affect the test case results then do try to provide an alternative without inaccuracies.
  • This is , so you may not approximate the result by sampling.


This is so the shortest submission (in bytes) for each language is the winner!


  2  ->  0/1            =  0
  3  ->  1/6            =  0.16666666666666
  5  ->  1/10           =  0.1
  7  ->  1/15           =  0.06666666666666
 11  ->  46/1155        =  0.03982683982684
 13  ->  44/1365        =  0.03223443223443
 17  ->  288/12155      =  0.02369395310571
 19  ->  33216/1616615  =  0.02054663602651

Inspired by this video

  • 3
    \$\begingroup\$ This is A342479 (numerators) and A342480 (denominators). \$\endgroup\$
    – Arnauld
    Commented Apr 3 at 8:06
  • \$\begingroup\$ I find the wording confusing, especially what is n and what is p. Do you mean this?: Given \$n\$, output the probability that the \$n\$-th prime is the second smallest distinct prime factor of an integer chosen uniformly on \$\{1,...,n\}\$ \$\endgroup\$
    – Luis Mendo
    Commented Apr 3 at 14:50
  • \$\begingroup\$ @LuisMendo, updated wording; you make take the input as the prime number itself \$p\$ or as \$n\$ to represent the n-th prime; should I perhaps remove the "uniformly random" part, the question asks for the probability over all positive integers \$\endgroup\$ Commented Apr 3 at 15:16
  • 1
    \$\begingroup\$ No, specifying that the randomness is "uniform" is necessary (as argued here). I still find the wording of the first paragraph confusing: uniformly chosen random positive integer on what set? The probability is not over all positive integers; in that case it could not be uniform \$\endgroup\$
    – Luis Mendo
    Commented Apr 3 at 15:32
  • 8
    \$\begingroup\$ Since uniform probability over a countably infinite set is problematic, you might instead say something like "the asymptotic density of the set of positive integers which have $p$ as their second-smallest distinct prime factor". en.wikipedia.org/wiki/Natural_density \$\endgroup\$ Commented Apr 3 at 15:43

12 Answers 12


Vyxal, 8 bytes


Try it Online! -2 thanks to tsh, -1 thanks to lyxal.

 ƛ   ;   # Map over range from 1 to...
¡        # n!
   ḣh    # Is the second element of
  Ǐ      # The prime factors of the number
      =  # Equal to the input
       ṁ # Take the mean of all those 1s/0s to get the probability
  • \$\begingroup\$ I don't speak Vyxal, but I'm wondering that is it possible to replace "Primes up to n" by "1 to n" (or simply \$n!\$) and saving some bytes here? \$\endgroup\$
    – tsh
    Commented Apr 3 at 5:31
  • \$\begingroup\$ @tsh Thanks! (filler) \$\endgroup\$
    – emanresu A
    Commented Apr 3 at 5:48
  • \$\begingroup\$ Try it Online! for 8 bytes by putting the = outside the map \$\endgroup\$
    – lyxal
    Commented Apr 3 at 7:09

Jelly, 11 bytes


Try it online!

Let \$P = p_1, p_2, \cdots, p_k\$ be the list of primes strictly less than the input prime \$p\$. Then, the chance that a number is a multiple of \$p_i\$ but not \$p_1, \cdots, p_{i-1}, p_{i+1}, \cdots, p_k\$ is equal to $$\frac{(p_1-1)(p_2-1) \cdots (p_{i-1}-1) (p_{i+1}-1) \cdots (p_k-1)}{p_1 p_2 \cdots p_k}$$

and since such cases are pairwise disjoint, we can sum them and divide by \$p\$ (meaning that the chosen number is a multiple of \$p\$). The denominators are all equal, so the code computes the sum of the numerators and divides by the product of all primes under \$p\$, including \$p\$.

ÆR            List of primes under p, including p
  µ           Set the list as the argument of the rest
   Ṗ          Remove p from the list
    ’         Decrement each number
     P-Ƥ      Products of outfixes
        S     Sum
         ÷P   Divide by product of all primes

05AB1E, 9 bytes


Port of @emanresuA's Vyxal answer.

Try it online or verify the smaller test cases (the larger ones time out).


!          # Take the factorial of the (implicit) input-integer
 L         # Pop and push a list in the range [1,n!]
  f        # Map each value to an inner list of its unique prime factors
   ζ       # Zip/transpose; swapping rows/columns,
           # using " " as filler for unequal length rows
    1è     # Pop and leave the second row of this matrix (0-based 1st row)
      Q    # Check for each item whether it's equal to the (implicit) input
           # (which will always be falsey for the spaces)
       ÅA  # Pop and take the arithmetic mean of these 0s and 1s
           # (which is output implicitly as result)

Jelly, 11 bytes


Try It Online!

Same idea as the Vyxal answer, which is the observation that we can just look at 1 through N!. Unfortunately, this is much longer in Jelly due to chaining rules and certain built-ins. For example, Jelly has no "unique prime factors" built-in so we have to get the prime factorization and then deduplicate it. Getting the prime factorization itself takes an extra byte as well. Finally, since arithmetic mean is two bytes, it's less efficient (because of chaining rules) to do it the way the Vyxal answer did it at the end.


!ÆfQḊḢƲ€ċ÷!    Main Link
!              Factorial
       €       For each from 1 .. N!
 Æf                Get its prime factorization
   Q               Remove duplicates
    Ḋ              Remove the first element
     Ḣ             Get the first element (0 if the list is empty)
        ċ      Count occurrences of N
         ÷     Divide by
          !    N!

Other Ideas

!ÆfQiɗ€ċ2÷! is also 11 bytes which instead checks for the index of N in each array and counts the number of 2s.

This also lets us do !ÆfQiɗ€=2Æm for 11 bytes which checks for equality with 2 to get a bit array and then takes the mean like the Vyxal answer does. Having =2 instead of just = saves us from needing to un-chain the = and Æm.


K (ngn/k), 25 24 bytes


Try it online!

Uses the formula from Bubbler's Jelly answer.

{               }.1`pri\    y := primes<x
          */\y-1            (y-1;prod(y-1))
      x,y,                  (x;y...;y-1;prod(y-1))
   %/|                      prod(y-1)/(x*prod(y)*(y-1))
 +/                         sum

PARI/GP, 40 bytes


Attempt This Online!

Based on the formula in @Bubbler's Jelly answer. Takes the derivative of the polynomial \$\frac{1}{p}\prod_{q}(1+\frac{x}{q})\$, where \$q\$ ranges over primes less than the input \$p\$, and evaluates it at \$x=-1\$.


05AB1E, 13 bytes


Attempt This Online!

Directly computes the answer. This seems way too long.

Roughly equivalent to

p = input()
a = primes_lt(p)
return sum(product((a-1)/a)/(a-1))/p

JavaScript (ES6), 66 bytes

Takes a prime number and returns the probability as a floating point number.


Try it online!


This is derived from the formula given on OEIS:

$$a(n)=\frac{1}{n}\prod_{\text{prime }q<n}\left(1-\frac{1}{q}\right)\times\sum_{\text{prime }q<n}\frac{1}{q-1}$$


n =>               // n = input
1 / n * (          // divide the final result by n
  h = q =>         // h is a recursive function taking q
  --n > 1 ?        //   decrement n; if it's greater than 1:
    ( g = _ =>     //     g is a recursive function that
      n % d-- ?    //     decrements d until it divides n
        g()        //     and returns ...
      :            //     ... 1/n - 1 if n is prime
        !d / n - 1 //     or -1 if n is composite (*)
    )(d = n - 1) * //     initial call to g with d = n - 1
    h(             //     multiply by the result of
      q - !d / ~-n //     a recursive call with q - 1/(n - 1)
    )              //     if n is prime or q if n is composite
  :                //   else:
    q              //     stop and return q
)``                // initial call to h with q zero'ish

* Because we need to coerce \$q\$ to a number for obscure golfing reasons, we multiply by a negative number at each call to the recursive function h and eventually multiply by \$q<0\$. This still returns a non-negative number as expected because the number of iterations is guaranteed to be odd (except for \$n=2\$, which immediately returns \$0\$).


Charcoal, 28 bytes


Attempt This Online! Link is to verbose version of code. Explanation: Port of @Bubbler's formula.


Input n.


Get the primes below n.


Vectorised divide the product of the decremented primes by the decremented primes, take the sum, and divide by the product of n and the primes. (The sum is taken using base 1 as this works for zero elements.)


Wolfram Language (Mathematica), 51 bytes


Try it online! An unnamed function taking a 1-indexed integer # as input and returning the fraction for the #th prime. Not super golfy, but demonstrates a generating function technique that folks could consider for their languages:

If \$p_k\$ denotes the \$k\$th prime, then the polynomial \$\prod_{k=1}^{n-1} \bigl( 1 + \frac{x-1}{p_k} \bigr)\$ expands into \$c_0 + c_1x + \cdots + c_{n-1}x^{n-1}\$ where \$c_j\$ is the asymptotic density of the set of integers divisible by exactly \$j\$ of the primes \$p_1,\ldots,p_{n-1}\$. (That's actually true for any \$n-1\$ primes, not just the first \$n-1\$.) We extract the coefficient \$c_1\$ by taking the derivative of the polynomial and then setting \$x=0\$, and then divide by \$p_n\$ manually.


Ruby, 75 bytes


Try it online!


Wait for Arnauld to post a Javascript solution, then write the equivalent in Ruby.


PARI/GP, 102 bytes

use the formula in A342479 (numerators) and A342480 (denominators)

$$ f(n) = \left\{ \begin{array}{l} p = \text{prime}(n), \\ q = \{ q \in \mathbb{N} \,|\, q < p \wedge \text{isprime}(q) \}, \\ \sum_{i=1}^{|\{q\}|} \frac{1}{q_i - 1} \cdot \left( \prod_{j=1}^{|\{q\}|} \frac{q_j - 1}{q_j} \right) \frac{1}{p} \end{array} \right. $$

102 bytes, it can be golfe more.

Golfed version. Attempt This Online!


Ungolfed version. Attempt This Online!

f(n) = {
  my(p = prime(n));
  my(q = select(p -> isprime(p), [1..p-1]));
  sum(i=1, #q, 1/(q[i] - 1) * prod(j=1, #q, (q[j] - 1)/q[j]) / p)

print(vector(8, n, f(n)));

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.