anchor.diff:          integer             :: i, j, k
   
           ps%length = 0
  -        allocate(ps%val(size(p, 1), size(p, 2), 4**3))
  +        allocate(ps%s(size(p, 1), size(p, 2), 4**3))
           allocate(ps%d(3, 3, 3, 4**3))
   
           p1 = p
  
positive.code:         end forall
      end function
  
      function contains_p (a, p) result(b)
          integer, intent(in) :: a(:, :, :, :)
          integer, intent(in) :: p(:, :)
          logical             :: b
          integer             :: i, j
  
          loop: do i = 1, size(a, 4)
              do j = 1, size(p, 2)
                  if (a(p(1, j), p(2, j), p(3, j), i) == 0) then
                      cycle loop
                  end if
              end do
  
negative.code:         end do
      end function
  
      subroutine place(c1, p, id, c2, ok)
          integer, intent(in)  :: c1(3, 3, 3), p(:, :), id
          integer, intent(out) :: c2(3, 3, 3)
          logical, intent(out) :: ok
          integer              :: i, x, y, z
  
          c2 = c1
          do i = 1, size(p, 2)
              x = p(1, i)
              y = p(2, i)
              z = p(3, i)
              if (c1(x, y, z) /= 0) then
                  ok = .false.
                  return
              end if
              c2(x, y, z) = id
          end do
          ok = .true.
      end subroutine
  
      function all_puts(cube, piece, id) result(cubes)
          integer, intent(in)                   :: cube(3, 3, 3), piece(:, :), id
          type(list)                            :: cubes, rots
          integer                               :: i, x, y, z, rmax(3)
          integer, dimension(3, size(piece, 2)) :: rot, px, py, pz
          integer                               :: c1(3, 3, 3)
          logical                               :: ok
  
          cubes%length = 0
          !allocate(cubes%s(size(piece, 1), size(piece, 2), 3**3))
          allocate(cubes%d(3, 3, 3, 3**3))
  
          rots = all_rots(piece)
          do i = 1, rots%length
              rot = rots%s(:, :, i)
              rmax = maxval(rot, 2)
              do x = 0, 3 - rmax(1)
                  px = rot
                  px(1, :) = px(1, :) + x
                  do y = 0, 3 - rmax(2)
                      py = px
                      py(2, :) = py(2, :) + y
                      do z = 0, 3 - rmax(3)
                          pz = py
                          pz(3, :) = pz(3, :) + z
                          call place(cube, pz, id, c1, ok)
                          if (ok) then
                              cubes%length = cubes%length + 1
                              !cubes%s(:, :, cubes%length) = pz
                              cubes%d(:, :, :, cubes%length) = c1
                          end if
                      end do
                  end do
              end do
          end do
      end function
  
      recursive subroutine search(ps, n, cube, sols)
          integer, intent(in)       :: ps(:, :, :), n, cube(3, 3, 3)
          type(list), intent(inout) :: sols
          type(list)                :: puts
          integer                   :: i
  
          if (n == 0) then
              sols%length = sols%length + 1
              sols%d(:, :, :, sols%length) = cube
              return
          end if
  
          ! TODO Remove searchs with tiny islands.
          puts = all_puts(cube, ps(:, :, n), n)
          ! TODO Remove puts that are rotations of each other.
          do i = 1, puts%length
              call search(ps, n - 1, puts%d(:, :, :, i), sols)
          end do
      end subroutine
  
      subroutine read_pieces(vecs)
          integer, allocatable, intent(inout) :: vecs(:, :, :)
          character(len=20) :: name_
  