Vyxal, 17 bytes

L½D₀$↔↔'∩JƛĠṠꜝ;?⁼

Try it Online!

Average \$ O(2^{n^2}) \$ bruteforcer. I have not managed to get this to do a 4x4 testcase, so a 3x3 one's in the link above. Outputs all possible solutions, but as uniqueness is guaranteed in the test cases this is effectively outputting one solution in a list.

     ↔            # All combinations of
   ₀$             # [0, 1]
L½                # Of length input/2
      ↔           # All combinations of that
  D               # Of length input/2
                  # (effectively this gives all n/2-by-n/2 boolean grids
                  # Each cut into chunks of length 
       '          # Keep only those where
        ∩J        # Appending the transpose
          ƛ   ;   # And over each 
           ĠṠꜝ    # summing consecutive groups and removing zeroes
               ?⁼ # results in the input?

Haskell, 109 bytes

Disclaimer: this is derived from @ThreeFx's answer. I helped him to golf down his answer but he seems to have lost interest to include my last substantial improvements, so I post them as new answer.

import Data.List
n=mapM id
a#b=[x|x<-n$(n$" #"<$a)<$b,g x==a,g(transpose x)==b]
g=map$max[0].map length.words

Usage example: [[2],[3],[3],[3,1],[1]] # [[2],[3],[4],[2],[2]] -> [[" ## "," ### ","### ","### #"," #"]].

Brute force. Try all combinations of and #, split int chunks of #, count the lengths and compare to the input.

PHP, 751 833 (720) 753 724 726 710 691 680 682 bytes

I was eager to build a specialised sequence increment and try my cartesian generator once more;
but dropped the cartesian in favour of backtracking to solve the large puzzle faster.

$p=[];foreach($r as$y=>$h){for($d=[2-($n=count($h)+1)+$u=-array_sum($h)+$w=count($r)]+array_fill($i=0,$n,1),$d[$n-1]=0;$i<1;$d[0]+=$u-array_sum($d)){$o=$x=0;foreach($d as$i=>$v)for($x+=$v,$k=$h[$i];$k--;)$o+=1<<$x++;if(($s[$y]|$o)==$o){$p[$y][]=$o;$q[$y]++;}for($i=0;$i<$n-1&$d[$i]==($i?1:0);$i++);if(++$i<$n)for($d[$i]++;$i--;)$d[$i]=1;}}
function s($i,$m){global$c,$w,$p;for(;!$k&&$i[$m]--;$k=$k&$m<$w-1?s($i,$m+1):$k){for($k=1,$x=$w;$k&&$x--;){$h=$c[$x];for($v=$n=$z=$y=0;$k&&$y<=$m;$y++)$n=$n*($f=($p[$y][$i[$y]]>>$x&1)==$v)+$k=$f?:($v=!$v)||$n==$h[$z++];if($k&$v)$k=$n<=$h[$z];}}return$k?is_array($k)?$k:$i:0;}
foreach(s($q,0)as$y=>$o)echo strrev(sprintf("\n%0{$w}b",$p[$y][$o]));

• expects hints in arrays $r for row hints, $c for column hints and $s for square hints.
• throws invalid argument supplied for foreach if it finds no solution.
• to get the correct byte count, use a physical \n and remove the other two line breaks.

description

1) from row hints
generate possible rows that satisfy the square hints
and remember their counts for each row index.

2) backtrack over the row combinations:
If the combination satisfies the column hints, seek deeper or return successful combination,
else try next possibility for this row

3) print solution

The last golfing had severe impact on the performance;
but I removed profiling assignments for the final benchmarks.

Replace $n=$n*($f=($p[$y][$i[$y]]>>$x&1)==$v)+$k=$f?:($v=!$v)||$n==$h[$z++];
with if(($p[$y][$i[$y]]>>$x&1)-$v){$k=($v=!$v)||$n==$h[$z++];$n=1;}else$n++;
to undo the last golfing step.

examples

For the small example (17 to 21 around 12 8 7 6.7 5.3 ms), use

$r=[[2],[3],[3],[3,1],[1]];$c=[[2],[3],[4],[2],[2]];$s=[0,0,0,0,0];

for the christmas puzzle:

• killed my small home server with the old solution
• killed the browser with test outputs
• now solved in 50 37.8 45.5 around 36 seconds

put data from question to a file christmas.nonogram and use this code to import:

$t=r;foreach(file('christmas.nonogram')as$h)if('-'==$h=trim($h))$t=c;elseif('#'==$h){$t=s;$f=count($h).b;}else
{$v=explode(' ',$h);if(s==$t)for($h=$v,$v=0,$b=1;count($h);$b*=2)$v+=$b*array_shift($h);${$t}[]=$v;}

breakdown

$p=[];  // must init $p to array or `$p[$y][]=$o;` will fail
foreach($r as$y=>$h)
{
    // walk $d through all combinations of $n=`hint count+1` numbers that sum up to $u=`width-hint sum`
    // (possible `0` hints for $h) - first and last number can be 0, all others are >0
    for(
        $d=[2-
            ($n=count($h)+1)+               // count(0 hint)=count(1 hint)+1
            $u=-array_sum($h)+$w=count($r)  // sum(0 hint) = width-sum(1 hint)
        ]                           // index 0 to max value $u-$n+2
        +array_fill($i=0,$n,1)      // other indexes to 1
        ,$d[$n-1]=0;                // last index to 0
                                    // --> first combination (little endian)
        $i<1;   // $i:0 before loop; -1 after increment; >=$n after the last combination
        $d[0]+=$u-array_sum($d) // (see below)
    )
    {
        // A: create row (binary value) from 1-hints $h and 0-hints $d
        $o=$x=0;
        foreach($d as$i=>$v)
            for($x+=$v,$k=$h[$i];$k--;)
                $o+=1<<$x++;
        // B: if $o satisfies the square hints
        if(($s[$y]|$o)==$o)
        {
            $p[$y][]=$o;    // add to possible combinations
            $q[$y]++;       // increase possibility counter
        }
        // C: increase $d
            // find lowest index with a value>min
                // this loop doesn´t need to go to the last index:
                // if all previous values are min, there is nothing left to increase
        for($i=0;$i<$n-1&$d[$i]==($i?1:0);$i++);
        if(++$i<$n)             // index one up; increase $d if possible
            for($d[$i]++        // increase this value
            ;$i--;)$d[$i]=1;    // reset everything below to 1
            // adjust $d[0] to have the correct sum (loop post condition)
    }
}

// search solution: with backtracking on the row combinations ...
function s($i,$m)
{
    global $c,$w,$p;
    for(;
        !$k // solution not yet found
        &&$i[$m]    // if $i[$m]==0, the previous iteration was the last one on this row: no solution
            --;     // decrease possibility index for row $m
        $k=$k&$m<$w-1? s($i,$m+1) : $k      // if ok, seek deeper while last row not reached ($m<$w-1)
    )
    {
        // test if the field so far satisfies the column hints: loop $x through columns
        for($k=1,$x=$w;$k&&$x--;)   // ok while $k is true
        {
            $h=$c[$x];
            // test column hints on the current combination: loop $y through rows up to $m
            for($v=$n=$z=   // $v=temporary value, $n=temporary hint, $z=hint index
                $y=0;$k&&$y<=$m;$y++)
                // if value has not changed, increase $n. if not, reset $n to 1
                // (or 0 for $k=false; in that case $n is irrelevant)
                $n=$n*  
                    // $f=false (int 0) when value has changed, true (1) if not
                    ($f=($p[$y][$i[$y]]>>$x&1)==$v)
                    +$k=$f?:    // ok if value has NOT changed, else
                        ($v=!$v)        // invert value. ok if value was 0
                        || $n==$h[$z    // value was 1: ok if temp hint equals current sub-hint
                        ++]             // next sub-hint
                ;
            // if there is a possibly incomplete hint ($v==1)
            // the incomplete hint ($n) must be <= the next sub-hint ($c[x][$z])
            // if $n was <$h[$z] in the last row, the previous column hints would not have matched
            if($k&$v)$k=$n<=$h[$z];
        }
        // ok: seek deeper (loop post condition)
        // not ok: try next possibility (loop pre condition)
    }
    return$k?is_array($k)?$k:$i:0;  // return solution if solved, 0 if not
}

// print solution
foreach(s($q,0)as$y=>$o)echo strrev(sprintf("\n%0{$w}b",$p[$y][$o]));

Brachylog, 70 69 bytes

[R:C]hlL~l:L:1f=.:3aR,.z:3aC,
tL,?he##ElL,E:2a
.<2,_1<
@b:4f:la
e.h1,

This takes a list of two lists (first the rows indicators, then the column ones). Each indicator is itself a list (for situtions like [3,1] on one row).

This version takes about 3 minutes to solve the 5 by 5 example of the challenge.

Way more efficient version, 91 bytes

[R:C]hlL~l:L:1f:Cz:3az:Rz:3a=.:4aR,.z:4aC,
tL,?he##ElL,E:2a
.<2,_1<
:+a#=,?h
@b:5f:la
e.h1,

Try it online!

This one is not complete brute force: the only difference is that this one imposes constraints on the values of cells such that the number of 1s in each row and column matches the numbers given as indicators in the input. The only brute force part is then in finding the one grid with those constraints for which the "blocks" of 1s match what is given as indication.

This one takes about 0.05 seconds on the 5 by 5 example of the challenge. This is still way too slow for the bonus case, as I have no idea how to express the blocks of ones separated by one or more zeroes in terms of constraints.

Explanation

I will explain below the 93 bytes version. The only difference between the two is the call to predicate 3 which doesn't exist in the 70 bytes version, and the numbering of the predicates (since there is one less).

• Main predicate:

  [R:C]     Input = [R, C]
  hlL       The length of R is L
  ~l        Create a list of length L
  :L:1f     Each element of that list is a sublist of length L with cells 0 or 1 (Pred 1)
            %%% Part unique to the 93 bytes version
  :Cz       Zip the rows of the list of lists with C
  :3a       The sum of 1s in each row is equal to the sum of the indicators (Pred 3)
  z         Transpose
  :Rz       Zip the columns of the list of lists with R
  :3a       The sum of 1s in each column is equal to the sum of the indicators (Pred 3)
            %%%
  =.        Assign values to the cells of the list of lists which satisfy the constraints
  :4aR,     The blocks of 1s must match the indicators on rows
  .z        Transpose
  :4aC,     The blocks of 1s must match the indicators on columns
  

• Predicate 1: Forces the rows to have a specific length, and that each cell is 0 or 1.

  tL,       L is the length given as second element of the input
  ?he       Take an element from the list
  ##ElL,    That element E is itself a list of length L
  E:2a      The elements of E are 0s and 1s (Pred 2)
  

• Predicate 2: Constrain a variable to be either 0 or 1

  .<2,      Input = Output < 2
  _1<       Output > -1
  

• Predicate 3: The sum of 1s in a list must be equal to the sum of indicators (e.g. if the indicator is [3:1] then the list must have sum 4)

  :+a       Sum the elements of the list and sum the indicator
  #=,       Both sums must be equal
  ?h        Output is the list
  

• Predicate 4: Check that blocks of 1s match the indicator

  @b        Split the list in blocks of the same value
  :5f       Find all blocks of 1s (Pred 5)
  :la       The list of lengths of the blocks results in the indicator (given as output)
  

• Predicate 5: True for blocks of 1s, false otherwise

  e.        Output is an element of the input
    h1,     Its first value is 1
  

Pyth, 91 72 71 bytes

D:GHdRq@@QdG.nCf.)TrH8V^,01^hQ2=TcNhQ=Y1VhQ=*Y*:H@TH1:H@CTH2)IYjbmjkdTb

A program that takes input of a list of the form [size, [horizontal clues], [vertical clues]] where each clue is a list of integers (empty clues are the empty list, []), and prints every solution, newline separated, in the form of a binary grid where 1 is shaded and 0 is unshaded.

This is a brute force, so is roughly O(2^n^2). It starts taking a very long time for larger puzzles, but will solve any arbritrary size given sufficient time.

Try it online

How it works

The program generates every possible layout by taking the repeated Cartesian product of [0, 1] with a length equal to size^2. This is then split into chunks, giving a list for each horizontal line. Each line is run-length encoded, filtered by the presence of 1 and flattened, leaving the clue for that line. This is then checked against the input. The above process is repeated for the transpose of the chunks, checking the vertical lines. If there is a hit, each chunk is concatenated, and the concatenated chunks are joined on newlines and implicitly printed, with a trailing newline.

D:GHdRq@@QdG.nCf.)TrH8V^,01^hQ2=TcNhQ=Y1VhQ=*Y*:H@TH1:H@CTH2)IYjbmjkdTb  Program. Input: Q
                            hQ                                           Q[0], size
                           ^  2                                          Square
                        ,01                                              [0, 1]
                       ^                                                 Cartesian product
                      V                                     )            For N in the Cartesian product:
                                 cNhQ                                    Split N into Q[0] chunks
                               =T                                        Assign that to T
                                     =Y1                                 Y=1
                                        VhQ                              For H in range [0, Q[0]-1]:
D:GHd                                                                     def :(G, H, d)
                   rH8                                                     Run-length-encode(H)
               f.)T                                                        Filter by presence of 1 in character part
            .nC                                                            Transpose and flatten, giving the clue
       @@QdG                                                               Q[d][G], the relevant input clue
     Rq                                                                    Return clue==input clue
                                               :H@TH1                     :(H, T, 1)
                                                     :H@CTH2              :(H, transpose(T), 2)
                                           =*Y*                           Y=Y*product of above two
                                                             IY           If Y:
                                                                 mjkdT     Conacatenate each element of T
                                                               jb          Join on newlines
                                                                      b    Add a newline and implicitly print

Thanks to @Pietu1998 for some tips

Haskell, 242 230 201 199 177 163 160 149 131 bytes

import Data.Lists
m=map
a#b=[x|x<-m(chunk$length b).mapM id$[0,1]<$(a>>b),g x==a,g(transpose x)==b]
g=m$list[0]id.m sum.wordsBy(<1)

Finally under 200 bytes, credit to @Bergi. Huge thanks to @nimi for helping almost halving the size.

Wow. Almost at half size now, partly because of me but mainly because of @nimi.

The magic function is (#). It finds all solutions of a given nonogram.

This is able to solve all cases, but may be super slow, since it's complexity is about O(2^(len a * len b)). A quick benchmark revealed 86GB allocated for a 5x5 nonogram.

Fun fact: It works for all nonograms, not only square ones.

How it works:

• a#b: Given lists of lists of integers which represent the number of squares, generate all grids (map(chunk$length b).mapM id$a>>b>>[[0,1]]) and filter the results to keep only the valid ones.
• g: Given a potential nonogram it sums the runs of 1's horizontally.

Javascript (ES6), 401 386 333 bytes

This is an early attempt. It's not very efficient but I was curious to test a solution using regular expressions on the binary representation of the rows and columns.

For example, it will translate the clue [3,1] into the following regular expression:

/^0*1{3}0+1{1}0*$/

Right now, this version is not taking square clues into account. I will probably add this later.

Code

(c,r)=>{W=c.length;w=[];S=0;M=(n,p)=>eval(`/^0*${p.map(v=>`1{${v}}`).join`0+`}0*$/`).exec(n);R=(y,i=0)=>S||(w[y]=r[y][i],y+1<W?R(y+1):c.every((c,y)=>(n=0,w.map((_,x)=>n+=w[W-1-x][y]),M(n,c)))&&(S=w.join`
`),r[y][i+1]&&R(y,i+1));r=r.map(r=>[...Array(1<<W)].map((_,n)=>((1<<30)|n).toString(2).slice(-W)).filter(n=>M(n,r)));return R(0)}

Output

The solution is displayed in binary format. Such as:

00110
01110
11100
11101
00001

Test

This is a simple test on the example grid.

let f =
(c,r)=>{W=c.length;w=[];S=0;M=(n,p)=>eval(`/^0*${p.map(v=>`1{${v}}`).join`0+`}0*$/`).exec(n);R=(y,i=0)=>S||(w[y]=r[y][i],y+1<W?R(y+1):c.every((c,y)=>(n=0,w.map((_,x)=>n+=w[W-1-x][y]),M(n,c)))&&(S=w.join`
`),r[y][i+1]&&R(y,i+1));r=r.map(r=>[...Array(1<<W)].map((_,n)=>((1<<30)|n).toString(2).slice(-W)).filter(n=>M(n,r)));return R(0)}

console.log(f(
  [[2],[3],[4],[2],[2]],
  [[2],[3],[3],[3,1],[1]]
));