#!/usr/bin/python
# -*- coding: utf-8 -*-
"""Subdividing an interval according to the golden ratio.
Uses a binary tree to represent divisions of the interval to see if
the multiples of the golden ratio, mod 1, always fall in the largest
undivided interval in small experiments. Up to 32000, they do seem to.
This program is somewhat inadequate to answer the question, because of
the problem of equal subintervals. The first apparent counterexample
is that phi * 5 % 1 is about 0.0902, which is in the interval
(0, phi*2%1 = 0.2361), leaving untouched the interval
(phi*1%1 = 0.6180, phi*3%1=0.8541) of length, uh, 0.2361. In fact
these intervals are exactly equal in size; phi*3 - phi = phi*2,
exactly, even outside the modulo-1 system.
However, this program is not adequate to deal with multiple “largest”
intervals, especially in the presence of arithmetic errors. The best
it can do is to compare the length of the interval being cut with the
length of the largest interval. Doing so, it doesn’t find a
counterexample with an error as large as 1e-12 within the first 4000
multiples of the golden ratio, or as large as 1e-11 within the first
32000.
One potential improvement to the program would be to balance the tree
upon point insertions. However, for the specific case of multiples of
the golden ratio, the tree remains reasonably well balanced without
any explicit balancing logic; it reaches a depth of only 17 for 4000
nodes, where 12 would be the theoretical minimum and 4001 the worst
case. At 32000 nodes it reaches a depth of 22.
I think a 2-3-4 tree would probably be the simplest variant to
implement.
"""
import math, sys
class Subdivision:
"""Represents an interval possibly subdivided into intervals.
Can efficiently perform the following operations:
- Return a subdivision with an additional point subdividing one of
the intervals;
- Return its longest interval;
- Return its ends;
- Return its length;
- Return all its division points;
- Return the interval within it containing a given point.
At least if you get lucky.
"""
def length(self):
return self.right_end() - self.left_end()
def points(self):
return [self.left_end()] + self.cuts() + [self.right_end()]
def __repr__(self):
return 'subdivision(%r)' % self.points()
class Interval(Subdivision):
"""Represents an interval not subdivided into intervals.
"""
def __init__(self, left, right):
self.left, self.right = left, right
def left_end(self): return self.left
def right_end(self): return self.right
def longest(self):
return self
def cut(self, where):
return Broken(Interval(self.left, where), Interval(where, self.right))
def find(self, point):
assert self.left <= point < self.right
return self
def cuts(self):
return []
def depth(self):
return 1
class Broken(Subdivision):
"""Represents an interval divided into intervals."""
def __init__(self, lefthalf, righthalf):
assert lefthalf.right_end() == righthalf.left_end()
self.lefthalf, self.righthalf = lefthalf, righthalf
ll = lefthalf.longest()
rl = righthalf.longest()
if ll.length() > rl.length():
self.longest_interval = ll
else:
self.longest_interval = rl
self.stored_depth = 1 + max(lefthalf.depth(), righthalf.depth())
def left_end(self): return self.lefthalf.left_end()
def right_end(self): return self.righthalf.right_end()
def longest(self): return self.longest_interval
def cut(self, where):
assert self.left_end() <= where < self.right_end()
if where < self.lefthalf.right_end():
return Broken(self.lefthalf.cut(where), self.righthalf)
else:
return Broken(self.lefthalf, self.righthalf.cut(where))
def find(self, point):
assert self.left_end() <= point < self.right_end()
if point < self.lefthalf.right_end():
return self.lefthalf.find(point)
else:
return self.righthalf.find(point)
def cuts(self):
return self.lefthalf.cuts() + [self.lefthalf.right_end()] + self.righthalf.cuts()
def depth(self):
return self.stored_depth
def phidemo(n=40):
phi = (1+math.sqrt(5))/2
interval = Interval(0, 1)
max_error = 0
for ii in range(1, n):
print "longest:", interval.longest()
cut = (phi * ii) % 1
print "cutting at phi * %s = %s" % (ii, cut)
found = interval.find(cut)
if found != interval.longest():
error = interval.longest().length() - found.length()
print "shorter? by", error
if error > max_error: max_error = error
interval = interval.cut(cut)
print "finally:", interval
print "depth:", interval.depth()
print "max_error:", max_error
if __name__ == '__main__':
if len(sys.argv) > 1:
phidemo(int(sys.argv[1]))
else:
phidemo()