code: def hessian_tensor_product(fun, argnum=0):
      
      fun_grad = grad(fun, argnum)
      def vector_dot_grad(*args, **kwargs):
          args, vector = args[:-1], args[-1]
          return np.tensordot(fun_grad(*args, **kwargs), vector, np.ndim(vector))
      return grad(vector_dot_grad, argnum)
docs: Builds a function that returns the exact Hessian-tensor product.
      The returned function has arguments (*args, tensor, **kwargs), and for
      vectors takes roughly 4x as long to evaluate as the original function.