With some "random" changes the time for computation on my computer went down from 1534s to 67s. What I did is exclude some pieces I found suspicious. Therefore I introduced basically the function exclude and had to do some bookkeeping of already places pieces.

  ```haskell
  len :: V3 Int -> Float
  len (V3 x y z) = sqrt (x'*x' + y'*y' + z'*z')
                where
                    x' = fromIntegral x :: Float
                    y' = fromIntegral y :: Float
                    z' = fromIntegral z :: Float
  exclude :: Shape -> Shape -> Bool
  exclude [a1, _, _, d1, _]  [a2, _, _, d2, _] = len (a2 - a1) < 1.1
                                                && len (d2 -   d1) < 1.1
  ```

  ```haskell
  conflict :: Shape -> Shape -> [Shape] -> Bool
  conflict occs piece oldPieces = any (\voxel -> elem voxel piece) occs ||
                any (\o -> exclude o piece) oldPieces
  ```

  ```haskell
  subsolutionsSmart :: Int   -- ^ How many pieces should be in the subsolution we're looking for?
                  -> Shape -- ^ Unavailable voxels, already occupied by some piece
                  -> [Shape]
                  -> [[Shape]]
  subsolutionsSmart 0 _ _ = [[]]
  subsolutionsSmart n occupiedVoxels oldPieces = do
    let freeVoxel = pickFreeVoxel occupiedVoxels
    newPiece <- filter
                  (\piece -> elem freeVoxel piece &&
                             not (conflict occupiedVoxels piece oldPieces)
                  )
                  allValidPieces
    let updatedOccupiedVoxels = newPiece <> occupiedVoxels
    let updatedOldPieces = (newPiece:oldPieces)
    otherPieces <- subsolutionsSmart (n - 1) updatedOccupiedVoxels updatedOldPieces
    return $ newPiece : otherPieces
  ```

  ```haskell
  allSolutionsSmart :: [[Shape]]
  allSolutionsSmart = subsolutionsSmart 25 [] []
  ```