chiark / gitweb /
1 #!/usr/bin/env python
3 import sys, subprocess, math, getopt
4 from fractions import Fraction
6 import bounds
8 verbose = False
9 def vprint(*args):
10     # In verbose mode, print diagnostics.
11     if verbose:
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
20 specified."""
21     if ret is None:
22         ret = [] # new array
23     maxfrag = min(maxfrag, n) if maxfrag is not None else n
24     if n == 0:
25         ret.append(prefix)
26     else:
27         for frag in xrange(minfrag, maxfrag+1):
28             find_partitions(n-frag, frag, maxfrag, prefix+(frag,), ret)
29     return 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
37     # the other.
38     while True:
39         vprint("Partitions of %d:" % n)
40         for np in nparts:
41             vprint("   ", np)
42         vprint("Partitions of %d:" % m)
43         for mp in mparts:
44             vprint("   ", mp)
46         oldlens = len(nparts), len(mparts)
47         nlengths = set(sum(nparts, ()))
48         mlengths = set(sum(mparts, ()))
49         new_nparts = []
50         for np in nparts:
51             s = set([k for k in np if k not in mlengths])
52             if len(s) == 0:
53                 new_nparts.append(np)
54             else:
55                 vprint("Winnowing %s (can't use %s)" % (
56                     np, ", ".join(map(str,sorted(s)))))
57         new_mparts = []
58         for mp in mparts:
59             s = set([k for k in mp if k not in nlengths])
60             if len(s) == 0:
61                 new_mparts.append(mp)
62             else:
63                 vprint("Winnowing %s (can't use %s)" % (
64                     mp, ", ".join(map(str,sorted(s)))))
65         nparts = new_nparts
66         mparts = new_mparts
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.")
72         return None
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
76     # ILP solver :-)
77     ilp_file = ""
78     nvarnames = {}
79     for np in nparts:
80         nvarnames[np] = "_".join(["n"] + map(str,np))
81     mvarnames = {}
82     for mp in mparts:
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)])
87     for var in 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
92     for k in nlengths:
93         ns = []
94         for np in nparts:
95             count = len([x for x in np if x == k])
96             if count > 0:
97                 ns.append((count, nvarnames[np]))
98         ms = []
99         for mp in mparts:
100             count = len([x for x in mp if x == k])
101             if count > 0:
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)
111     if p.wait() != 0:
112         vprint("ILP solver failed")
113         return None
114     else:
115         ncounts = {}
116         mcounts = {}
117         for line in stdout.splitlines():
118             words = line.split()
119             if len(words) == 0:
120                 pass # rule out for future elifs
121             elif words[:2] == "n_" and words != "0":
122                 ncounts[tuple(map(int, words[2:].split("_")))] = int(words)
123             elif words[:2] == "m_" and words != "0":
124                 mcounts[tuple(map(int, words[2:].split("_")))] = int(words)
125         return ncounts, mcounts
127 def search(n, m):
128     if n % m == 0:
129         # Trivial special case.
130         return (m, 1, {(m,):n}, {(m,)*(n/m):m})
131     best = (0,)
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*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
138                 break
139         if best == bound:
140             break # terminate early if we know it's the best answer
141     return best
143 def search_and_report(n, m):
144     best = search(n, m)
145     d = best
146     print "%d into %d best min fragment found: %s" % (n, m, best)
147     print "  Cut up %d sticks of length %d like this:" % (n, m)
148     for row, count in sorted(best.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.items(), reverse=True):
152         print "    %d x (%s)" % (count, " + ".join([str(Fraction(k,d)) for k in col]))
154 def main():
155     global verbose
156     opts, args = getopt.gnu_getopt(sys.argv[1:], "v")
157     for opt, val in opts:
158         if opt == "-v":
159             verbose = True
160         else:
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__":
166     main()