The Lumber Room

"Consign them to dust and damp by way of preserving them"

Posts Tagged ‘programming

Some playing with Python

with 2 comments

A long time ago, Diophantus (sort of) discussed integer solutions to the equation

\displaystyle  x^2 + y^2 = z^2

(solutions to this equation are called Pythagorean triples).

Centuries later, in 1637, Fermat made a conjecture (now called Fermat’s Last Theorem, not because he uttered it in his dying breath, but because it was the last one to be proved — in ~1995) that

\displaystyle x^n + y^n = z^n

has no positive integer solutions for n \ge 3. In other words, his conjecture was that none of the following equations has a solution:

\displaystyle x^3 + y^3 = z^3

\displaystyle x^4 + y^4 = z^4

\displaystyle x^5 + y^5 = z^5

\displaystyle x^6 + y^6 = z^6

… and so on. An nth power cannot be partitioned into two nth powers.

About a century later, Euler proved the n = 3 case of Fermat’s conjecture, but generalized it in a different direction: he conjectured in 1769 that an nth power cannot be partitioned into fewer than n nth powers, namely

\displaystyle z^n = \sum_{i = 1}^k x_i^k

has no solutions with k < n. So his conjecture was that (among others) none of the following equations has a solution:

\displaystyle z^3 = a^3 + b^3

\displaystyle z^4 = a^4 + b^4 + c^4

\displaystyle z^5 = a^5 + b^5 + c^5 + d^5

\displaystyle z^6 = a^6 + b^6 + c^6 + d^6 + e^6

… and so on.

This conjecture stood for about two centuries, until abruptly it was found to be false, by Lander and Parkin who in 1966 simply did a direct search on the fastest (super)computer at the time, and found this counterexample:

\displaystyle 27^5 + 84^5 + 110^5 + 133^5 = 144^5

(It is still one of only three examples known, according to Wikipedia.)

Now, how might you find this solution on a computer today?

In his wonderful (as always) post at bit-player, Brian Hayes showed the following code:

import itertools as it

def four_fifths(n):
    """Return smallest positive integers ((a,b,c,d),e) such that
       a^5 + b^5 + c^5 + d^5 = e^5; if no such tuple exists
       with e < n, return the string 'Failed'."""
    fifths = [x**5 for x in range(n)]
    combos = it.combinations_with_replacement(range(1,n), 4)
    while True:
        try:
            cc = combos.next()
            cc_sum = sum([fifths[i] for i in cc])
            if cc_sum in fifths:
                return(cc, fifths.index(cc_sum))
        except StopIteration:
            return('Failed')

to which, if you add (say) print four_fifths(150) and run it, it returns the correct answer in about 47 seconds on my laptop.

The if cc_sum in fifths: line inside the loop is an O(n) cost each time it’s run, so with a simple improvement to the code (using a set instead) and rewriting it a bit, we can write the following full program:

import itertools
 
def find_counterexample(n):
  fifth_powers = [x**5 for x in range(n)]
  fifth_powers_set = set(fifth_powers)
  for xs in itertools.combinations_with_replacement(range(1, n), 4):
    xs_sum = sum([fifth_powers[i] for i in xs])
    if xs_sum in fifth_powers_set:
      return (xs, fifth_powers.index(xs_sum))
  return 'Failed'

print find_counterexample(150)

which finishes in about 8.5 seconds.

Great!

But there’s something unsatisfying about this solution, which is that it assumes there’s a solution with all four numbers on the LHS less than 150. After all, changing the function invocation to find_counterexample(145) makes it run a second faster even, but how could we know to do without already knowing the solution? Besides, we don’t have a fixed 8- or 10-second budget; what we’d really like is a program that keeps searching till it finds a solution or we abort it (or it runs out of memory or something), with no other fixed termination condition.

The above program used the given “n” as an upper bound to generate the combinations of 4 numbers; is there a way to generate all combinations when we don’t know an upper bound on them?

Yes! One of the things I learned from Knuth volume 4 is that if you simply write down each combination in descending order and order them lexicographically, the combinations you get for each upper bound are a prefix of the list of the next bigger one, i.e., for any upper bound, all the combinations form a prefix of the same infinite list, which starts as follows (line breaks for clarity):

1111,
2111, 2211, 2221, 2222,
3111, 3211, 3221, 3222, 3311, 3321, 3322, 3331, 3332, 3333,
4111, ...
     ... 9541, 9542, 9543, 9544, 9551, ... 9555, 9611, ...

There doesn’t seem to be a library function in Python to generate these though, so we can write our own. If we stare at the above list, we can figure out how to generate the next combination from a given one:

  1. Walk backwards from the end, till you reach the beginning or find one that’s less than the previous element.
  2. Increase that element, set all the following ones to 1s, and continue.

We could write, say, the following code for it:

def all_combinations(r):
  xs = [1] * r
  while True:
    yield xs
    for i in range(r - 1, 0, -1):
      if xs[i] < xs[i - 1]:
        break
    else:
      i = 0
    xs[i] += 1
    xs[i + 1:] = [1] * (r - i - 1)

(The else block on a for loop is an interesting Python feature: it is executed if the loop wasn’t terminated with break.) We could even hard-code the r=4 case, as we’ll see later below.

For testing whether a given number is a fifth power, we can no longer simply lookup in a fixed precomputed set. We can do a binary search instead:

def is_fifth_power(n):
  assert n > 0
  lo = 0
  hi = n
  # Invariant: lo^5 < n <= hi^5
  while hi - lo > 1:
    mid = lo + (hi - lo) / 2
    if mid ** 5 < n:
      lo = mid
    else:
      hi = mid
  return hi ** 5 == n

but it turns out that this is slower than one based on looking up in a growing set (as below).

Putting everything together, we can write the following (very C-like) code:

largest_known_fifth_power = (0, 0)
known_fifth_powers = set()
def is_fifth_power(n):
  global largest_known_fifth_power
  while n > largest_known_fifth_power[0]:
    m = largest_known_fifth_power[1] + 1
    m5 = m ** 5
    largest_known_fifth_power = (m5, m)
    known_fifth_powers.add(m5)
  return n in known_fifth_powers

def fournums_with_replacement():
  (x0, x1, x2, x3) = (1, 1, 1, 1)
  while True:
    yield (x0, x1, x2, x3)
    if x3 < x2:
      x3 += 1
      continue
    x3 = 1
    if x2 < x1:
      x2 += 1
      continue
    x2 = 1
    if x1 < x0:
      x1 += 1
      continue
    x1 = 1
    x0 += 1
    continue

if __name__ == '__main__':
  tried = 0
  for get in fournums_with_replacement():
    tried += 1
    if (tried % 1000000 == 0):
      print tried, 'Trying:', get
    rhs = get[0]**5 + get[1]**5 + get[2]**5 + get[3]**5
    if is_fifth_power(rhs):
      print 'Found:', get, rhs
      break

which is both longer and slower (takes about 20 seconds) than the original program, but at least we have the satisfaction that it doesn’t depend on any externally known upper bound.

I originally started writing this post because I wanted to describe some experiments I did with profiling, but it’s late and I’m sleepy so I’ll just mention it.

python -m cProfile euler_conjecture.py

will print relevant output in the terminal:

         26916504 function calls in 26.991 seconds

   Ordered by: standard name

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1   18.555   18.555   26.991   26.991 euler_conjecture.py:1()
 13458164    4.145    0.000    4.145    0.000 euler_conjecture.py:12(fournums_with_replacement)
 13458163    4.292    0.000    4.292    0.000 euler_conjecture.py:3(is_fifth_power)
      175    0.000    0.000    0.000    0.000 {method 'add' of 'set' objects}
        1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}

Another way to view the same thing is to write the profile output to a file and read it with cprofilev:

python -m cProfile -o euler_profile.out euler_conjecture.py
cprofilev euler_profile.out

and visit http://localhost:4000 to view it.

Of course, simply translating this code to C++ makes it run much faster:

#include <array>
#include <iostream>
#include <map>
#include <utility>

typedef long long Int;
constexpr Int fifth_power(Int x) { return x * x * x * x * x; }

std::map<Int, int> known_fifth_powers = {{0, 0}};
bool is_fifth_power(Int n) {
  while (n > known_fifth_powers.rbegin()->first) {
    int m = known_fifth_powers.rbegin()->second  + 1;
    known_fifth_powers[fifth_power(m)] = m;
  }
  return known_fifth_powers.count(n);
}

std::array<Int, 4> four_nums() {
  static std::array<Int, 4> x = {1, 1, 1, 0};
  int i = 3;
  while (i > 0 && x[i] == x[i - 1]) --i;
  x[i] += 1;
  while (++i < 4) x[i] = 1;
  return x;
}

std::ostream& operator<<(std::ostream& os, std::array<Int, 4> x) {
  os << "(" << x[0] << ", " << x[1] << ", " << x[2] << ", " << x[3] << ")";
  return os;
}

int main() {
  while (true) {
    std::array<Int, 4> get = four_nums();
    Int rhs = fifth_power(get[0]) + fifth_power(get[1]) + fifth_power(get[2]) + fifth_power(get[3]);
    if (is_fifth_power(rhs)) {
      std::cout << "Found: " << get << " " << known_fifth_powers[rhs] << std::endl;
      break;
    }
  }
}

and


clang++ -std=c++11 euler_conjecture.cc && time ./a.out

runs in 2.43s, or 0.36s if compiled with -O2.

But I don’t have a satisfactory answer to how to make our Python program which takes 20 seconds as fast as the 8.5-second known-upper-bound version.


Edit [2015-05-08]: I wrote some benchmarking code to compare all the different “combination” functions.

import itertools

# Copied from the Python documentation
def itertools_equivalent(iterable, r):
    pool = tuple(iterable)
    n = len(pool)
    if not n and r:
        return
    indices = [0] * r
    yield tuple(pool[i] for i in indices)
    while True:
        for i in reversed(range(r)):
            if indices[i] != n - 1:
                break
        else:
            return
        indices[i:] = [indices[i] + 1] * (r - i)
        yield tuple(pool[i] for i in indices)

# Above function, specialized to first argument being range(1, n)
def itertools_equivalent_specialized(n, r):
  indices = [1] * r
  yield indices
  while True:
    for i in reversed(range(r)):
      if indices[i] != n - 1:
        break
    else:
      return
    indices[i:] = [indices[i] + 1] * (r - i)
    yield indices

# Function to generate all combinations of 4 elements
def all_combinations_pythonic(r):
  xs = [1] * r
  while True:
    yield xs
    for i in range(r - 1, 0, -1):
      if xs[i] < xs[i - 1]:
        break
    else:
      i = 0
    xs[i] += 1
    xs[i + 1:] = [1] * (r - i - 1)

# Above function, written in a more explicit C-like way
def all_combinations_clike(r):
  xs = [1] * r
  while True:
    yield xs
    i = r - 1
    while i > 0 and xs[i] == xs[i - 1]:
      i -= 1
    xs[i] += 1
    while i < r - 1:
      i += 1
      xs[i] = 1

# Above two functions, specialized to r = 4, using tuple over list.
def fournums():
  (x0, x1, x2, x3) = (1, 1, 1, 1)
  while True:
    yield (x0, x1, x2, x3)
    if x3 < x2:
      x3 += 1
      continue
    x3 = 1
    if x2 < x1:
      x2 += 1
      continue
    x2 = 1
    if x1 < x0:
      x1 += 1
      continue
    x1 = 1
    x0 += 1
    continue

# Benchmarks for all functions defined above (and the library function)
def benchmark_itertools(n):
  for xs in itertools.combinations_with_replacement(range(1, n), 4):
    if xs[0] >= n:
      break
def benchmark_itertools_try(n):
  combinations = itertools.combinations_with_replacement(range(1, n), 4)
  while True:
    try:
      xs = combinations.next()
      if xs[0] >= n:
        break
    except StopIteration:
      return
def benchmark_itertools_equivalent(n):
  for xs in itertools_equivalent(range(1, n), 4):
    if xs[0] >= n:
      break
def benchmark_itertools_equivalent_specialized(n):
  for xs in itertools_equivalent_specialized(n, 4):
    if xs[0] >= n:
      break
def benchmark_all_combinations_pythonic(n):
  for xs in all_combinations_pythonic(4):
    if xs[0] >= n:
      break
def benchmark_all_combinations_clike(n):
  for xs in all_combinations_clike(4):
    if xs[0] >= n:
      break
def benchmark_fournums(n):
  for xs in fournums():
    if xs[0] >= n:
      break

if __name__ == '__main__':
  benchmark_itertools(150)
  benchmark_itertools_try(150)
  benchmark_itertools_equivalent(150)
  benchmark_itertools_equivalent_specialized(150)
  benchmark_all_combinations_pythonic(150)
  benchmark_all_combinations_clike(150)
  benchmark_fournums(150)

As you can see, I chose inside the benchmarking function the same statement that would cause all_combinations to terminate, and have no effect for the other combination functions.
When run with

python -m cProfile benchmark_combinations.py

the results include:

    2.817 benchmark_combinations.py:80(benchmark_itertools)
    8.583 benchmark_combinations.py:84(benchmark_itertools_try)
  126.980 benchmark_combinations.py:93(benchmark_itertools_equivalent)
   46.635 benchmark_combinations.py:97(benchmark_itertools_equivalent_specialized)
   44.032 benchmark_combinations.py:101(benchmark_all_combinations_pythonic)
   18.049 benchmark_combinations.py:105(benchmark_all_combinations_clike)
   10.923 benchmark_combinations.py:109(benchmark_fournums)

Lessons:

  • Calling itertools.combinations_with_replacement is by far the fastest, taking about 2.7 seconds. It turns out that it’s written in C, so this would be hard to beat. (Still, writing it in a try block is seriously bad.)
  • The “equivalent” Python code from the itertools documentation (benchmark_itertools_combinations_with_replacment) is about 50x slower.
  • Gets slightly better when specialized to numbers.
  • Simply generating all combinations without an upper bound is actually faster.
  • It can be made even faster by writing it in a more C-like way.
  • The tuples version with the loop unrolled manually is rather fast when seen in this light, less than 4x slower than the library version.

Written by S

Sun, 2015-02-08 at 00:03:38

Boolean blindness

leave a comment »

(Just digesting the first page of Google search results.)

One of the lessons from functional programming is to encode as much information as possible into the types. Almost all programmers understand to some extent that types are helpful: they know not to store everything as void* (in C/C++) or as Object (in Java). They even know not to use say double for all numbers, or string for everything (as shell scripts / Tcl do). (This pejoratively called “stringly typed”.)

A corollary, from slide 19 here, is that (when your type system supports richer types) a boolean type is almost always the wrong choice, as it carries too little information.

The name “Boolean blindness” for this seems to have been coined by Dan Licata when taught a course at CMU as a PhD student.

From here (blog post by Robert Harper, his advisor at CMU):

There is no information carried by a Boolean beyond its value, and that’s the rub. As Conor McBride puts it, to make use of a Boolean you have to know its provenance so that you can know what it means.
[…]
Keeping track of this information (or attempting to recover it using any number of program analysis techniques) is notoriously difficult. The only thing you can do with a bit is to branch on it, and pretty soon you’re lost in a thicket of if-then-else’s, and you lose track of what’s what.
[…]

The problem is computing the bit in the first place. Having done so, you have blinded yourself by reducing the information you have at hand to a bit, and then trying to recover that information later by remembering the provenance of that bit. An illustrative example was considered in my article on equality:

fun plus x y = if x=Z then y else S(plus (pred x) y)

Here we’ve crushed the information we have about x down to one bit, then branched on it, then were forced to recover the information we lost to justify the call to pred, which typically cannot recover the fact that its argument is non-zero and must check for it to be sure. What a mess! Instead, we should write

fun plus x y = case x of Z => y | S(x') => S(plus x' y)

No Boolean necessary, and the code is improved as well! In particular, we obtain the predecessor en passant, and have no need to keep track of the provenance of any bit.

Some commenter there says

To the best of my knowledge, Ted Codd was the first to point out, in his relational model, that there is no place for Boolean data types in entity modeling. It is a basic design principle to avoid characterizing data in terms of Boolean values, since there is usually some other piece of information you are forgetting to model, and once you slam a Boolean into your overall data model, it becomes very hard to version towards a more exact model (information loss).

An example from Hacker News (on an unrelated post):

Basically, the idea is that when you branch on a conditional, information is gained. This information may be represented in the type system and used by the compiler to verify safety, or it can be ignored. If it is ignored, the language is said to have “boolean blindness”.
Example:


  if (ptr == NULL) {
    ... a ...
  } else {
    ... b ...
  }

In branch a and branch b, different invariants about ptr hold. But the language/compiler are not verifying any of these invariants.
Instead, consider:

  data Maybe a = Nothing | Just a


This defines a type “Maybe”, parameterized by a type variable “a”, and it defines two “data constructors” that can make a Maybe value: “Nothing” which contains nothing, and “Just” which contains a value of type “a” inside it.
This is known as a “sum” type, because the possible values of the type are the sum of all data constructor values.
We could still use this sum data-type in a boolean-blind way:


  if isJust mx then
    .. use fromJust mx .. -- UNSAFE!
  else
    .. can't use content of mx ..

However, using pattern-matching, we can use it in a safe way. Assume “mx” is of type “Maybe a”:

  case mx of
    Nothing -> ... There is no value of type "a" in our scope
    Just x -> ... "x" of type "a" is now available in our scope!


So when we branch on the two possible cases for the “mx” type, we gain new type information that gets put into our scope.
“Maybe” is of course a simple example, but this is also applicable for any sum type at all.

Another from notes of someone called HXA7241:

A nullable pointer has two ‘parts’: null, and all-other-values. These can be got at with a simple if-else conditional: if p is not null then ... else ... end. And there is still nothing wrong with that. The problem arrives when you want to handle the parts – and you lack a good type system. What do you do with the non-null pointer value? You just have to put it back into a pointer again – a nullable pointer: which is what you started with. So where has what you did with the test been captured? Nowhere. When you handle that intended part elsewhere you have to do the test again and again.

A Reddit discussion: Don Stewart says

Not strongly typed, richly typed. Where evidence isn’t thrown away. Agda is about the closest thing to this we have.

Relation to “object-oriented programming”

There’s also a good (and unintentionally revealing) example there by user munificent. Consider Harper’s example, which in more explicit Haskell could be:


data Num = Zero | Succ Num

plus :: Num -> Num -> Num
plus Zero     y = y
plus (Succ x) y = Succ (plus x y)

We might write it in Java as the following:

interface Num {}
class Zero implements Num {}
class Succ implements Num {
  public Succ(Num pred) {
    this.pred = pred;
  }
  public final Num pred;
}

Num plus(Num x, Num y) {
  if (x instanceof Succ) {      // (1)
    Succ xAsSucc = (Succ)x;       // (2)
    return new Succ(plus(xAsSucc.pred, y));
  } else {
    return y;
  }
}

Here instanceof returns a boolean (comment (1)), but doesn’t carry with it any information about what that boolean represents (namely that x is an instance of Succ) so when we get to the next line (2) we’re forced to do an unsafe cast. As programmers, we know it’s safe from context, but the compiler doesn’t.

There is a way in Java of avoiding this situation (where the programmer has context the compiler doesn’t):


interface Num {
  Num plus(Num other);
}

class Zero implements Num {
  public Num plus(Num other) {
    return other;
  }
}

class Succ implements Num {
  public Succ(Num pred) {
    this.pred = pred;
  }

  public Num plus(Num other) {
    return new Succ(pred.plus(y));
  }

  public final Num pred;
}

But see what has happened (by user aaronia):

There’s a rub though — in your first code snippet, “plus” was written as a non-member function; anyone can write it as a stand alone library. It was modular. But, as you demonstrated, it had the redundant type check.
However, your second code snippet has lost this modularity. You had to add new methods to the Num data type.

I think this is the pitfall of Object-Oriented Programming (which is like “Pants-Oriented Clothing”): objects can talk about themselves, but it’s harder to talk about objects.

If you want to write a function that does similar things for a bunch of types, you have to write similar function definitions in all those different classes. These function definitions do not stay together, so there is no way of knowing or ensuring that they are similar.

Written by S

Sat, 2015-01-31 at 14:55:00

Posted in computer science

Tagged with

C++ vs Java terminology glossary

leave a comment »

E.g. there is no such thing as a “method” in C++. As far as the C++ standard is concerned, functions inside classes are always called “member functions”.

Here is a handout from this course (taught by, AFAICT, Jeffrey S. Leon) that helps me understand the Java people.

comparable-terms

Written by S

Sat, 2015-01-31 at 11:53:06

Posted in programming

Tagged with

Visualizing product of permutations

with 4 comments

A simple pedagogical trick that may come in handy: represent a permutation \sigma using arrows (curved lines) from k to \sigma(k) for each k. Then, the product of two permutations can be represented by just putting the two corresponding figures (sets of arrows) one below the other, and following the arrows.

Representing permutations and products of permutations.

Representing permutations and products of permutations.

The figure is from an article called Symmetries by Alain Connes, found via the Wikipedia article on Morley’s trisector theorem (something entirely unrelated to permutations, but the article covers both of them and more).

I’m thinking how one might write a program to actually draw these: if we decide that the “height” of the figure is some h, then each arrow needs to go from some (k, 0) to (\sigma(k), h) (using here the usual screen convention of x coordinate increasing from left to right, and y coordinate increasing from top to bottom). Further, each curve needs to have vertical slope at its two endpoints, so that successive curves can line up smoothly. The constraint on starting point, ending point, and directions at the endpoints defines almost a quadratic Bezier curve, except that here the two directions are parallel. So it’s somewhere between a quadratic and the (usual) cubic Bezier curve, which is given by the start point, end point, and derivatives at the start and end point. (Here we only care about the direction of the derivative; we can pick some arbitrary magnitude to fix the curve: the larger we pick, the more smooth it will look at the ends, at the cost of smoothness in the interior.)

Even knowing the curve, how do we generate an image?

Written by S

Thu, 2014-03-06 at 23:15:44

The potions puzzle by Snape in Harry Potter and the Philosopher’s Stone

with one comment

Sometime this week, I reread the first Harry Potter book (after at least 10 years… wow, has it been that long?), just for contrast after reading Rowling’s adult novel The Casual Vacancy (on which more later). Anyway, in the penultimate chapter there is a puzzle:

He pulled open the next door, both of them hardly daring to look at what came next — but there was nothing very frightening in here, just a table with seven differently shaped bottles standing on it in a line.
“Snape’s,” said Harry. “What do we have to do?”
They stepped over the threshold, and immediately a fire sprang up behind them in the doorway. It wasn’t ordinary fire either; it was purple. At the same instant, black flames shot up in the doorway leading onward. They were trapped.
“Look!” Hermione seized a roll of paper lying next to the bottles. Harry looked over her shoulder to read it:

Danger lies before you, while safety lies behind,
Two of us will help you, whichever you would find,
One among us seven will let you move ahead,
Another will transport the drinker back instead,
Two among our number hold only nettle wine,
Three of us are killers, waiting hidden in line.
Choose, unless you wish to stay here forevermore,
To help you in your choice, we give you these clues four:
First, however slyly the poison tries to hide
You will always find some on nettle wine’s left side;
Second, different are those who stand at either end,
But if you would move onward, neither is your friend;
Third, as you see clearly, all are different size,
Neither dwarf nor giant holds death in their insides;
Fourth, the second left and the second on the right
Are twins once you taste them, though different at first sight.

I became curious about whether this is just a ditty Rowling made up, or the puzzle actually makes sense and is consistent. It turns out she has constructed it well. Let’s take a look. This investigation can be carried out by hand, but we’ll be lazy and use a computer, specifically Python. The code examples below are all to be typed in an interactive Python shell, the one that you get by typing “python” in a terminal.

So what we have are seven bottles, of which one will take you forward (F), one will let you go back (B), two are just nettle wine (N), and three are poison (P).

>>> bottles = ['F', 'B', 'N', 'N', 'P', 'P', 'P']

The actual ordering of these 7 bottles is some ordering (permutation) of them. As 7 is a very small number, we can afford to be inefficient and resort to brute-force enumeration.

>>> import itertools
>>> perms = [''.join(s) for s in set(itertools.permutations(bottles))]
>>> len(perms)
420

The set is needed to remove duplicates, because otherwise itertools.permutations will print 7! “permutations”. So already the number of all possible orderings is rather small (it is \frac{7!}{2!3!} = 420). We can look at a sample to check whether things look fine.

>>> perms[:10]
['PNFNPBP', 'NPPBNFP', 'FNNPBPP', 'PPPFNBN', 'NPPNBFP', 'PFNNBPP', 'NPBPPFN', 'NBNPPFP', 'NPPFBNP', 'BNPFNPP']

Now let us try to solve the puzzle. We can start with the first clue, which says that wherever a nettle-wine bottle occurs, on its left is always a poison bottle (and in particular therefore, a nettle-wine bottle cannot be in the leftmost position). So we must restrict the orderings to just those that satisfy this condition.

>>> def clue1(s): return all(i > 0 and s[i-1] == 'P' for i in range(len(s)) if s[i]=='N')
...
>>> len([s for s in perms if clue1(s)])
60

(In the code, the 7 positions are 0 to 6, as array indices in code generally start at 0.)
Then the second clue says that the bottles at the end are different, and neither contains the potion that lets you go forward.

>>> def clue2(s): return s[0] != s[6] and s[0] != 'F' and s[6] != 'F'
...
>>> len([s for s in perms if clue1(s) and clue2(s)])
30

The third clue says that the smallest and largest bottles don’t contain poison, and this would be of help to Harry and Hermione who can see the sizes of the bottles. But as we readers are not told the sizes of the bottles, this doesn’t seem of any help to us; let us return to this later.

The fourth clue says that the second-left and second-right bottles have the same contents.

>>> def clue4(s): return s[1] == s[5]
...
>>> len([s for s in perms if clue1(s) and clue2(s) and clue4(s)])
8

There are now just 8 possibilities, finally small enough to print them all.

>>> [s for s in perms if clue1(s) and clue2(s) and clue4(s)]
['PPNBFPN', 'BPNPFPN', 'BPFPNPN', 'BPPNFPN', 'PNPFPNB', 'BPNFPPN', 'PPNFBPN', 'PNFPPNB']

Alas, without knowing which the “dwarf” and “giant” bottles are, we cannot use the third clue, and this seems as far as we can go. We seem to have exhausted all the information available…

Almost. It is reasonable to assume that the puzzle is meant to have a solution. So even without knowing where exactly the “dwarf” and “giant” bottles are, we can say that they are in some pair of locations that ensure a unique solution.

>>> def clue3(d, g, s): return s[d]!='P' and s[g]!='P'
...
>>> for d in range(7):
...   for g in range(7):
...     if d == g: continue
...     poss = [s for s in perms if clue1(s) and clue2(s) and clue4(s) and clue3(d,g,s)]
...     if len(poss) == 1:
...       print d, g, poss[0]
... 
1 2 PNFPPNB
1 3 PNPFPNB
2 1 PNFPPNB
2 5 PNFPPNB
3 1 PNPFPNB
3 5 PNPFPNB
5 2 PNFPPNB
5 3 PNPFPNB

Aha! If you look at the possible orderings closely, you will see that we are down to just two possibilities for the ordering of the bottles.
Actually there is some scope for quibbling in what we did above: perhaps we cannot say that there is a unique solution determining the entire configuration; perhaps all we can say is that the puzzle should let us uniquely determine the positions of just the two useful bottles. Fortunately, that gives exactly the same set of possibilities, so this distinction happens to be inconsequential.

>>> for d in range(7):
...   for g in range(7):
...     if d == g: continue
...     poss = [(s.index('F'),s.index('B')) for s in perms if clue1(s) and clue2(s) and clue4(s) and clue3(d,g,s)]
...     if len(set(poss)) == 1:
...       print d, g, [s for s in perms if clue1(s) and clue2(s) and clue4(s) and clue3(d,g,s)][0]
... 
1 2 PNFPPNB
1 3 PNPFPNB
2 1 PNFPPNB
2 5 PNFPPNB
3 1 PNPFPNB
3 5 PNPFPNB
5 2 PNFPPNB
5 3 PNPFPNB

Good. Note that there are only two configurations above. So with only the clues in the poem, and the assumption that the puzzle can be solved, we can narrow down the possibilities to two configurations, and be sure of the contents of all the bottles except the third and fourth. We know that the potion that lets us go forward is in either the third or the fourth bottle.

In particular we see that the last bottle lets us go back, and indeed this is confirmed by the story later:

“Which one will get you back through the purple flames?”
Hermione pointed at a rounded bottle at the right end of the line.
[…]
She took a long drink from the round bottle at the end, and shuddered.

But we don’t know which of the two it is, as we can’t reconstruct all the relevant details of the configuration. Perhaps we can reconstruct something with the remaining piece of information from the story?

“Got it,” she said. “The smallest bottle will get us through the black fire — toward the Stone.”
Harry looked at the tiny bottle.
[…]
Harry took a deep breath and picked up the smallest bottle.

So we know that the bottle that lets one move forward is in fact in the smallest one, the “dwarf”.

>>> for d in range(7):
...   for g in range(7):
...     poss = [s for s in perms if clue1(s) and clue2(s) and clue4(s) and clue3(d,g,s)]
...     if len(poss) == 1 and poss[0][d] == 'F':
...       print d, g, poss[0]
... 
2 1 PNFPPNB
2 5 PNFPPNB
3 1 PNPFPNB
3 5 PNPFPNB

This narrows the possible positions of the smallest and largest bottles (note that it says the largest bottle is one that contains nettle wine), but still leaves the same two possibilities for the complete configuration. So we can stop here.

Conclusion
What we can conclude is the following: apart from the clues mentioned in the poem, the “dwarf” (the smallest bottle) was in either position 2 (third from the left) or 3 (fourth from the left). The biggest bottle was in either position 1 (second from the left) or 5 (sixth from the left). With this information about the location of the smallest bottle (and without necessarily assuming the puzzle has a unique solution!), Hermione could determine the contents of all the bottles. In particular she could determine the location of the two useful bottles: namely that the bottle that lets you go back was the last one, and that the one that lets you go forward was the smallest bottle.

>>> for (d,g) in [(2,1), (2,5), (3,1), (3,5)]:
...   poss = [s for s in perms if clue1(s) and clue2(s) and clue4(s) and clue3(d, g, s)]
...   assert len(poss) == 1
...   s = poss[0]
...   assert s.index('B') == 6
...   assert s.index('F') == d
...   print (d,g), s
... 
(2, 1) PNFPPNB
(2, 5) PNFPPNB
(3, 1) PNPFPNB
(3, 5) PNPFPNB

It is not clear why she went to the effort to create a meaningful puzzle, then withheld details that would let the reader solve it fully. Perhaps some details were removed during editing. As far as making up stuff for the sake of a story goes, though, this is nothing; consider for instance the language created for Avatar which viewers have learned.

Added:
See also http://www.zhasea.com/logic/snape.html which does it by hand, and has a perfectly clear exposition (it doesn’t try the trick of guessing that solution is unique before reaching for the additional information from the story).

Written by S

Mon, 2012-10-22 at 02:00:12

How does Tupper’s self-referential formula work?

with 16 comments

[I write this post with a certain degree of embarrassment, because in the end it turns out (1) to be more simple than I anticipated, and (2) already done before, as I could have found if I had internet access when I did this. :-)]

The so-called “Tupper’s self-referential formula” is the following, due to Jeff Tupper.

Graph the set of all points {(x,y)} such that

\displaystyle  \frac12 < \left\lfloor \mathrm{mod} \left( \left\lfloor{\frac{y}{17}}\right\rfloor 2^{-17\lfloor x \rfloor - \mathrm{mod}(\lfloor y \rfloor, 17)}, 2 \right) \right\rfloor

in the region

\displaystyle  0 < x < 106

\displaystyle  N < y < N+17

where N is the following 544-digit integer:
48584506361897134235820959624942020445814005879832445494830930850619
34704708809928450644769865524364849997247024915119110411605739177407
85691975432657185544205721044573588368182982375413963433822519945219
16512843483329051311931999535024137587652392648746133949068701305622
95813219481113685339535565290850023875092856892694555974281546386510
73004910672305893358605254409666435126534936364395712556569593681518
43348576052669401612512669514215505395545191537854575257565907405401
57929001765967965480064427829131488548259914721248506352686630476300

The result is the following graph:

Figure 1: The graph of the formula, in some obscure region, is a picture of the formula itself.

Whoa. How does this work?

At first sight this is rather too incredible for words.

But after a few moments we can begin to guess what is going on, and see that—while clever—this is perhaps not so extraordinary after all. So let us calmly try to reverse-engineer this feat.

Read the rest of this entry »

Written by S

Tue, 2011-04-12 at 13:05:20

Posted in mathematics

Tagged with , , ,

Firebug “console is undefined”

with 2 comments

If you’re using Firebug and don’t want to bother removing or commenting out all the console.log debug messages for users who aren’t, put this at the top:

if(typeof(console) === "undefined" || typeof(console.log) === "undefined") 
    var console = { log: function() { } };

This reminds me of the trick I use for C/C++ of putting debugging statements inside a macro:

D(cerr<<"The number is "<<n<<endl;);

where the macro D is defined as

#define D(A) A

when you want the debugging code to run, and

#define D(A) 

when you don’t. :-)

(The final semicolon after the D() is for Emacs to indent it reasonably.)

Update: Of course, the above are just hacky hacks to save typing. The “right” way for conditional C code is usually to use #ifdef DEBUG and the like, and the right way around the Firebug problem is to use something more general like this code at metamalcolm.

Written by S

Wed, 2009-08-05 at 14:32:15

Follow

Get every new post delivered to your Inbox.

Join 99 other followers