Here's a method I made for detecting matching columns. It checks each pair of columns row by row, and does an early out as soon as a mismatch is found (which knocks out 12 of the 28 iterations after the first row). Stack hackage is winning!

; Returns PE if there is a match
; Does 28 iterations, interrupts must be disabled
ld d,matrix >> 8
ld hl,0
add hl,sp
ld sp,masktable
loop:
pop bc
ld e,b
ld a,(de) \ and c \ ret po
inc e
ld a,(de) \ and c \ ret po
inc e
ld a,(de) \ and c \ ret po
inc e
ld a,(de) \ and c \ ret po
inc e
ld a,(de) \ and c \ ret po
inc e
ld a,(de) \ and c \ ret po
inc e
ld a,(de) \ and c \ ret po
inc e
ld a,(de) \ and c \ ret po
end:
ld sp,hl
ret

masktable:
.dw (matrix&$FF)*256+%00000011, loop
.dw (matrix&$FF)*256+%00000101, loop
.dw (matrix&$FF)*256+%00000110, loop
.dw (matrix&$FF)*256+%00001001, loop
.dw (matrix&$FF)*256+%00001010, loop
.dw (matrix&$FF)*256+%00001100, loop
.dw (matrix&$FF)*256+%00010001, loop
.dw (matrix&$FF)*256+%00010010, loop
.dw (matrix&$FF)*256+%00010100, loop
.dw (matrix&$FF)*256+%00011000, loop
.dw (matrix&$FF)*256+%00100001, loop
.dw (matrix&$FF)*256+%00100010, loop
.dw (matrix&$FF)*256+%00100100, loop
.dw (matrix&$FF)*256+%00101000, loop
.dw (matrix&$FF)*256+%00110000, loop
.dw (matrix&$FF)*256+%01000001, loop
.dw (matrix&$FF)*256+%01000010, loop
.dw (matrix&$FF)*256+%01000100, loop
.dw (matrix&$FF)*256+%01001000, loop
.dw (matrix&$FF)*256+%01010000, loop
.dw (matrix&$FF)*256+%01100000, loop
.dw (matrix&$FF)*256+%10000001, loop
.dw (matrix&$FF)*256+%10000010, loop
.dw (matrix&$FF)*256+%10000100, loop
.dw (matrix&$FF)*256+%10001000, loop
.dw (matrix&$FF)*256+%10010000, loop
.dw (matrix&$FF)*256+%10100000, loop
.dw (matrix&$FF)*256+%11000000, end
Edit: Changed return value to be the P/V flag, saving a byte

Edit2: Here's a more size-optimized version that is interrupt-compatible.
; Returns NC if there is a match
; Does C(8,2) iterations
ld b,$80
ld d,matrix >> 8
ld hl,innerloop
outerloop:
srl b
jr c,end
ld c,b
scf
innerloop:
rl c
jr c,outerloop
ld e,matrix & $FF
push hl
ld a,(de) \ and c \ ret po
inc e
ld a,(de) \ and c \ ret po
inc e
ld a,(de) \ and c \ ret po
inc e
ld a,(de) \ and c \ ret po
inc e
ld a,(de) \ and c \ ret po
inc e
ld a,(de) \ and c \ ret po
inc e
ld a,(de) \ and c \ ret po
inc e
ld a,(de) \ and c \ ret po
pop hl
end:
ret
Edit3: Here's an even faster hardcore version!
; Returns PE if there is a match
; Does 28 iterations, interrupts must be disabled
ld d,matrix >> 8
ld l,matrix & $FF
ld (saveSP),sp
ld sp,masktable
loop1:
pop bc
ld e,l
ld a,(de) \ and c \ ret po
inc e
ld a,(de) \ and c \ ret po
inc e
ld a,(de) \ and c \ ret po
inc e
ld a,(de) \ and c \ ret po
inc e
ld a,(de) \ and c \ ret po
inc e
ld a,(de) \ and c \ ret po
inc e
ld a,(de) \ and c \ ret po
inc e
ld a,(de) \ and c \ ret po
end:
ld sp,(saveSP)
ret

loop2:
ld e,l
ld a,(de) \ and b \ ret po
inc e
ld a,(de) \ and b \ ret po
inc e
ld a,(de) \ and b \ ret po
inc e
ld a,(de) \ and b \ ret po
inc e
ld a,(de) \ and b \ ret po
inc e
ld a,(de) \ and b \ ret po
inc e
ld a,(de) \ and b \ ret po
inc e
ld a,(de) \ and b \ ret po
ld sp,(saveSP)
ret

masktable:
.dw %0000001100000101, loop2, loop1
.dw %0000011000001001, loop2, loop1
.dw %0000101000001100, loop2, loop1
.dw %0001000100010010, loop2, loop1
.dw %0001010000011000, loop2, loop1
.dw %0010000100100010, loop2, loop1
.dw %0010010000101000, loop2, loop1
.dw %0011000001000001, loop2, loop1
.dw %0100001001000100, loop2, loop1
.dw %0100100001010000, loop2, loop1
.dw %0110000010000001, loop2, loop1
.dw %1000001010000100, loop2, loop1
.dw %1000100010010000, loop2, loop1
.dw %1010000011000000, loop2, end

Ah yes this thread, I sort of forget to ever update it.

Ok, in a sort of high level overview, the best algorithm I know for generating instances is this:
 - generate a fully filled in board (using that pseudocode algorithm from the opening post, there are some issues)
 - try some random subsets of those bits as instances, try to solve them to find out whether they're solvable instances
 - from all solvable instances, take the hardest one

So, how to select a fully filled board.
In 8x8 puzzles, there are only 4111116 such configurations (much fewer if you don't count symmetric solutions, there are many sorts of symmetry). For other sizes, the number is similarly "unexpectedly low" (compared to the full 2^size)
Select a random number up to the number of possible configurations, generate boards in lexicographical order until you reach that index (alternatively: store them all in the program, that array is pretty small). That's surprisingly fast if you use the following tricks:
 - consider only valid rows. There's no point in doing it bit by bit. Do a whole row at the time. There are only 34 possible rows of length 8.
 - consider only rows that don't immediately violate the no-3-in-a-row rule vertically.
 - consider only rows that don't immediately violate the 4-zeroes-per-column rule, using the column-counting approach from the pseudocode in my first post. The last row is implied by the column counts.
Using those 3 tricks, the only constraint you can violate is no-two-equal-columns, that's why I concentrated on that. That's a pretty annoying check to do.
Tip: you can use (797 * row) >> 8 to map the 34 valid length-8-rows to non-overlapping indexes in [0, 64]. Using that, you can check for duplicate columns by transposing the board and then ORing together 1 << row_to_index(rows in the transposed board). If that has fewer than 8 set bits, at least one column appears twice.
Still it's probably not fast enough (may take a second, or several on slower hardware) and there a big speed difference between high indexes and low indexes. You can cut that time down enormously by storing an array of how many boards there are if you start with a certain row. That array is only 34 (for 8x8) ints long, and is in fact perfectly symmetrical around the middle, so you can store only half of it if you want (not that it matters..)
Changing the indexes of the boards slightly, leaving perfect lexicographical order behind, the symmetry around the middle is due to the fact that the last half of the boards is the same as the first half, but both reversed in order an inverted in its bits. You can cut it in half there, and when generating a random board, choose which half and which index within that half.
Similar trickery can be done with the second row, cutting the time it takes down so much that you can afford to do a binary search over all the boards (you don't need that, I just needed it because I wanted to find the indexes of the rotated and mirrored boards (the index of the inverted board can be found trivially due to the symmetry that I already mentioned) - I was building a database with Very Hard Instances and thought I might as well put all symmetrically equivalent boards in it as well).

Then, how to solve an instance.
Filling the cells implied by the no-3-in-a-row rule is easily done with some bitmath.
The rest is more complicated, but it can still be done iteratively. No need to do a recursive search.
For this, I use a lookup table to handle all cases where cell(s) can be filled in just to do knowledge about the row it's in. That handles everything you can fill in using any combination of no-3-in-a-row and 4-zeroes-per-row.
That doesn't do everything, obviously. First, it disregards columns, but that's easy: just transpose the board and try again.
Second, it disregards things you can know from the no-equal-rows rule. For that I enumerate all possible rows, test them against the currently filled in cells in the row, then against any surrounding rows for the no-3-in-a-row rule, then against the currently already completed rows, and then against the column-counts. If it passes all those tests, then it's possible to fill in this row here. All possible rows for one position are combined to see if they have any bits that are always the same, those bits are "filled in".
Then the board is transposed and the same is done for the columns. Until either there is no change or the board is completely filled. If it's completely filled, then the instance is solvable.

Ok, so now, how to test how hard an instance is.
For this, I actually generate the specific moves that you would play. A "move" is when you have some reason to fill in a cell. Generating actual moves is harder and especially more verbose than implicitly having moves as in the solver. And there's an other problem: when there are multiple moves, not all of them are created equal - some lead to easier total paths to the end. To avoid scoring puzzles too high, I use A* to find the easiest path to the solution (A* yes so it has a heuristic - a very bad one though, the number of cells currently not filled in. That only works if you score the hardness of moves in a certain way. Be very careful that you don't accidentally make a non-admissible heuristic, I spent several days debugging that).