'''
Factorials with BIG numbers, i'm talking huge. Ever wonder what (2**128)! is like?

Uses python's Decimal module. Stirling approximation for the factorial calculation.

The bigger the numbers - the higher Emax needs to be. The closer to zero - the
farther down Emin needs to be pushed.

built for python 3
'''

import time
import pprint
import math
import decimal
from decimal import Decimal
import profile

ACCURACY_NEEDED = Decimal.from_float(0.1)
TESTS = 10, 100, 1000, 10000

# generic constants
PI = Decimal.from_float(math.pi)
E = Decimal.from_float(math.e)

def stirling_factorial(n):
    # for some odd reason, python's Decimal.from_float throws on a Decimal
    if not isinstance(n, Decimal):
        n = Decimal.from_float(n)
    root = (PI * (Decimal('2.0') * n + Decimal.from_float(1 / 3)))**Decimal('0.5')
    return root * n**n * E**(-n)

def compare_to_factorial():
    print('---how close is stirling(i) to i!')
    results = []
    for i in TESTS:
        fact = math.factorial(i)
        stir = stirling_factorial(i)
        diff = abs(fact - stir) / fact
        results.append(diff < ACCURACY_NEEDED)
        print(i, diff)

    return results

def check_consecutive_numbers():
    results = []
    print('---how good is n!/(n-1)! = n')

    for i in range(64, 148):
        orig = 2**i
        x = stirling_factorial(orig)
        y = stirling_factorial(orig - 1)
        maybe_orig = x / y

        # this turns out to be negative,
        # I don't know what that means...
        diff = abs(orig - maybe_orig) / orig

        print('orig vs maybe: %d!' % i, diff)
        
        results.append(diff < ACCURACY_NEEDED)

    return results
    
    
def quality_test():
    results = []
    
    results += check_consecutive_numbers()
    results += compare_to_factorial()
    
    if all(results):
        print(':)')
        return True
    else:
        # something went wrong, check Emax, Emin, and prec
        print(':(')
        return False

def speed_test():
    # timeit only takes small bytes of code
    # and profile doesn't work for some odd reason.
    num = 10000
    t1 = time.time()
    x = math.factorial(num)
    t2 = time.time()
    y = stirling_factorial(num)
    t3 = time.time()
    print('---the speed of the function')
    print('factorial: %g' % (t2 - t1))
    print('stirling factorial: %g' % (t3 - t2))
    

def bigger_numbers():
    print('---now checking what these giant numbers look like')
    x = stirling_factorial(2**128)
    print('big number:', str(x))
    s = Decimal.from_float(2**128)
    n = Decimal.from_float(65000000)
    print((s**n))
    print(stirling_factorial(s - n))
    print((s**n) * stirling_factorial(s - n))
    print(stirling_factorial(s))
    num = stirling_factorial(s) / ((s**n) * stirling_factorial(s - n))
    print('another number:', str(num))


decimal.getcontext().Emax = 10**600
decimal.getcontext().Emin = -10**600

# the default prec 28 means stirling works only up to 2**93
decimal.getcontext().prec = 50

#profile.run('quality_test()')
quality_test()
speed_test()

bigger_numbers()