Skip to content
48 changes: 32 additions & 16 deletions ad_gcd_lcm.py
Original file line number Diff line number Diff line change
@@ -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
Expand All @@ -7,19 +7,35 @@
# rational numbers.

# author: ali dasdan
# updated 2024 onwards by Tim Abdiukov.

import sys
import getopt
import re
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)

Expand Down Expand Up @@ -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:
Expand All @@ -96,15 +112,15 @@ def main():
if l < 1: raise Exception()
else:
raise Exception()
except Exception, msg:
except Exception as msg:
at_exit(msg)

if s == None:
at_exit(None)

# 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")

Expand All @@ -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 = []
Expand All @@ -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)
Expand All @@ -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()
Expand Down
56 changes: 38 additions & 18 deletions ad_rat_by_cont_frac.py
Original file line number Diff line number Diff line change
@@ -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
Expand All @@ -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)

Expand Down Expand Up @@ -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
Expand All @@ -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,
Expand Down Expand Up @@ -110,28 +130,28 @@ 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:
try:
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:
Expand All @@ -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()
Expand Down
43 changes: 32 additions & 11 deletions ad_rat_by_exhaustive.py
Original file line number Diff line number Diff line change
@@ -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
Expand All @@ -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)

Expand All @@ -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
Expand Down Expand Up @@ -73,28 +94,28 @@ 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:
try:
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:
Expand All @@ -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()
Expand Down
Loading