from mpmath import nsum, inf, nprod
import math
def fgh(iter, num):
    if iter == 0:
        return num + 1
    elif iter == 1:
        return num * 2
    elif iter == 2:
        return 2 ** num
    else:
        result = fgh(iter - 1, num)
        for i in range(num - 1):
            result = fgh(iter - 1, result)
        return result

def sgh(iter, num):
    if iter == 0:
        return 0
    else:
        return sgh(iter - 1, num) + 1

def mgh(iter, num):
    if iter == 0:
        return num + 1
    else:
        return mgh(iter - 1, mgh(iter - 1, num))

def hh(iter, num):
    if iter == 0:
        return num
    else:
        return hh(iter - 1, num + 1)

def fact(num):
    if num == 1:
        return 1
    else:
        return num * fact(num - 1)

def pochhammer_symbol(num, iter):
    if num == 1 or iter == 0:
        return num
    else:
        result = num
        for it in range(1, iter):
            result *= it + num
        return result

def qfact(k, q):
    result = 1
    chapter = 1
    for iter in range(1, k):
        chapter += q ** iter
        result *= chapter
    return result

def qseries(a, q, n):
    result = 1 - a
    for _ in range(1, n):
        result *= 1 - a * (q ** _)
    return result

def qps(a, q, k=inf):
    if k > 0:
        return nprod(lambda j: 1 - a * (q ** j), [0, k - 1])
    elif k == 0:
        return 1
    elif k < 0:
        return nprod(lambda j: (1 - a * (q ** -j)) ** -1, [0, math.abs(k - 1)])
def qbc(n, k):
    if 0 <= k and k < n:
        return fact(n) / (fact(k) * fact(n - k))
    else:
        return 0

def qhypgeom(asl, bsl, q, z):
    return nsum(lambda x: (nprod(lambda a: qps(asl[int(a)], q, x), [0, len(asl) - 1]) / nprod(lambda b: qps(bsl[int(b)], q, x), [0, len(bsl) - 1]) * ((z ** x) / qps(q, q, x)) * ((-1 ** x) * (q ** (x * (x - 1)/2))) ** (1 + len(bsl) - len(asl))), [0, inf])

def qe(q, x):
    return nsum(lambda x: (z ** x) / qps(q, q, x), [0, inf])
def qE(q, x):
    return qps(-x, q, inf)
def gamma(z):
    #Euler's limit
    return (1.0 / z) * nprod(lambda x: ((1 + (1.0 / x)) ** z) * ((1 + (z / x)) ** -1), [1, inf])
def qgamma(q, z):
    return (qps(q, q, inf) / qps(q ** z, q, inf)) * ((1 - q) ** (1 - x))
def beta(q, a, b):
    return (qgamma(q, a) * qgamma(q, b)) / qgamma(a + b)