diff --git a/ad_gcd_lcm.py b/ad_gcd_lcm.py index 03ac328..66cf3dc 100755 --- a/ad_gcd_lcm.py +++ b/ad_gcd_lcm.py @@ -1,4 +1,4 @@ -#!/usr/bin/python +#!/usr/bin/python3 # find the gcd and lcm of a given list of floating-point numbers. this # code uses a best rational approximation algorithm to solve the @@ -7,6 +7,7 @@ # rational numbers. # author: ali dasdan +# updated 2024 onwards by Tim Abdiukov. import sys import getopt @@ -14,12 +15,27 @@ from math import * from ad_rat_by_fast_farey import * +USE_NUMPY_HI_PRECISION = True +if(USE_NUMPY_HI_PRECISION): + try: + import numpy as np + float2 = lambda x: np.longdouble(x) + except ModuleNotFoundError: + print("Please disable USE_NUMPY_HI_PRECISION or install prerequisite,") + print("```") + print("pip install numpy") + print("```") + exit() +else: + print("WARN: High precision NumPy logic is not used") + float2 = lambda x: float(x) + def show_usage(): - print "Usage: " + sys.argv[0] + " -l/--limit=int>=1 -n/--nums=quoted list of at least 2 numbers or math expression returning float (w/o spaces)" + print("Usage: " + sys.argv[0] + " -l/--limit=int>=1 -n/--nums=quoted list of at least 2 numbers or math expression returning float (w/o spaces)") def at_exit(msg): if msg != '' and msg != None: - print "Error:", msg + print("Error:", msg) show_usage() sys.exit(0) @@ -84,7 +100,7 @@ def main(): try: opts, args = getopt.getopt(sys.argv[1:], 'l:n:', ['limit=', 'nums']) - except getopt.GetoptError, msg: + except getopt.GetoptError as msg: at_exit(msg) for o, a in opts: @@ -96,7 +112,7 @@ def main(): if l < 1: raise Exception() else: raise Exception() - except Exception, msg: + except Exception as msg: at_exit(msg) if s == None: @@ -104,7 +120,7 @@ def main(): # get the list from the input string lst = re.split("[\s,;|:]+", s) - print "input:", lst + print("input:", lst) if len(lst) <= 1 or l == None: at_exit("Bad args") @@ -113,10 +129,10 @@ def main(): for x in lst: try: nums.append(eval(x)) - except Exception, msg: + except Exception as msg: at_exit(msg) - print "input evaluated:", nums + print("input evaluated:", nums) # convert to best rational approximations ns = [] @@ -130,17 +146,17 @@ def main(): for n, d in zip(ns, ds): rats.append(str(n) + "/" + str(d)) srats = ", ".join(rats) - print "best_rats= [" + srats + "]" + print("best_rats= [" + srats + "]") # compute gcd and lcm of the elts in the list n_gcd, d_gcd = rgcd_lst(ns, ds) n_lcm, d_lcm = rlcm_lst(ns, ds) - gcd = float(n_gcd) / d_gcd - lcm = float(n_lcm) / d_lcm + gcd = float2(n_gcd) / d_gcd + lcm = float2(n_lcm) / d_lcm print("gcd: rat= %d / %d val= %g" % (n_gcd, d_gcd, gcd)) - print("lcm: rat= %d / %d val= %g lcm/gcd= %g" % (n_lcm, d_lcm, lcm, float(lcm)/gcd)) + print("lcm: rat= %d / %d val= %g lcm/gcd= %g" % (n_lcm, d_lcm, lcm, float2(lcm)/gcd)) # compute the ratios and measure error (difference from the # closest int) @@ -150,21 +166,21 @@ def main(): n_div_gcd = [] lcm_div_n = [] for x in nums: - r = float(x) / gcd + r = float2(x) / gcd r_int = int(round(r)) err = abs(r - r_int) if err > max_err_gcd: max_err_gcd = err n_div_gcd.append(r_int) - r = float(lcm) / x + r = float2(lcm) / x r_int = int(round(r)) err = abs(r - r_int) if err > max_err_lcm: max_err_lcm = err lcm_div_n.append(r_int) - print "nums_div_gcd= ", n_div_gcd, "max_err= %g" % (max_err_gcd) - print "lcm_div_nums= ", lcm_div_n, "max_err= %g" % (max_err_lcm) + print("nums_div_gcd= ", n_div_gcd, "max_err= %g" % (max_err_gcd)) + print("lcm_div_nums= ", lcm_div_n, "max_err= %g" % (max_err_lcm)) if __name__ == "__main__": main() diff --git a/ad_rat_by_cont_frac.py b/ad_rat_by_cont_frac.py index 5e3daf7..b3a1080 100755 --- a/ad_rat_by_cont_frac.py +++ b/ad_rat_by_cont_frac.py @@ -1,4 +1,4 @@ -#!/usr/bin/python +#!/usr/bin/python3 # find the best rational approximation n/d where d <= l to the given # target number t. here the best means the one with the smallest @@ -15,17 +15,33 @@ # just keep the last partial product of these matrics. # author: ali dasdan (based on the C code by d. eppstein) +# updated 2024 onwards by Tim Abdiukov. import sys import getopt from math import * +USE_NUMPY_HI_PRECISION = True +if(USE_NUMPY_HI_PRECISION): + try: + import numpy as np + float2 = lambda x: np.longdouble(x) + except ModuleNotFoundError: + print("Please disable USE_NUMPY_HI_PRECISION or install prerequisite,") + print("```") + print("pip install numpy") + print("```") + exit() +else: + print("WARN: High precision NumPy logic is not used") + float2 = lambda x: float(x) + def show_usage(): - print sys.argv[0] + " -h/--help [-e/--error=float>=0] -l/--limit=int>1 -t/--target=float or quoted math expr returning float" + print(sys.argv[0] + " -h/--help [-e/--error=float>=0] -l/--limit=int>1 -t/--target=float or quoted math expr returning float") def at_exit(msg): if msg != "" and msg != None: - print "Error:", msg + print("Error:", msg) show_usage() sys.exit(0) @@ -61,8 +77,8 @@ def find_best_rat(l, t): tmp = m10 * ai + m11 m11 = m10 m10 = tmp - if (x == float(ai)): break - x = float(1.0) / (x - float(ai)) + if (x == float2(ai)): break + x = float2(1.0) / (x - float2(ai)) ai = int(x) # now remaining x is between 0 and 1/ai. approx as either o or 1/m @@ -71,18 +87,22 @@ def find_best_rat(l, t): # first try zero n1 = m00 d1 = m10 - err1 = (t - float(n1) / d1) + err1 = (t - float2(n1) / d1) # try the other possibility - ai = int(float(l - m11) / m10) + ai = int(float2(l - m11) / m10) n2 = m00 * ai + m01 d2 = m10 * ai + m11 - err2 = (t - float(n2) / d2) + err2 = (t - float2(n2) / d2) - if abs(err1) <= abs(err2): - return err1, n1, d1, niter - else: - return err2, n2, d2, niter + # optimization + best_err = min(abs(err1), abs(err2)) + + # assert precision wasn't lost + assert(type(best_err) == type(float2(pi))) + + # return data + return best_err, n2, d2, niter # this function takes in an error bound err_in, an int limit l, and # the target fraction to approximate t and returns err_out, n, d, @@ -110,7 +130,7 @@ def main(): opts, args = getopt.getopt(sys.argv[1:], 'he:l:t:', ['help', 'error=', 'limit=', 'target=']) - except getopt.GetoptError, msg: + except getopt.GetoptError as msg: at_exit(msg) for o, a in opts: @@ -118,20 +138,20 @@ def main(): if o in ('-h', '--help'): raise Exception() elif o in ('-e', '--error'): - eps = float(a) + eps = float2(a) if eps <= 0: raise Exception() elif o in ('-l', '--limit'): l = int(a) if l < 1: raise Exception() elif o in ('-t', '--target'): try: - t = float(eval(a)) - except Exception, msg: + t = float2(eval(a)) + except Exception as msg: at_exit(msg) if t <= 0: raise Exception() else: raise Exception() - except Exception, msg: + except Exception as msg: at_exit(msg) if t == None or l==None: @@ -147,7 +167,7 @@ def main(): if eps == None: print("target= %f best_rat= %d / %d max_denom= %d err= %g abs_err= %g niter= %d" % (t, n, d, l, err, abs(err), niter)) else: - print("target= %f best_rat= %d / %d max_denom= %d err= %g abs_err= %g abs_err/error= %g niter= %d" % (t, n, d, l, err, abs(err), float(abs(err)) / eps, niter)) + print("target= %f best_rat= %d / %d max_denom= %d err= %g abs_err= %g abs_err/error= %g niter= %d" % (t, n, d, l, err, abs(err), float2(abs(err)) / eps, niter)) if __name__ == '__main__': main() diff --git a/ad_rat_by_exhaustive.py b/ad_rat_by_exhaustive.py index 922e328..5798434 100755 --- a/ad_rat_by_exhaustive.py +++ b/ad_rat_by_exhaustive.py @@ -1,4 +1,4 @@ -#!/usr/bin/python +#!/usr/bin/python3 # find the best rational approximation n/d where d <= l to the given # target number t. here the best means the one with the smallest @@ -8,17 +8,34 @@ # up to and including l. # author: ali dasdan +# updated 2024 onwards by Tim Abdiukov. + import sys import getopt from math import * +USE_NUMPY_HI_PRECISION = True +if(USE_NUMPY_HI_PRECISION): + try: + import numpy as np + float2 = lambda x: np.longdouble(x) + except ModuleNotFoundError: + print("Please disable USE_NUMPY_HI_PRECISION or install prerequisite,") + print("```") + print("pip install numpy") + print("```") + exit() +else: + print("WARN: High precision NumPy logic is not used") + float2 = lambda x: float(x) + def show_usage(): - print sys.argv[0] + " -h/--help [-e/--error=float>=0] -l/--limit=int=>1 -t/--target=float or quoted math expr returning float" + print(sys.argv[0] + " -h/--help [-e/--error=float>=0] -l/--limit=int=>1 -t/--target=float or quoted math expr returning float") def at_exit(msg): if msg != "" and msg != None: - print "Error:", msg + print("Error:", msg) show_usage() sys.exit(0) @@ -40,11 +57,15 @@ def find_best_rat(l, t): for n in range(1, l + 1): niter += 1 - err = abs(t - float(n) / d) + err = abs(t - float2(n) / d) if err < best_err: best_err = err best_n, best_d = n, d + # assert precision wasn't lost + assert(type(best_err) == type(float2(pi))) + + # return data return best_err, best_n, best_d, niter # this function takes in an error bound err_in, an int limit l, and @@ -73,7 +94,7 @@ def main(): opts, args = getopt.getopt(sys.argv[1:], 'he:l:t:', ['help', 'error=', 'limit=', 'target=']) - except getopt.GetoptError, msg: + except getopt.GetoptError as msg: at_exit(msg) for o, a in opts: @@ -81,20 +102,20 @@ def main(): if o in ('-h', '--help'): raise Exception() elif o in ('-e', '--error'): - eps = float(a) + eps = float2(a) if eps <= 0: raise Exception() elif o in ('-l', '--limit'): l = int(a) if l < 1: raise Exception() elif o in ('-t', '--target'): try: - t = float(eval(a)) - except Exception, msg: + t = float2(eval(a)) + except Exception as msg: at_exit(msg) if t <= 0: raise Exception() else: raise Exception() - except Exception, msg: + except Exception as msg: at_exit(msg) if t == None or l==None: @@ -108,11 +129,11 @@ def main(): else: err, n, d, niter = find_best_rat_with_err_bound(eps, l, t_frac) n += t_int * d - err = (t - float(n) / d) + err = (t - float2(n) / d) if eps == None: print("target= %f best_rat= %d / %d max_denom= %d err= %g abs_err= %g niter= %d" % (t, n, d, l, err, abs(err), niter)) else: - print("target= %f best_rat= %d / %d max_denom= %d err= %g abs_err= %g abs_err/error= %g niter= %d" % (t, n, d, l, err, abs(err), float(abs(err)) / eps, niter)) + print("target= %f best_rat= %d / %d max_denom= %d err= %g abs_err= %g abs_err/error= %g niter= %d" % (t, n, d, l, err, abs(err), float2(abs(err)) / eps, niter)) if __name__ == '__main__': main() diff --git a/ad_rat_by_farey.py b/ad_rat_by_farey.py index ec307d8..5843a64 100755 --- a/ad_rat_by_farey.py +++ b/ad_rat_by_farey.py @@ -1,4 +1,4 @@ -#!/usr/bin/python +#!/usr/bin/python3 # find the best rational approximation n/d where d <= l to the given # target number t. here the best means the one with the smallest @@ -12,17 +12,34 @@ # naming convention: n=numerator, d=denominator, l=left, r=right # author: ali dasdan +# updated 2024 onwards by Tim Abdiukov. import sys import getopt from math import * +USE_NUMPY_HI_PRECISION = True +if(USE_NUMPY_HI_PRECISION): + try: + import numpy as np + float2 = lambda x: np.longdouble(x) + except ModuleNotFoundError: + print("Please disable USE_NUMPY_HI_PRECISION or install prerequisite,") + print("```") + print("pip install numpy") + print("```") + exit() +else: + print("WARN: High precision NumPy logic is not used") + float2 = lambda x: float2(x) + + def show_usage(): - print sys.argv[0] + " -h/--help [-e/--error=float>=0] -l/--limit=int=>1 -t/--target=float or quoted math expr returning float" + print(sys.argv[0] + " -h/--help [-e/--error=float>=0] -l/--limit=int=>1 -t/--target=float or quoted math expr returning float") def at_exit(msg): if msg != "" and msg != None: - print "Error:", msg + print("Error:", msg) show_usage() sys.exit(0) @@ -55,7 +72,7 @@ def find_best_rat(l, t): # find the mediant med=nm/dm with dm <= l nm, dm = nl + nr, dl + dr if dm > l: break - med = float(nm) / dm + med = float2(nm) / dm # branch based on t's position in nl/dr < med < nr/dr if t_frac == med: @@ -75,8 +92,8 @@ def find_best_rat(l, t): n, d = nm, dm else: # find out the endpoint closest to t_frac - errl = abs(t_frac - float(nl) / dl) - errr = abs(t_frac - float(nr) / dr) + errl = abs(t_frac - float2(nl) / dl) + errr = abs(t_frac - float2(nr) / dr) if errl <= errr: n, d = nl, dl else: @@ -84,7 +101,12 @@ def find_best_rat(l, t): # convert the solution to a rat for t n += t_int * d - err = (t - float(n) / d) + err = (t - float2(n) / d) + + # assert precision wasn't lost + assert(type(err) == type(float2(pi))) + + # return data return err, n, d, niter # this function takes in an error bound err_in, an int limit l, and @@ -113,7 +135,7 @@ def main(): opts, args = getopt.getopt(sys.argv[1:], 'he:l:t:', ['help', 'error=', 'limit=', 'target=']) - except getopt.GetoptError, msg: + except getopt.GetoptError as msg: at_exit(msg) for o, a in opts: @@ -121,20 +143,20 @@ def main(): if o in ('-h', '--help'): raise Exception() elif o in ('-e', '--error'): - eps = float(a) + eps = float2(a) if eps <= 0: raise Exception() elif o in ('-l', '--limit'): l = int(a) if l < 1: raise Exception() elif o in ('-t', '--target'): try: - t = float(eval(a)) - except Exception, msg: + t = float2(eval(a)) + except Exception as msg: at_exit(msg) if t <= 0: raise Exception() else: raise Exception() - except Exception, msg: + except Exception as msg: at_exit(msg) if t == None or l==None: @@ -150,7 +172,7 @@ def main(): if eps == None: print("target= %f best_rat= %d / %d max_denom= %d err= %g abs_err= %g niter= %d" % (t, n, d, l, err, abs(err), niter)) else: - print("target= %f best_rat= %d / %d max_denom= %d err= %g abs_err= %g abs_err/error= %g niter= %d" % (t, n, d, l, err, abs(err), float(abs(err)) / eps, niter)) + print("target= %f best_rat= %d / %d max_denom= %d err= %g abs_err= %g abs_err/error= %g niter= %d" % (t, n, d, l, err, abs(err), float2(abs(err)) / eps, niter)) if __name__ == '__main__': main() diff --git a/ad_rat_by_fast_farey.py b/ad_rat_by_fast_farey.py index 9398915..00526ab 100755 --- a/ad_rat_by_fast_farey.py +++ b/ad_rat_by_fast_farey.py @@ -1,4 +1,4 @@ -#!/usr/bin/python +#!/usr/bin/python3 # find the best rational approximation n/d where d <= l to the given # target number t. here the best means the one with the smallest @@ -13,17 +13,33 @@ # naming convention: n=numerator, d=denominator, l=left, r=right # author: ali dasdan +# updated 2024 onwards by Tim Abdiukov. import sys import getopt from math import * +USE_NUMPY_HI_PRECISION = True +if(USE_NUMPY_HI_PRECISION): + try: + import numpy as np + float2 = lambda x: np.longdouble(x) + except ModuleNotFoundError: + print("Please disable USE_NUMPY_HI_PRECISION or install prerequisite,") + print("```") + print("pip install numpy") + print("```") + exit() +else: + print("WARN: High precision NumPy logic is not used") + float2 = lambda x: float(x) + def show_usage(): - print sys.argv[0] + " -h/--help [-e/--error=float>=0] -l/--limit=int>1 -t/--target=float or quoted math expr returning float" + print(sys.argv[0] + " -h/--help [-e/--error=float>=0] -l/--limit=int>1 -t/--target=float or quoted math expr returning float") def at_exit(msg): if msg != "" and msg != None: - print "Error:", msg + print("Error:", msg) show_usage() sys.exit(0) @@ -61,8 +77,8 @@ def find_best_rat(l, t): # also find the largest bl such that x < (bl * nl + nr) / (bl * dl + dr). ntmp = nl - t_frac * dl dtmp = t_frac * dr - nr - br = int(floor(float(ntmp) / dtmp)) - bl = int(floor(float(dtmp) / ntmp)) + br = int(floor(float2(ntmp) / dtmp)) + bl = int(floor(float2(dtmp) / ntmp)) side = 0 # left:-1, init:0, right:1 if bl == 0: bl = 1 @@ -80,15 +96,15 @@ def find_best_rat(l, t): if dm > l: if side == -1: - br = max(1, int(floor(float(l - dl) / dr))) + br = max(1, int(floor(float2(l - dl) / dr))) elif side == 1: - bl = max(1, int(floor(float(l - dr) / dl))) + bl = max(1, int(floor(float2(l - dr) / dl))) # recalculate nm, dm, and med after the jump nm = nl * bl + nr * br dm = dl * bl + dr * br if dm > l: break - med = float(nm) / dm + med = float2(nm) / dm # branch based on t's position in nl/dr < med < nr/dr if t_frac == med: @@ -108,8 +124,8 @@ def find_best_rat(l, t): n, d = nm, dm else: # find out the endpoint closest to t_frac - errl = abs(t_frac - float(nl) / dl) - errr = abs(t_frac - float(nr) / dr) + errl = abs(t_frac - float2(nl) / dl) + errr = abs(t_frac - float2(nr) / dr) if errl <= errr: n, d = nl, dl else: @@ -117,7 +133,12 @@ def find_best_rat(l, t): # convert the solution to a rat for t n += t_int * d - err = (t - float(n) / d) + err = (t - float2(n) / d) + + # assert precision wasn't lost + assert(type(err) == type(float2(pi))) + + # return data return err, n, d, niter # this function takes in an error bound err_in, an int limit l, and @@ -146,7 +167,7 @@ def main(): opts, args = getopt.getopt(sys.argv[1:], 'he:l:t:', ['help', 'error=', 'limit=', 'target=']) - except getopt.GetoptError, msg: + except getopt.GetoptError as msg: at_exit(msg) for o, a in opts: @@ -154,20 +175,20 @@ def main(): if o in ('-h', '--help'): raise Exception() elif o in ('-e', '--error'): - eps = float(a) + eps = float2(a) if eps <= 0: raise Exception() elif o in ('-l', '--limit'): l = int(a) if l < 1: raise Exception() elif o in ('-t', '--target'): try: - t = float(eval(a)) - except Exception, msg: + t = float2(eval(a)) + except Exception as msg: at_exit(msg) if t <= 0: raise Exception() else: raise Exception() - except Exception, msg: + except Exception as msg: at_exit(msg) if t == None or l==None: @@ -183,7 +204,7 @@ def main(): if eps == None: print("target= %f best_rat= %d / %d max_denom= %d err= %g abs_err= %g niter= %d" % (t, n, d, l, err, abs(err), niter)) else: - print("target= %f best_rat= %d / %d max_denom= %d err= %g abs_err= %g abs_err/error= %g niter= %d" % (t, n, d, l, err, abs(err), float(abs(err)) / eps, niter)) + print("target= %f best_rat= %d / %d max_denom= %d err= %g abs_err= %g abs_err/error= %g niter= %d" % (t, n, d, l, err, abs(err), float2(abs(err)) / eps, niter)) if __name__ == '__main__': main()