code: def reduce_var_square(
      var_square, avg_diff, count, op, axis, sum_func
  ):
      moment = op.moment
      dtype = op.dtype
      kw = dict(axis=axis, dtype=dtype, keepdims=bool(op.keepdims))
  
      reduced_var_square = var_square[..., moment - 2].sum(
          **kw
      ) + sum_func(count * avg_diff**moment, **kw)
      for i in range(1, moment - 1):
          coeff = factorial(moment) / float(
              factorial(i) * factorial(moment - i)
          )
          reduced_var_square += coeff * sum_func(
              var_square[..., moment - i - 2] * avg_diff**moment, **kw
          )
      return reduced_var_square
  
imports: import numpy
  from math import factorial
  from ...None import opcodes
  from ...serialization.serializables import Int32Field
  from ..datasource import tensor
  from ..array_utils import as_same_device, device, get_array_module
  from .core import TensorReduction, TensorReductionMixin, numel
label_desc: Use image node_id 1 for calling a global function with example usage: reduce_var_square(var_square, avg_diff, count, op, axis, sum_func) and returns: reduced_var_square