Edit2: I tried reading the Haskell code. (I don't know Haskell. I don't even know how many L's it has):
table =
[ ( 65536 / (2^32) , 1 - 121.0 ), ( 749731840 / (2^32) , 1 - 0.5 )
[...]
fun x = sum [ p * log (1 + b * x) | (p,b) <- table ]
It looks like you're using newton's method to maximize the sum of the products - but you're using the payout multipliers from the player's point of view. You need to look at it from the other side. The house never wins 121x. The player's big wins are the house's big losses. We need to calculate each house profit or loss as a multiplier of the amount the house is risking.
No, it's from the house's point of view. Notice the 1 - 121.0 in the second component of the first entry? It takes the stake of 1 into account and subtracts the payout to get the profit (-120) for the house in that case. Sorry I should've normalized the table a bit more.
Edit: Oh I see what you're doing. You're normalizing the stake for the house. Let me ponder about that some more.
Ok, I now believe that you are correct and my initial calculation is invalid. The payouts need indeed be normalized by the house's risk, i.e. the biggest amount the house can lose, which is 120 in this case. So the code isn't faulty per se, but the payout table was.