This Python program, inspired by an ASCII-art generator by Darius
Bacon, generates a PostScript graph showing how the powers of φ, the
golden ratio (1+√5)/2, are distributed such that the fractional part
of each power falls into the largest remaining hole between the
fractional parts of the earlier powers — although in some cases there
is more than one such “largest remaining hole”.
#!/usr/bin/python
from __future__ import division
phi = (1 + 5**.5) / 2
n = 200
print '%!'
print '36 dup translate 540 dup scale 0 0 moveto 0 1 lineto 1 1 lineto 1 0
lineto closepath gsave 0.9 setgray fill grestore clip 0 setlinewidth'
for i in range(1, n):
print 'newpath', i * phi % 1, i/n, 'moveto 0 -2 rlineto stroke'
print 'showpage'
(End of `golden.py`.)
Example usage:
: kra...@inexorable:~/devel/inexorable-misc ; ./golden.py > golden.ps
: kra...@inexorable:~/devel/inexorable-misc ; evince golden.ps
(`evince` is GNOME’s viewer for PostScript files.)
The X-axis inside the gray square is the interval from 0 to 1 within
which the fractional parts fall; the Y-axis is the power to which φ is
raised, starting from 0 at the top of the square and continuing to n =
200 at the bottom.
Unfortunately, although the result is visually fairly striking, it
isn’t very convincing. So I wrote an overly elaborate program to
confirm this property experimentally on some examples, although of
course that can be misleading. This program is called `intervals.py`.
#!/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()
(End of `intervals.py`.)
After writing 150 lines of code to test this property, I tried
rewriting it in Lua to see how it came out. It came out 4× as fast,
and roughly equally clear, I think. The log of the problems I ran into
is perhaps interesting when considering human interface questions
related to programming languages. This is `intervals.lua`:
#!/usr/bin/lua
-- Analyze the distribution of the fractional parts of the multiples
-- of the golden ratio.
-- The original was in Python; I rewrote it in Lua to see if LuaJIT
-- could run it faster than Python. The answer:
-- YES, running it with 32000 points takes
-- 56.5 seconds with CPython 2.6.2,
-- 34.3 seconds with Lua 5.1.4, and
-- 15.8 seconds with LuaJIT 2.0.0-beta2.
-- A 4x speedup isn’t bad, but 40x is probably doable.
-- Log of my problems, as a Lua novice:
-- * Didn’t have lua-mode installed for Emacs.
-- * Tried to use “load('intervals.lua')” instead of
-- “loadfile('intervals.lua')()” to test in the REPL.
-- * Forgot to return a value from the “interval” function.
-- * Couldn’t remember that “does not equal” is spelled “~=”.
-- * Couldn't remember that “not” is spelled “not”, not “!”.
-- * Tried to use “loadfile('intervals.lua')” instead of
-- “print(loadfile('intervals.lua'))” to see syntax errors.
-- * Left out parens on a method call: “if foo:bar() < baz:bar then”,
-- getting “function arguments expected near 'then'”.
-- * Tried to declare a local variable with “var cuts”
-- instead of “local cuts”, getting “'=' expected near 'cuts'”.
-- * Forgot that I had translated “length” as a function instead of
-- a method, getting “attempt to call method 'length' (a nil value)”.
-- * Was initializing an instance variable as “rv.longest = foo” instead
-- of “rv.longest_interval = foo”, so my "longest" method returned nil,
-- so I was getting “attempt to index local 'i' (a nil value)” in a
-- totally different function.
function concat(t, items)
for _, item in ipairs(items) do table.insert(t, item) end
end
function interval_length(i) return i:right_end() - i:left_end() end
function points(i)
-- I was concatenating tables with unpack() and then I got this with
32000 items:
-- /home/kragen/pkgs/luajit: intervals.lua:9: too many results to unpack
-- So now I’m using a concat() function.
local p = {i:left_end()}
concat(p, i:cuts())
table.insert(p, i:right_end())
return p
end
function repr(i) return 'subdivision(['..table.concat(points(i), ',
')..'])' end
-- a subdivision not subdivided into intervals; a leaf node.
function interval(left, right)
local rv = {left = left, right = right}
function rv:left_end() return self.left end
function rv:right_end() return self.right end
function rv:longest() return self end
function rv:cut(where) return broken(interval(self.left, where),
interval(where, self.right)) end
function rv:find(point)
assert(self.left <= point)
assert(point < self.right)
return self
end
function rv:cuts() return {} end
function rv:depth() return 1 end
return rv
end
-- a subdivision broken down further into intervals
function broken(lefthalf, righthalf)
assert(lefthalf:right_end() == righthalf:left_end())
local rv = { lefthalf = lefthalf, righthalf = righthalf }
local ll = lefthalf:longest()
local rl = righthalf:longest()
if interval_length(ll) > interval_length(rl)
then rv.longest_interval = ll
else rv.longest_interval = rl
end
rv.stored_depth = 1 + math.max(lefthalf:depth(), righthalf:depth())
function rv:left_end() return self.lefthalf:left_end() end
function rv:right_end() return self.righthalf:right_end() end
function rv:longest() return self.longest_interval end
function rv:cut(where)
assert(self:left_end() <= where)
assert(where < self:right_end())
if where < self.lefthalf:right_end() then
return broken(self.lefthalf:cut(where), self.righthalf)
else
return broken(self.lefthalf, self.righthalf:cut(where))
end
end
function rv:find(point)
assert(self:left_end() <= point)
assert(point < self:right_end())
if point < self.lefthalf:right_end() then
return self.lefthalf:find(point)
else
return self.righthalf:find(point)
end
end
function rv:cuts()
local cuts = self.lefthalf:cuts()
table.insert(cuts, self.lefthalf:right_end())
concat(cuts, self.righthalf:cuts())
return cuts
end
function rv:depth() return self.stored_depth end
return rv
end
function phidemo(n)
if not n then n = 40 end
local phi = (1+math.sqrt(5))/2
local iv = interval(0, 1)
local max_error = 0
for ii = 1, n do
print("longest: "..repr(iv:longest()))
local cut = (phi * ii) % 1
print("cutting at phi * "..ii.." = "..cut)
local found = iv:find(cut)
if found ~= iv:longest() then
local err = interval_length(iv:longest()) - interval_length(found)
print("shorter? by "..err)
if err > max_error then max_error = err end
end
iv = iv:cut(cut)
end
print("finally: "..repr(iv))
print("depth: "..iv:depth())
print("max_error: "..max_error)
end
if #{...} > 0 then phidemo(...) end
(End of `intervals.lua`.)
Some of this software is available via
git clone http://canonical.org/~kragen/sw/inexorable-misc.git
(or in <http://canonical.org/~kragen/sw/inexorable-misc>) in the file
`golden.py`.
Like everything else posted to kragen-hacks without a notice to the
contrary, this software is in the public domain.
--
To unsubscribe: http://lists.canonical.org/mailman/listinfo/kragen-hacks
