Parity of $\lfloor 1/(x y) \rfloor$ not equally distributed

You're dividing the square $S = \{ (x,y) \colon 0 < x < 1, 0 < y < 1\}$ into two regions according to the parity of $\lfloor 1/(xy) \rfloor$, separated by the segments of the hyperbolas $xy = 1/n$ ($n=2,3,4,\ldots$) contained in $S$. There's no reason to expect that the two regions have the same area. If I did this right, the area between the $n$-th hyperbola and the top right corner of the square is $$ A(n) := 1 - \frac{1 + \log n}{n} $$ so the discrepancy between odd and even values of $\lfloor 1/(xy) \rfloor$ is $$ (A(2)-A(1)) - (A(3)-A(2)) + (A(4)-A(3)) - + \cdots $$ which is numerically $0.066556553635\ldots$ according to the gp calculation

A(n) = 1 - (log(n)+1)/n
sumalt(n=1, (-1)^n*(A(n)-A(n+1)))

So we expect about 53.33% odd and 46.67% even values, which seems consistent with your experiment.

P.S. Using a formula I found in MO Question 140547, I gather that this number $0.066556553635\ldots$ has the closed form $$ (\log 2)^2 + \bigl(2 (1 - \gamma) \log 2\bigr) - 1, $$ where $\gamma$ is Euler's constant $0.5772156649\ldots$.

P.P.S. I see that I didn't address the end of the original question: "If I use the ceiling function instead of the floor, the bias reverses, [...] if I round to the nearest integer instead, then about 48% are odd." The first part is clear because changing $\lfloor 1/(xy) \rfloor$ to $\lceil 1/(xy) \rceil$ switches even and odd values (except in the negligible case that $1/(xy)$ is an exact integer). For the nearest-integer function, the discrepancy between odd and even values is $$\bigl(A(3/2)-A(1)\bigr) - \bigl(A(5/2)-A(3/2)\bigr) + \bigl(A(7/2)-A(5/2)\bigr) - + \cdots $$ which evaluates numerically to $-0.03500998166\ldots$ (using sumalt in gp as before), which again is consistent with observation (48.25% odd, 51.75% even). There's still a "closed form" for this discrepancy, but more complicated: $$ -3 + 4 \log(2) + \pi \bigl(1 + \log(\pi/2) - \gamma - 4 \log\Gamma(3/4) \bigr). $$ This requires evaluation of $\log(1)/1 - \log(3)/3 + \log(5)/5 - \log(7)/7 + - \cdots$, which can be achieved by differentiating the functional equation for the Dirichlet L-function $L(s,\chi_4) = 1 - 3^{-s} + 5^{-s} - 7^{-s} + - \cdots$ and evaluating at $s=1$.

I suspect that it has to do with the fact that the most likely outcome is $\lfloor 1/(x y) \rfloor=1$ (which happens to be odd), followed with an application of the Strong Law of Small Numbers. Here are some more details. Let $Z$ be the random variable $xy$. Note that $Z$ takes values in $[0,1]$ and $\lfloor 1/Z \rfloor$ is odd if and only if

$$ Z \in (1/2, 1] \cup (1/4, 1/3] \cup (1/6, 1/5] \cup \dots $$

Note that this set has measure more than $1/2$. However, the distribution of $Z$ is of course not uniform on $[0,1]$ (it is actually skewed towards $0$ instead of $1$). I suspect what ends up happening is that the distribution of $\lfloor 1/Z \rfloor$ is almost perfectly split between even and odd for say $Z < 1/5$, and the discrepancy is thus a result of what happens for $Z \geq 1/5$ (where one easily sees that odd wins out).