Python mpmath.exp() Examples

def gammainc(a, x, dps=50, maxterms=10**8):
    """Compute gammainc exactly like mpmath does but allow for more
    summands in hypercomb. See

    mpmath/functions/expintegrals.py#L134
    
    in the mpmath github repository.

    """
    with mp.workdps(dps):
        z, a, b = mp.mpf(a), mp.mpf(x), mp.mpf(x)
        G = [z]
        negb = mp.fneg(b, exact=True)

        def h(z):
            T1 = [mp.exp(negb), b, z], [1, z, -1], [], G, [1], [1+z], b
            return (T1,)

        res = mp.hypercomb(h, [z], maxterms=maxterms)
        return mpf2float(res) 

def test_spence_circle():
    # The trickiest region for spence is around the circle |z - 1| = 1,
    # so test that region carefully.

    def spence(z):
        return complex(mpmath.polylog(2, 1 - z))

    r = np.linspace(0.5, 1.5)
    theta = np.linspace(0, 2*pi)
    z = (1 + np.outer(r, np.exp(1j*theta))).flatten()
    dataset = []
    for z0 in z:
        dataset.append((z0, spence(z0)))

    dataset = np.array(dataset)
    FuncData(sc.spence, dataset, 0, 1, rtol=1e-14).check()


# ------------------------------------------------------------------------------
# sinpi and cospi
# ------------------------------------------------------------------------------ 

def compute_b_mp(self, sigma, q, order, verbose=False):
    """Compute B_lambda for arbitrary lambda by numerical integration."""
    mu0, _, mu = self._distributions_mp(sigma, q)
    b_lambda_fn = lambda z: mu0(z) * (mu0(z) / mu(z)) ** order
    b_numeric = self._integral_inf_mp(b_lambda_fn)

    if verbose:
      _, z1 = rdp_accountant._compute_zs(sigma, q)
      print("z1 = ", z1)
      print("x in the Taylor series = ", q / (1 - q) * np.exp(
          (2 * z1 - 1) / (2 * sigma ** 2)))

      b0_numeric = self._integral_bounded_mp(b_lambda_fn, -np.inf, z1)
      b1_numeric = self._integral_bounded_mp(b_lambda_fn, z1, +np.inf)

      print("B: numerically {} = {} + {}".format(b_numeric, b0_numeric,
                                                 b1_numeric))
    return float(b_numeric) 

def gammainc(a, x, dps=50, maxterms=10**8):
    """Compute gammainc exactly like mpmath does but allow for more
    summands in hypercomb. See

    mpmath/functions/expintegrals.py#L134
    
    in the mpmath github repository.

    """
    with mp.workdps(dps):
        z, a, b = mp.mpf(a), mp.mpf(x), mp.mpf(x)
        G = [z]
        negb = mp.fneg(b, exact=True)

        def h(z):
            T1 = [mp.exp(negb), b, z], [1, z, -1], [], G, [1], [1+z], b
            return (T1,)

        res = mp.hypercomb(h, [z], maxterms=maxterms)
        return mpf2float(res) 

def _compute_delta(log_moments, eps):
  """Compute delta for given log_moments and eps.
  Args:
    log_moments: the log moments of privacy loss, in the form of pairs
      of (moment_order, log_moment)
    eps: the target epsilon.
  Returns:
    delta
  """
  min_delta = 1.0
  for moment_order, log_moment in log_moments:
    if moment_order == 0:
      continue
    if math.isinf(log_moment) or math.isnan(log_moment):
      sys.stderr.write("The %d-th order is inf or Nan\n" % moment_order)
      continue
    if log_moment < moment_order * eps:
      min_delta = min(min_delta,
                      math.exp(log_moment - moment_order * eps))
  return min_delta 

def compute_a(sigma, q, lmbd, verbose=False):
  lmbd_int = int(math.ceil(lmbd))
  if lmbd_int == 0:
    return 1.0

  a_lambda_first_term_exact = 0
  a_lambda_second_term_exact = 0
  for i in xrange(lmbd_int + 1):
    coef_i = scipy.special.binom(lmbd_int, i) * (q ** i)
    s1, s2 = 0, 0
    for j in xrange(i + 1):
      coef_j = scipy.special.binom(i, j) * (-1) ** (i - j)
      s1 += coef_j * np.exp((j * j - j) / (2.0 * (sigma ** 2)))
      s2 += coef_j * np.exp((j * j + j) / (2.0 * (sigma ** 2)))
    a_lambda_first_term_exact += coef_i * s1
    a_lambda_second_term_exact += coef_i * s2

  a_lambda_exact = ((1.0 - q) * a_lambda_first_term_exact +
                    q * a_lambda_second_term_exact)
  if verbose:
    print "A: by binomial expansion    {} = {} + {}".format(
        a_lambda_exact,
        (1.0 - q) * a_lambda_first_term_exact,
        q * a_lambda_second_term_exact)
  return _to_np_float64(a_lambda_exact) 

def test_lanczos_sum_expg_scaled(self):
        maxgamma = 171.624376956302725
        e = np.exp(1)
        g = 6.024680040776729583740234375

        def gamma(x):
            with np.errstate(over='ignore'):
                fac = ((x + g - 0.5)/e)**(x - 0.5)
                if fac != np.inf:
                    res = fac*_lanczos_sum_expg_scaled(x)
                else:
                    fac = ((x + g - 0.5)/e)**(0.5*(x - 0.5))
                    res = fac*_lanczos_sum_expg_scaled(x)
                    res *= fac
            return res

        assert_mpmath_equal(gamma,
                            mpmath.gamma,
                            [Arg(0, maxgamma, inclusive_a=False)],
                            rtol=1e-13) 

def gammainc(a, x, dps=50, maxterms=10**8):
    """Compute gammainc exactly like mpmath does but allow for more
    summands in hypercomb. See

    mpmath/functions/expintegrals.py#L134
    
    in the mpmath github repository.

    """
    with mp.workdps(dps):
        z, a, b = mp.mpf(a), mp.mpf(x), mp.mpf(x)
        G = [z]
        negb = mp.fneg(b, exact=True)

        def h(z):
            T1 = [mp.exp(negb), b, z], [1, z, -1], [], G, [1], [1+z], b
            return (T1,)

        res = mp.hypercomb(h, [z], maxterms=maxterms)
        return mpf2float(res) 

def test_loggamma_taylor():
    # Test around the zeros at z = 1, 2.

    r = np.logspace(-16, np.log10(LOGGAMMA_TAYLOR_RADIUS), 10)
    theta = np.linspace(0, 2*np.pi, 20)
    r, theta = np.meshgrid(r, theta)
    dz = r*np.exp(1j*theta)
    z = np.r_[1 + dz, 2 + dz].flatten()

    dataset = []
    for z0 in z:
        dataset.append((z0, complex(mpmath.loggamma(z0))))
    dataset = np.array(dataset)

    FuncData(sc.loggamma, dataset, 0, 1, rtol=5e-14).check()


# ------------------------------------------------------------------------------
# rgamma
# ------------------------------------------------------------------------------ 

def test_loggamma_taylor_transition():
    # Make sure there isn't a big jump in accuracy when we move from
    # using the Taylor series to using the recurrence relation.

    r = LOGGAMMA_TAYLOR_RADIUS + np.array([-0.1, -0.01, 0, 0.01, 0.1])
    theta = np.linspace(0, 2*np.pi, 20)
    r, theta = np.meshgrid(r, theta)
    dz = r*np.exp(1j*theta)
    z = np.r_[1 + dz, 2 + dz].flatten()

    dataset = []
    for z0 in z:
        dataset.append((z0, complex(mpmath.loggamma(z0))))
    dataset = np.array(dataset)

    FuncData(sc.loggamma, dataset, 0, 1, rtol=5e-14).check() 

def zpkfreqz(z, p, k, worN=None):
    """
    Frequency response of a filter in zpk format, using mpmath.

    This is the same calculation as scipy.signal.freqz, but the input is in
    zpk format, the calculation is performed using mpath, and the results are
    returned in lists instead of numpy arrays.
    """
    if worN is None or isinstance(worN, int):
        N = worN or 512
        ws = [mpmath.pi * mpmath.mpf(j) / N for j in range(N)]
    else:
        ws = worN

    h = []
    for wk in ws:
        zm1 = mpmath.exp(1j * wk)
        numer = _prod([zm1 - t for t in z])
        denom = _prod([zm1 - t for t in p])
        hk = k * numer / denom
        h.append(hk)
    return ws, h 

def test_expi_complex():
    dataset = []
    for r in np.logspace(-99, 2, 10):
        for p in np.linspace(0, 2*np.pi, 30):
            z = r*np.exp(1j*p)
            dataset.append((z, complex(mpmath.ei(z))))
    dataset = np.array(dataset, dtype=np.complex_)

    FuncData(sc.expi, dataset, 0, 1).check()


# ------------------------------------------------------------------------------
# expn
# ------------------------------------------------------------------------------ 

def test_log(self):
        with mp.workdps(30):
            logcoeffs = mp.taylor(lambda x: mp.log(1 + x), 0, 10)
            expcoeffs = mp.taylor(lambda x: mp.exp(x) - 1, 0, 10)
            invlogcoeffs = lagrange_inversion(logcoeffs)
            mp_assert_allclose(invlogcoeffs, expcoeffs) 

def test_igam_fac(self):
        def mp_igam_fac(a, x):
            return mpmath.power(x, a)*mpmath.exp(-x)/mpmath.gamma(a)

        assert_mpmath_equal(_igam_fac,
                            mp_igam_fac,
                            [Arg(0, 1e14, inclusive_a=False), Arg(0, 1e14)],
                            rtol=1e-10) 

def _noncentral_chi_pdf(t, df, nc):
    res = mpmath.besseli(df/2 - 1, mpmath.sqrt(nc*t))
    res *= mpmath.exp(-(t + nc)/2)*(t/nc)**(df/4 - 1/2)/2
    return res 

def gammaincc(a, x, dps=50, maxterms=10**8):
    """Compute gammaincc exactly like mpmath does but allow for more
    terms in hypercomb. See

    mpmath/functions/expintegrals.py#L187

    in the mpmath github repository.

    """
    with mp.workdps(dps):
        z, a = a, x
        
        if mp.isint(z):
            try:
                # mpmath has a fast integer path
                return mpf2float(mp.gammainc(z, a=a, regularized=True))
            except mp.libmp.NoConvergence:
                pass
        nega = mp.fneg(a, exact=True)
        G = [z]
        # Use 2F0 series when possible; fall back to lower gamma representation
        try:
            def h(z):
                r = z-1
                return [([mp.exp(nega), a], [1, r], [], G, [1, -r], [], 1/nega)]
            return mpf2float(mp.hypercomb(h, [z], force_series=True))
        except mp.libmp.NoConvergence:
            def h(z):
                T1 = [], [1, z-1], [z], G, [], [], 0
                T2 = [-mp.exp(nega), a, z], [1, z, -1], [], G, [1], [1+z], a
                return T1, T2
            return mpf2float(mp.hypercomb(h, [z], maxterms=maxterms)) 

def pdf_gauss_mp(x, sigma, mean):
  return mp.mpf(1.) / mp.sqrt(mp.mpf("2.") * sigma ** 2 * mp.pi) * mp.exp(
      - (x - mean) ** 2 / (mp.mpf("2.") * sigma ** 2)) 

def _mu1_over_mu0(self, x, sigma):
    # Closed-form expression for N(1, sigma^2) / N(0, sigma^2) at x.
    return exp((2 * x - 1) / (2 * sigma**2)) 

def _butter_analog_poles(n):
    """
    Poles of an analog Butterworth lowpass filter.

    This is the same calculation as scipy.signal.buttap(n) or
    scipy.signal.butter(n, 1, analog=True, output='zpk'), but mpmath is used,
    and only the poles are returned.
    """
    poles = []
    for k in range(-n+1, n, 2):
        poles.append(-mpmath.exp(1j*mpmath.pi*k/(2*n)))
    return poles 

def gammaincc(a, x, dps=50, maxterms=10**8):
    """Compute gammaincc exactly like mpmath does but allow for more
    terms in hypercomb. See

    mpmath/functions/expintegrals.py#L187

    in the mpmath github repository.

    """
    with mp.workdps(dps):
        z, a = a, x
        
        if mp.isint(z):
            try:
                # mpmath has a fast integer path
                return mpf2float(mp.gammainc(z, a=a, regularized=True))
            except mp.libmp.NoConvergence:
                pass
        nega = mp.fneg(a, exact=True)
        G = [z]
        # Use 2F0 series when possible; fall back to lower gamma representation
        try:
            def h(z):
                r = z-1
                return [([mp.exp(nega), a], [1, r], [], G, [1, -r], [], 1/nega)]
            return mpf2float(mp.hypercomb(h, [z], force_series=True))
        except mp.libmp.NoConvergence:
            def h(z):
                T1 = [], [1, z-1], [z], G, [], [], 0
                T2 = [-mp.exp(nega), a, z], [1, z, -1], [], G, [1], [1+z], a
                return T1, T2
            return mpf2float(mp.hypercomb(h, [z], maxterms=maxterms)) 

def _mu1_over_mu0(self, x, sigma):
    # Closed-form expression for N(1, sigma^2) / N(0, sigma^2) at x.
    return mp.exp((2 * x - 1) / (2 * sigma**2)) 

def _pdf_gauss_mp(self, x, sigma, mean):
    return 1. / mp.sqrt(2. * sigma ** 2 * mp.pi) * mp.exp(
        -(x - mean) ** 2 / (2. * sigma ** 2)) 

def gammaincc(a, x, dps=50, maxterms=10**8):
    """Compute gammaincc exactly like mpmath does but allow for more
    terms in hypercomb. See

    mpmath/functions/expintegrals.py#L187

    in the mpmath github repository.

    """
    with mp.workdps(dps):
        z, a = a, x
        
        if mp.isint(z):
            try:
                # mpmath has a fast integer path
                return mpf2float(mp.gammainc(z, a=a, regularized=True))
            except mp.libmp.NoConvergence:
                pass
        nega = mp.fneg(a, exact=True)
        G = [z]
        # Use 2F0 series when possible; fall back to lower gamma representation
        try:
            def h(z):
                r = z-1
                return [([mp.exp(nega), a], [1, r], [], G, [1, -r], [], 1/nega)]
            return mpf2float(mp.hypercomb(h, [z], force_series=True))
        except mp.libmp.NoConvergence:
            def h(z):
                T1 = [], [1, z-1], [z], G, [], [], 0
                T2 = [-mp.exp(nega), a, z], [1, z, -1], [], G, [1], [1+z], a
                return T1, T2
            return mpf2float(mp.hypercomb(h, [z], maxterms=maxterms))