#!/usr/bin/python
# -*- coding: utf-8 -*-
"""Optimal addition chains.
In Stepanov's Three Algorithmic Journeys, he poses the problem
(Exercise 1.4 on p. 18, the 27th page, of draft 0.3) of finding
optimal addition chains for multiplying by N, for any N under 100. He
illustrates the problem by providing a multiply-by-15 algorithm:
b = a + a + a
c = b + b
c + c + b
This uses five additions, while the standard Rhind-papyrus
multiplication algorithm (aka mediation and duplation, Ethiopian
multiplication, or Russian peasant multiplication, and used more or
less universally in modern computers) uses six:
d = a + a + a
e = d + d + a
e + e + a
This seems like a problem that's eminently suited to dynamic
programming. But it turns out that it doesn't yield quite so easily;
in the above case, c is 6a, and is achieved in only 4 additions, but
in a way that, it turns out, pleasantly computes 2a and 3a rather than
the usual 2a and 4a as its components. Computing 6a as 2a + 4a rather
than 3a + 3a is just fine as an optimal addition chain for 6a, but the
second one is better when you want to compute 15a, because if you
computed 12a in a way using four additions that gives you 3a for free,
then you only need one more addition to get to 15a.
Worse, the best addition chain to use to compute some subpart of your
final sum may not even be optimal for that subpart, I think. Two
suboptimal chains for different subparts might have enough overlap
that they turn out to jointly comprise the optimal addition chain.
An "addition chain for K" is some set of numbers such that:
- K is in the set.
- Every number in the set, except 1, is the sum of two other members
of the set.
This represents a procedure for calculating K times a number x:
calculate each of the multiples nK, except for 1K, by adding two of
the multiples you've already calculated. A smallest such set
represents an optimal addition chain for K; if it has N members, it
represents a chain needing N-1 additions.
A brute-force approach needing no worse than 2**(K/2) time is simply
to search. For any partial addition chain, you can search through all
the possible ways to complete it, recursively, starting by trying to
complete {1, K-1}, unless it's already complete, and then trying
{2, K-2}, and so on, until those things cross. Then you just take the
shortest one of those K/2 possibilities and add K to it.
In that form, it's pretty much exactly 2**(K/2) (it needs 2872
recursions to find the optimal chain for 15), but with memoization,
it's probably better (needing only 302). But with memoization, it
needs potentially 2**(K/2) space.
An optimal chain for 63, needing 8 additions, is [1, 2, 3, 6, 12, 48,
24, 60, 63]. The Rhind algorithm instead uses the chain [1, 2, 4, 8,
16, 32, 48, 56, 60, 62, 63], needing 10 additions. I expect these
cases of 2^N-1 to be a kind of worst case for the Rhind algorithm,
while 2^N should be the best case.
This program is able to calculate optimal addition chains for numbers
up to 40 almost instantly, and up to 75 in about a minute. It gets to
90 in about 5 minutes, and 100 in 143 minutes (just over 17 minutes of
CPU time).
What is this useful for? Probably nothing, in practice, because it
counts doubling a number as being just as expensive as adding two
arbitrary numbers. But, in circuitry, doubling a number — or actually
multiplying it by any power of 2 — is just a matter of hooking up the
wires a little to the left. And for manual decimal calculation,
multiplying 2830 by 63 is probably easier by using the multiplication
table:
2830
* 63
----
8490
+ 16980
--------
178290
than by using the addition chain:
2830 + 2830 = 5660 (2x)
5660 + 2830 = 8490 (3x)
8490 + 8490 = 16980 (6x)
16980 + 16980 = 33960 (12x)
33960 + 33960 = 67920 (24x)
67920 + 67920 = 135840 (48x)
135840 + 33960 = 169800 (60x) (note that in base ten this is 6x shifted 1)
169800 + 8490 = 178290 (63x)
But I thought it was an interesting problem, anyway.
"""
import sys
def main(n):
for ii in range(1, n):
print sorted(optimal_chain(set([ii])))
chains = {}
def optimal_chain(start):
start = frozenset(start)
# Avoid uselessly duplicating work (and storage); 1x will always
# be known, so there's no use considering, and memoizing, the case
# where we don't compute it.
if 1 not in start:
start = frozenset(start | set([1]))
#print "trying to reach", start
if start in chains:
return chains[start]
if is_chain(start):
# Empirically, memoizing this case to avoid calls to is_chain
# slows the program down.
return start
K = max(start)
candidates = []
reduced = start - set([K])
for ii in range(1, K):
Ki = K - ii
if Ki < ii:
break
candidates.append(optimal_chain(reduced | set([Ki, ii])))
result = min(candidates, key=len) | set([K])
chains[start] = result
return result
def is_chain(xs):
# This is an O(N^2) algorithm, and since N gets up to about 10,
# you might conceivably be able to get a speedup of up to around
# 10 by optimizing it.
return all(c == 1 or any(c - a in xs for a in xs) for c in xs)
if __name__ == '__main__':
main(int(sys.argv[1]) if len(sys.argv) > 1 else 101)