The basic idea is to fix the margins and simulate from the set of tables with those margins.
Consider the 2x2 table:
2 1
0 4
The margins are:
x x 3
x x 4
2 5
The possible tables with those margins are:
0 3 1 2 2 1
2 2 1 3 0 4
and their probabilities under the hypothesis of no association can be computed (consider the top left cell, say, and it will be hypergeometric).
Consequently we can simulate from that distribution over $2\times 2$ tables under the null and compute the distribution of any statistic we wish, and so obtain p-values in the usual fashion when sampling the null distribution*. The case for $r\times c$ tables is an extension of this but hopefully this is enough to get the idea. [There's some discussion of ways that $r \times c$ tables might be simulated here How to simulate effectiveness of treatment in R? and gung gives discussion of the situation if you don't have both margins fixed]
Note that in R this simulation is done using the algorithm of Patefield (1981) [1], as is explained in the help. (the function r2dtable will simulate for you, if you wanted to check the performance of the chi-square against some other statistic under fixed margins).
* it's also possible to generate all tables (if your tables aren't too big) and get an exact permutation test. Clever algorithms exist for just looking at some statistic for tables in the tail -- and doing so quite efficiently -- which makes it feasible to do exact tests for surprisingly large tables (to me at least, considering the scale of the combinatorial explosion).
[1] Patefield, W. M. (1981)
"Algorithm AS159. An efficient method of generating r x c tables with given row and column totals"
Applied Statistics, 30, 91-97.
chisq.test()you have problems with. As it stands, your question reads like "I don't want to read the paper, someone please read it and explain it to me", which is not a good fit to CV. – Stephan Kolassa May 24 '17 at 06:29