[Question]
  [
Write a function (using as few bytes as possible) that takes a bi-dimensional array of any number of columns and rows in which:
* `0` represents empty block,
* `1` represents snake block.
The function must return the number of possible paths the snake traveled.
**Example 1:**
Input:
```
[
  [1,1,1,1,1],
  [0,0,0,0,1],
  [0,0,0,0,1],
]
```
Output: `2`
In the example above, the function will return `2` because the answer is either one of:
[![enter image description here](https://i.stack.imgur.com/hzXbC.png)](https://i.stack.imgur.com/hzXbC.png)
**Example 2:**
Input:
```
[
  [1,1,1,1],
  [0,0,1,1],
  [0,0,1,1],
]
```
Output: `6`
In this example the function will return `6` because the answer is either one of:
[![enter image description here](https://i.stack.imgur.com/Myx2r.png)](https://i.stack.imgur.com/Myx2r.png)
**Note:** 
When assessing the input, you can assume that:
* The arrays representing columns will always have the same sizes (so the arrays are rectangular);
* There exists at least 1 valid path;
* The snake cannot walk through the edges (as can happen in some versions of snake);
* The snake will always have at least 2 blocks;
* The snake cannot move diagonally;
* The paths are directed. (so, two paths ending at different positions but otherwise looking exactly the same are not the same path, it’ll add up to the total)
  
[Answer]
# [Wolfram Language (Mathematica)](https://www.wolfram.com/wolframscript/), 16+83=99 bytes
Library import statement (16 bytes):
```
<<Combinatorica`
```
Actual function body (83 bytes):
```
Length@HamiltonianCycle[MakeGraph[#~Position~1~Join~{1>0},##||Norm[#-#2]==1&],All]&
```
[Try it online!](https://tio.run/##bcixCsIwEIDhvU9RCHSKYJxbqXRQRKV7KXiW1B4mF4m3SNu8erTgKP/0/RZ40BYYO4gxzytnb0jAzn/HNemLk6Y7D@UBLBp2hEDVuzO6OcND7z08h0aE2r2Q0VFQ4eiQwqi261kKMU0X520jVmLTFoXKWrkzps1i7ZG47MsxSdNRyV@zXLiWS3@YzPED "Wolfram Language (Mathematica) – Try It Online")
---
Note that the question just ask for the number of Hamiltonian path in the graph.
However, (for some reason) the `HamiltonianPath` function doesn't really work with directed graph ([example](https://tio.run/##fY6xCsIwFEX3/kbBKUWaro1UHOqgUqVOIeCzfaTBJC0xLYj47THg7nrOgXsN@AENeNVBCOdZoedluRvNXVnwo4v4JhKZs9qpvnYwDZySQiQLa@dJ45NfwEpcV@9IPzGk7AgP/IULOY3O8DRLqWAsX0Vd/NGkfU2Yba62Vw47j71IGqes5we00g/VHozSfrQKbBNP85RstRaruC1zIimR8UAIXw)). So, I used the workaround described in [this Mathematica.SE](https://mathematica.stackexchange.com/questions/105196) question:
* Add an vertex (called `True`) that is connected to all other vertices.
* Count the number of Hamiltonian cycle on the resulting graph.
The graph is constructed using `MakeGraph` (annoyingly there are no directly equivalent built-in), using the boolean function `##||Norm[#-#2]==1&`, which returns `True` if and only if one of the arguments is `True` or the distance between the two vertices are `1`.
---
`Tr[1^x]` cannot be used instead of `Length@x`, and `<2` cannot be used instead of `==1`.
---
`HamiltonianPath` can be used if the graph is undirected, with the function body takes **84** bytes (exactly 1 byte more than the current submission):
```
Length@HamiltonianPath[MakeGraph[#~Position~1,Norm[#-#2]==1&,Type->Undirected],All]&
```
[Try it online!](https://tio.run/##hY3BSsNAEIbv@xQDgaJlC8ZzIxEP9qASQU8x4NhMm6HZ2bA7QaSkrx6j9SIW5L/M9898jENtyKHyGsfHnknL5fLGuzcWVB@m9rWCszl8r4Aj9JFqUA@x77pAMcKkwzsGYdnC/NxsspfxjmSrTb5Cx616YZRi@lLe445uA3ZNmRwKH1nZyyG1Dz64Mlkkl1WWpTP79NHR4upZag60Vqore9221WwsAovmm3xvAPap/clgv/DCHvMHzWDM/@IJPCUeb34NZhg/AQ "Wolfram Language (Mathematica) – Try It Online")
[Answer]
# JavaScript (ES6), ~~154~~ 134 bytes
```
m=>m.map((r,Y)=>r.map(g=(_,x,y,r=m[y=1/y?y:Y])=>r&&r[x]&&[-1,0,1,2].map(d=>r[r[x]=0,/1/.test(m)?g(_,x+d%2,y+~-d%2):++n,x]=1)),n=0)|n/4
```
[Try it online!](https://tio.run/##jY3LCoMwEEX3foWbSoKjJtKVEP0OCaGIL1pMlCjFQOmvp8ZdaQWZxXC498w8qmc11/o@LZEam9Z2zEqWy1hWE0IaSsxyvUPP0A1WMKCZ5IbRxBQmK4XLg0DzVQQBjygQoJCK3Wi2iLuEEUhoEi/tvCCJi94dCptLCiZ8R9vGWRgq2HoUY1CM4JdKrrYe1TwObTyMPeoQ93yfU9hHgAMCbn5AYOydUf/AgXrQPqN@fT2ATbUf "JavaScript (Node.js) – Try It Online")
## How?
### Method
Starting from each possible cell, we flood-fill the matrix, clearing all cells on our way. Whenever the matrix contains no more **1**'s, we increment the number **n** of possible paths.
Each valid path is counted 4 times because of the direction chosen on the last cell, which actually doesn't matter. Therefore, the final result is **n / 4**.
### Recursive function
Instead of calling the recursive function **g()** from the callback of the second **map()** like this...
```
m=>m.map((r,y)=>r.map((_,x)=>(g=(x,y,r=m[y])=>...g(x+dx,y+dy)...)(x,y)))
```
...we define the recursive function **g()** *directly as the callback* of **map()**:
```
m=>m.map((r,Y)=>r.map(g=(_,x,y,r=m[y=1/y?y:Y])=>...g(_,x+dx,y+dy)...))
```
Despite the rather long formula `y=1/y?y:Y` which is needed to set the initial value of **y**, this saves 2 bytes overall.
### Commented code
```
m =>                           // given the input matrix m[][]
  m.map((r, Y) =>              // for each row r[] at position Y in m[][]:
    r.map(g = (                //   for each entry in r[], use g() taking:
      _,                       //     - the value of the cell (ignored)
      x,                       //     - the x coord. of this cell
      y,                       //     - either the y coord. or an array (1st iteration),
                               //       in which case we'll set y to Y instead
      r = m[y = 1 / y ? y : Y] //     - r = the row we're currently located in
    ) =>                       //       (and update y if necessary)
      r && r[x] &&             //     do nothing if this cell doesn't exist or is 0
      [-1, 0, 1, 2].map(d =>   //     otherwise, for each direction d,
        r[                     //     with -1 = West, 0 = North, 1 = East, 2 = South:
          r[x] = 0,            //       clear the current cell
          /1/.test(m) ?        //       if the matrix still contains at least one '1':
            g(                 //         do a recursive call to g() with:
              _,               //           a dummy first parameter (ignored)
              x + d % 2,       //           the new value of x
              y + ~-d % 2      //           the new value of y
            )                  //         end of recursive call
          :                    //       else (we've found a valid path):
            ++n,               //         increment n
          x                    //       \_ either way,
        ] = 1                  //       /  do r[x] = 1 to restore the current cell to 1
      )                        //     end of map() over directions
    ),                         //   end of map() over the cells of the current row
    n = 0                      //   start with n = 0
  ) | n / 4                    // end of map() over the rows; return n / 4
```
[Answer]
# [Jelly](https://github.com/DennisMitchell/jelly), ~~12~~ 11 bytes
```
ŒṪŒ!ạƝ€§ÐṂL
```
[Try it online!](https://tio.run/##y0rNyan8///opIc7Vx2dpPhw18Jjcx81rTm0/PCEhzubfP4fbgfy/v@P5ormUlCINtSBwlgdENdABwKxcIEIWQdCARZOLFcsAA "Jelly – Try It Online")
---
## Explanation.
```
ŒṪ               Positions of snake blocks.
  Œ!             All permutations.
                 For each permutation:
    ạƝ€             Calculate the absolute difference for each neighbor pair
       §            Vectorized sum.
                 Now we have a list of Manhattan distance between snake
                    blocks. Each one is at least 1.
        ÐṂL      Count the number of minimum values.
                    Because it's guaranteed that there exists a valid snake,
                    the minimum value is [1,1,1,...,1].
```
[Answer]
# [Python 2](https://docs.python.org/2/), ~~257~~ ~~246~~ ~~241~~ ~~234~~ ~~233~~ ~~227~~ ~~214~~ 210 bytes
```
lambda b:sum(g(b,i,j)for j,l in e(b)for i,_ in e(l))
e=enumerate
def g(b,x,y):d=len(b[0])>x>-1<y<len(b);c=eval(`b`);c[d*y][d*x]=0;return d and b[y][x]and('1'not in`c`or sum(g(c,x+a,y)+g(c,x,y+a)for a in(1,-1)))
```
[Try it online!](https://tio.run/##dY7BbsMgEETvfMXeAvVGMj0mIT9CUQ0Gt44wjiiuzNe7xK5UqXW1Epq3zGjnntP7GJ6XTrwsXg/GajCnj2mgb9RgjzfWjRFu6KEP4KhZscfXDT1jxAkXpsFFnRyxroNHbsbMTlZ4F6iRtWLX@Xrkl3xZF@zcCvepPW1MU7S0T1mVZ1aiPkeXphjAgg4WjCwfsyqSHvghjKkcbdqmFNj6tThXulyqVom50ms7XWyU45EzxpZ77EOCjkoCIDl@j8IH1rjNDipGyF70x7kDf0L/@H6Fli8 "Python 2 – Try It Online")
---
Saved
* -8 bytes, thanks to Kevin Cruijssen
* -14 bytes, thanks to user202729
[Answer]
# Python 2, 158 bytes
```
E=enumerate
g=lambda P,x,y:sum(g(P-{o},*o)for o in P if x<0 or abs(x-o[0])+abs(y-o[1])<2)+0**len(P)
lambda L:g({(x,y)for y,r in E(L)for x,e in E(r)if e},-1,0)
```
[Try it online!](https://tio.run/##bYw9C4MwEIb3/Iob7/SE2KGD2LFbB3dxsDRaQRNJLUTE325jSqGUcje8z3084zzdjT5s2/mk9HNQtp6UaE99PVxvNRTseM4ezwFbLJLFrBwZaowFA52GAroGXC7BD@rrA11iSllRvOfZ57Si/ECxjKJeaSxINB/vJWtxQe8OspntrjvjJaBj9UZL3q9WTlKWtI220xM2WAqAMuVQFe8gWf6DihiOJH7@YO/3qU/yG8ImvKVH2l4) 
[Answer]
# [Haskell](https://www.haskell.org/), 187 155 bytes
```
r=filter
l=length
(a,b)?(x,y)=abs(a-x)+abs(b-y)==1
l#x=sum[p!r(/=p)l|p<-x]
p![]=1
p!l=l#r(p?)l
f x|l<-[(i,j)|i<-[0..l x-1],j<-[0..l(x!!0)-1],x!!i!!j>0]=l#l
```
[Try it online!](https://tio.run/##dY/BbsMgEETvfAUoPYAKrp1cQ/whiANW7AZn4yLsSETyvzubOFIjlbKXmbc7I3Fy47kFWJaoOw9TGwloaIfv6US4k42oeZI3oV0zcqeS@HyIRiHRFYFN0uP1YgKL/EsHAXPYq2RJYMbiOjCs2kQeagGko2mGvTLcy17MHlVZFECTqqzsX44nxkrxICg8Y/2htNgAS@dbOFZUU0MoPlPJ11i5glKukwHoLSHPhu2fhvfzjH3L7n6z/xxnsxfnB0wef57LEP0w0Q/a0fVHObjNwd1yBw "Haskell – Try It Online")
[Answer]
# [APL (Dyalog Extended)](https://github.com/abrudz/dyalog-apl-extended), 28 bytes
```
+/×/1=2(1⊥∘|-)/x[⌂pmat≢x←⍸⎕]
```
[Try it online!](https://tio.run/##SyzI0U2pTMzJT9dNrShJzUtJTfn/qKM97b@2/uHp@oa2RhqGj7qWPuqYUaOrqV8R/ainqSA3seRR56KKR20THvXueNQ3NRak4X8a16O2iUCeV7C/n4J6dLShDhTG6kQb6EAgKjtWnQuHJqgyVFasOgA "APL (Dyalog Extended) – Try It Online")
A full program that accepts a matrix.
]