3 import sys, subprocess, math, getopt
4 from fractions import Fraction
10 # In verbose mode, print diagnostics.
12 sys.stdout.write(" ".join(map(str,args)) + "\n")
14 def find_partitions(n, minfrag, maxfrag=None, prefix=(), ret=None):
15 """Find all partitions of n into integer pieces at least minfrag.
17 Returns a list of tuples.
19 For recursive calls: appends to an existing 'ret' if none is
23 maxfrag = min(maxfrag, n) if maxfrag is not None else n
27 for frag in xrange(minfrag, maxfrag+1):
28 find_partitions(n-frag, frag, maxfrag, prefix+(frag,), ret)
31 def try_one_minfrag(n, m, minfrag, d):
32 # Find all possible partitions of the two stick lengths.
33 nparts = find_partitions(n*d, minfrag, m*d)
34 mparts = find_partitions(m*d, minfrag)
36 # Winnow by discarding any partition using a length not present in
39 vprint("Partitions of %d:" % n)
42 vprint("Partitions of %d:" % m)
46 oldlens = len(nparts), len(mparts)
47 nlengths = set(sum(nparts, ()))
48 mlengths = set(sum(mparts, ()))
51 s = set([k for k in np if k not in mlengths])
55 vprint("Winnowing %s (can't use %s)" % (
56 np, ", ".join(map(str,sorted(s)))))
59 s = set([k for k in mp if k not in nlengths])
63 vprint("Winnowing %s (can't use %s)" % (
64 mp, ", ".join(map(str,sorted(s)))))
67 if oldlens == (len(nparts), len(mparts)):
68 break # we have converged
70 if len(nparts) == 0 or len(mparts) == 0:
71 vprint("No partitions available.")
73 # Now we need to look for an integer occurrence count of each
74 # nparts row, summing to m, and one for each mparts row, summing
75 # to n, with the right constraints. We do this by appealing to an
80 nvarnames[np] = "_".join(["n"] + map(str,np))
83 mvarnames[mp] = "_".join(["m"] + map(str,mp))
84 varlist = sorted(nvarnames.values()) + sorted(mvarnames.values())
85 # Have to try to minimise _something_!
86 ilp_file += "min: %s;\n" % " + ".join(["%d * %s" % (i+1, name) for i, name in enumerate(varlist)])
88 ilp_file += "%s >= 0;\n" % var
89 ilp_file += " + ".join(sorted(nvarnames.values())) + " = %d;\n" % m
90 ilp_file += " + ".join(sorted(mvarnames.values())) + " = %d;\n" % n
91 assert nlengths == mlengths
95 count = len([x for x in np if x == k])
97 ns.append((count, nvarnames[np]))
100 count = len([x for x in mp if x == k])
102 ms.append((count, mvarnames[mp]))
103 ilp_file += " + ".join(["%d * %s" % t for t in ns]) + " = "
104 ilp_file += " + ".join(["%d * %s" % t for t in ms]) + ";\n"
105 for var in sorted(nvarnames.values()) + sorted(mvarnames.values()):
106 ilp_file += "int %s;\n" % var
108 p = subprocess.Popen(["lp_solve", "-lp"], stdin=subprocess.PIPE,
109 stdout=subprocess.PIPE)
110 stdout, stderr = p.communicate(ilp_file)
112 vprint("ILP solver failed")
117 for line in stdout.splitlines():
120 pass # rule out for future elifs
121 elif words[0][:2] == "n_" and words[1] != "0":
122 ncounts[tuple(map(int, words[0][2:].split("_")))] = int(words[1])
123 elif words[0][:2] == "m_" and words[1] != "0":
124 mcounts[tuple(map(int, words[0][2:].split("_")))] = int(words[1])
125 return ncounts, mcounts
129 # Trivial special case.
130 return (m, 1, {(m,):n}, {(m,)*(n/m):m})
132 bound, _ = bounds.upper_bound(n, m)
133 for d in xrange(1, n+1):
134 for k in xrange(int(bound*d), int(math.ceil(best[0]*d))-1, -1):
135 result = try_one_minfrag(n, m, k, d)
136 if result is not None:
137 best = (Fraction(k, d), d) + result
140 break # terminate early if we know it's the best answer
143 def search_and_report(n, m):
146 print "%d into %d best min fragment found: %s" % (n, m, best[0])
147 print " Cut up %d sticks of length %d like this:" % (n, m)
148 for row, count in sorted(best[3].items(), reverse=True):
149 print " %d x (%s)" % (count, " + ".join([str(Fraction(k,d)) for k in row]))
150 print " Reassemble as %d sticks of length %d like this:" % (m, n)
151 for col, count in sorted(best[2].items(), reverse=True):
152 print " %d x (%s)" % (count, " + ".join([str(Fraction(k,d)) for k in col]))
156 opts, args = getopt.gnu_getopt(sys.argv[1:], "v")
157 for opt, val in opts:
161 assert False, "unrecognised option '%s'" % opt
162 m, n = sorted(map(int,args[:2]))
163 search_and_report(n, m)
165 if __name__ == "__main__":