Print a specific value in this generated binary matrix

J - 42 38 char

Pretty fast, 100% accurate, and well within the memory constraints.

([{+2&(~:/@#:@#@],~:/\,(p:>:)&#)0:)&<:

The strategy is as follows: we will calculate successive antidiagonals of this matrix, performing a pairwise XOR to move along and adding the current Thue-Morse and prime bits to the ends. We then pull the required digit out of the antidiagonal when we get there.

Explanation by explosion:

(                                 )&<:  NB. decrement each of x and y
     &(                        )        NB. apply the following function...
   +                                    NB. ... (x-1)+(y-1) times...
                                0:      NB. ... starting with a zero:
    2             ~:/\                  NB.   pairwise XOR on the argument
                      ,(p:>:)&#         NB.   append prime bit (is 1+length prime?)
       ~:/@#:@#@],                      NB.   prepend TM bit (XOR of binary)
 [{                                     NB. take the x-th bit (at index x-1)

Usage of this verb is x m y for M(x, y) as specified in the question, where m is the verb.

   5 ([{+2&(~:/@#:@#@],~:/\,(p:>:)&#)0:)&<: 8
0
   m =: ([{+2&(~:/@#:@#@],~:/\,(p:>:)&#)0:)&<:
   1+i.16
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
   m/~ 1+i.16
0 1 1 0 1 0 0 1 1 0 0 1 0 1 1 0
1 0 1 1 0 0 0 1 0 0 0 1 1 0 1 1
1 1 0 1 1 1 1 0 0 0 0 1 0 0 1 0
0 1 1 0 1 0 1 1 1 1 1 0 0 0 1 1
1 0 1 1 0 0 1 0 1 0 1 1 1 1 0 1
0 0 1 0 0 0 1 1 0 0 1 0 1 0 0 1
1 1 0 0 0 0 1 0 0 0 1 1 0 0 0 1
0 1 1 1 1 1 0 0 0 0 1 0 0 0 0 1
0 1 0 1 0 1 1 1 1 1 0 0 0 0 0 1
0 1 1 0 0 1 0 1 0 1 1 1 1 1 1 0
1 0 1 1 1 0 0 1 1 0 1 0 1 0 1 1
0 0 1 0 1 1 1 0 1 1 0 0 1 1 0 1
1 1 0 0 1 0 1 1 0 1 1 1 0 1 1 0
0 1 1 1 0 0 1 0 0 1 0 1 1 0 1 1
0 1 0 1 1 1 0 0 0 1 1 0 1 1 0 1
0 1 1 0 1 0 0 0 0 1 0 0 1 0 0 1

To save keystrokes we don't try to tell if we still need more prime or Thue-Morse bits, so we compute the entire antidiagonal to get the bit we want. However, 8192 m 8192 still runs in less than 0.07 s and about 100 KiB on my modest laptop.

Mathematica – 100% accuracy, 223 193 189 bytes

f=(r=Array[0&,Max@##];For[s=2,s<=#+#2,++s,For[i=Max[1,s-#2],i<=Min[s-1,#],++i,j=s-i;r[[j]]=Which[i==1,PrimeQ@j,j==1,OddQ@Total@IntegerDigits[i-1,2],0<1,Xor@@r[[j-1;;j]]]]];If[r[[#2]],1,0])&

Here is a legible version:

f[x_,y_] := (
   r = Array[0 &, Max[x,y]];
   For[s = 2, s <= x + y, ++s,
    For[
     i = Max[1, s - y],
     i <= Min[s - 1, x],
     ++i,

     j = s - i;
     r[[j]] = Which[
       i == 1,
       PrimeQ@j,
       j == 1,
       OddQ@Total@IntegerDigits[i - 1, 2],
       0 < 1,
       r[[j - 1]]~Xor~r[[j]]
       ]
     ]
    ];
   If[r[[y]], 1, 0]
   );

I basically precompute along diagonals of constant x+y.

Features:

• It's accurate.
• It runs in O(x*y).
• f[8192,8192] takes about 400 seconds. I suppose there is room for improvement (maybe RotateLeft could replace the inner loop).

• It only uses one array of up to max(x,y) intermediate results in memory. So there is no necessity to use more than about 32k (assuming 32-bit integers) for the algorithm itself (plus, whatever Mathematica uses). In fact, Mathematica uses 31M by itself on my system, but this works without an issue:

  MemoryConstrained[f[8192, 8192], 2^21]
  

Matlab: 100% accuracy, 120 characters, unreasonable execution time

function z=M(x,y)
if y==1 z=(x>1)*isprime(x);elseif x==1 z=mod(sum(dec2bin(y-1)-48),2);else z=xor(M(x,y-1),M(x-1,y));end

To use:

> M(4,4)
ans =
      0
> M(1, 9)
ans =
      1