'Compute 0 to t-th differential moments for Gaussian variable.
E[(P(x+s)/P(x+s-1)-1)^t]
= sum_{i=0}^t (t choose i) (-1)^{t-i} E[(P(x+s)/P(x+s-1))^i]
= sum_{i=0}^t (t choose i) (-1)^{t-i} E[exp(-i*(2*x+2*s-1)/(2*sigma^2))]
= sum_{i=0}^t (t choose i) (-1)^{t-i} exp(i(i+1-2*s)/(2 sigma^2))
Args:
sigma: the noise sigma, in the multiples of the sensitivity.
s: the shift.
t: 0 to t-th moment.
Returns:
0 to t-th moment as a tensor of shape [t+1].'

def _differential_moments(self, sigma, s, t):

    assert (t <= self._max_moment_order), ('The    order    of    %d    is    out    of    the    upper    bound    %d.' % (t, self._max_moment_order))
    binomial = tf.slice(self._binomial_table, [0, 0], [(t + 1), (t + 1)])
    signs = numpy.zeros(((t + 1), (t + 1)), dtype=numpy.float64)
    for i in range((t + 1)):
        for j in range((t + 1)):
            signs[(i, j)] = (1.0 - (2 * ((i - j) % 2)))
    exponents = tf.constant([((j * ((j + 1.0) - (2.0 * s))) / ((2.0 * sigma) * sigma)) for j in range((t + 1))], dtype=tf.float64)
    x = tf.multiply(binomial, signs)
    y = tf.multiply(x, tf.exp(exponents))
    z = tf.reduce_sum(y, 1)
    return z
