fast_float/script/mushtak_lemire.py

76 lines
2.2 KiB
Python
Raw Permalink Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

#
# Reference :
# Noble Mushtak and Daniel Lemire, Fast Number Parsing Without Fallback (to appear)
#
all_tqs = []
# Generates all possible values of T[q]
# Appendix B of Number parsing at a gigabyte per second.
# Software: Practice and Experience 2021;51(8):17001727.
for q in range(-342, -27):
power5 = 5**-q
z = 0
while (1 << z) < power5:
z += 1
b = 2 * z + 2 * 64
c = 2**b // power5 + 1
while c >= (1 << 128):
c //= 2
all_tqs.append(c)
for q in range(-27, 0):
power5 = 5**-q
z = 0
while (1 << z) < power5:
z += 1
b = z + 127
c = 2**b // power5 + 1
all_tqs.append(c)
for q in range(0, 308 + 1):
power5 = 5**q
while power5 < (1 << 127):
power5 *= 2
while power5 >= (1 << 128):
power5 //= 2
all_tqs.append(power5)
# Returns the continued fraction of numer/denom as a list [a0; a1, a2, ..., an]
def continued_fraction(numer, denom):
# (look at page numbers in top-left, not PDF page numbers)
cf = []
while denom != 0:
quot, rem = divmod(numer, denom)
cf.append(quot)
numer, denom = denom, rem
return cf
# Given a continued fraction [a0; a1, a2, ..., an], returns
# all the convergents of that continued fraction
# as pairs of the form (numer, denom), where numer/denom is
# a convergent of the continued fraction in simple form.
def convergents(cf):
p_n_minus_2 = 0
q_n_minus_2 = 1
p_n_minus_1 = 1
q_n_minus_1 = 0
convergents = []
for a_n in cf:
p_n = a_n * p_n_minus_1 + p_n_minus_2
q_n = a_n * q_n_minus_1 + q_n_minus_2
convergents.append((p_n, q_n))
p_n_minus_2, q_n_minus_2, p_n_minus_1, q_n_minus_1 = p_n_minus_1, q_n_minus_1, p_n, q_n
return convergents
# Enumerate through all the convergents of T[q] / 2^137 with denominators < 2^64
found_solution = False
for j, tq in enumerate(all_tqs):
for _, w in convergents(continued_fraction(tq, 2**137)):
if w >= 2**64:
break
if (tq*w) % 2**137 > 2**137 - 2**64:
print(f"SOLUTION: q={j-342} T[q]={tq} w={w}")
found_solution = True
if not found_solution:
print("No solutions!")