code: def dmp_zz_heu_gcd(f, g, u, K):
      """
      Heuristic polynomial GCD in `Z[X]`.
  
      Given univariate polynomials `f` and `g` in `Z[X]`, returns
      their GCD and cofactors, i.e. polynomials ``h``, ``cff`` and ``cfg``
      such that::
  
            h = gcd(f, g), cff = quo(f, h) and cfg = quo(g, h)
  
      The algorithm is purely heuristic which means it may fail to compute
      the GCD. This will be signaled by raising an exception. In this case
      you will need to switch to another GCD method.
  
      The algorithm computes the polynomial GCD by evaluating polynomials
      f and g at certain points and computing (fast) integer GCD of those
      evaluations. The polynomial GCD is recovered from the integer image
      by interpolation. The evaluation process reduces f and g variable by
      variable into a large integer.  The final step is to verify if the
      interpolated polynomial is the correct GCD. This gives cofactors of
      the input polynomials as a side effect.
  
      Examples
      ========
  
      >>> from sympy.polys import ring, ZZ
      >>> R, x,y, = ring("x,y", ZZ)
  
      >>> f = x**2 + 2*x*y + y**2
      >>> g = x**2 + x*y
  
      >>> R.dmp_zz_heu_gcd(f, g)
      (x + y, x + y, x)
  
      References
      ==========
  
      .. [1] [Liao95]_
  
      """
      if not u:
          return dup_zz_heu_gcd(f, g, K)
  
      result = _dmp_rr_trivial_gcd(f, g, u, K)
  
      if result is not None:
          return result
  
      gcd, f, g = dmp_ground_extract(f, g, u, K)
  
      f_norm = dmp_max_norm(f, u, K)
      g_norm = dmp_max_norm(g, u, K)
  
      B = K(2 * min(f_norm, g_norm) + 29)
  
      x = max(
          min(B, 99 * K.sqrt(B)),
          2
          * min(
              f_norm // abs(dmp_ground_LC(f, u, K)),
              g_norm // abs(dmp_ground_LC(g, u, K)),
          )
          + 4,
      )
  
      for i in range(0, HEU_GCD_MAX):
          ff = dmp_eval(f, x, u, K)
          gg = dmp_eval(g, x, u, K)
  
          v = u - 1
  
          if not (dmp_zero_p(ff, v) or dmp_zero_p(gg, v)):
              h, cff, cfg = dmp_zz_heu_gcd(ff, gg, v, K)
  
              h = _dmp_zz_gcd_interpolate(h, x, v, K)
              h = dmp_ground_primitive(h, u, K)[1]
  
              cff_, r = dmp_div(f, h, u, K)
  
              if dmp_zero_p(r, u):
                  cfg_, r = dmp_div(g, h, u, K)
  
                  if dmp_zero_p(r, u):
                      h = dmp_mul_ground(h, gcd, u, K)
                      return h, cff_, cfg_
  
              cff = _dmp_zz_gcd_interpolate(cff, x, v, K)
  
              h, r = dmp_div(f, cff, u, K)
  
              if dmp_zero_p(r, u):
                  cfg_, r = dmp_div(g, h, u, K)
  
                  if dmp_zero_p(r, u):
                      h = dmp_mul_ground(h, gcd, u, K)
                      return h, cff, cfg_
  
              cfg = _dmp_zz_gcd_interpolate(cfg, x, v, K)
  
              h, r = dmp_div(g, cfg, u, K)
  
              if dmp_zero_p(r, u):
                  cff_, r = dmp_div(f, h, u, K)
  
                  if dmp_zero_p(r, u):
                      h = dmp_mul_ground(h, gcd, u, K)
                      return h, cff_, cfg
  
          x = 73794 * x * K.sqrt(K.sqrt(x)) // 27011
  
      raise HeuristicGCDFailed("no luck")
  
imports: from sympy.polys.densearith import dup_sub_mul, dup_neg, dmp_neg, dmp_add, dmp_sub, dup_mul, dmp_mul, dmp_pow, dup_div, dmp_div, dup_rem, dup_quo, dmp_quo, dup_prem, dmp_prem, dup_mul_ground, dmp_mul_ground, dmp_mul_term, dup_quo_ground, dmp_quo_ground, dup_max_norm, dmp_max_norm
  from sympy.polys.densebasic import dup_strip, dmp_raise, dmp_zero, dmp_one, dmp_ground, dmp_one_p, dmp_zero_p, dmp_zeros, dup_degree, dmp_degree, dmp_degree_in, dup_LC, dmp_LC, dmp_ground_LC, dmp_multi_deflate, dmp_inflate, dup_convert, dmp_convert, dmp_apply_pairs
  from sympy.polys.densetools import dup_clear_denoms, dmp_clear_denoms, dup_diff, dmp_diff, dup_eval, dmp_eval, dmp_eval_in, dup_trunc, dmp_ground_trunc, dup_monic, dmp_ground_monic, dup_primitive, dmp_ground_primitive, dup_extract, dmp_ground_extract
  from sympy.polys.galoistools import gf_int, gf_crt
  from sympy.polys.polyconfig import query
  from sympy.polys.polyerrors import MultivariatePolynomialError, HeuristicGCDFailed, HomomorphismFailed, NotInvertible, DomainError
label_desc: Use image node_id 37 for calling a global function with example usage: dmp_zz_heu_gcd(f, g, u, K) and returns: dup_zz_heu_gcd, result, h,  cff_,  cfg_, h,  cff,  cfg_, h,  cff_,  cfg