[Question]
  [
Over at our friends at [Puzzling.SE](http://puzzling.stackexchange.com), the following puzzle was posted: [Is this chromatic puzzle always solvable?](https://puzzling.stackexchange.com/questions/36479/is-this-chromatic-puzzle-always-solvable) by Edgar G. You can play it [here](http://www.puzzlopia.com/puzzles/pixelea/play).
## Puzzle explanation
Given a `m x n` grid with tiles of three different colours, you may select **any two adjacent tiles**, if their colours are **different**. These two tiles are then converted to the third colour, i.e., the one colour not represented by these two tiles. The puzzle is solved if **all tiles have the same colour**. Apparently, one can [prove](https://puzzling.stackexchange.com/a/36491/10180) that this puzzle is always solvable, if neither `m` nor `n` are divisible by 3.
[![8x8 three colour puzzle](https://i.stack.imgur.com/W7Jv2.jpg)](https://i.stack.imgur.com/W7Jv2.jpg)
Of course, this begs for a solving algorithm. You will write a function or program that solves this puzzle. Note that functions with 'side effects' (i.e., the output is on `stdout` rather than in some awkward data type return value) are explicitly allowed.
## Input & Output
The input will be an `m x n` matrix consisting of the integers `1`, `2` and `3` (or `0`, `1`, `2` if convenient). You may take this input in any sane format. Both `m` and `n` are `>1` and not divisible by 3. You may assume the puzzle is not solved 
You will then solve the puzzle. This will involve a repeated selection of two adjacent tiles to be 'converted' (see above). You will output the two coordinates of these tiles for each step your solving algorithm took. This may also be in any sane output format. You are free to choose between 0-based and 1-based indexing of your coordinates, and whether rows or columns are indexed first. Please mention this in your answer, however.
Your algorithm should run within reasonable time on the original 8x8 case. Brute-forcing it completely is explicitly disallowed, i.e. your algorithm should run under `O(k^[m*(n-1)+(m-1)*n])` with `k` the number of steps needed for the solution. The solution however is **not** required to be optimal. The proof given in the linked question may give you an idea as to how to do this (e.g., first do all columns using only vertically adjacent tiles, and then do all rows)
## Test cases
In these test cases, the coordinates are 1-based and rows are indexed first (like MATLAB/Octave and probably many others).
```
Input: 
[1 2]
Output: (result: all 3's)
[1 1],[1,2]
Input:
[ 1 2
  3 1 ]
Output: (result: all 1's)
[1 1],[2 1]        (turn left column into 2's)
[2 1],[2 2]        (turn right column into 3's)
[1 1],[1 2]        (turn top row into 1's)
[2 1],[2 2]        (turn bottom row into 1's)
Input:
[1 2 3 2
 3 2 1 1]
Output: (result: all 3's)
[1 1],[1 2] 
[1 3],[1 4] 
[1 2],[1 3] 
[1 1],[1 2] 
[1 2],[1 3] 
[1 1],[1 2]
[1 3],[1 4]
[2 1],[2 2]
[1 1],[2 1]
[1 2],[2 2]
[1 3],[2 3]
[1 4],[2 4]
```
If desired, I may post a pastebin of larger test cases, but I think this should be sufficient.
  
[Answer]
## Octave, ~~334~~ 313 bytes
Since the challenge may seem a bit daunting, I present my own solution. I did not formally prove that this method works (I guess that will come down to proving that the algorithm will never get stuck in a loop), but so far it works perfectly, doing 100x100 testcases within 15 seconds. Note that I chose to use a function with side effects rather than one that returns all the coordinates since that saved me a few bytes. Coordinates are row-major, 1-based, and formatted as `row1 col1 row2 col2`. Input colours are `0,1,2` since this works better with `mod`, at the cost of having to use `numel` rather than `nnz`. Golfed version: *Edit: saved another few bytes by using a technique from Kevin Lau's answer.*
```
function P(p)
k=0;[m,n]=size(p);t=@(v)mod(sum(v)*numel(v),3);for c=1:n
p(:,c)=V(p(:,c));end
k=1;for r=1:m
p(r,:)=V(p(r,:));end
function a=V(a)
while any(a-(g=t(a)))
isempty(q=find(diff(a)&a(1:end-1)-g,1))&&(q=find(diff(a),1));
a([q q+1])=t(a([q q+1]));if k
disp([r q r q+1])
else
disp([q c q+1 c])
end;end;end;end
```
Example GIF of the solving algorithm:
[![enter image description here](https://i.stack.imgur.com/njo6G.gif)](https://i.stack.imgur.com/njo6G.gif)
Ungolfed version:
```
function solveChromaticPuzzle(p)
[m,n]=size(p);                           % Get size
t=@(v)mod(sum(v)*numel(v),3);            % Target colour function
for c=1:n                                % Loop over columns
    p(:,c)=solveVec(p(:,c));             % Solve vector
end
for r=1:m                                % Loop over rows
    p(r,:)=solveVec(p(r,:));
end
    function a=solveVec(a)               % Nested function to get globals
        g=t(a);                          % Determine target colour
        while(any(a~=g))                 % While any is diff from target...
            % Do the finding magic. Working left-to-right, we find the
            % first pair that can be flipped (nonzero diff) that does not
            % have the left colour different from our goal colour
            q=min(intersect(find(diff(a)),find(a~=g)));
            if(isempty(q))               % In case we get stuck...
                q=find(diff(a),1);       % ... just flip the first possible
            end;
            a([q q+1])=t(a([q q+1]));    % Do the actual flipping.
            if(exist('r'))               % If we're working per row
                disp([r q r q+1])        % Print the pair, using global row
            else
                disp([q c q+1 c])        % Print the pari, using global col
            end
        end
    end
end
```
[Answer]
## Lua, ~~594~~ ~~575~~ 559 Bytes
**Warning** There's still lots of work before I'm done with this golfing! I should be able to take that under 500 Bytes, at the very least. For the moment, it's the first solution that worked, and I'm still working on it.
I will provide a full explanation once I'm done.
```
function f(t)s=#t a=","for i=1,s do p=t[i]for i=1,s
do p.Q=p.Q and p.Q+p[i]or p[i]end p.Q=(p.Q*#p)%3 for i=1,s do for j=1,#p-1 do
x=p[j]y=p[j+1]u=x~=y and(p.S and p.R==p.S or x~=p.Q)v=(x+y)*2p[j]=u and v%3or x
p[j+1]=u and v%3or y print(i..a..j,i..a..j+1)end
p.R=p.S p.S=table.concat(p)end end
for i=1,s do Q=Q and Q+t[i][1]or t[i][1]end Q=(Q*s)%3 for i=1,s
do for j=1,s-1 do p=t[j]q=t[j+1]x=p[1]y=q[1]u=x~=y and(S and R==S or x~=Q)v=(x+y)*2
for k=1,#p do p[k]=u and v%3or x q[k]=u and v%3or y
print(j..a..k,j+1..a..k)end Y=Y and Y..x or x end
R=S S=Y end end
```
[Answer]
# Ruby, 266 bytes
More-or-less just a port of the Octave solution, except it solves rows first instead of columns. Input is an array of arrays, with the inner arrays being the rows. Output moves are `[row, column, row, column]`. [Test suite](https://repl.it/C9jd)
```
->m{t=->v{v.size*v.inject(:+)%3}
s=->a,x,r{g=t[a]
(q=(r=0..a.size-2).find{|i|a[i]!=a[i+1]&&g!=a[i]}||r.find{|i|a[i]!=a[i+1]}
a[q,2]=[t[a[q,2]]]*2
p r ?[x,q,x,q+1]:[q,x,q+1,x])while[]!=a-[g]}
m.size.times{|i|s[m[i],i,1]}
m=m.shift.zip *m
m.size.times{|i|s[m[i],i,p]}}
```
### Ungolfed with explanation
```
->m{                                  # Start lambda function, argument `m`
  t=->v{v.size*v.inject(:+)%3}        # Target color function
  s=->a,x,r{                          # Function to determine proper moves
                                      #   a = row array, x = row index, r = horizontal
    g=t[a]                            # Obtain target color
    (
      q=(r=0..a.size-2).find{|i|      # Find the first index `i` from 0 to a.size-2 where...
        a[i]!=a[i+1]                  # ...that element at `i` is different from the next...
        &&g!=a[i]                     # ...and it is not the same as the target color
      } || r.find{|i|a[i]!=a[i+1]}    # If none found just find for different colors
      a[q,2]=[t[a[q,2]]]*2            # Do the color flipping operation
      p r ?[x,q,x,q+1]:[q,x,q+1,x]    # Print the next move based on if `r` is truthy
    ) while[]!=a-[g]}                 # While row is not all the same target color, repeat
m.size.times{|i|                      # For each index `i` within the matrix's rows...
  s[m[i],i,1]                         # ...run the solving function on that row
                                      #   (`r` is truthy so the moves printed are for rows)
}
m=m.shift.zip *m                      # Dark magic to transpose the matrix
m.size.times{|i|s[m[i],i,p]}}         # Run the solving function on all the columns now
                                      #   (`r` is falsy so the moves printed are for columns)
```
[Answer]
# Rust, ~~496~~ 495 bytes
Sadly I can't beat OP, but for a Rust answer I am quite satisfied with the bytecount.
```
let s=|mut v:Vec<_>,c|{
let p=|v:&mut[_],f,t|{
let x=|v:&mut[_],i|{
let n=v[i]^v[i+1];v[i]=n;v[i+1]=n;
for k in f..t+1{print!("{:?}",if f==t{(k,i,k,i+1)}else{(i,k,i+1,k)});}};
let l=v.len();let t=(1..4).find(|x|l*x)%3==v.iter().fold(0,|a,b|a+b)%3).unwrap();
let mut i=0;while i<l{let c=v[i];if c==t{i+=1;}else if c==v[i+1]{
let j=if let Some(x)=(i+1..l).find(|j|v[j+1]!=c){x}else{i-=1;i};x(v,j);}else{x(v,i);}}t};
p(&mut (0..).zip(v.chunks_mut(c)).map(|(i,x)|{p(x,i,i)}).collect::<Vec<_>>(),0,c-1usize)};
```
**Input:** a vector of numbers as well as the number of columns. E.g.
```
s(vec!{1,2,1,3},2);
```
outputs
```
 (row1,col1,row2,col2)
```
to the command line.
I first solve every row and then solve the resulting column only once ,but print the steps for all columns. So it is actually quite efficient :-P.
**With formating:**
```
let s=|mut v:Vec<_>,c|{  
    let p=|v:&mut[_],f,t|{     // solves a single row/column
        let x=|v:&mut[_],i|{   // makes a move and prints it 
            let n=v[i]^v[i+1]; // use xor to calculate the new color
            v[i]=n;
            v[i+1]=n;
            for k in f..t{
                print!("{:?}",if f==t{(k,i,k,i+1)}else{(i,k,i+1,k)});
            }
        };
        let l=v.len();
        // find target color
        // oh man i am so looking forward to sum() being stabilized
        let t=(1..4).find(|x|(l*x)%3==v.iter().fold(0,|a,b|a+b)%3).unwrap();
        let mut i=0;
        while i<l{
            let c=v[i];
            if c==t{             // if the color is target color move on
                i+=1;
            }else if c==v[i+1]{ // if the next color is the same
                                // find the next possible move
                let j=if let Some(x)=(i+1..l).find(|j|v[j+1]!=c){x}else{i-=1;i};
                x(v,j);
            }else{              // colors are different so we can make a move
                x(v,i);         
            }
        }
        t
    };
    // first solve all rows and than sovle the resulting column c times 
    p(&mut (0..).zip(v.chunks_mut(c)).map(|(i,x)|p(x,i,i)).collect::<Vec<_>>(),0,c-1usize)
};
```
**Edit:**
now returns the color of the solution which saves me a semicolon^^
[Answer]
# [Befunge](https://fr.wikipedia.org/wiki/Befunge), ~~197~~ ~~368~~ ~~696~~ 754 Bytes
---
(yes, I'm doing some reverse code golf, the more bytes the better)
  
I was thinking it could be a challenge to write this algorithm in Befunge and that it could be fun
~~I would like it to be a community program, so if anyone wants to work on it, please, do so.~~
In the end, I made everything alone so far, so I'll finish on my own (it's almost done)
---
What's done yet: a troll-shaped code
```
&:19p&:09p:v:p94g92g90  <
 v94+1:g94&_$59g1+59p1-:|
 >p59gp1-: ^    vp95g93$<
v94\-1<v:g<     >  v
>g:1+v^_$v^90p94g92<
v5p94<   3>1+59p   ^
>9gg+\:^ %g   v93:g95<           v3_2         v
v1pg95g94<^95<>g-v>$v^           v ^-%3*2\gg9<
>9g39g+59g-1-|v-1_^ #^pg95+g92g90<1_09g:29g+5^
       ;  >  >  09g:29g+59gg\3%-# !^v         <
          ^p95<                  ^  <
     v  p96-1+g90g92<
     v                  p95_@
            >59g1+:39g-19g-^
     v    >p 79g:59gg\1+59gp$$$$$29g49pv
    > p59g^ |<<<<<<<<<<<<<<<<<<!-g96g94<
>:79^>29g49p>69g1+59gg49g:59gg\1+49p- v
^\-\6+gg95+1\g< v         !-g96:<-1g94_^
>"]",52*,::59g^v_::1+59gg\59gg-v^ <
^ .-g93g95,.-g<>:69g- v  v-g96:_1+^
>+" [,]",,,\29^       >#v_$:49g2-v
^1:.-g93g95,.-g92\,"[ ":<        <
```
(yes, it's a troll, believe me)
  
Basically, it reads an array and computes the move to do to solve the rows, given an input as
```
(number of rows) (number of columns) 1 2 3 1 1 3 2 1 2 ....
```
(the whole array is passed as a list [row1, row2, row3,…])
output is
```
[col row],[col',row']
[col2 row2],[col2',row2']
...
```
rows and cols both start at 0.
  
Now that rows are solved, it's almost done ! Hooray !
---
Explanation: (will be updated later on)
![image](https://i.stack.imgur.com/CxCxS.png)
So there are 5 main parts :
* The first one, in green, reads the input line and writes one row of the array
* The second one, in orange, passes to the next row of the array
* The third, in the blue, sums a row
* The fourth one, in hot pink, takes the modulus 3 of the sum, saves it at the right of the concerned row, and go to the next row
* Finally, in the red, the part that computes the target color from the previously computed number. This part is really dumb and should probably be rewritten, but I didn't figure out how I could do that in a nice way (passed from 197 bytes to 368 with just that part)
Grey parts are initialisations
---
Here is a a deeper explanation of the module that finds the to boxes to combine (which is coded here, by the way)
```
                                       B
            @                          v
            |                  !-g96g94<
ENTRY>29g49p>69g1+59gg49g:59gg\1+49p- v
                v         !-g96:<-1g94_^
               v_::1+59gg\59gg-v^ <
               >:69g- v  v-g96:_1+^
                      >#v_$:49g2-v
                    CALL<        <
```
The CALL part is when the instruction pointer is going to another module, to combine to boxes. It comes back to this module through the 'B' entry
Here is some pseudo code:
(currentx is related to the array reading)
 For:
```
    69g1+59gg  // read target color
    49g:59gg\1+49p // read current color and THEN shift the currentx to the next box
    if ( top != top ){  // if the current color is bad
        49g1-          //  put the current place  (or currentx - 1) in the stack
        While:
            if ( :top != 69g ){   // if you didn't reach the end of the list
                ::1+              // copy twice, add 1
                if ( 59gg == \59gg ){ // if the next color is the same than current
                   1+                // iterate
                   break While;
                }
            }
        : // copies j (the previous raw iterator)
        if ( top == 69g ){  // if you reached the end of the row (which mean you can't swap with the color after that point)
            $:              // discard j's copy and copy target
            49g2-           // put the place just before the color change on the stack
            combine_with_next;
        } else {
            combine_with_next;
        }
        29g49p   // back to the beginning of the row (there was some changes int the row)
    }
    if ( 49g != 69g ) // if you didn't reach the end of the list
        break For:
```
---
Note that if you want to test it, you will have to put some trailing space and trailing new lines so that there is enough space to store the array, if you wish to use [the interpret](http://www.quirkster.com/iano/js/befunge.html) I linked.
22 + the number of rows in input as trailing lines, and 34 + the number of columns as trailing spaces on one line should be ok.
[Answer]
# C, 404 bytes
My first code golf, I'm pretty happy with how it turned out. Took way too long, though. It's not fully standard C, just whatever will compile under gcc with no special flags (and ignoring warnings). So there's a nested function in there. The function `f` takes the dimensions `m` and `n` as its first arguments, and as it's third argument takes takes an (int pointer) to an array of size `m`×`n` (indexed by rows first). The other arguments are dummy arguments, you don't need to pass them in, they're just there to save bytes on declaring variables. It writes each changed pair to STDOUT in the format `row1,col1:row1,col1;`, with the semicolon separating the pairs. Uses 0-based indexing.
```
#define A a[i*o+p
#define B A*c
f(m,n,a,i,j,o,p,b,c,d)int*a;{
int t(x){printf("%d,%d:%d,%d;",b?i:c+x,b?c+x:i,b?i:c+1+x,b?c+1+x:i);}
o=n;p=b=1;for(;~b;b--){
for(i=0;i<m;i++){c=j=0;
for(;j<n;)c+=A*j++];d=c*n%3;
for(j=0;j<n-1;j++) 
while(A*j]^d|A*j+p]^d){
for(c=j;B]==B+p];c++);
if(c<n-2&B]==d&2*(B+p]+A*(c+2)])%3==d)
B+p]=A*(c+2)]=d,t(1);else
B]=B+p]=2*(B]+B+p])%3,
t(0);}}o=m;p=m=n;n=o;o=1;}}
```
I used a slightly different algorithm than OP for solving individual rows/columns. It goes something like this(pseudocode):
```
for j in range [0, rowlength):
    while row[j] != targetCol or row[j+1] != targetCol:
        e=j
        while row[e] == row[e+1]:
            e++             //e will never go out of bounds
        if e<=rowLength-3 and row[e] == targetCol 
                and (row[e+1] != row[e+2] != targetCol):
            row.changeColor(e+1, e+2)
        else:
            row.changeColor(e, e+1)
```
The `for(;~b;b--)` loop executes exactly twice, on the second pass it solves columns instead of rows. This is done by swapping `n` and `m`, and changing the values of `o` and `p` which are used in pointer arithmetic to address the array.
Here's a version that's ungolfed, with a test main, and prints the whole array after every move (press enter to step 1 turn):
```
#define s(x,y)b?x:y,b?y:x
#define A a[i*o+p
#define B A*e
f(m,n,a,i,j,o,p,b,c,d,e)int*a;{
    int t(x){
        printf("%d,%d:%d,%d;\n",s(i,e+x),s(i,e+1+x));
        getchar();
        printf("\n");
        for(int i2=0;i2<(b?m:n);i2++){
            for(int j2=0;j2<(b?n:m);j2++){
                printf("%d ",a[i2*(b?n:m)+j2]);
            }
            printf("\n");
        }
        printf("\n");
    }
    printf("\n");
    b=1;
    for(int i2=0;i2<(b?m:n);i2++){
        for(int j2=0;j2<(b?n:m);j2++){
            printf("%d ",a[i2*(b?n:m)+j2]);
        }
        printf("\n");
    }
    printf("\n");
    o=n;p=1;
    for(b=1;~b;b--){
        for(i=0;i<m;i++){
            c=0;
            for(j=0;j<n;j++) c+= a[i*o+p*j];
            d=0;
            d = (c*n)%3;
            for(j=0;j<n-1;j++) {
                while(a[i*o+p*j]!=d||a[i*o+p*j+p]!=d){
                    for(e=j;a[i*o+p*e]==a[i*o+p*e+p];e++);
                    if(e<=n-3 && a[i*o+p*e]==d 
                            && 2*(a[i*o+p*e+p]+a[i*o+p*(e+2)])%3==d){
                        a[i*o+p*e+p]=a[i*o+p*(e+2)]=d;
                        t(1);
                    }else{
                        a[i*o+p*e]=a[i*o+p*e+p] = 2*(a[i*o+p*e]+a[i*o+p*e+p])%3;
                        t(0);
                    }
                }
            }
        }
        o=m;p=m=n;n=o;o=1;
    }
}
main(){
    int g[11][11] = 
    {
        {0,2,1,2,2,1,0,1,1,0,2},
        {2,1,1,0,1,1,2,0,2,1,0},
        {1,0,2,1,0,1,0,2,1,2,0},
        {0,0,2,1,2,0,1,2,0,0,1},
        {0,2,1,2,2,1,0,0,0,2,1},
        {2,1,1,0,1,1,2,1,0,0,2},
        {1,0,2,1,0,1,0,2,2,1,2},
        {0,0,2,1,2,0,1,0,1,2,0},
        {1,2,0,1,2,0,0,2,1,2,0},
        {2,1,1,0,1,1,2,1,0,0,2},
        {0,2,1,0,1,0,2,1,0,0,2},
    };
    #define M 4
    #define N 7
    int grid[M][N];
    for(int i=0;i<M;i++) {
        for(int j=0;j<N;j++) {
            grid[i][j] = g[i][j];
        }
    }
    f(M,N,grid[0],0,0,0,0,0,0,0,0);
};
```
]