code: def cyclotomic(ctx, n, z):
      n = int(n)
      if n < 0:
          raise ValueError("n cannot be negative")
      p = ctx.one
      if n == 0:
          return p
      if n == 1:
          return z - p
      if n == 2:
          return z + p
      # Use divisor product representation. Unfortunately, this sometimes
      # includes singularities for roots of unity, which we have to cancel out.
      # Matching zeros/poles pairwise, we have (1-z^a)/(1-z^b) ~ a/b + O(z-1).
      a_prod = 1
      b_prod = 1
      num_zeros = 0
      num_poles = 0
      for d in range(1, n + 1):
          if not n % d:
              w = ctx.moebius(n // d)
              # Use powm1 because it is important that we get 0 only
              # if it really is exactly 0
              b = -ctx.powm1(z, d)
              if b:
                  p *= b**w
              else:
                  if w == 1:
                      a_prod *= d
                      num_zeros += 1
                  elif w == -1:
                      b_prod *= d
                      num_poles += 1
      # print n, num_zeros, num_poles
      if num_zeros:
          if num_zeros > num_poles:
              p *= 0
          else:
              p *= a_prod
              p /= b_prod
      return p
  
imports: from ..libmp.backend import xrange
  import math
  import cmath
label_desc: Use image node_id 45 for calling a global function with example usage: cyclotomic(ctx, n, z) and returns: p, p, unknown, unknown