Official xkcd puzzle [SOLUTIONS]

[BaldAdonis]

Right, so, here's what I've got:
To help visualize the function, construct a mx4 matrix, labelling the rows {0, 1, 2,...} and columns {0, 1, 2, 3}. Then F(N) is the unique (i, j) entry, where 4i+j=N. The matrix entries are determined as follows:

Partition the natural numbers, including zero, into "houses" (my written solution has a little pentagon around the number, I'll use n here) of size 64. That is, {0,...,63} are in house n=0, {64,...,127} are in n=1, and so on.
Choose a house n as follows: to start, take houses (0 2 3 1) in that order. Given a previous set of 4 houses, choose the next four by taking the inverse permutation, staring at 0, of the last four. I make this sound more complicated than it is: if the last four were (0 2 3 1), take (0 1 3 2). If the last four were (0 1 3 2), take (0 2 3 1). Place the houses in the matrix; the order of elements within each house will be determined later.

Partition each house into four blocks (this time it's a square around the number, or k) of size 16. That is, {0,...,15 (mod 64)} are in k=0, and so on, as above. For more on (mod x), see below.
Choose a block k as follows: to start, take blocks (0 1 3 2) to form the first house. Given a previous block, invert as above.

Partition each block into four groups (it's a circle! or i, and you thought I was going for the triangle) of size 4. Each group i takes one of two forms: if it is in the first or second block k, in n=0, {0, 1, 2, 3 (mod 16)} belong in i=0, and so on. If it is in the third or fourth block, in n=0, then {0, 2, 5, 7 (mod 16)} are in i=0, {1, 3, 4, 6 (mod 16)} are in i=1, add 8 to i=0 to get i=2, add 8 to i=1 to get i=3. Given a previous house n, the i groups in the next house are opposite, that is, the first and second blocks contain the alternating i groups, the third are fourth contain the standard i groups.
Choose a group i as follows: to start, take (0 2 3 1). Invert as above.

If you have trouble following the i groups, write down a 4x4 block of 16 integers. The standard i group is any row of 4 integers. The alternating i group takes the first and third of one row, and second and fourth of the adjacent row, and puts them in numerical order. i=0 starts at the first row and goes to the second, i=1 starts at the second and goes to the first. i=2 and 3 are similar.

Each group i has four elements, denoted (0, 1, 2, 3) in numerical order.
Place each element from i as follows: in block k=0, to start, for the first i (always i=0) take (0 1 3 2). For the next i, invert. For the following i, use the inverse (in this case, (0 2 3 1)) again. For the final i, invert the previous, but BEGIN the cycle at 3, that is (3 2 0 1).
Given any block, the arrangement of elements of i in the next block in the inverse of the previous block, with the final i always beginning at 3.

Yikes! If you follow all of this notation, here are the next 10 (and for good measure, six more), starting at N=36:
56 61 63 58 (n=0, k=3, i=2)
57 60 62 59 (n=0, k=3, i=3)
54 52 51 49 (n=0, k=3, i=1).
32 37 39 34 (n=0, k=2, i=0)

To find a specific value F(N), write N=x_0 + (4^1)x_1 +(4^2)x_2 + ...., where each x is chosen from {0, 1, 2, 3}.
We must then find the (x_0)th entry of the (x_1)th group of the (x_2)th block of the (x_3)th house.... Greater values can be defined inductively following the pattern of blocks and houses. So 4 houses form a street (n_2, with a hexagonal stop sign!), and the streets are chosen in the same manner as blocks: first (0 1 3 2), then the inverse. 4 streets make a neighbourhood (this neighbourhood is heptagonal shaped, n_3), with neighbourhoods chosen just like houses.
So F(N) = y_0 + (4^1)y_1 +(4^2)y_2 + ...., where each y is found by applying the recursions above.

As an example, we find the 31337th number, N=31336= 0 +2(4^1)+2(4^2)+1(4^3)+2(4^4)+2(4^5)+3(4^6)+1(4^7)

x_8 =0
x_7=1, so F(N) is in the second choice for n_5, since x_8 is even, it is the initial order, which is the same as the initial order of houses, and the second chosen is 2. So y_7=2.
x_6 =3, since x_7 is odd, it is in the inverted order, the same as inverted order of blocks, so y_6=1.
x_5 =2, since x_6 is odd, it is in the inverted order, so y_5=3.
x_4 =2, since x_5 is even, it is in the initial order, so y_4=3.
x_3 =1, since x_4 is even, it is in the initial order, so y_3=2.
x_2 =2, since x_3 is odd, it is in the inverted order, so y_2=3.
x_1 =2. Now, since x_3=1, x_2=2, i is a standard group. Since x_2 is even, the groups are in initial order, so y_1=3
x_0 =0, so y_0=0.
Hence F(31336)=40892.
Similarly, F(31337)=40893.

This function is clearly bijective, since every N takes a value, and the integers {0, 1,...,(4^n)-1} are mapped to {0, 1,...,(4^n)-1}.

To find "a (mod b)", write a=bx + y, where x, y>=0, y<b. Then y=a (mod b). So {0, 2, 5, 7} = {16, 18, 21, 23} = {32, 34, 37, 39} (mod 16).

Snow Day! No exams. Cleaned up typos.

*Edit*
Just reread the question, and it asks for a general algorithm, which i didn't state explicitly. Write N=x_0 + (4^1)x_1 +(4^2)x_2 + .... as above. For each x_k, k>1, evaluate y_k by reading down. Invert x_k if x_(k+1) is odd.

x_k = 0 1 2 3
y_k = 0 2 3 1 if k is odd and in initial order, or even and inverted
y_k = 0 1 3 2 if k is even and in initial order, or odd and inverted

For x_1, determine if the i group is standard or alternating: standard if x_3 is even and x_2=0, 1 or x_3 is odd and x_2=2, 3, alternating otherwise. Find the group and element position as above, using x_1 and x_0. If it is standard, the values are y_1 and y_0. If it is alternating, consider the following 4x4 array, with terms (x_1, x_0), and use the row and column as y_1, y_0 respectively.

(0,0) (1,0) (0,1) (1,1)
(1,2) (0,2) (1,3) (0,3)
(2,0) (3,0) (2,1) (3,1)
(3,2) (2,2) (3,3) (2,3)

**Edit**
To find y_0 in standard group:

x_0 = 0 1 2 3
y_0 = 0 2 3 1 if x_1=0 and x_2 is odd, or x_1=1, 2 and x_2 is even
y_0 = 0 1 3 2 if x_1=1, 2 and x_2 is odd, or x_1=0 and x_2 is even
y_0 = 3 2 0 1 if x_1=3, and x_2 is even
y_0 = 3 1 0 2 if x_1=3, and x_2 is odd.

I've decided I hate it when exercises are left to the reader.
71=3+1(4)+0(4^2)+1(4^3)
x_3=1 => y_3=2
x_2=0 => y_2=0
x_3 odd, x_2=0 => alternating group
x_1=1, x_0=3 => y_1=1, y_0=2 by consulting the array.
F(71)=2+1(4)+0(4^2)+2(4^3)=134.


***It occurs to me that the alternating groups may want to rotate in the (0132) manner, more on this to come as I imagineer it.

[xkcd]

In retrospect I shouldn't have said "post your answer", but I didn't expect one so quick :)

This is an interesting approach, and I'm so busy with holiday packing that I can't actually check all the way through it for the next 24 hours. But it looks pretty good!

Now, there's still a shirt open for 'most elegant algorithm' -- I already have a conceptual approach that lets me generate the numbers. That's how I got the original list. What I'd like to see is a couple lines of code or a simple equation that takes in the argument, does a trick, and spits out the result. I feel like there should be a simple iterative solution with shuffling bits, but I don't have it.

[Gelsamel]

What BaldAdonis said except I came up with it 10 years ago


xD

[phlip]

So... it doesn't actually have anything to do with the Gray code?

Man, I was on the wrong track...

[Token]

Now, there's still a shirt open for 'most elegant algorithm' -- I already have a conceptual approach that lets me generate the numbers. That's how I got the original list. What I'd like to see is a couple lines of code or a simple equation that takes in the argument, does a trick, and spits out the result. I feel like there should be a simple iterative solution with shuffling bits, but I don't have it.

I suspect you shall only be giving out one t-shirt, depending on your definition of "simple".

[BaldAdonis]

I'm not much for programming, but I reckon this should work if someone wants to try it out.

Operations:
A) (1 --> 2--> 3), 0 fixed
B) (2--> 3), 0, 1 fixed
C) (1--> 2--> 0--> 3)
D) (0--> 3--> 2), 1 fixed

Write N in base 4 notation, x_k...x_3.x_2.x_1.x_0.
Apply operations on x_k to get y_k, k>1.
If k + x_(k+1) is odd, apply A.
If k + x_(k+1) is even, apply B.

--> Apply A to x_1 if 1+x_2+x_3+... is odd
--> Apply B to x_1 if 1+x_2+x_3+... is even

-->If ((x_1=0, x_2 + x_3 is odd) OR (x_1=1, 2, x_2 + x_3 is even) apply A to x_0.
--> If ((x_1=0, x_2 + x_3 is even) OR (x_1=1, 2, x_2 + x_3 is odd) apply B to x_0.
--> If (x_1=3, x_1 + x_2 is odd) apply C to x_0.
--> If (x_1=3, x_1 + x_2 is even) apply D to x_0.
Output y_k....y_3.y_2.y_1.y_0.


*Edited to make it nicer to look at*
**Edited to remove alternating groups**

[xkcd]

So... it doesn't actually have anything to do with the Gray code?

Man, I was on the wrong track...

It still might. This is something of a different approach from the one I used to generate the list :)

[BaldAdonis]

Looking at Wikipedia's entry on the Gray code, it seems like it could do the trick, if you reversed it in appropriate places. I came at this ignoring everything binary, which was clearly intended by the problem, so don't take my solution to be the best. In fact, I think my solution may slip at a few points, but we only got the first 36 to test, so I won't ever really know for sure.

[Cosmologicon]

Well, I give up. I've been staring at this sequence for hours now. I'll just throw up what I've found and hope that helps someone else.

The glaring pattern I found is that F(N) XOR F(N+1) is always one of {1, 2, 5, 10, 21, 42, ... }. These numbers are, in binary, {1, 10, 101, 1010, 10101, 101010, ... } This is highly analagous to a Gray code, in which XORed consecutive numbers are always one of {1, 10, 100, 1000, 10000, ... }.

The analogy with a standard Gray code is not complete, though, because the ordering of the XORed values does not follow the same pattern as for a Gray code. If you did it like that, you would get the following sequence:

Code: Select all

   N   G(N) G(N) bin      G(N) ^ G(N-1)  F(N) ^ F(N-1)
   0     0  00000000           0             0
   1     1  00000001           1             1
   2     3  00000011           2             2
   3     2  00000010           1             1
   4     7  00000111           5            10
   5     6  00000110           1             2
   6     4  00000100           2             1
   7     5  00000101           1             2
   8    15  00001111          10             5
   9    14  00001110           1             2
  10    12  00001100           2             1
  11    13  00001101           1             2
  12     8  00001000           5            10
  13     9  00001001           1             1
  14    11  00001011           2             2
  15    10  00001010           1             1
  16    31  00011111          21            21
  17    30  00011110           1             2
  18    28  00011100           2             1
  19    29  00011101           1             2
  20    24  00011000           5             5
  21    25  00011001           1             1
  22    27  00011011           2             2
  23    26  00011010           1             1
  24    16  00010000          10            10
  25    17  00010001           1             1
  26    19  00010011           2             2
  27    18  00010010           1             1
  28    23  00010111           5             5
  29    22  00010110           1             2
  30    20  00010100           2             1
  31    21  00010101           1             2
  32    63  00111111          42            42
  33    62  00111110           1             2
  34    60  00111100           2             5
  35    61  00111101           1             2


This function shares several properties with F(N). It's a bijection, or at least it appears to be, and consecutive numbers differ in either even bits or odd bits only. The only difference is that the sequencing of bit swaps is slightly different. Looking at the right two columns (the ^ indicates exclusive-or) they're tantalizingly similar, but the first one is easy to generate, and I just can't come up with a pattern for the rightmost one.

Well... good luck!

[Cosmologicon]

BaldAdonis, I'm trying to implement the algorithm in your post, but I don't understand what you mean by "apply A to x_0", for instance. Does that mean you change the value of y_0 or x_0 or both?

[Cosmologicon]

Nevermind, I think I get the general idea. But if so, I'm not getting the right answers. Just to begin with, I get F(1) = 2 rather than F(1) = 1. Here are the relevant operations:

A) (1 --> 2--> 3), 0 fixed

-->If ((x_1=0, x_2 + x_3 is even) OR (x_1=1, 2, x_2 + x_3 is odd) apply A to x_0.

Since x_1 = 0 and x_2 + x_3 = 0 is even, I apply A to x_0, converting it from a 1 to a 2. Is that not right?

Originally I accidentally swapped the A and the B in this step and the next one, and it looked much closer to correct. Should they be swapped by any chance?

[GreedyAlgorithm]

I must admit I cannot rationalize F(34)=55 and F(35)=53. If those really are the correct values, I submit that my sequence F' is nicer, and more in keeping with the newly revealed comic (which is now outside my office). For 0<=x<=33 and x=71, F'(x)=F(x), but F'(34)=51 and F'(35)=49.

F'(x) for 36<=x<=46:
52, 53, 55, 54, 60, 61, 63, 62, 59, 57, 56

F'(31337) = 36994.

To generate F', though not elegantly:

Code: Select all

int N = 42;
int[][][] state =
{
   { {1,0,0}, {0,1,0}, {0,1,1}, {3,0,1} },
   { {0,0,0}, {1,0,1}, {1,1,1}, {2,1,0} },
   { {3,1,1}, {2,0,1}, {2,0,0}, {1,1,0} },
   { {2,1,1}, {3,1,0}, {3,0,0}, {0,0,1} }
};
int M = 15; //datatype has at least M*2 bits
int x = 0; //(x,y) coordinates in space, N units along space-filling curve
int y = 0;
int m = (1<<(M-1)); //current scale
int s = (M-1)%2==0 ? 0 : 1; //starting state
for (int i = M-1; i >= 0; i--) { //for each decreasing scale
   int r = (N>>>(i*2))%4; //value of N at that place in base 4
   x += state[s][r][1]*m; //add the offset (0 or 1) times the scale
   y += state[s][r][2]*m;
   s = state[s][r][0]; //transition to a new curve direction
   m = m>>>1; //decrease scale
}

//convert (x,y) to F(N)
int f = 0; //start value at corner
m = 1; //current scale
while(x>0) { //each bit of x
   f += (x%2)*m; //add the offset (0 or 1) times the scale
   x = x>>>1; //delete lowest bit
   m = m<<2; //increase scale
}
m = 2;
while(y>0) {
   f += (y%2)*m;
   y = y>>>1;
   m = m<<2;
}
System.out.println(f);


[eridius]

Here's a program I wrote in logo (using the ucblogo interpreter, I don't know if this program is portable to other logo dialects). It figures out the number by actually drawing the hilbert fractal. Unfortunately I think the prettiest picture shows up when drawing using degree 8, but to preserve the mapping present in the original problem the degree needs to be odd.

Update: I've made it use the appropriate degree, and then swap the coordinates if the degree is even. This preserves the mapping, but lets the pretty degree 8 draw

Code: Select all

to R
    if :n = 0 [stop]
    make "n :n - 1
    rt 90 L F
    lt 90 R F R
    lt 90 F L rt 90
    make "n :n + 1 ; restore :n to value before procedure was called
end

to L
    if :n = 0 [stop]
    make "n :n - 1
    lt 90 R F
    rt 90 L F L
    rt 90 F R lt 90
    make "n :n + 1
end

to F
    if :count = 0 [
        make "coord (map [abs round ? / :dist] list  (first pos) + 250 (last pos) - 250)
        ifelse :swap = 1 [print reverse :coord] [print :coord]
        throw "hilbert.result
    ]
    make "count :count - 1
    FD :dist
end

to abs :n
    ifelse :n < 0 [op 0 - :n] [op :n]
end

to HILBERT :n :count
    hideturtle
    fence
    penup
    setpos [-250 250]
    setheading 180
    clean
    pendown
   
    local "dist
    make "dist 500 / (power 2 :n - 1 / power 2 :n)
   
    catch "hilbert.result [L]
end

type "|Enter a number to calculate: |
make "num readword
make "degree ((int (ln :num + 1) / (ln 4)) + 1)
; if degree is even, we need to flip the coordinates from the output
; because every other degree is flipped
ifelse (modulo :degree 2) = 0 [make "swap 1] [make "swap 0]
HILBERT :degree :num

bye


Note: I couldn't figure out any way to find the bounds of the window, but ucblogo, at least on OS X in X11, defaults to a window of 500x500, hence the [-250 250] starting position. If you try and run this and the default bounds on your system are different you may need to change this.

[Strilanc]

This is a hard one.

Maybe a useful property: If you use the function on the output repeatedly, you seem to usually get your input back eventually. For example, 4-8-12-7-9-14-4, or 6-11-13-6. The only exception I found to this was 20-21-22-23-22. But I'm sure that's just because of 20's attitude towards making this problem as hard as possible.

[Torn Apart By Dingos]

I must admit I cannot rationalize F(34)=55 and F(35)=53. (...) For 0<=x<=33 and x=71, F'(x)=F(x), but F'(34)=51 and F'(35)=49.

I don't get it either. With your F' values, the sequence is just a slight modification of the Hilbert curve.

[Strilanc]

Alright, I have an iterative algorithm that matches all the output you gave. The only part that feels "wrong" is the exception for F(34) and F(35).
EDIT: removed the F(34) and F(35) since they've been corrected now. Also added F(31337).

Algorithm (Java):

Code: Select all

/** An iterative solution to the puzzle
 * Essentially loops through [0, 1, 3, 2] back and forth to get the answer
 * @Author: Craig Gidney */
public class SolDigit {
   private int val=0, //current value
     inc=1, //direction
     index=0, //digit position
     count=0; //number of increments of this digit
   private SolDigit next=null, prev=null; //for linking this list

   /** When run, display output of F on [0, 35] U [71, 71] */
   public static void main(String args[]) {
      SolDigit d = new SolDigit();
      int a;
      //Print out 0-36
      for (a = 0; a < 36; a++) {
         System.out.println(d);
         d.increment();
      }
      //Print out F(71)
      for (; a < 71; a++)
         d.increment();
      System.out.println(d);
      //Print out F(31337)
      for (; a < 31337; a++)
         d.increment();
      System.out.println(d);      
   }   
   
   /** turns our values into the digits they actually represent (just swaps 2 and 3) */
   private int digit() {
      return (val == 2) ? 3 : (val == 3) ? 2 : val;
   }
   
   /** Increment to the next output of F */
   public void increment() {
      
       //increment
      val = (val + inc + 4) % 4;
       count++;
      
       //switch directions at these values
       if (count % 16 == 4 || count  % 16 == 8 || count  % 16 == 11 || count  % 16 == 15)
           inc *= -1;
       //switch direction of previous digits at val == 2
       if (val == 2)
           for (SolDigit p = prev; p != null; p = p.prev)
               p.inc *= -1;
      
       //increment next digit when we roll over
       if (count % 4 == 0) {
           //create it if it doesn't exist
          if (next == null) {
              next = new SolDigit();
              next.index = index + 1;
              next.inc = (int)Math.pow(-1, next.index); //initial direction
               next.prev = this;
           }
           next.increment();
      }
   }
   
   /** Sum the digits' values */
   public int actualVal()  {
      //Leftshift (base 4) next digit's value and add ours
      return ((next != null) ? next.actualVal()*4 : 0) + digit();
   }
   
   public String toString() {
      return "F(" + count + ") = " + actualVal();
   }
}

Output:

Code: Select all

F(0) = 0
F(1) = 1
F(2) = 3
F(3) = 2
F(4) = 8
F(5) = 10
F(6) = 11
F(7) = 9
F(8) = 12
F(9) = 14
F(10) = 15
F(11) = 13
F(12) = 7
F(13) = 6
F(14) = 4
F(15) = 5
F(16) = 16
F(17) = 18
F(18) = 19
F(19) = 17
F(20) = 20
F(21) = 21
F(22) = 23
F(23) = 22
F(24) = 28
F(25) = 29
F(26) = 31
F(27) = 30
F(28) = 27
F(29) = 25
F(30) = 24
F(31) = 26
F(32) = 48
F(33) = 50
F(34) = 51
F(35) = 49
F(71) = 134
F(31337) = 37757

[xkcd]

Boy, is there egg on my face. The last two numbers were indeed wrong (I did this by hand, but that should be no excuse). The last four should read:

32 48
33 50
34 51
35 49

More to come.

[xkcd]

BaldAdonis, your answer for F(31337) in your main post is off, but GreedyAlgorithm got it right.

[Strilanc]

BaldAdonis, your answer for F(31337) in your main post is off, but GreedyAlgorithm got it right.

I got a completely different value: F(31337) = 37757. But can you really say an answer is "wrong" if it matches all the given input?

Ok, maybe an exception for algorithms that just go "If x = y then return z, if x = y2 then return z2, etc." Those are wrong.

[WAHa.06x36]

Code: Select all

sub fb($)
{
   my @table=([0,1,2,3],[0,2,1,3],[3,2,1,0],[3,1,2,0]);

   my $n=shift;
   if($n==0) { return (0,0) }
   elsif($n==1) { return (1,1) }
   elsif($n==2) { return (3,1) }
   elsif($n==3) { return (2,2) }
   {
      my ($b1,$trans1)=fb($n&3);
      my ($b2,$trans2)=fb($n>>2);

      return ((($b2<<3)&0xaaaaaaa8)|(($b2<<1)&0x55555554)|$table[$trans2][$b1],$trans1^$trans2);
   }
}


I have a version that makes sense, too, but it's not as funny, so I posted this instead.

[DeadCatX2]

I took GreedyAlgorithm's algorithm and ported it to C++, and then proceeded to look for neat patterns.

Check out this list of powers of 2

Code: Select all

N       f(N)    base2                   f(N)^f(N-1)     base2
1       1       000000000000000000001   1               000000000000000000001
2       3       000000000000000000011   2               000000000000000000010
4       8       000000000000000001000   10              000000000000000001010
8       12      000000000000000001100   5               000000000000000000101
16      16      000000000000000010000   21              000000000000000010101
32      48      000000000000000110000   42              000000000000000101010
64      128     000000000000010000000   170             000000000000010101010
128     192     000000000000011000000   85              000000000000001010101
256     256     000000000000100000000   341             000000000000101010101
512     768     000000000001100000000   682             000000000001010101010
1024    2048    000000000100000000000   2730            000000000101010101010
2048    3072    000000000110000000000   1365            000000000010101010101
4096    4096    000000001000000000000   5461            000000001010101010101
8192    12288   000000011000000000000   10922           000000010101010101010
16384   32768   000001000000000000000   43690           000001010101010101010
32768   49152   000001100000000000000   21845           000000101010101010101
65536   65536   000010000000000000000   87381           000010101010101010101
131072  196608  000110000000000000000   174762          000101010101010101010
262144  524288  010000000000000000000   699050          010101010101010101010
524288  786432  011000000000000000000   349525          001010101010101010101
1048576 1048576 100000000000000000000   1398101         101010101010101010101

[Locoluis]

Hi.

This is not what you've been asking for in the puzzle. It's just a function that would take an integer and return the pair of coordinates.

Code: Select all

/**
 * Author: Luis Alejandro González Miranda.
 *
 * Parameters:
 * number - The number you want to lookup
 * scale - The scale of the number, the total of elements in the curve. The number should be in the range [0 .. scale[
 * cx, cy - The returned coordinates in terms of qs (see below)
 */
void hilpoint(long number, long scale, double *cx, double *cy)
{
   /* This is a description of the shape of each fundamental pattern */
   int patterns[4][4]={{0,1,3,2},{0,2,3,1},{3,1,0,2},{3,2,0,1}};

   /* This is a description of what happens to each pattern the next iteration */
   int transforms[4][4]={{1,0,2,0},{0,3,1,1},{2,2,0,3},{3,1,3,2}};

   /* Quadrant coordinates in pixels; should be initialized to the width
    * and height of your arena, which must be square.
    * The coordinates will be expressed in numbers in the range [0..qs[
    */
   double qs=1.0;

   /* This variable should be initialized according to the scale of the
    * numbers at the Hilbert curve. We will need to divide by this number,
    * and the integer part of the result should be always between 0 and 3.
    *
    * For example, if you're making a map of Class A IP addresses, your
    * range is [0..256[; therefore, the value of divisor is 256/4, or 64;
    */
   long divisor=scale/4;

   /* A copy of your number */
   long tmp=number;

   /* The start pattern. Don't modify */
   long pattern=1;

   (*cx)=(*cy)=0.0;

   while(divisor>0) {
      int cell;
      long mantissa;
      
      mantissa=tmp/divisor;
      tmp=tmp%divisor;
      divisor/=4;
      cell=patterns[pattern][mantissa];
      pattern=transforms[pattern][cell];
      qs/=2.0;

      if(cell & 1) (*cx)+=qs;
      if(cell & 2) (*cy)+=qs;

   }
}


[BaldAdonis]

BaldAdonis, I'm trying to implement the algorithm in your post, but I don't understand what you mean by "apply A to x_0", for instance. Does that mean you change the value of y_0 or x_0 or both?

Sorry about that, A and B should be swapped at that point. I wrote it out calling them "even" and "odd" permutations, when I transferred it from paper to desktop it got garbled.
What I meant was, input x_0, cycle through (1 2 3), get y_0.

Boy, is there egg on my face. The last two numbers were indeed wrong (I did this by hand, but that should be no excuse). The last four should read:

32 48
33 50
34 51
35 49

This changes everything! Alternating groups aren't needed, which means the entire half of the new algorithm (if x_2=2,3...) is unnecessary.

[WAHa.06x36]

The code I posted earlier has the nice property that it is "easily" invertible. So here is a complete example with both the normal and inverse function:

Code: Select all

#!/usr/bin/perl

use strict;

for(0..100,31337)
{
   my ($bit)=fb($_);
   my ($inv)=fbi($bit);
   print "$_: $bit: $inv\n";
}

sub fb($)
{
   my @table=([0,1,2,3],[0,2,1,3],[3,2,1,0],[3,1,2,0]);
   my $n=shift;
   if($n==0) { return (0,0) }
   elsif($n==1) { return (1,1) }
   elsif($n==2) { return (3,1) }
   elsif($n==3) { return (2,2) }
   else
   {
      my ($b1,$trans1)=fb($n&3);
      my ($b2,$trans2)=fb($n>>2);
      return ((($b2<<3)&0xaaaaaaa8)|(($b2<<1)&0x55555554)|$table[$trans2][$b1],$trans1^$trans2);
   }
}

sub fbi($)
{
   my @table=([0,1,2,3],[0,2,1,3],[3,2,1,0],[3,1,2,0]);
   my $n=shift;
   if($n==0) { return (0,0) }
   elsif($n==1) { return (1,1) }
   elsif($n==2) { return (3,2) }
   elsif($n==3) { return (2,1) }
   else
   {
      my ($b2,$trans2)=fbi( (($n&0xaaaaaaa8)>>3) | (($n&0x55555554)>>1) );
      my ($b1,$trans1)=fbi($table[$trans2][$n&3]);
      return (($b2<<2)|$b1,$trans1^$trans2);
   }
}


Edit: Oops, fixed a last-minute bug. Should work now.

[BaldAdonis]

Right, so with the (new) values, I've touched up my solution, and generated a list. xkcd, how did you generate the list in the first place?

Operations:
A) (1 --> 2--> 3 --> 1), 0 fixed
B) (2--> 3 --> 2), 0, 1 fixed
C) (1--> 2--> 0--> 3 --> 1)
D) (0--> 3--> 2 --> 0), 1 fixed
To apply an operation, input x_k, permute, and output y_k. So for A,
x=0 => y=0
x=1 => y=2
x=2 => y=3
x=3 => y=1

Write N in base 4 notation, x_k...x_3.x_2.x_1.x_0.
Apply operations on x_k to get y_k, k>1.
If k + x_(k+1) is odd, apply A.
If k + x_(k+1) is even, apply B.

--> Apply A to x_1 if 1+x_2+x_3+... is odd
--> Apply B to x_1 if 1+x_2+x_3+... is even

-->If x_2=x_3=...=x_k=0
----> If x_1=0, apply B to x_0.
----> If x_1=1,2, apply A to x_0.
----> If x_1=3, apply C to x_0.
Else
----> If x_1=0, apply A to x_0.
----> If x_1=1,2, apply B to x_0.
----> If x_1=3, apply D to x_0.
----> Output y_k....y_3.y_2.y_1.y_0.

0 = 000 --> 000 = 0
1 = 001 --> 001 = 1
2 = 002 --> 003 = 3
3 = 003 --> 002 = 2
4 = 010 --> 020 = 8
5 = 011 --> 022 = 10
6 = 012 --> 023 = 11
7 = 013 --> 021 = 9
8 = 020 --> 030 = 12
9 = 021 --> 032 = 14
10 = 022 --> 033 = 15
11 = 023 --> 031 = 13
12 = 030 --> 013 = 7
13 = 031 --> 012 = 6
14 = 032 --> 010 = 4
15 = 033 --> 011 = 5

16 = 100 --> 100 = 16
17 = 101 --> 102 = 18
18 = 102 --> 103 = 19
19 = 103 --> 101 = 17
20 = 110 --> 110 = 20
21 = 111 --> 111 = 21
22 = 112 --> 113 = 23
23 = 113 --> 112 = 22
24 = 120 --> 130 = 28
25 = 121 --> 131 = 29
26 = 122 --> 133 = 31
27 = 123 --> 132 = 30
28 = 130 --> 123 = 27
29 = 131 --> 121 = 25
30 = 132 --> 120 = 24
31 = 133 --> 122 = 26

32 = 200 --> 300 = 48
33 = 201 --> 302 = 50
34 = 202 --> 303 = 51
35 = 203 --> 301 = 49
36 = 210 --> 320 = 56
37 = 211 --> 321 = 57
38 = 212 --> 323 = 59
39 = 213 --> 322 = 58
40 = 220 --> 330 = 60
41 = 221 --> 331 = 61
42 = 222 --> 333 = 63
43 = 223 --> 332 = 62
44 = 230 --> 313 = 55
45 = 231 --> 311 = 53
46 = 232 --> 310 = 52
47 = 233 --> 312 = 54

71 = 1013 --> 2012 = 134

31337 = 13221221 --> 21332331 = 40893

[rektide]

^-- dont answer that q... yet.

[WAHa.06x36]

Further refinement: A completely iterative version in C99, using only bitwise operations and table lookups, and using more useful coordinates instead of the wacky interleaved numbers. Both normal and inverse functions:

Code: Select all

#include <stdio.h>

void f(int n,int *xp,int *yp);
int fi(int x,int y);

int main(int argc,char **argv)
{
   for(int i=0;i<102;i++)
   {
      int x,y;

      if(i==101) i=31337;
      f(i,&x,&y);

      int il=0;
      for(int j=0;j<16;j++) il|=( (x<<j)&(1<<(2*j)) )|( (y<<(j+1))&(1<<(2*j+1)) );

      printf("%d: %d: %d %s\n",i,il,fi(x,y),i!=fi(x,y)?"!!!":"");
   }
   return 0;
}

static int transform_table[4][4]={{0,1,2,3},{0,2,1,3},{3,2,1,0},{3,1,2,0}};
static int locations[4]={0,1,3,2};

void f(int n,int *xp,int *yp)
{
   static int transforms[4]={1,0,0,3};
   int x=0,y=0;
   int trans=0;
   for(int i=30;i>=0;i-=2)
   {
      int m=(n>>i)&3;
      int bits=transform_table[trans][locations[m]];

      x=(x<<1)|((bits>>1)&1);
      y=(y<<1)|(bits&1);

      trans^=transforms[m];
   }
   *xp=x; *yp=y;
}

int fi(int x,int y)
{
   static int transforms[4]={1,0,3,0};
   int n=0;
   int trans=0;
   for(int i=15;i>=0;i--)
   {
      int m=transform_table[trans][((y>>i)&1)|(((x>>i)&1)<<1)];
      int bits=locations[m];

      n=(n<<2)|bits;

      trans^=transforms[m];
   }
   return n;
}


[Cosmologicon]

I really wanted to find a recursive solution, but in the end, I couldn't get anything better than this boring old non-recurvise one. The only thing I really like about this is that there are no lookup tables for the basic curve. Well, try to find them in the code anyway, if you feel up to it....

Code: Select all

unsigned F (const unsigned N) {
  unsigned m = 1, n = N, r = 0, p, a = 0;
  while (n >>= 2) ++m, a ^= 1;
  while (m)
    (r <<= 2) |= (p = ((n = N >> (--m << 1) & 3) > 1) ^ n) ^
      3 * (a & p << 1 >> p ^ a >> 1),
    a ^= n * (n & 1) ^ (n < 2);
  return r;
}


EDIT: combined a couple of lines into one. No need to be excessively verbose.

[BaldAdonis]

Drawing out the points N in a plane, with F(N) as the binary interleaved coordinates, results in a Hilbert curve. Turns out that this binary coordinate graph is called the Z-curve, and what the function does is map the Hilbert curve position to the Z-curve, giving coordinates for it. Then it hits me that this last comic used the Hilbert curve, and everything starts to make sense.

The Hilbert curve has some interesting symmetries about it. Consider a block of 4^n connected points: The Hilbert curve starts and ends on the same side of the square around those points, so that there are only 4 combinations of different curves. Also, every curve goes in only one direction, based on the construction of the Hilbert curve.
Curve A begins at the top left, and goes to the bottom left. (0-->0, 1-->1, 2-->3, 3-->2). Looks like a backwards C.
Curve B goes from top left to top right. (0-->0, 1-->2, 2-->3, 3-->1). Looks like a U.
Curve C goes from bottom right to top right. (0-->3, 1-->2, 2-->0, 3-->1). Looks like a C.
Curve D goes from bottom right to botom left. (0-->3, 1-->1, 2-->0, 3-->2). Looks like an inverted U.

This holds for any set of 4^n connected points. The first 16 points connect along ABBC, resulting in (draw it!) a B curve, order 2 (where the order of 4^n connected points is n). This block of 16 connects to the next block of 16 (BAAD), which is similar to A, order 2. This connects to another A curve (BAAD, there's only one way to do it), and finally to a D curve, order 2 (CDDA), giving 64 connected points in an A formation, order 3. This connects to another block of 64, a B curve, order 3, and from here, the self symmetry of the Hilbert curve takes care of the rest.
The symmetries are:
BAAD = A
ABBC = B
DCCB = C
CDDA = D

Back to the problem at hand: finding the coordinates for a given point on the Hilbert curve. Write the number in base 4, x_k.x_(k-1)...x_2.x_1.x_0. Each digit x_i describes to which block of 4^i the point belongs, within a block of 4^(i+1)
The curve begins in the top left corner, so the first 4^(k+1) points must be in either A or B formation. It's easy to see that if k is even, it is A, if odd, it is B (consider k=1, and note that the formation alternates by adding 1).

Now for a spicy new algorithm:
Input two values: F from {A, B, C, D} for the formation
and x_i from {0, 1, 2, 3} for the point within the fromation.
Output new formation and y_i.

Each formation contains two pieces of information: where the point x_i will end up in terms of a Z-curve, and which formation will result from focusing on the point x_i within the curve.

A=BAAD, (0-->0, 1-->1, 2-->3, 3-->2).
B=ABBC, (0-->0, 1-->2, 2-->3, 3-->1).
C=DCCB, (0-->3, 1-->2, 2-->0, 3-->1).
D=CDDA, (0-->3, 1-->1, 2-->0, 3-->2).

Since there are two pieces of information, we have two functions. F(x_i) returns the value for y_i, and is found by taking x_i to be the left element and tracing to the right element in the sets described above. So B(0)=0, B(1)=2, B(2)=3, B(3)=1.
[F, x_i] returns a formation, and is found by taking the (x_i)th letter equivalent to the formation F, beginning at 0. So [C, 0]=D, [C, 1]=[C, 2]=C, [C, 3]=B.

The y_i value output should be saved, while the formation output should be used in the next iteration as the new formation. I'm sure there's a slick way of programming that.

Iterate this k+1 times, then output y_k.y_(k-1)...y_1.y_0 (or convert to decimal, and output).

For example: 71=1.0.1.3 in base 4. The maximum k value is 3, which is odd, so we input B, 1.
First output is then B(1)=2, since B takes 1 to 2, and [B, 1]=B, since B=ABBC, and the second value is B. So y_3=2.
Input B, 0. Output B(0)=0, [B, 0]=A. So y_2=0.
Input A, 1. Output A(1)=1, [A, 1]=A. So y_1=1
Input A, 3. Output A(3)=2, [A, 3]=D. So y_0=2, thus F(71)=2012=134.

Another example: 31337=13221221. k=7, so we input B, 1.
First output is B(1)=2, [B, 1]=B. y_7=2.
Input B, 3. Output B(3)=1, [B, 3]=C. y_6=1
Input C, 2. Output C(2)=0, [C, 2]=C. y_5=0
Input C, 2. Output C(2)=0, [C, 2]=C. y_4=0
Input C, 1. Output C(1)=2, [C, 1]=C. y_3=2
Input C, 2. Output C(2)=0, [C, 2]=C. y_2=0
Input C, 2. Output C(2)=0, [C, 2]=C. y_1=0
Input C, 1. Output C(2)=2, [C, 2]=C. y_5=2

F(31337)=21002002=36994.
You can check by drawing this out as a section of the Hilbert curve. Start by drawing a B curve, shaped like a U. You'll probably want to cover a full sheet of paper with it. Starting at the top left corner, point zero, move to point 1, since x_7=1. Note that this corner is number 2 in a Z-curve (or using binary coordinates), and hence y_7=2. Circle this corner, and draw the corresponding curve around that point (it's a B curve again). Starting at point 0, move to point 3, as x_6=3. Note that this is 1 in a Z-curve, so y_6=1. Draw the curve around this point, this time a C. Circle point 2, (0 in a Z-curve) remembering that a C curve moves from bottom right to top right, then draw the curve around that. Continue until you finish with 6 embedded C curves, and a very pretty section of Hilbert curves of increasing complexity.

[Cosmologicon]

Wow, BaldAdonis, I approached this problem exactly like you did. I even called the four shapes A, B, C, and D just like you! My variable a encodes whether you're currently in an A, B, C, or D subcurve as 0, 1, 2, or 3. I called it m instead of k, but still. I think if you wrote out a program to implement your algorithm and then played around with the binary operations to compactify it, you'd get the same thing I have. Good explanation!

[BaldAdonis]

I wrote a program with the (diminished) C++ I picked up a few years ago, but I don't have a compiler so I have no idea if it works. Basically just a function that takes F={A, B, C, D} and x={0, 1, 2, 3}, iterates and spits out values for y. The worst part was converting to base 4, and back again.

[Cosmologicon]

When I don't have a compiler handy I use Javascript. The basic syntax is very similar. It should probably work okay if you take out the ints. If you're interested, just drop it into a <script> tag in the header of an HTML file, then make the body tag say:

<body onload="alert(F(31337))">

If you've gotten it 100% right to begin with (a miracle in programming) it should show you the right answer when you load the document. You can pm me too if you want help with the Javascript syntax.

[WAHa.06x36]

Here's how to find the recursive solution: Each level of blocks follows the curve, but alternatingly transposed. To find the coordinates of a point, find the coordinates for the number divided by 4, flip the X and Y axes, multiply by two, and then add the position within the 4-point block. This requires the combined transformations from the earlier steps, so that needs to be calculated at the same time as the position.

This is what my first solution does. When expressed as such, it is fairly easy to find the reverse operation, too, just by studying the steps the code takes.

My final solution replaces the recursion with iteration, and does away with the coordinate swap step by instead swapping the coordinates in the subblock being studied, and generates two coordinates instead of a single interleaved number. It ends up being pretty much equivalent to Cosmologicon's solution, except that he uses the interleaved number, and that he's expressed the transforms in code instead of lookup tables.

[WAHa.06x36]

Additionally: The transformations form the group C2*C2, which means they can be expressed as the set of integerrs 0,1,2,3 and the combining operation as XOR, which makes calculating combined transforms very easy.

[Cosmologicon]

I do like your recursive perl version, WAHa.06x36, but why does it return two numbers? I assumed at first the second number was an iteration count, which is less than ideal in my opinion, but now I have no clue.

[DeadCatX2]

Code: Select all

unsigned F (const unsigned N) {
  unsigned m = 1, n = N, r = 0, p, a = 0;
  while (n >>= 2) ++m, a ^= 1;
  while (m)
    (r <<= 2) |= (p = ((n = N >> (--m << 1) & 3) > 1) ^ n) ^
      3 * (a & p << 1 >> p ^ a >> 1),
    a ^= n * (n & 1) ^ (n < 2);
  return r;
}


...wow. Just...wow.

Could you explain it a little?

I mean, I'm looking at the symbols, and I see the trees (xors, shifts, multiplies, etc..), but no forest.

Again, kudos for coming up with an amazing bit manipulation.

[WAHa.06x36]

As I said, when doing it recursively, you need to both know the position on the curve at the next level up, and also the orientation of that block. The first number is the position, the second the orientation. For the final answert you can throw away the orientation, but the intermediate steps need it.

[Cosmologicon]

Yeah.... that makes sense. I was sort of hoping there'd be a way around that, though. You know, some way to transform the results of the deeper level of recursion to account for the orientation shift. Then you could base F(N) solely on lower F values and not extra information. That's really the beauty of a recursive solution over an iterative one, to me. But I spent hours looking for it and I didn't come up with anything as good as yours, so maybe there's nothing to be found. Ah well.

[Cosmologicon]

An explanation for DeadCatX2.

Really most of it's not cleverness; it's just obfuscation. In writing this post, though, I realized how to do it even more efficiently, without the variable p. So here's my improved version, with spaces removed for even less readability:

Code: Select all

unsigned F (const unsigned N) {
  unsigned m = 1, n = N, r = 0, a = 0;
  while (n>>=2) a=1&m++;
  while (m--) r=r*4|(n=N>>m+m&3)^n>>1^3*(((a^=n<3?!n:3)+(n&1)&2)>>1);
  return r;
}

I saved a couple characters in the main loop by replacing "<<2" with "*4", etc, but usually I prefer bit operations when they're possible. Now, the things like the nested assignments are just silliness. So if you cleaned it up and made it readable you'd get the following code:

Code: Select all

unsigned F (const unsigned N) {
  unsigned m = 1, n = N, r = 0, a = 0;
  while (n >>= 2) ++m, a ^= 1;
  while (m--) {
    n = N >> (m << 1) & 3;
    unsigned d = n ^ n >> 1;         // lookup: (0,1,3,2)
    if (a & n & 1 ^ a >> 1) d ^= 3;
    r = r << 2 | d;
    a ^= n * (n & 1) ^ (n < 2);      // lookup: (1,0,0,3)
  }
  return r;
}

Read BaldAdonis's post from earlier today to understand the inner workings. n, d, and m (which BaldAdonis calls x_i, y_i, and k) are respectively the base-4 digits of the input, the base-4 digits of the output, and the number of base-4 digits in each. r, of course, is the running total of the output. a keeps track of the orientation information. In BaldAdonis's notation, a is 0, 1, 2, or 3, for curve types A, B, C, and D, respectively. The only tricky part, I guess, is the two lines where you get d as a function of n and a.

It occurs to me that if you were doing this with a stack-based language, you could just push all the base-4 digits of N onto the stack in the first loop and pop them off in the second loop. This would look even neater, and you wouldn't need the variable m. In C++, though, you don't have access to "the" stack. You can use a stack object, but it's not as pretty:

Code: Select all

#include <stack>
unsigned F (const unsigned N) {
  unsigned n = N, r = 0, a = 1;
  std::stack<unsigned> s;
  do { s.push(n & 3); a ^= 1; } while (n /= 4);
  while (!s.empty()) {
    n = s.top(); s.pop();
    unsigned d = n & 2 ? n ^ 1 : n;
    if (a & n & 1 ^ a >> 1) d ^= 3;
    r = r * 4 + d;
    a ^= n < 3 ? !n : 3;
  }
  return r;
}

[valexodamium]

In python:

Code: Select all

def Hilbert(input, exp, list, result=0, start=0):
    posit = [0, 1, 3, 2]
    for i in range(len(list)):
        if ((start+list[i]*(4**exp)) <= input) & ((start+(list[i]+1)*(4**exp)) > input):
            start = start + list[i]*(4**exp)
            result = result*4 + posit[i]
            if list[i] == 3:
                b = list[0]
                list[0] = list[2]
                list[2] = b
            if list[i] == 0:
                b = list[1]
                list[1] = list[3]
                list[3] = b
            break
    if exp>0:
        Hilbert(input, exp-1, list, result, start)
    else:
        print result


def FindExp(input):
    exp = 0
    while(4**(exp+1))<=input:exp+=2
    return exp


def Hilb(input):
    Hilbert(input, FindExp(input), [0, 1, 2, 3])


print "Enter the number to be converted, q to quit, or l to loop"
while True:
   i = raw_input("> ")
   if i == "q":
      break
   elif i == "l":
      print "Up to what?"
      u = int(raw_input("> "))
      for x in range(u):
         Hilb(x)
   else:
      i = int(i)
      Hilb(i)


there. It estimates how big the grid should be, then recursively decides in which quadrant the number will be found, until the number is found. Flips are added when needed...
Hopefuly this isn't to bad...