This is now fixed thanks to a discussion with Reinhard Burger and a more aggressive (albeit simplistic) classification of cases.

One pattern is that the totally recombinant offspring are all between zero and one. Perhaps it is only the cases that include nonrecombinant offspring that needs "adjustment."
As we can see, the values over 1 are in the row and column that produces nonrecombinant targets. This pattern occurs for every target offspring.
Here are the exact words from Burger:
"...where $R_I (j_Ij_J, k_Ik_J > i_Ii_J)$ is the probability that a randomly chosen gamete produced by a jk individual is i and each of the sets I and J is passed to the next generation without recombination." In other words, there is only recombination between some locus in I and some locus in J.

One solution is to view the R(jk>i) as a relative frequency and divide by the total summed across possible offspring. Then we'd get a probability. We could only apply this to values that have been calculated (e.g. the diagonal intersection value 1.0, is not calculated). I predict this would be computationally expensive, but it would be worth it. There may be a way to hold the values in cache and execute the division at the end instead of looping twice.
The test for this idea is this: does it still yield valid results under the tlta model?

As I suspected:
The formula is wrong somehow (or my interpretation of it). Unless the delta's are calculated differently from how I think they are, the formula is wrong.

Here's the output from the new algorithm:
promising, but I had already uncovered by hand that Bürger's formula appears to be wrong (or I completely misunderstand it).

Partitions are easily achieved by iterating over odd natural numbers from 1 (0001) to
(((1<<geno)  1) & ~2) & ((1 << geno)  1) (e.g. 1101).
It's so simple! No need for a separate function; no need for an extra data structure.

The algorithm is not guaranteed to work in any way. What we need is to follow an established algorithm instead of making wild guesses.
Now the task is to find an algorithm that can generate the partitions needed to use the formula in Reinhard Burger's book: R_{ij \rightarrow k} = \sum_I R_I
the algorithm we need is to generate the sets I, which are all the nonempty disjoint pairwise partitions. In other words, we need all the possible ways to split up the genome. Once we find that we can generate the deltas easily (differences between the genomes of the parents and the offspring).

I've made a few more changes to rec_total(), meaning we need to carefully check the threelocus result. The twolocus result is fine.

New algorithm identifies cases that need to call rec_total() completely; values appear correct for three loci. We cannot close this bug until we have a systematic check (such as an externallygenerated table that agrees).

A new algorithm is in place; however we cannot officially close this until we produce some way to check the threelocus case. The twolocus case is easy to check by hand. I will attempt to check it by hand tomorrow.

We need to verify the output of rec_test3; or we need to prove this rather than through calculation (e.g. Reiersol algebra).

rec_total: TLTA problem solved; we now need to reduce larger genome problems to the two locus case. That is, we need to loop over pairs of loci while properly multiplying the total probability

Fixed logic so that rec_gen_table no longer produces impossible offspring; however, inspection of table for offspring 0 from rec_test3 shows that probabilities are incorrect

Definition: impossible offspring  offspring that two parents cannot
produce through recombination
While developing my age dependent sexual selection simulation, I have
found that rec_gen_table produces nonzero offspring probabilities for
parents that cannot produce that particular offspring. For example,
for offspring 010 (2):
The entry at position (0,1) (and its transpose) is incorrect; it
should be zero. I discovered this from seeing that rec_mating was
producing offspring frequencies for offspring genotypes that should
not be possible: I entered allele frequencies of zero, and was getting
out nonzero allele frequencies (impossble, again).
The tabular output was produced by copying sparse.h into my source
directory and using the same procedure as in tests/prtable.c of the
haploid source tree.
