Here is my optimized in-place Rust implementation [1].

It is a very tricky algorithm to get right. There are many edge cases that only happen for ~2^-64 fractions of the input space, so hard to find even with fuzz-testing. Best strategy is to implement it for small limbs first, and fuzz that well.

[1] https://github.com/0xProject/OpenZKP/blob/master/algebra/u25...