Inspired by this post & TF comment I tried symbollic regression [0] Basically it uses genetic algorithm to find a formula that matches known input and output vectors with minimal loss I tried to force it to use pi constant but was unable I don't have much expreience with this library but I'm sure with more tweaks you'll get the right result

  from pysr import PySRRegressor

  def f(n):
      if n % 15 == 0:
          return 3
      elif n%5 == 0:
          return 2
      elif n%3 == 0:
          return 1
      return 0

  n = 500
  X = np.array(range(1,n)).reshape(-1,1)
  Y = np.array([f(n) for n in range(1,n)]).reshape(-1,1)
  model = PySRRegressor(
          maxsize=25,
          niterations=200,  # < Increase me for better results
          binary_operators=["+", "*"],
          unary_operators=["cos", "sin", "exp"],
          elementwise_loss="loss(prediction, target) = (prediction - target)^2",

)

  model.fit(X,Y)

Result I got is this:

((cos((x0 + x0) * 1.0471969) * 0.66784626) + ((cos(sin(x0 * 0.628323) * -4.0887628) + 0.06374673) * 1.1508249)) + 1.1086457

with compleixty 22 loss: 0.000015800686 The first term is close to 2/3 * cos(2pi*n/3) which is featured in the actual formula in the article. the constant doesn't compare to 11/15 though

[0] https://github.com/MilesCranmer/PySR