Number of coprimes of $n$ divisible by 3

Question

I'm looking for a formula for $C(n)$ := the number of coprimes of $n$ in the range $[1, n]$ divisible by 3, where $n$ is any positive integer. The formula should be quick to compute, preferably at similar speed as computing $\varphi(n)$, the totient of $n$, so it can use the factorization of numbers of similar size as $n$.

I know that if $n$ is divisible by 3, then $C(n)=0$.

I know that if $n$ is not divisible by 3, then $\varphi(n)-C(n)$ equals to the number of coprimes of $n$ in the range $[1, n]$ not divisible by 3, which equals to the number of coprimes of $3n$ in the range $[1, n]$. For this quantity there is a formula at https://math.stackexchange.com/a/983287/5758: $\varphi(3n) / 3$, but it's proven true only if $\varphi(3n)/3$ is an integer. I need a general formula which works for any positive integer $n$.

Answer 1

I've found an algorithm to compute $C(n)$ which is a bit slower than I wanted: it first factorizes $n$ and then it uses the inclusion-exclusion principle on $c$ sets, where $c$ is the number of prime factors of $n$. So the time needed is the factorization time of $n$ plus $O(c\cdot 2^c)$ for computing all the terms in the inclusion-exclusion principle.

Let $A_i$ := the set of integers in $[1,n]$ divisible by 3 and divisible by $p_i$, where $p_i$ is the $i$th prime factor of $n$. In this case $C(n) = n - |\cup A_i|$.

It's easy to compute the size of the intersections of any number of sets in $A_i$, for example if $p_1=3$, $p_2=7$, $p_3=11$, then $|A_1\cap A_2\cap A_3|=\lfloor n / 3 / 7 / 11\rfloor$. Once we have the sizes of all $2^c$ set intersections, we can use the inclusion-exclusion principle to sum them with alternating signs to get $|\cup A_i|$. And from that we get $C(n)$.

Here is a Python function which computes $C(n)$:

def C(n):
  """Returns the count of numbers in xrange(3, 1 + n, 3) coprime to n."""
  ps = []
  for p in factorize(n):
    if not ps or ps[-1] != p:  # Skip duplicate primes.
      ps.append(p)
  c = len(ps)
  result = 0
  for ii in xrange(1 << c):  
    sign = 1  
    product = 3
    for j in xrange(c):
      if ii & (1 << j):
        sign = -sign
        if ps[j] != 3:
          product *= ps[j]
    result += sign * (n // product)
  return result