code: def heugcd(f, g):
      """
      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.heuristicgcd import heugcd
      >>> 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
  
      >>> h, cff, cfg = heugcd(f, g)
      >>> h, cff, cfg
      (x + y, x + y, x)
  
      >>> cff*h == f
      True
      >>> cfg*h == g
      True
  
      References
      ==========
  
      .. [1] [Liao95]_
  
      """
      assert f.ring == g.ring and f.ring.domain.is_ZZ
  
      ring = f.ring
      x0 = ring.gens[0]
      domain = ring.domain
  
      gcd, f, g = f.extract_ground(g)
  
      f_norm = f.max_norm()
      g_norm = g.max_norm()
  
      B = domain(2 * min(f_norm, g_norm) + 29)
  
      x = max(
          min(B, 99 * domain.sqrt(B)),
          2 * min(f_norm // abs(f.LC), g_norm // abs(g.LC)) + 4,
      )
  
      for i in range(0, HEU_GCD_MAX):
          ff = f.evaluate(x0, x)
          gg = g.evaluate(x0, x)
  
          if ff and gg:
              if ring.ngens == 1:
                  h, cff, cfg = domain.cofactors(ff, gg)
              else:
                  h, cff, cfg = heugcd(ff, gg)
  
              h = _gcd_interpolate(h, x, ring)
              h = h.primitive()[1]
  
              cff_, r = f.div(h)
  
              if not r:
                  cfg_, r = g.div(h)
  
                  if not r:
                      h = h.mul_ground(gcd)
                      return h, cff_, cfg_
  
              cff = _gcd_interpolate(cff, x, ring)
  
              h, r = f.div(cff)
  
              if not r:
                  cfg_, r = g.div(h)
  
                  if not r:
                      h = h.mul_ground(gcd)
                      return h, cff, cfg_
  
              cfg = _gcd_interpolate(cfg, x, ring)
  
              h, r = g.div(cfg)
  
              if not r:
                  cff_, r = f.div(h)
  
                  if not r:
                      h = h.mul_ground(gcd)
                      return h, cff_, cfg
  
          x = 73794 * x * domain.sqrt(domain.sqrt(x)) // 27011
  
      raise HeuristicGCDFailed("no luck")
  
label_desc: Use image node_id 1 for calling a global function with example usage: heugcd(f, g) and returns: h,  cff_,  cfg_, h,  cff,  cfg_, h,  cff_,  cfg