Created
October 18, 2017 18:30
-
-
Save jaishirb/e7ff934090ce96ef620aaeeac6701ddc to your computer and use it in GitHub Desktop.
Get the minimum n and k values in order to get a value greater than or equal to a constant c.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
from scipy.special import lambertw | |
from scipy.special import gamma | |
import math | |
import time | |
e = math.e | |
def l_func(x): | |
c = 0.036534 | |
return math.log((x + c) / math.sqrt(2 * math.pi)) | |
def inv_gamma(x): | |
global e | |
return round(l_func(x) / lambertw([l_func(x) / e]).real[0] - 0.5) | |
def fact(n): | |
return gamma(n + 1) | |
def combine(n, k): | |
return fact(n) / (fact(n - k) * fact(k)) | |
def get_lim(c, t): | |
while combine(t, t / 2) < c: | |
t += 1 | |
return t | |
def digits(n): | |
return math.log10(n) + 1 | |
def get_t(t, c, d=2): | |
prev = d | |
while combine(t, t / 2) <= c: | |
prev = t | |
d += 1 | |
t = inv_gamma(c ** d + 1) | |
return prev | |
def get_min_k(t, max_k, c): | |
min_k = max_k - 1 | |
while combine(t, min_k) >= c: | |
min_k -= 1 | |
return min_k + 1 | |
if __name__ == '__main__': | |
print("Write the result of the combinatorial: ", end='') | |
c = int(input()) | |
if c > 0: | |
start = time.time() | |
t = inv_gamma(c ** 2 + 1) | |
t = int(get_t(t, c)) | |
print("start at: {}".format(t)) # print where brute force start | |
t = int(get_lim(c, t)) | |
print("Real val: {}".format(t)) | |
max_k = int(t / 2) if t % 2 == 0 else int((t + 1) / 2) | |
min_k = get_min_k(t, max_k, c) | |
found = int(combine(t, min_k)) | |
print("the min value of n is: {}, max k value: {}, min k value: {}.".format(t, max_k, min_k)) | |
print("Value found: {}".format(found)) | |
print("Distance between values: {}".format(found - c)) | |
print("Elapsed time: {}".format(time.time() - start)) | |
else: | |
print("Indeterminate") |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment