Better way to get Fisher Exact?

As you have found out yourself, the hard part is generating all possible matrices that have a specific row- and column-sum. I don't think that using Reduce for this is a good way to go if you are targeting an algorithm that always works. Consider this simple example with your approach

Clearly, the p-value should be 1. in this example. I have tried several random matrices and sometimes, Reduce doesn't give the result you are expecting in your setting.

This leads me to point number two: How could one tackle the problem of generating all matrices? The naive solution of generating all matrices of a certain dimension and filter out the ones we need is clearly too complex and will blow up. I haven't searched the net carefully, but on a first glance, I did not come over some clever implementation.

Although, what we have is the information about the col- and row-sums. Therefore, one way is to build all possible matrices row by row. We start with the first row: Its row-sum gives us the predicate if it is a correct row and the col-sums gives us boundaries for the iteration. Example, if we know the row-sum is 10, and the col-sums are say {5,4,6} then we could go through all possibilities

Do[
  (* check if sum of {i1,i2,i3} is 10 and return valid row *),
  {i1, 0, 5}, {i2, 0, 4}, {i3, 0, 6}
]

If we found a valid row r1, we know that the possibilies for r2 are not as many, because since the columns need to add up too, our current column sum changed to colSum - r1. We can do this step by step until we processed all rows and in the end, our current column sum needs to be 0.

Since we don't know upfront, which dimension and what column and row sums we get, we need to create the loop dynamically:

allRows[sum_Integer, colSums_List] := 
  Module[{createRows, $e, $iterators, $all},
   If[Min[colSums] < 0, Return[{}]];
   $iterators = 
    Table[{$e[i], 0, Min[colSums[[i]], sum]}, {i, Length[colSums]}];
   $all = Last@Reap@Do[With[{col = First /@ $iterators},
        If[Total[col] === sum,
         Sow[col]
         ]], Evaluate[Sequence @@ $iterators]];
   If[$all === {},
    {},
    First[$all]
    ]
   ];

Here are all valid rows of the example above

Building now all possible matrices (in a nested list) can be done by recursively calling allRows:

build[{{rowSum1_, rest___}, colSums_}, result_] := Module[{$rows},
   $rows = allRows[rowSum1, colSums];
   build[{{rest}, colSums - #}, {result, #}] & /@ $rows;
   ];
build[{{}, colSums_}, result_] := 
  If[Total[colSums] === 0, Sow[result]];

The thing that is left are some definitions you have yourself and creating the fisherTest that uses the created possible matrices to sum up all p-values

sumMatrix[mat_] := Total /@ {Transpose[mat], mat};

getP[mat_] := Module[{$sumMat = sumMatrix[mat], $tot},
  $tot = Total[First[$sumMat]];
  Times @@ (Factorial /@ 
       Flatten[$sumMat]) /(Factorial[$tot] (Times @@ (Factorial /@ 
          Flatten[mat]))) // N
  ]

fisherTest[mat_] := Module[{$poss, $pCut, $pVals},
  $poss = First@Last@Reap[build[sumMatrix[mat], {}]];
  $pCut = getP[mat];
  $pVals = With[{m = ArrayReshape[Flatten[#], Dimensions[mat]]},
      getP[m]] & /@ $poss;
  Total[Select[$pVals, # <= $pCut &]]
  ]

My fisherTest will be slower than yours, but if the possibilies are not too large, it will always return the correct result.

#Advanced

For the truly brave, I will add a version of allRows that uses in-time compiling which means that it creates dynamically a compiled function for a certain vector length n which can be called to find all valid rows. This compiling is only done once, so on the next run the function will be at hand and there is no compile-time.

With this function, I'm getting similar timings as the Reduce implementation:

ClearAll[allRows];
allRows[n_] := allRows[n] = Block[{iters, body, cs, e, sum, vec},
    iters = 
     Table[With[{i = i}, {Unique[e], 0, 
        Min[sum, Compile`GetElement[cs, i]]}], {i, n}];
    vec = List @@ iters[[All, 1]];
    (Compile[{{sum, _Integer, 0}, {cs, _Integer, 1}},
        Module[{result = Internal`Bag[Most[{0}], 1]},
         Do[If[#2 == sum, Internal`StuffBag[result, #1, 1]],
          ##3
          ];
         Internal`BagPart[result, All]
         ]] &)[vec, Total[vec], Sequence @@ iters]
    ];
allRows[sum_Integer, colSums_List] := 
 Partition[allRows[Length[colSums]][sum, colSums], Length[colSums]]

Quick test

m = RandomInteger[10, {3, 3}];

FisherExact[m] // AbsoluteTiming

During evaluation of {cutoff,7.67668*10^-6}
(* {0.319039, 0.00172301} *)

fisherTest[m] // AbsoluteTiming
(* {1.12482, 0.00172301} *)